Search code examples
rggplot2geospatialr-sf

Longitude and latitude data not generating correct polygons using maps/mapdata/ggplot2 packages


I am trying to draw a map for all province within Egypt using maps and ggplot2 package in R, but I did not find any sub-region (provinces) for Egypt available in the databases of these packages. So, I downloaded coordinates (longitude and latitude) for all provinces in Egypt from google (and other website), and tried to generate the map, but it gave me a totally crumbled picture which is definitely not Egypt' provinces. This made me think that the coordinates from a searching engine (e.g. google) is not like the one used by R. I do understand that R is dealing with coordinates in +ve and negative values, whereas public databases is not. Do I need to convert it somehow to a "R" format ? Could any of you suggest a way around this, especially I did not find Egypt's province long/lat in the databases of R packages ?

Thanks

I have tried this code successfully for Egypt as a whole country

ggplot() + 
  geom_polygon(data = mapdata_Egypt, aes(x=long, y = lat, group = group), fill = NA, color = "red") + 
  coord_fixed(1.3)

But I could not find Egypt's province coordinates to draw them.


Solution

  • Whilst you have not provided links to where you sourced your data, or the code you are using to generate your mapdata_Egypt object, here is a repex that will achieve your goal.

    When dealing with point, line, or polygon geographical data, you should use the sf package. Manipulating and plotting geographical data is relatively easy with sf, and ggplot2 handles sf objects really well. Your issue is most likely a result of R not knowing in what order to plot your raw coordinates to construct your polygon features, hence your "totally crumbled picture". Using sf::st_read() avoids this issue as it directly translates geographical data, in this case a shapfile, to an object R understands e.g. and sf object.

    The data used in this example can be found here.

    Load packages and data into R

    library(sf)
    library(ggplot2)
    
    # Load shapefile previously unzipped into your working directory
    egypt_sf <- st_read("egy_admbnda_adm2_capmas_20170421.shp")
    
    # View subset of columns and rows in egypt_sf, misalignment of column values
    # due to mix of Arabic/English text
    egypt_sf[,1:4]
    # Simple feature collection with 365 features and 4 fields
    # Geometry type: MULTIPOLYGON
    # Dimension:     XY
    # Bounding box:  xmin: 24.69451 ymin: 21.99913 xmax: 36.90871 ymax: 31.67028
    # Geodetic CRS:  WGS 84
    # First 10 features:
    #           ADM2_EN                ADM2_AR ADM2_PCODE       ADM2_REF                      geometry
    # 1    10 Ramadan 1        قسم اول مدينة العاشر من رمض     EG1309   10 Ramadan 1  MULTIPOLYGON (((31.70323 30...
    # 2    10 Ramadan 2        قسم ثان مدينه العاشر من رمض     EG1321   10 Ramadan 2  MULTIPOLYGON (((31.72752 30...
    # 3         15 Mayu                قسم  15 مايو     EG0103        15 Mayu MULTIPOLYGON (((31.37739 29...
    # 4     6 October-1              قسم أول 6 أكتوبر      EG2107    6 October-1 MULTIPOLYGON (((30.98266 30...
    # 5     6 October-2              قسم ثان 6 أكتوبر      EG2121    6 October-2 MULTIPOLYGON (((30.954 29.9...
    # 6  A-Dakhla Oasis             مركز الوحات الداخلة      EG3202 A-Dakhla Oasis MULTIPOLYGON (((30.53004 27...
    # 7      A L Labban                  قسم  اللبان     EG0209     A L Labban MULTIPOLYGON (((29.89994 31...
    # 8           Abdin                  قسم عابدين     EG0109          Abdin MULTIPOLYGON (((31.24768 30...
    # 9           Abnub                  مركز أبنوب     EG2504          Abnub MULTIPOLYGON (((31.20966 27...
    # 10          Abour                   قسم العبو     EG1413          Abour MULTIPOLYGON (((31.4632 30...
    

    To plot, simply use the geom_sf() function from ggplot2

    ggplot() +
      geom_sf(data = egypt_sf)
    

    result

    Update based on OP's comment

    It is not entirely clear what it is you want exactly, but here is some code that hopefully reflects some of the things I believe you want:

    # Subset egypt_sf based on value in a column in egypt_sf
    cairo_sf <- filter(egypt_sf, ADM1_EN == "Cairo")
    
    cairo_sf[,c(1, 3:4, 9)]
    # Simple feature collection with 42 features and 4 fields
    # Geometry type: MULTIPOLYGON
    # Dimension:     XY
    # Bounding box:  xmin: 31.21676 ymin: 29.73371 xmax: 31.90697 ymax: 30.19757
    # Geodetic CRS:  WGS 84
    # First 10 features:
    #             ADM2_EN ADM2_PCODE         ADM2_REF ADM1_EN                       geometry
    # 1           15 Mayu     EG0103          15 Mayu   Cairo MULTIPOLYGON (((31.37739 29...
    # 2             Abdin     EG0109            Abdin   Cairo MULTIPOLYGON (((31.24768 30...
    # 3         Ain Shams     EG0130        Ain Shams   Cairo MULTIPOLYGON (((31.36298 30...
    # 4      Al Azbakiyya     EG0113     Al Azbakiyya   Cairo MULTIPOLYGON (((31.25389 30...
    # 5  Al Darb al-Ahmar     EG0114 Al Darb al-Ahmar   Cairo MULTIPOLYGON (((31.25604 30...
    # 6        Al Khalifa     EG0108       Al Khalifa   Cairo MULTIPOLYGON (((31.30736 30...
    # 7      Al Matariyya     EG0125     Al Matariyya   Cairo MULTIPOLYGON (((31.31662 30...
    # 8          Al Sahil     EG0121         Al Sahil   Cairo MULTIPOLYGON (((31.25724 30...
    # 9          Al Salam     EG0132         Al Salam   Cairo MULTIPOLYGON (((31.43298 30...
    # 10    Al Sharabiyya     EG0118    Al Sharabiyya   Cairo MULTIPOLYGON (((31.27195 30...
    
    # Plot subsetted sf
    ggplot() +
      geom_sf(data = cairo_sf)
    
    # Plot both sf datasets at full extent, colour cairo_sf
    ggplot() +
      geom_sf(data = egypt_sf) +
      geom_sf(data = cairo_sf, aes(fill = ADM1_EN)) +
      scale_fill_discrete(name = "Area of\nInterest")
    
    # Plot both sf datasets, colour cairo_sf, use coord_sf() to zoom in
    ggplot() +
      geom_sf(data = egypt_sf) +
      geom_sf(data = cairo_sf, aes(fill = ADM1_EN)) +
      coord_sf(xlim = c(31, 32),
               ylim = c(29.6, 30.4)) +
      scale_fill_discrete(name = "Area of\nInterest")