Spatial analysis with R

1 Libraries and data

library(tidyverse)
library(magrittr)
library(sf)

1.1 Out toy data

• Bogota: Population by UPZ
• Bogota (with SES)
• Bogota: mass transit stations
• Bogota: mass transit trunk lines
• Bogota: Bicycle parking spots

Disclaimer: These files might be wrong or outdated. Do not use them for research or serious projects.

2 Warm up

2.1 Read (import)

Zones <- st_read(dsn = "Data/Bogota_UPZ_Population.shp")

Reading layer `Bogota_UPZ_Population' from data source 
  `C:\Users\Orlan\Dropbox\Teaching\SpatialAnalysis\Tutorials\02_R_Spatial\Data\Bogota_UPZ_Population.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 113 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -74.22358 ymin: 4.435827 xmax: -74.00978 ymax: 4.83066
Geodetic CRS:  WGS 84

Zones_SES <- st_read("Data/Bogota_NSE.shp")

Reading layer `Bogota_NSE' from data source 
  `C:\Users\Orlan\Dropbox\Teaching\SpatialAnalysis\Tutorials\02_R_Spatial\Data\Bogota_NSE.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 45051 features and 13 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -74.22339 ymin: 4.466349 xmax: -74.00643 ymax: 4.832934
Geodetic CRS:  MAGNA-SIRGAS

BRT_Stations <- st_read("Data/Bogota_BRT_Stations.shp")

Reading layer `Bogota_BRT_Stations' from data source 
  `C:\Users\Orlan\Dropbox\Teaching\SpatialAnalysis\Tutorials\02_R_Spatial\Data\Bogota_BRT_Stations.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 149 features and 16 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -74.20546 ymin: 4.531715 xmax: -74.04357 ymax: 4.768818
Geodetic CRS:  WGS 84

BRT_Trunk <- st_read("Data/Bogota_BRT_TrunckLines.shp")

Reading layer `Bogota_BRT_TrunckLines' from data source 
  `C:\Users\Orlan\Dropbox\Teaching\SpatialAnalysis\Tutorials\02_R_Spatial\Data\Bogota_BRT_TrunckLines.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 13 features and 8 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: -74.20615 ymin: 4.55645 xmax: -74.04325 ymax: 4.770763
Geodetic CRS:  WGS 84

Bicycle_ParkingSpots <- st_read("Data/Bogota_BicycleParkingSpots.shp")

Reading layer `Bogota_BicycleParkingSpots' from data source 
  `C:\Users\Orlan\Dropbox\Teaching\SpatialAnalysis\Tutorials\02_R_Spatial\Data\Bogota_BicycleParkingSpots.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 255 features and 11 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -8257660 ymin: 506987.2 xmax: -8240328 ymax: 535725.3
Projected CRS: WGS 84 / Pseudo-Mercator

class(Zones)

[1] "sf"         "data.frame"

class(Zones_SES)

[1] "sf"         "data.frame"

class(BRT_Stations)

[1] "sf"         "data.frame"

class(BRT_Trunk)

[1] "sf"         "data.frame"

class(Bicycle_ParkingSpots)

[1] "sf"         "data.frame"

dim(Zones)

[1] 113   8

dim(Zones_SES)

[1] 45051    14

dim(BRT_Stations)

[1] 149  17

dim(BRT_Trunk)

[1] 13  9

dim(Bicycle_ParkingSpots)

[1] 255  12

2.2 Write (export)

# st_write(obj = BRT_Stations, dsn = "ChooseTheName.shp")

2.3 Basic mapping

The following is a mess!

plot(Zones)

If you want to see the geometry, you do this:

plot(st_geometry(Zones))

For a choroplet map:

plot(Zones["poblacion_"])

The intuitive way fails here:

plot(Zones$poblacion_)

But the tidy intuitive way does not

plot(Zones %>% select(poblacion_))

You can improve your map in many ways. Neverheless, we will use the tmap library for all the map-making and to produce images. Another alternative is to use ggplot2. You can also produce interactive maps with the leaflet library.

2.4 CSR

Planet earth is not a flat surface and is not perfectly round. Nevertheless, when you make maps you have to put a certain geographic area in two dimensions. This gets even more complicated if you consider mountains and different sea levels.

Projections are the way we transform the shape of the earth (or an area) to only two dimensions. There are many projections and all have some level of distortion. You can see the section 2.4 from the Geocomputation With R book for a clear introduction to the topic. So far, you have to know that every sf object has a Coordinate Referense System (CRS) that can be configured using EPSG codes. There are geographic projections where the notion of distance is not related to a regular understanding of distance (it is not in meter or icnhes). In Projected projections on the other hand, the distance is in meters.

In the past, many sf functions did not distinguish geographic projections from projected when running some operation sus as buffers or distances calculations. I understand this is not longer an isse, but I recommend you to doublecheck every output.

How to get the CRS?

st_crs(x = Zones)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

st_crs(BRT_Stations)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

st_crs(Zones) == st_crs(BRT_Stations)

[1] TRUE

The function st_transform() changes the CRS (re-project) of any sf object. On the other hand, the st_set_crs() function assigns (sets) a CRS, but is not projecting the object. The st_set_crs() function is useful, for example, when you manually collect data or when you have a file with missing CRS

Zones_Projected <- st_transform(Zones, 3116)

st_crs(Zones) == st_crs(Zones_Projected) 

[1] FALSE

st_is_longlat(Zones)

[1] TRUE

st_is_longlat(Zones_Projected)

[1] FALSE

3 Basic tidy operations

An sf objects works just like as data.frame or tibble. It is a table (databe) with one additional column for the geometry that is “ignored” for certain operations.

3.1 summary

summary(Zones$poblacion_)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    167   24532   50262   64671   93076  302006       1 

glimpse(Zones)

Rows: 113
Columns: 8
$ cod_loc    <dbl> 2, 3, 7, 8, 8, 8, 8, 11, 11, 11, 13, 15, 19, 2, 3, 4, 9, 10, 11, …
$ nomb_loc   <chr> "CHAPINERO", "SANTA FE", "BOSA", "KENNEDY", "KENNEDY", "KENNEDY",…
$ cod_upz    <chr> "88", "91", "86", "44", "48", "82", "113", "3", "23", "28", "107"…
$ nom_upz    <chr> "EL REFUGIO", "SAGRADO CORAZON", "EL PORVENIR", "AMERICAS", "TIMI…
$ area_urban <dbl> 335.97849, 146.18927, 461.03212, 380.96940, 430.43831, 317.32166,…
$ poblacion_ <dbl> 30763, 5879, 73629, 84584, 147298, 174145, 20993, 167, 36521, 302…
$ densidad_u <dbl> 91.562410, 40.214989, 159.704707, 222.023080, 342.204667, 548.796…
$ geometry   <MULTIPOLYGON [°]> MULTIPOLYGON (((-74.03579 4..., MULTIPOLYGON (((-74.…

3.2 Filter

Zones_A <- Zones %>% filter(poblacion_ < 64000)
dim(Zones)

[1] 113   8

dim(Zones_A)

[1] 64  8

plot(st_geometry(Zones_A))

3.3 Mutate

Zones %<>% 
  mutate(Density = poblacion_/area_urban) %>% 
  mutate(Classification = if_else(poblacion_ < 64000, "small", "large"))

3.4 group_by()

Zones %>% group_by(nomb_loc) %>% 
  summarise(Total_UPZ = n())

Simple feature collection with 19 features and 2 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: -74.22358 ymin: 4.435827 xmax: -74.00978 ymax: 4.83066
Geodetic CRS:  WGS 84

Zones %>% group_by(nomb_loc) %>% 
  summarise(Total_UPZ = n(),
            TotalPopulation = sum(poblacion_))

Simple feature collection with 19 features and 3 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: -74.22358 ymin: 4.435827 xmax: -74.00978 ymax: 4.83066
Geodetic CRS:  WGS 84

3.5 No geometry

Zones_Data <- Zones %>% st_set_geometry(NULL)

class(Zones_Data)

[1] "data.frame"

Zones_Data

Zones_Data <- Zones %>% st_drop_geometry()
Zones_Data

Zones %>% 
  st_set_geometry(NULL) %>% 
  group_by(nomb_loc) %>% 
  summarise(Total_UPZ = n(),
            TotalPopulation = sum(poblacion_))

4 Spatial Operations

4.1 Union

Bogota <- st_union(Zones)
plot(st_geometry(Bogota))

# st_combine() is not the best way to go

Bogota_Localidades <- Zones %>% 
  group_by(cod_loc) %>% 
  summarise(Total_UPZ = n(),
            TotalPopulation = sum(poblacion_))

plot(st_geometry(Bogota_Localidades))

Bogota_Localidades_Name <- Zones %>% 
  group_by(nomb_loc) %>% 
  summarise(Total_UPZ = n(),
            TotalPopulation = sum(poblacion_))

4.2 Joining

Localidad_SantaFe <- Bogota_Localidades_Name %>% filter(nomb_loc == "SANTA FE")

Localidad_Chapinero <- Bogota_Localidades_Name %>% filter(nomb_loc == "CHAPINERO")


SantaFe_Chapinero <- bind_rows(Localidad_SantaFe, Localidad_Chapinero)
# SantaFe_Chapinero <- rbind(Localidad_SantaFe, Localidad_Chapinero)

plot(st_geometry(SantaFe_Chapinero))

4.3 Centroids

Zones_Centroids <- st_centroid(Zones)

Warning: st_centroid assumes attributes are constant over geometries of x

plot(st_geometry(Zones_Centroids))

plot(st_geometry(Zones))
plot(st_geometry(Zones_Centroids), add = TRUE)

4.4 Points in polygons

plot(st_geometry(Bicycle_ParkingSpots))

plot(st_geometry(Bogota_Localidades))
plot(st_geometry(Bicycle_ParkingSpots), add = TRUE)

st_crs(Bogota_Localidades) == st_crs(Bicycle_ParkingSpots)

[1] FALSE

st_crs(Bogota_Localidades) #EPSG: 4326

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

st_crs(Bicycle_ParkingSpots) # EPSG: 3857

Coordinate Reference System:
  User input: WGS 84 / Pseudo-Mercator 
  wkt:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06°S and 85.06°N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]

Bicycle_ParkingSpots <- st_transform(Bicycle_ParkingSpots, 4326)
st_crs(Bogota_Localidades) == st_crs(Bicycle_ParkingSpots)

[1] TRUE

plot(st_geometry(Bogota_Localidades))
plot(st_geometry(Bicycle_ParkingSpots), add = TRUE)

SpotsInLocalities <- st_join(Bicycle_ParkingSpots, Bogota_Localidades)

class(SpotsInLocalities)

[1] "sf"         "data.frame"

SpotsInLocalities

Simple feature collection with 255 features and 14 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -74.17982 ymin: 4.549555 xmax: -74.02413 ymax: 4.806854
Geodetic CRS:  WGS 84
First 10 features:
   OBJECTID                            NOMBRE                            HORARIO
1        58 Universidad Central (Sede Centro) L-S 06,00 - 22:30  D 06:00 - 16:00
2        59        Tribunales - Rama Judicial                  L-V 07:00 - 20:00
3        60                       Portal Suba                  L-D 05:00 - 24:00
4        61                   Portal Américas                  L-D 05:00 - 24:00
5        62           Estación Transversal 86                  L-D 05:00 - 24:00
6        63                 Estación Banderas                  L-D 05:00 - 24:00
7        64           Estación Mundo Aventura                  L-D 05:00 - 24:00
8        65                 Estación Marsella                  L-D 05:00 - 24:00
9        66                  Estación Pradera                  L-D 05:00 - 24:00
10       67                    Portal del Sur                  L-D 05:00 - 24:00
          DIRECCION CUPOS    TIPOLOGIA  AÑO         X        Y LOCALIDAD SELLO
1        CL 21 4 40   126      Privado 2017 -74.06866 4.605609         3     1
2       KR 57 43 91    60      Privado 2017 -74.09393 4.646599        13     1
3      Troncal Suba   710 Transmilenio 2018 -74.09431 4.746746        11     1
4  Troncal Americas   785 Transmilenio 2018 -74.17294 4.629054         8     1
5  Troncal Americas    84 Transmilenio 2018 -74.15295 4.635053         8     1
6  Troncal Americas   101 Transmilenio 2018 -74.14555 4.631256         8     1
7  Troncal Americas    32 Transmilenio 2018 -74.13522 4.630140         8     1
8  Troncal Americas    32 Transmilenio 2018 -74.12860 4.629488         8     1
9  Troncal Americas    32 Transmilenio    0 -74.11969 4.628609        16     1
10  Troncal NQS Sur   409 Transmilenio    0 -74.16935 4.596727         7     1
   cod_loc Total_UPZ TotalPopulation                   geometry
1        3         5          103544 POINT (-74.06866 4.605609)
2       13         6          143891 POINT (-74.09393 4.646599)
3       11        12         1018450 POINT (-74.09431 4.746746)
4        8        12          997693 POINT (-74.17294 4.629054)
5        8        12          997693 POINT (-74.15295 4.635053)
6        8        12          997693 POINT (-74.14555 4.631256)
7        8        12          997693  POINT (-74.13522 4.63014)
8        8        12          997693  POINT (-74.1286 4.629488)
9       16         5          258368 POINT (-74.11969 4.628609)
10       7         5          554389 POINT (-74.16935 4.596727)

names(Bicycle_ParkingSpots)

 [1] "OBJECTID"  "NOMBRE"    "HORARIO"   "DIRECCION" "CUPOS"     "TIPOLOGIA"
 [7] "AÑO"       "X"         "Y"         "LOCALIDAD" "SELLO"     "geometry" 

names(SpotsInLocalities)

 [1] "OBJECTID"        "NOMBRE"          "HORARIO"         "DIRECCION"      
 [5] "CUPOS"           "TIPOLOGIA"       "AÑO"             "X"              
 [9] "Y"               "LOCALIDAD"       "SELLO"           "cod_loc"        
[13] "Total_UPZ"       "TotalPopulation" "geometry"       

SpotsInLocalities %>% 
  st_drop_geometry() %>% 
  group_by(cod_loc) %>% summarise(ParkingSpots = n())

Bogota_Localidades %<>%
  left_join(
    SpotsInLocalities %>% 
      st_drop_geometry() %>% 
      group_by(cod_loc) %>% summarise(ParkingSpots = n())
  )

Joining, by = "cod_loc"

plot(Bogota_Localidades["ParkingSpots"])

4.5 Buffers (points)

Temp <- st_centroid(SantaFe_Chapinero)

Warning: st_centroid assumes attributes are constant over geometries of x

Temp <- st_transform(Temp, 3116)

st_is_longlat(Temp)

[1] FALSE

Point1 <- Temp %>% slice(1)
Point2 <- Temp %>% slice(2)

Buffer1 <- st_buffer(Point1, 2000)  #The units are un meters
Buffer2 <- st_buffer(Point2, 2000)

plot(st_geometry(st_transform(SantaFe_Chapinero, 3116)))
plot(st_geometry(Buffer1), add = TRUE)

plot(st_geometry(Buffer2), add = TRUE)

Buffers <- bind_rows(Buffer1, Buffer2)
# st_write(Buffers, "Data/Buffers_Tests.shp")

rm(Buffers)

# st_write(Temp, "Data/Centroids_Test.shp")

You probably want to change the EPSG before exporing the files.

Take a look at the buffers in QGIS, what do you think? Here is the answer

4.6 Buffers (lines)

plot(st_geometry(BRT_Trunk))

st_is_longlat(BRT_Trunk)

[1] TRUE

BRT_Trunk_Proj <- st_transform(BRT_Trunk, 3116)

Lines_Buffer <- st_buffer(BRT_Trunk, 600)  #The units are un meters

plot(st_geometry(Lines_Buffer))

Lines_Buffer_Union_1 <- st_combine(Lines_Buffer)
plot(st_geometry(Lines_Buffer_Union_1, singleSide = FALSE))

Lines_Buffer_Union_2 <- st_union(Lines_Buffer)
plot(st_geometry(Lines_Buffer_Union_2, singleSide = FALSE))

4.7 Intersection

TwoPolygons <- Bogota_Localidades %>% 
  st_transform(3116) %>% st_centroid() %>% slice(1,2) %>% st_buffer(6000) %>% 
  st_transform(4326)

Warning: st_centroid assumes attributes are constant over geometries of x

plot(st_geometry(TwoPolygons))

Polygon1 <- TwoPolygons %>% slice(1)
Polygon2 <- TwoPolygons %>% slice(2)

Intersection <- st_intersection(Polygon1,Polygon2)

Warning: attribute variables are assumed to be spatially constant throughout all geometries

#Try st_intersects(Polygon1, Polygon2)

plot(st_geometry(Intersection))

plot(st_geometry(TwoPolygons))
plot(st_geometry(Intersection), col = "darkgreen", add = TRUE)

4.8 Difference

Difference1 <- st_difference(Polygon1, Polygon2)

Warning: attribute variables are assumed to be spatially constant throughout all geometries

plot(st_geometry(TwoPolygons))
plot(st_geometry(Difference1), col = "darkgreen", add = TRUE)

Difference2 <- st_difference(Polygon1, Polygon2)

Warning: attribute variables are assumed to be spatially constant throughout all geometries

plot(st_geometry(TwoPolygons))
plot(st_geometry(Difference2), col = "darkgreen", add = TRUE)

4.9 Distances

Points <- Bogota_Localidades %>% st_centroid()

Warning: st_centroid assumes attributes are constant over geometries of x

st_distance(Points[1,], Points[2,])

Units: [m]
        [,1]
[1,] 8633.54

st_distance(Points[1,], Points %>% slice(1:10))

Units: [m]
     [,1]    [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]     [,9]
[1,]    0 8633.54 15602.66 21003.34 27513.17 21173.11 21667.24 17458.33 13225.83
        [,10]
[1,] 9430.646

4.10 Areas (of polygons)

Bogota_Localidades

Simple feature collection with 19 features and 4 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: -74.22358 ymin: 4.435827 xmax: -74.00978 ymax: 4.83066
Geodetic CRS:  WGS 84

Bogota_Localidades %>% st_area()

Units: [m^2]
 [1] 38206434 13150137  7453510 16142036 32611697 10525532 24023801 38704812 33374184
[10] 35685857 67015948 11938006 14261459  6569395  4953664 17307013  1820002 18571827
[19] 35153482

Bogota_Localidades %>% 
  mutate(AreaTotal = st_area(.))

Simple feature collection with 19 features and 5 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: -74.22358 ymin: 4.435827 xmax: -74.00978 ymax: 4.83066
Geodetic CRS:  WGS 84

4.11 Length (of lines)

BRT_Trunk %>% st_length()

Units: [m]
 [1] 11818.012  7299.874  9727.726  9720.276 12107.235  1559.293  6593.208  2284.442
 [9]  8249.366 11498.975 10042.141  8886.344  2755.109

5 Recommended readings

• Spatial Manipulation with sf: CHEAT SHEET
• Chapters 2 to 7 in the Geocomputation with R book
• Chapters 1 to 6 in tne Spatial Data Science book
• All articles in the sf website