Spatial and Geographical data

While single geographical locations can be identified with a single pair of coordinates, the geographical extent of trade routes or countries are series of such coordinate values and cannot easily accommodated with the observation-by-variable format of data frames. Instead such geographical data constitutes its own data type, which is discussed in this chapter. The chapter also discusses how such spatial data can be connected with other data about geographical entities, e.g. the population density or GDP per capita of countries. Further it discusses the import of geographical data from shape files and OpenStreetMap OSM files and the definition of and conversion between cartographic projections.

Below is the supporting material for the various sections of the chapter.

The Structure of Spatial Data

  • The structure of objects of classes defined in the sp package

    • Script file: spatial-datastructures-sp.R

      The script makes use of the two packages:

      The sp package is a requirement of cshapes and provides the necessary data class definitions.

    • Interactive notebook:

      In [1]:
      options(jupyter.rich_display=FALSE) # Create output as usual in R
      
      In [2]:
      library(cshapes)
      
      Loading required package: sp
      
      Loading required package: maptools
      
      Checking rgeos availability: TRUE
      
      Loading required package: plyr
      
      
      In [3]:
      cshape.sp <- cshp()
      
      Warning message:
      “readShapePoly is deprecated; use rgdal::readOGR or sf::st_read”
      
      In [4]:
      # Obtaining the data for the Gambia from the Cshapes data base
      Gambia <- subset(cshape.sp,CNTRY_NAME=="The Gambia")
      
      In [5]:
      plot(Gambia)
      
      In [6]:
      class(Gambia)
      
      [1] "SpatialPolygonsDataFrame"
      attr(,"package")
      [1] "sp"
  • The structure of objects of classes defined in the sf package

    • Script file: spatial-datastructures-sf.R

      The script makes use of the sf package, which is available from https://cran.r-project.org/package=sf

      Data file used in the script: ged191.RData which is contained in the ZIP archive ged191-RData.zip downloaded from https://ucdp.uu.se/downloads/index.html#ged_global. At https://ucdp.uu.se/downloads/index.html#ged_global you will now find a newer version of these data, which however have a different format.

    • Interactive notebook:

      In [1]:
      options(jupyter.rich_display=FALSE) # Create output as usual in R
      
      In [2]:
      # The data used here is available from https://ucdp.uu.se/downloads/ged/ged191-RData.zip
      # download.file("https://ucdp.uu.se/downloads/ged/ged191-RData.zip",
      #              "ged191-RData.zip")
      
      In [3]:
      unzip("ged191-RData.zip")
      load("ged191.RData")
      
      In [4]:
      class(ged191)
      
      [1] "sf"         "data.frame"
      In [5]:
      # Number of observations
      nrow(ged191)
      
      [1] 152616
      In [6]:
      # Number of variables
      ncol(ged191)
      
      [1] 43
      In [7]:
      library(sf)
      
      Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0
      
      
      In [8]:
      options(width=140)
      # Using the adapted 'print()' method from the 'sf' package
      print(ged191)
      
      Simple feature collection with 152616 features and 42 fields
      geometry type:  POINT
      dimension:      XY
      bbox:           xmin: -117.3 ymin: -37.81361 xmax: 155.8967 ymax: 61.25
      CRS:            EPSG:4326
      First 10 features:
         year             side_b active_year                         side_a priogrid_gid conflict_new_id country_id      region
      1  2013              MUJAO           1             Government of Mali       150833           11347        432      Africa
      2  2004          Civilians           1                            FNI       132182             583        490      Africa
      3  2007             PARECO           0                           CNDP       127138            4600        490      Africa
      4  2008 Kashmir insurgents           1            Government of India       179070             364        750        Asia
      5  2008 Kashmir insurgents           1            Government of India       179069             364        750        Asia
      6  2011                PKK           1           Government of Turkey       184761             354        640 Middle East
      7  1997          Civilians           0                           ULFA       167587             523        750        Asia
      8  2013              APCLS           1 Government of DR Congo (Zaire)       127859             283        490      Africa
      9  2010         Al-Shabaab           1          Government of Somalia       132931             337        520      Africa
      10 2003               MDJT           1             Government of Chad       160235             288        483      Africa
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       source_article
      1        "Agence France Presse,2013-01-12,French pilot killed in Mali helicopter raid: defence ministry";"Agence France Presse,2013-01-28,Mali town counts dead after French strikes on Islamists";"Agence France Presse,2013-01-12,French gunships stop Mali Islamist advance";"Reuters News,2013-01-12,More than 100 killed in French air strikes and fighting in Mali";"Agence France Presse,2013-01-15,Eighty-six wounded in Mali fighting in Mopti, Gao: Red Cross";"Agence France Presse,2013-01-12,10 civilians including 3 children killed in Mali unrest: HRW";"Jeune Afrique Blog Défense,2013-01-30,Mali : retour sur la bataille décisive de Konna"
      2                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                Reuters 2004-01-23\nS/2004/573
      3                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   http://www.monuc.org/News.aspx?newsId=16672
      4                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                SATP, 2008-3-5
      5                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 AFP, 2008-7-1
      6                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  R 20111219, BBC Eur 20111219
      7                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          SATP
      8  "BBC Monitoring Africa,2013-03-10,UN mission backs DRCongo army to halt militia advance in east";"Business Day,2013-03-08,Fresh violence in Congo worries Red Cross";"All Africa,2013-03-06,Briefing - Militias in Masisi";"Voice of America Press Releases and Documents,2013-03-05,Risk of 'Ethnic War' in Eastern Congo Town";"Reuters News,2013-03-05,At least 70 dead, thousands flee fighting in east Congo";"Agence France Presse,2013-03-04,DRCongo rebels seize Kitchanga, 80 dead: UN";"All Africa,2013-03-06,Briefing - Militias in Masisi";"MSF,2013-03-07,Democratic Republic of Congo: Thousands have fled as violence continues in Kitchanga"
      9                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          Xinhua 1 July 2010 "Eleven die in Mogadishu fighting
      10                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          AFP 200303027, 29, 31, WMRC Daily Analysis 20030401
         type_of_violence date_prec
      1                 1         2
      2                 3         2
      3                 2         2
      4                 1         2
      5                 1         2
      6                 1         2
      7                 3         1
      8                 1         2
      9                 1         2
      10                1         2
                                                                                                                                                                                                                                                                                                                                                                                                    source_headline
      1  French pilot killed in Mali helicopter raid: defence ministry;Mali town counts dead after French strikes on Islamists;French gunships stop Mali Islamist advance;More than 100 killed in French air strikes and fighting in Mali;Eighty-six wounded in Mali fighting in Mopti, Gao: Red Cross;10 civilians including 3 children killed in Mali unrest: HRW;Mali : retour sur la bataille décisive de Konna
      2                                                                                                                                                                                                                                                                                                                                                                                                            
      3                                                                                                                                                                                                                                                                                                                                                                                                            
      4                                                                                                                                                                                                                                                                                                                                                                                                            
      5                                                                                                                                                                                                                                                                                                                                                                                                            
      6                                                                                                                                                                                                                                                                                                                                                                                                            
      7                                                                                                                                                                                                                                                                                                                                                                                                            
      8     UN mission backs DRCongo army to halt militia advance in east;Fresh violence in Congo worries Red Cross;Briefing - Militias in Masisi;Risk of 'Ethnic War' in Eastern Congo Town;At least 70 dead, thousands flee fighting in east Congo;DRCongo rebels seize Kitchanga, 80 dead: UN;Briefing - Militias in Masisi;Democratic Republic of Congo: Thousands have fled as violence continues in Kitchanga
      9                                                                                                                                                                                                                                                                                                                                                                                                            
      10                                                                                                                                                                                                                                                                                                                                                                                                           
                        conflict_name                                dyad_name event_clarity    where_coordinates   date_end side_b_new_id
      1              Mali: Government               Government of Mali - MUJAO             1           Konna town 2013/01/12          1161
      2               FNI - Civilians                          FNI - Civilians             1         Gobu village 2004/01/16             1
      3                 CNDP - PARECO                            CNDP - PARECO             1       Ngungu village 2007/12/18           896
      4                India: Kashmir Government of India - Kashmir insurgents             1 Chithi Bandi village 2008/03/05           325
      5                India: Kashmir Government of India - Kashmir insurgents             1     Kupwara district 2008/07/01           325
      6             Turkey: Kurdistan               Government of Turkey - PKK             1       Dicle district 2011/12/19           323
      7              ULFA - Civilians                         ULFA - Civilians             1          Assam State 1997/07/04             1
      8  DR Congo (Zaire): Government   Government of DR Congo (Zaire) - APCLS             1       Kitchanga town 2013/03/05          1200
      9           Somalia: Government       Government of Somalia - Al-Shabaab             1       Mogadishu city 2010/07/01           717
      10             Chad: Government                Government of Chad - MDJT             1          Bardai town 2003/03/27           452
         dyad_new_id high
      1        12571   50
      2         1050  200
      3         5210    4
      4          792    2
      5          792   12
      6          781   30
      7          990    1
      8        10509   80
      9          750   11
      10         616   72
                                                                                                                                            source_office
      1  Agence France Presse;Agence France Presse;Agence France Presse;Reuters News;Agence France Presse;Agence France Presse;Jeune Afrique Blog Défense
      2                                                                                                                                                  
      3                                                                                                                                                  
      4                                                                                                                                                  
      5                                                                                                                                                  
      6                                                                                                                                                  
      7                                                                                                                                                  
      8      BBC Monitoring Africa;Business Day;All Africa;Voice of America Press Releases and Documents;Reuters News;Agence France Presse;All Africa;MSF
      9                                                                                                                                                  
      10                                                                                                                                                 
         gwnob deaths_civilians best number_of_sources                                                                             source_date
      1   <NA>                3   31                 7            2013-01-12;2013-01-28;2013-01-12;2013-01-12;2013-01-15;2013-01-12;2013-01-30
      2   <NA>              200  200                -1                                                                                        
      3   <NA>                4    4                -1                                                                                        
      4   <NA>                0    2                -1                                                                                        
      5   <NA>                0   12                -1                                                                                        
      6   <NA>                0   13                -1                                                                                        
      7   <NA>                1    1                -1                                                                                        
      8   <NA>                0   80                 8 2013-03-10;2013-03-08;2013-03-06;2013-03-05;2013-03-05;2013-03-04;2013-03-06;2013-03-07
      9   <NA>                0   11                -1                                                                                        
      10  <NA>                0    2                -1                                                                                        
         longitude where_prec side_a_new_id     id                                        source_original                        adm_1 gwnoa
      1   -3.89474          1            72  67972 French Defence Minister Jean-Yves Le Drian;  residents                 Mopti region   432
      2   30.77583          1           606  23385                                         local official               Ituri province  <NA>
      3   28.87694          2           426  24255                                                                  Nord Kivu province  <NA>
      4   74.72549          1           141  82612                                                   SATP      Jammu and Kashmir State   750
      5   74.16729          3           141  82645                                                   army      Jammu and Kashmir State   750
      6   40.08422          3           115 158618            senior military official, Cihan news agency          Diyarbakır province   640
      7   93.00000          4           326  87396                                                   SATP                  Assam State  <NA>
      8   29.08000          1            89  60774                      deputy UN spokesman Eduardo Buey.           Nord Kivu province   490
      9   45.36667          1            95  25360                                     emergency services              Banaadir region   520
      10  17.00055          1            87  11995                                                        Borkou-Ennedi-Tibesti region   483
                  country              adm_2  latitude                    geom_wkt date_start deaths_a deaths_b deaths_unknown low
      1              Mali       Mopti cercle 14.943290 POINT (-3.894740 14.943290) 2013/01/11       12       16              0  31
      2  DR Congo (Zaire)    Djugu territory  1.769722  POINT (30.775833 1.769722) 2004/01/14        0        0              0 100
      3  DR Congo (Zaire)   Masisi territory -1.651944 POINT (28.876944 -1.651944) 2007/12/13        0        0              0   4
      4             India Bandipore district 34.378934 POINT (74.725490 34.378934) 2008/03/04        0        2              0   2
      5             India   Kupwara district 34.372601 POINT (74.167293 34.372601) 2008/06/30        1       11              0  12
      6            Turkey     Dicle district 38.382704 POINT (40.084215 38.382704) 2011/12/17        0       13              0  13
      7             India                    26.000000 POINT (93.000000 26.000000) 1997/07/04        0        0              0   1
      8  DR Congo (Zaire) Rutshuru territory -1.250000 POINT (29.080000 -1.250000) 2013/03/03        0        0             80  70
      9           Somalia Mogadishu district  2.066667  POINT (45.366667 2.066667) 2010/06/30        0        0             11  11
      10             Chad Tibesti department 21.358333 POINT (17.000550 21.358333) 2003/03/25        0        2              0   2
                           geometry
      1   POINT (-3.89474 14.94329)
      2   POINT (30.77583 1.769722)
      3  POINT (28.87694 -1.651944)
      4   POINT (74.72549 34.37893)
      5    POINT (74.16729 34.3726)
      6    POINT (40.08422 38.3827)
      7               POINT (93 26)
      8         POINT (29.08 -1.25)
      9   POINT (45.36667 2.066667)
      10  POINT (17.00055 21.35833)
      
      In [9]:
      # Examining the structure of an "sf" object with 'str()'
      str(ged191)
      
      Classes ‘sf’ and 'data.frame':	152616 obs. of  43 variables:
       $ year             : num  2013 2004 2007 2008 2008 ...
       $ side_b           : chr  "MUJAO" "Civilians" "PARECO" "Kashmir insurgents" ...
       $ active_year      : num  1 1 0 1 1 1 0 1 1 1 ...
       $ side_a           : chr  "Government of Mali" "FNI" "CNDP" "Government of India" ...
       $ priogrid_gid     : num  150833 132182 127138 179070 179069 ...
       $ conflict_new_id  : num  11347 583 4600 364 364 ...
       $ country_id       : num  432 490 490 750 750 640 750 490 520 483 ...
       $ region           : chr  "Africa" "Africa" "Africa" "Asia" ...
       $ source_article   : chr  "\"Agence France Presse,2013-01-12,French pilot killed in Mali helicopter raid: defence ministry\";\"Agence Fran"| __truncated__ "Reuters 2004-01-23\nS/2004/573" "http://www.monuc.org/News.aspx?newsId=16672" "SATP, 2008-3-5" ...
       $ type_of_violence : num  1 3 2 1 1 1 3 1 1 1 ...
       $ date_prec        : num  2 2 2 2 2 2 1 2 2 2 ...
       $ source_headline  : chr  "French pilot killed in Mali helicopter raid: defence ministry;Mali town counts dead after French strikes on Isl"| __truncated__ "" "" "" ...
       $ conflict_name    : chr  "Mali: Government" "FNI - Civilians" "CNDP - PARECO" "India: Kashmir" ...
       $ dyad_name        : chr  "Government of Mali - MUJAO" "FNI - Civilians" "CNDP - PARECO" "Government of India - Kashmir insurgents" ...
       $ event_clarity    : num  1 1 1 1 1 1 1 1 1 1 ...
       $ where_coordinates: chr  "Konna town" "Gobu village" "Ngungu village" "Chithi Bandi village" ...
       $ date_end         : chr  "2013/01/12" "2004/01/16" "2007/12/18" "2008/03/05" ...
       $ side_b_new_id    : num  1161 1 896 325 325 ...
       $ dyad_new_id      : num  12571 1050 5210 792 792 ...
       $ high             : num  50 200 4 2 12 30 1 80 11 72 ...
       $ source_office    : chr  "Agence France Presse;Agence France Presse;Agence France Presse;Reuters News;Agence France Presse;Agence France "| __truncated__ "" "" "" ...
       $ gwnob            : chr  NA NA NA NA ...
       $ deaths_civilians : num  3 200 4 0 0 0 1 0 0 0 ...
       $ best             : num  31 200 4 2 12 13 1 80 11 2 ...
       $ number_of_sources: num  7 -1 -1 -1 -1 -1 -1 8 -1 -1 ...
       $ source_date      : chr  "2013-01-12;2013-01-28;2013-01-12;2013-01-12;2013-01-15;2013-01-12;2013-01-30" "" "" "" ...
       $ longitude        : num  -3.89 30.78 28.88 74.73 74.17 ...
       $ where_prec       : num  1 1 2 1 3 3 4 1 1 1 ...
       $ side_a_new_id    : num  72 606 426 141 141 115 326 89 95 87 ...
       $ id               : num  67972 23385 24255 82612 82645 ...
       $ source_original  : chr  "French Defence Minister Jean-Yves Le Drian;  residents" "local official" "" "SATP" ...
       $ adm_1            : chr  "Mopti region" "Ituri province" "Nord Kivu province" "Jammu and Kashmir State" ...
       $ gwnoa            : chr  "432" NA NA "750" ...
       $ country          : chr  "Mali" "DR Congo (Zaire)" "DR Congo (Zaire)" "India" ...
       $ adm_2            : chr  "Mopti cercle" "Djugu territory" "Masisi territory" "Bandipore district" ...
       $ latitude         : num  14.94 1.77 -1.65 34.38 34.37 ...
       $ geom_wkt         : chr  "POINT (-3.894740 14.943290)" "POINT (30.775833 1.769722)" "POINT (28.876944 -1.651944)" "POINT (74.725490 34.378934)" ...
       $ date_start       : chr  "2013/01/11" "2004/01/14" "2007/12/13" "2008/03/04" ...
       $ deaths_a         : num  12 0 0 0 1 0 0 0 0 0 ...
       $ deaths_b         : num  16 0 0 2 11 13 0 0 0 2 ...
       $ deaths_unknown   : num  0 0 0 0 0 0 0 80 11 0 ...
       $ low              : num  31 100 4 2 12 13 1 70 11 2 ...
       $ geometry         :sfc_POINT of length 152616; first list element:  'XY' num  -3.89 14.94
       - attr(*, "sf_column")= chr "geometry"
      
      In [10]:
      # Obtaining the geometry contained in 'ged191'
      gged191 <- st_geometry(ged191)
      gged191
      
      POINT (-3.89474 14.94329)
      
      POINT (30.77583 1.769722)
      
      POINT (28.87694 -1.651944)
      
      POINT (74.72549 34.37893)
      
      POINT (74.16729 34.3726)
      
      
      Geometry set for 152616 features 
      geometry type:  POINT
      dimension:      XY
      bbox:           xmin: -117.3 ymin: -37.81361 xmax: 155.8967 ymax: 61.25
      CRS:            EPSG:4326
      First 5 geometries:
      In [11]:
      # Accessing an individual shape
      gged191[[1]]
      
      POINT (-3.89474 14.94329)
      
      
      In [12]:
      # The class of a shape object
      class(gged191[[1]])
      
      [1] "XY"    "POINT" "sfg"  
  • Converting “sp” objects to “sf” objects:

    • Script file: converting-sp-to-sf.R

      The script makes use of the two packages:

    • Interactive notebook:

      In [1]:
      options(jupyter.rich_display=FALSE) # Create output as usual in R
      
      In [2]:
      library(sf)
      
      Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0
      
      
      In [3]:
      library(cshapes)
      
      Loading required package: sp
      
      Loading required package: maptools
      
      Checking rgeos availability: TRUE
      
      Loading required package: plyr
      
      
      In [4]:
      cshapes.1990 <- cshp(as.Date("1990-01-01"))
      cshapes.1990 <- as(cshapes.1990,"sf")
      
      Warning message:
      “readShapePoly is deprecated; use rgdal::readOGR or sf::st_read”
      
      In [5]:
      options(width=200)
      print(cshapes.1990[1:10])
      
      Simple feature collection with 171 features and 10 fields
      geometry type:  MULTIPOLYGON
      dimension:      XY
      bbox:           xmin: -180 ymin: -55.90223 xmax: 180 ymax: 83.11387
      CRS:            +proj=longlat +ellps=WGS84
      First 10 features:
                 CNTRY_NAME         AREA       CAPNAME CAPLONG     CAPLAT FEATUREID COWCODE COWSYEAR COWSMONTH COWSDAY                       geometry
      0              Guyana  211982.0050    Georgetown   -58.2   6.800000         0     110     1966         5      26 MULTIPOLYGON (((-58.17262 6...
      1            Suriname  145952.2740    Paramaribo   -55.2   5.833333         1     115     1975        11      25 MULTIPOLYGON (((-55.12796 5...
      2 Trinidad and Tobago    5041.7290 Port-of-Spain   -61.5  10.650000         2      52     1962         8      31 MULTIPOLYGON (((-61.07945 1...
      3           Venezuela  916782.2172       Caracas   -66.9  10.500000         3     101     1946         1       1 MULTIPOLYGON (((-66.31029 1...
      4               Samoa    2955.2124          Apia  -172.0 -13.800000         4     990     1976        12      15 MULTIPOLYGON (((-172.5965 -...
      5               Tonga     464.7473    Nuku'alofa  -175.0 -21.100000         5     955     1999         9      14 MULTIPOLYGON (((-175.1453 -...
      6           Argentina 2787442.0977  Buenos Aires   -58.7 -34.600000         6     160     1946         1       1 MULTIPOLYGON (((-71.85916 -...
      7             Bolivia 1092697.4356        La Paz   -68.2 -16.500000         7     145     1946         1       1 MULTIPOLYGON (((-62.19884 -...
      8              Brazil 8523619.5715      Brasilia   -47.9 -15.800000         8     140     1960         4      21 MULTIPOLYGON (((-44.69501 -...
      9               Chile  745808.4936      Santiago   -70.7 -33.500000         9     155     1946         1       1 MULTIPOLYGON (((-73.0421 -4...
      
      In [6]:
      SthAmCntry.names <- c(
          "Argentina",
          "Bolivia",
          "Brazil",
          "Chile",
          "Colombia",
          "Ecuador",
          "Guyana",
          "Paraguay",
          "Peru",
          "Suriname",
          "Uruguay",
          "Venezuela")
      
      In [7]:
      SthAmCountries <-
          subset(cshapes.1990,
                 CNTRY_NAME %in% SthAmCntry.names)
      
      In [8]:
      Brazil <- subset(cshapes.1990,CNTRY_NAME=="Brazil")
      Chile <-  subset(cshapes.1990,CNTRY_NAME=="Chile")
      Colombia <-  subset(cshapes.1990,CNTRY_NAME=="Colombia")
      
      In [9]:
      cap.latlong <- with(cshapes.1990,cbind(CAPLONG,CAPLAT))
      
      In [10]:
      cap.latlong <- lapply(1:nrow(cap.latlong),
                            function(i)cap.latlong[i,])
      
      In [11]:
      cap.latlong <- lapply(cap.latlong,st_point)
      cap.latlong <- st_sfc(cap.latlong)
      
      In [12]:
      cshapes.capitals.1990 <- cshapes.1990
      st_geometry(cshapes.capitals.1990) <- cap.latlong
      
      In [13]:
      st_crs(cshapes.capitals.1990) <- st_crs(cshapes.1990)
      
      In [14]:
      print(cshapes.capitals.1990[1:10])
      
      Simple feature collection with 171 features and 10 fields
      geometry type:  POINT
      dimension:      XY
      bbox:           xmin: -175 ymin: -41.3 xmax: 179 ymax: 64.15
      CRS:            +proj=longlat +ellps=WGS84
      First 10 features:
                 CNTRY_NAME         AREA       CAPNAME CAPLONG     CAPLAT FEATUREID COWCODE COWSYEAR COWSMONTH COWSDAY               geometry
      0              Guyana  211982.0050    Georgetown   -58.2   6.800000         0     110     1966         5      26      POINT (-58.2 6.8)
      1            Suriname  145952.2740    Paramaribo   -55.2   5.833333         1     115     1975        11      25 POINT (-55.2 5.833333)
      2 Trinidad and Tobago    5041.7290 Port-of-Spain   -61.5  10.650000         2      52     1962         8      31    POINT (-61.5 10.65)
      3           Venezuela  916782.2172       Caracas   -66.9  10.500000         3     101     1946         1       1     POINT (-66.9 10.5)
      4               Samoa    2955.2124          Apia  -172.0 -13.800000         4     990     1976        12      15     POINT (-172 -13.8)
      5               Tonga     464.7473    Nuku'alofa  -175.0 -21.100000         5     955     1999         9      14     POINT (-175 -21.1)
      6           Argentina 2787442.0977  Buenos Aires   -58.7 -34.600000         6     160     1946         1       1    POINT (-58.7 -34.6)
      7             Bolivia 1092697.4356        La Paz   -68.2 -16.500000         7     145     1946         1       1    POINT (-68.2 -16.5)
      8              Brazil 8523619.5715      Brasilia   -47.9 -15.800000         8     140     1960         4      21    POINT (-47.9 -15.8)
      9               Chile  745808.4936      Santiago   -70.7 -33.500000         9     155     1946         1       1    POINT (-70.7 -33.5)
      
      In [15]:
      Brasilia <- subset(cshapes.capitals.1990,CNTRY_NAME=="Brazil")
      Santiago <-  subset(cshapes.capitals.1990,CNTRY_NAME=="Chile")
      Bogota <-  subset(cshapes.capitals.1990,CNTRY_NAME=="Colombia")
      
      In [16]:
      graypal <- function(n)gray.colors(n,start=.2,end=.9,alpha=.5)
      plot(SthAmCountries,pal=graypal)
      
      Warning message:
      “plotting the first 9 out of 24 attributes; use max.plot = 24 to plot all”
      
      In [17]:
      plot(st_geometry(SthAmCountries))
      
      In [18]:
      plot(st_geometry(Brazil))
      
      In [19]:
      save(cshapes.1990,cshapes.capitals.1990,file="cshapes-1990.RData")
      save(Brazil,Chile,Colombia,
           Brasilia,Santiago,Bogota,
           SthAmCountries,
           file="south-america-1990.RData")
      

Spatial Relations and Operations

  • Spatial relations

    • Script file: spatial-relations.R

      The script makes use of the two packages:

      Data file used in the script: south-america-1990.RData created by earlier script converting-sp-to-sf.R

    • Interactive notebook:

      In [1]:
      options(jupyter.rich_display=FALSE) # Create output as usual in R
      
      In [2]:
      library(sf)
      
      Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0
      
      
      In [3]:
      library(cshapes)
      
      Loading required package: sp
      
      Loading required package: maptools
      
      Checking rgeos availability: TRUE
      
      Loading required package: plyr
      
      
      In [4]:
      load("south-america-1990.RData")
      
      In [5]:
      st_contains(Brazil,Brasilia,sparse=FALSE)
      
      although coordinates are longitude/latitude, st_contains assumes that they are planar
      
      
           [,1]
      [1,] TRUE
      In [6]:
      st_contains(Brazil,Bogota,sparse=FALSE)
      
      although coordinates are longitude/latitude, st_contains assumes that they are planar
      
      
           [,1] 
      [1,] FALSE
      In [7]:
      load("cshapes-1990.RData")
      ThreeCountries <-
          subset(cshapes.1990,
                 CNTRY_NAME %in% c("Brazil","Chile","Colombia"))
      rownames(ThreeCountries) <- ThreeCountries$CNTRY_NAME
      
      In [8]:
      ThreeCapitals <-
          subset(cshapes.capitals.1990,
                 CNTRY_NAME %in% c("Brazil","Chile","Colombia"))
      rownames(ThreeCapitals) <- ThreeCapitals$CAPNAME
      
      In [9]:
      st_contains(ThreeCountries,ThreeCapitals,sparse=FALSE)
      
      although coordinates are longitude/latitude, st_contains assumes that they are planar
      
      
           [,1]  [,2]  [,3] 
      [1,]  TRUE FALSE FALSE
      [2,] FALSE  TRUE FALSE
      [3,] FALSE FALSE  TRUE
      In [10]:
      structure(
         st_contains(ThreeCountries,ThreeCapitals,sparse=FALSE),
         dimnames=list(rownames(ThreeCountries),rownames(ThreeCapitals))
      )
      
      although coordinates are longitude/latitude, st_contains assumes that they are planar
      
      
               Brasilia Santiago Bogota
      Brazil    TRUE    FALSE    FALSE 
      Chile    FALSE     TRUE    FALSE 
      Colombia FALSE    FALSE     TRUE 
      In [11]:
      st_touches(Brazil,Colombia,sparse=FALSE)
      
      although coordinates are longitude/latitude, st_touches assumes that they are planar
      
      
           [,1]
      [1,] TRUE
      In [12]:
      st_touches(Brazil,Chile,sparse=FALSE)
      
      although coordinates are longitude/latitude, st_touches assumes that they are planar
      
      
           [,1] 
      [1,] FALSE
      In [13]:
      structure(
         st_touches(ThreeCountries,sparse=FALSE),
         dimnames=list(rownames(ThreeCountries),rownames(ThreeCountries))
      )
      
      although coordinates are longitude/latitude, st_touches assumes that they are planar
      
      
               Brazil Chile Colombia
      Brazil   FALSE  FALSE  TRUE   
      Chile    FALSE  FALSE FALSE   
      Colombia  TRUE  FALSE FALSE   
      In [14]:
      save(ThreeCountries,ThreeCapitals,file="three-countries.RData")
      
  • Properties of spatial objects

    • Script file: spatial-proterties.R

      The script makes use of the two packages:

      Data files used in the script:

      created by earlier scripts spatial-datastructures-sf.R and spatial-relations.R

    • Interactive notebook:

      In [1]:
      options(jupyter.rich_display=FALSE) # Create output as usual in R
      
      In [2]:
      library(sf)
      
      Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0
      
      
      In [3]:
      load("south-america-1990.RData")
      load("three-countries.RData")
      
      In [4]:
      st_area(Brazil)
      
      8.472314e+12 [m^2]
      In [5]:
      in_km2 <- function(x) units::set_units(x,"km^2")
      in_km2(st_area(Brazil))
      
      8472314 [km^2]
      In [6]:
      in_km2(st_area(SthAmCountries))
      
      Units: [km^2]
       [1]  210585.6  144986.2  910860.6 2780991.8 1086612.6 8472314.3  744390.5
       [8]  255309.8  398803.1 1290857.6  177861.7 1135173.0
      In [7]:
      structure(in_km2(st_area(SthAmCountries)),
                names=as.character(SthAmCountries$CNTRY_NAME))
      
      Units: [km^2]
         Guyana  Suriname Venezuela Argentina   Bolivia    Brazil     Chile   Ecuador 
       210585.6  144986.2  910860.6 2780991.8 1086612.6 8472314.3  744390.5  255309.8 
       Paraguay      Peru   Uruguay  Colombia 
       398803.1 1290857.6  177861.7 1135173.0 
      In [8]:
      st_distance(Brasilia,Bogota)
      
      Units: [m]
              [,1]
      [1,] 3663768
      In [9]:
      st_distance(Chile,Bogota)
      
      Units: [m]
              [,1]
      [1,] 2496612
      In [10]:
      # This takes a while, because R needs to figure out which points of the borders are the closest to one another
      st_distance(Chile,Colombia)
      
      Units: [m]
              [,1]
      [1,] 1468577
      In [11]:
      in_km <- function(x) units::set_units(x,"km")
      in_km(st_distance(Brasilia,Bogota))
      
      Units: [km]
               [,1]
      [1,] 3663.768
      In [12]:
      in_km(st_distance(ThreeCapitals))
      
      Units: [km]
               [,1]     [,2]     [,3]
      [1,]    0.000 3014.942 3663.768
      [2,] 3014.942    0.000 4232.052
      [3,] 3663.768 4232.052    0.000
  • Unions and intersections (of artificial geometric objects)

    • Script file: artificial-unions-intersections.R

      The script makes use of the sf package, which is available from https://cran.r-project.org/package=sf.

    • Interactive notebook:

      In [1]:
      options(jupyter.rich_display=FALSE) # Create output as usual in R
      
      In [2]:
      library(sf)
      
      Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0
      
      
      In [3]:
      A <- rbind(
          c(0,0),
          c(2,0),
          c(0,2),
          c(0,0)
      )
      A <- st_polygon(list(A))
      B <- rbind(
          c( 1, 1),
          c( 1,-1),
          c(-1,-1),
          c(-1, 1),
          c( 1, 1)
      )
      B <- st_polygon(list(B))
      
      In [4]:
      A
      
      POLYGON ((0 0, 2 0, 0 2, 0 0))
      
      
      In [5]:
      B
      
      POLYGON ((1 1, 1 -1, -1 -1, -1 1, 1 1))
      
      
      In [6]:
      op <- par(mai=c(0,0,0,0),mfrow=c(1,3),xpd=NA)
      plot(A,xlim=c(-1,2),ylim=c(-1,2))
      plot(B,lty=2,add=TRUE)
      text(0,1.5,"A",pos=2)
      text(1,-.5,"B",pos=4)
      text(.5,-1.5,"Two shapes A and B",pos=1)
      plot(st_union(A,B),col="gray70",xlim=c(-1,2),ylim=c(-1,2))
      plot(A,lty=3,add=TRUE)
      plot(B,lty=3,add=TRUE)
      text(.5,-1.5,"st_union(A,B)",pos=1)
      plot(A,lty=3,xlim=c(-1,2),ylim=c(-1,2))
      plot(B,lty=3,add=TRUE)
      plot(st_intersection(A,B),add=TRUE,col="gray70")
      text(.5,-1.5,"st_intersection(A,B)",pos=1)
      par(op)
      
      In [7]:
      op <- par(mai=c(0,0,0,0),mfrow=c(1,3),xpd=NA)
      plot(st_difference(A,B),col="gray70",,xlim=c(-1,2),ylim=c(-1,2))
      plot(A,lty=3,add=TRUE)
      plot(B,lty=3,add=TRUE)
      text(.5,-1.5,"st_difference(A,B)",pos=1)
      plot(st_difference(B,A),col="gray70",,xlim=c(-1,2),ylim=c(-1,2))
      plot(A,lty=3,add=TRUE)
      plot(B,lty=3,add=TRUE)
      text(.5,-1.5,"st_difference(B,A)",pos=1)
      plot(st_sym_difference(A,B),col="gray70",,xlim=c(-1,2),ylim=c(-1,2))
      text(.5,-1.5,"st_sym_difference(A,B)",pos=1)
      par(op)
      
  • Subsetting spatial feature objects

    • Script file: spatial-subsets.R

      The script makes use of the sf package, which is available from https://cran.r-project.org/package=sf.

      Data files used in the script:

      created by earlier scripts spatial-datastructures-sf.R and converting-sp-to-sf.R

    • Interactive notebook:

      In [1]:
      options(jupyter.rich_display=FALSE) # Create output as usual in R
      
      In [2]:
      library(sf)
      
      Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0
      
      
      In [3]:
      load("south-america-1990.RData")
      load("ged191.RData")
      
      In [4]:
      # This fails due to different coordinate reference systems
      Colombia.conflicts <- ged191[Colombia,]
      
      Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared): st_crs(x) == st_crs(y) is not TRUE
      Traceback:
      
      1. ged191[Colombia, ]
      2. `[.sf`(ged191, Colombia, )
      3. lengths(op(x, i, ...))
      4. op(x, i, ...)
      5. st_intersects.sf(x, i, ...)
      6. st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared)
      7. stopifnot(st_crs(x) == st_crs(y))
      In [5]:
      st_crs(Colombia)
      
      Coordinate Reference System:
        User input: +proj=longlat +ellps=WGS84 
        wkt:
      GEOGCS["WGS 84",
          DATUM["unknown",
              SPHEROID["WGS84",6378137,298.257223563]],
          PRIMEM["Greenwich",0],
          UNIT["degree",0.0174532925199433]]
      In [6]:
      st_crs(ged191)
      
      Coordinate Reference System:
        User input: EPSG:4326 
        wkt:
      GEOGCS["WGS 84",
          DATUM["WGS_1984",
              SPHEROID["WGS 84",6378137,298.257223563,
                  AUTHORITY["EPSG","7030"]],
              AUTHORITY["EPSG","6326"]],
          PRIMEM["Greenwich",0,
              AUTHORITY["EPSG","8901"]],
          UNIT["degree",0.0174532925199433,
              AUTHORITY["EPSG","9122"]],
          AUTHORITY["EPSG","4326"]]
      In [7]:
      # Here we make sure that the coordinate reference systems match
      ged191_ellips <- st_transform(ged191,st_crs(Colombia))
      Colombia.conflicts <- ged191_ellips[Colombia,]
      
      although coordinates are longitude/latitude, st_intersects assumes that they are planar
      
      although coordinates are longitude/latitude, st_intersects assumes that they are planar
      
      
      In [8]:
      mypal <- function(n)gray.colors(n,start=.2,end=.9,alpha=.5)
      plot(st_geometry(Colombia))
      plot(Colombia.conflicts["deaths_civilians"],
           add=TRUE,pch=19,cex=.2,
           pal=mypal,
           nbreaks=30
           )
      
      In [9]:
      st_circ <- function(x,dist.km){
          dist.degr <- 360*dist.km/40007.863
          st_buffer(st_geometry(x),dist=dist.degr)
      }
      
      In [10]:
      Bogota.200km <- st_circ(Bogota,dist.km=200)
      
      Warning message in st_buffer.sfc(st_geometry(x), dist = dist.degr):
      “st_buffer does not correctly buffer longitude/latitude data”
      dist is assumed to be in decimal degrees (arc_degrees).
      
      
      In [11]:
      Bogota.conflicts <- ged191_ellips[Bogota.200km,]
      
      although coordinates are longitude/latitude, st_intersects assumes that they are planar
      
      although coordinates are longitude/latitude, st_intersects assumes that they are planar
      
      
      In [12]:
      plot(st_geometry(Colombia))
      plot(st_geometry(Colombia.conflicts),
           add=TRUE,pch=1,cex=.3,col="gray80"
           )
      plot(st_geometry(Bogota.conflicts),
           add=TRUE,pch=19,cex=.3)
      plot(Bogota.200km,lty=3,add=TRUE)
      
  • Aggregating spatial feature objects

    • Script file: spatial-aggregates.R

      The script makes use of the sf package, which is available from https://cran.r-project.org/package=sf.

      Data files used in the script:

      created by earlier scripts spatial-datastructures-sf.R and converting-sp-to-sf.R

    • Interactive notebook:

      In [1]:
      options(jupyter.rich_display=FALSE) # Create output as usual in R
      
      In [2]:
      library(sf)
      
      Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0
      
      
      In [3]:
      load("south-america-1990.RData")
      load("ged191.RData")
      load("cshapes-1990.RData")
      
      In [4]:
      ged191_ellips <- st_transform(ged191,st_crs(cshapes.1990))
      
      In [5]:
      # Civilian deaths per country
      aggregate(ged191_ellips["deaths_civilians"],by=SthAmCountries,sum)
      
      although coordinates are longitude/latitude, st_intersects assumes that they are planar
      
      
         deaths_civilians geometry                      
      0    20             MULTIPOLYGON (((-58.17262 6...
      1    NA             MULTIPOLYGON (((-55.12796 5...
      3   792             MULTIPOLYGON (((-66.31029 1...
      6     0             MULTIPOLYGON (((-71.85916 -...
      7     0             MULTIPOLYGON (((-62.19884 -...
      8   164             MULTIPOLYGON (((-44.69501 -...
      9    NA             MULTIPOLYGON (((-73.0421 -4...
      10   15             MULTIPOLYGON (((-79.45167 -...
      11    0             MULTIPOLYGON (((-57.67267 -...
      12 1021             MULTIPOLYGON (((-69.49973 -...
      13   NA             MULTIPOLYGON (((-58.38889 -...
      27 6015             MULTIPOLYGON (((-74.86081 1...
      In [6]:
      # Civilian deaths per country, with country names
      within(
          aggregate(ged191_ellips["deaths_civilians"],by=SthAmCountries,sum),
          country <- SthAmCountries$CNTRY_NAME)
      
      although coordinates are longitude/latitude, st_intersects assumes that they are planar
      
      
         deaths_civilians geometry                       country  
      0    20             MULTIPOLYGON (((-58.17262 6... Guyana   
      1    NA             MULTIPOLYGON (((-55.12796 5... Suriname 
      3   792             MULTIPOLYGON (((-66.31029 1... Venezuela
      6     0             MULTIPOLYGON (((-71.85916 -... Argentina
      7     0             MULTIPOLYGON (((-62.19884 -... Bolivia  
      8   164             MULTIPOLYGON (((-44.69501 -... Brazil   
      9    NA             MULTIPOLYGON (((-73.0421 -4... Chile    
      10   15             MULTIPOLYGON (((-79.45167 -... Ecuador  
      11    0             MULTIPOLYGON (((-57.67267 -... Paraguay 
      12 1021             MULTIPOLYGON (((-69.49973 -... Peru     
      13   NA             MULTIPOLYGON (((-58.38889 -... Uruguay  
      27 6015             MULTIPOLYGON (((-74.86081 1... Colombia 
      In [7]:
      st_circ <- function(x,dist.km){
          dist.degr <- 360*dist.km/40007.863
          st_buffer(st_geometry(x),dist=dist.degr)
      }
      
      In [8]:
      Bogota.region <- st_circ(Bogota,dist.km=200)
      Colombia.rest <- st_difference(st_geometry(Colombia),Bogota.region)
      
      Warning message in st_buffer.sfc(st_geometry(x), dist = dist.degr):
      “st_buffer does not correctly buffer longitude/latitude data”
      dist is assumed to be in decimal degrees (arc_degrees).
      
      although coordinates are longitude/latitude, st_difference assumes that they are planar
      
      
      In [9]:
      # Civilian deaths in the Bogota region and the rest of Colombia
      aggregate(ged191_ellips["deaths_civilians"],
                by=c(Bogota.region,Colombia.rest),
                sum)
      
      although coordinates are longitude/latitude, st_intersects assumes that they are planar
      
      
        deaths_civilians geometry                      
      1 1021             POLYGON ((-72.30035 4.6, -7...
      2 4994             MULTIPOLYGON (((-74.86081 1...

Cartographic Projections

  • Script file: cartographic-projections.R

    The script makes use of the sf package, which is available from https://cran.r-project.org/package=sf.

    Data file used in the script: US_flat.RData which contains prepared data originally from the maps package, available from https://cran.r-project.org/package=maps.

  • Interactive notebook:

    Geographical projections of a mainland USA map (with apologies to Alaskans and Hawaiians)

    In [1]:
    options(jupyter.rich_display=FALSE) # Create output as usual in R
    
    In [2]:
    library(sf)
    
    Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0
    
    
    In [3]:
    load("US_flat.RData")
    US_flat
    
       geometry                       State               
    1  MULTIPOLYGON (((-87.46201 3... Alabama             
    2  MULTIPOLYGON (((-114.6374 3... Arizona             
    3  MULTIPOLYGON (((-94.05103 3... Arkansas            
    4  MULTIPOLYGON (((-120.006 42... California          
    5  MULTIPOLYGON (((-102.0552 4... Colorado            
    6  MULTIPOLYGON (((-73.49902 4... Connecticut         
    7  MULTIPOLYGON (((-75.80231 3... Delaware            
    8  MULTIPOLYGON (((-77.13731 3... District of Columbia
    9  MULTIPOLYGON (((-85.01548 3... Florida             
    10 MULTIPOLYGON (((-80.89018 3... Georgia             
    11 MULTIPOLYGON (((-117.0266 4... Idaho               
    12 MULTIPOLYGON (((-90.64192 4... Illinois            
    13 MULTIPOLYGON (((-88.02351 3... Indiana             
    14 MULTIPOLYGON (((-91.22634 4... Iowa                
    15 MULTIPOLYGON (((-94.63544 3... Kansas              
    16 MULTIPOLYGON (((-82.56895 3... Kentucky            
    17 MULTIPOLYGON (((-94.05103 3... Louisiana           
    18 MULTIPOLYGON (((-70.73737 4... Maine               
    19 MULTIPOLYGON (((-77.716 39.... Maryland            
    20 MULTIPOLYGON (((-70.45089 4... Massachusetts       
    21 MULTIPOLYGON (((-90.41273 4... Michigan            
    22 MULTIPOLYGON (((-96.42879 4... Minnesota           
    23 MULTIPOLYGON (((-88.18966 3... Mississippi         
    24 MULTIPOLYGON (((-95.75271 4... Missouri            
    25 MULTIPOLYGON (((-104.0491 4... Montana             
    26 MULTIPOLYGON (((-104.0606 4... Nebraska            
    27 MULTIPOLYGON (((-114.0415 3... Nevada              
    28 MULTIPOLYGON (((-72.47343 4... New Hampshire       
    29 MULTIPOLYGON (((-74.72516 4... New Jersey          
    30 MULTIPOLYGON (((-103.0063 3... New Mexico          
    31 MULTIPOLYGON (((-73.92874 4... New York            
    32 MULTIPOLYGON (((-75.89399 3... North Carolina      
    33 MULTIPOLYGON (((-104.0491 4... North Dakota        
    34 MULTIPOLYGON (((-80.51776 4... Ohio                
    35 MULTIPOLYGON (((-102.0495 3... Oklahoma            
    36 MULTIPOLYGON (((-116.9292 4... Oregon              
    37 MULTIPOLYGON (((-79.76718 4... Pennsylvania        
    38 MULTIPOLYGON (((-71.84318 4... Rhode Island        
    39 MULTIPOLYGON (((-83.10753 3... South Carolina      
    40 MULTIPOLYGON (((-96.43452 4... South Dakota        
    41 MULTIPOLYGON (((-89.71946 3... Tennessee           
    42 MULTIPOLYGON (((-94.49792 3... Texas               
    43 MULTIPOLYGON (((-114.0472 4... Utah                
    44 MULTIPOLYGON (((-73.36152 4... Vermont             
    45 MULTIPOLYGON (((-75.64188 3... Virginia            
    46 MULTIPOLYGON (((-123.0198 4... Washington          
    47 MULTIPOLYGON (((-79.49789 3... West Virginia       
    48 MULTIPOLYGON (((-87.80006 4... Wisconsin           
    49 MULTIPOLYGON (((-109.0511 4... Wyoming             
    In [4]:
    plot(st_geometry(US_flat),
         graticule=TRUE,axes=TRUE)
    
    In [5]:
    bbox_US <- st_bbox(US_flat)
    
    In [6]:
    c(xcenter = mean(bbox_US[c("xmin","xmax")]),
      ycenter = mean(bbox_US[c("ymin","ymax")]))
    
      xcenter   ycenter 
    -95.84438  37.25658 
    In [7]:
    laea <- st_crs("+proj=laea +lon_0=-95.8 +lat_0=37.3")
    US_proj <- st_transform(US_flat,laea)
    plot(st_geometry(US_proj),
         graticule=TRUE,axes=TRUE)