Exercises for Chapter 11

In these exercises, we use again the database spatialdata with the two datasets we have used in the chapter. You need to have the PostGIS extension enabled for this database, as explained in the chapter. You can read these datasets and create the corresponding tables using the code from the chapter:

library(RPostgres)
library(tidyverse)
library(sf)

db <- dbConnect(Postgres(),
  dbname = "spatialdata",
  user = "postgres",
  password = "pgpasswd"
)

# dbExecute(db, "DROP EXTENSION postgis")
dbExecute(db, "CREATE EXTENSION postgis")

events <- read.csv(file.path("ch11", "ged.csv"))
dbWriteTable(db, "events", events)
dbExecute(db, "ALTER TABLE events ADD COLUMN geom geometry(point, 4326)")
dbExecute(db, "UPDATE events SET geom = st_setSRID(st_point(longitude, latitude), 4326)")
events <- st_as_sf(events, coords = c("longitude", "latitude"), crs = 4326)

municipalities <- st_read(file.path("ch11", "bosnia.shp"), crs = 4326) %>%
  rename(geom = geometry)
st_write(municipalities, dsn = db, layer = "municipalities")

Exercise 1: Retrieving the Location of the First Event in R/sf

Where did violence start in Bosnia? In this exercise, we use our two datasets to determine in which municipality the first event recorded by the UCDP Geo-referenced Event Dataset occurred. How can you use the joined dataset we created in the chapter for this purpose?

Solution

We need the dataset where events were joined with municipalities based on their location, as demonstrated in the chapter:

joined <- st_join(events, municipalities)

Sorting this dataset by event date puts the earliest event and the corresponding municipality at the top:

joined %>%
  arrange(date_start) %>%
  select(date_start, name) %>%
  slice(1:3)

## Simple feature collection with 3 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 18.31778 ymin: 43.505 xmax: 18.77861 ymax: 45.09
## Geodetic CRS:  WGS 84
##   date_start          name                geometry
## 1 1992-04-27       Srbinje POINT (18.77861 43.505)
## 2 1992-04-27 Novo Sarajevo  POINT (18.38333 43.85)
## 3 1992-04-28         Odzak  POINT (18.31778 45.09)

Exercise 2: Retrieving the Location of the First Event in PostGIS

This exercise repeats the task from the previous one, but using PostGIS. Can you solve the task with a single SQL statement?

Solution

In PostgreSQL/PostGIS, all we need to do is amend the spatial join from the chapter, such that it orders the joined dataset by event date:

dbGetQuery(db, "SELECT events.date_start, municipalities.name FROM municipalities JOIN events ON st_contains(municipalities.geom, events.geom) ORDER BY date_start LIMIT 3")

##   date_start          name
## 1 1992-04-27       Srbinje
## 2 1992-04-27 Novo Sarajevo
## 3 1992-04-28         Odzak

Exercise 3: Distance Calculation in R/sf

In this exercise, we calculate distances between municipalities. Distances are often used in spatial analysis to model the declining influence of one unit over another - for example, it is typically assumed that state strength declines the further we get away from a country’s capital. To incorporate this into our analysis, we need to amend our municipalities data frame with a new column that contains the distances from the capital of Bosnia, Sarajevo. The sf package has a useful function for this purpose: st_distance(). Using the package documentation, familiarize yourself with this function. What happens if you simply apply this function to the municipalities dataset? To compute the distance from Sarajevo, we use the “Centar” district in our dataset, which is the center district of Sarajevo. From the output of the function call, how can you extract a column with the distances from this district, and append it to the `municipalities’ dataset?

Solution

The st_distance() function can be applied to objects of type sf. By default, it takes two datasets, and computes a matrix with the distances between each pair of objects from the two lists. If you only provide one dataset, it combines the dataset with itself, and computes distances between each pair of objects from this dataset (output is not displayed here, since it is a huge matrix):

st_distance(municipalities)

In this matrix, cell i,j contains the distance between the municipality with index i and the one with index j. Therefore, all entries i,i along the diagonale are zero, and those between adjacent units too. The distances are provided in meters. All we need to do is find out which index the “Centar” districts has

municipalities[municipalities$name == "Centar", ]

## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 18.38905 ymin: 43.84151 xmax: 18.44457 ymax: 43.94167
## Geodetic CRS:  WGS 84
##     id   name                           geom
## 107 65 Centar POLYGON ((18.38905 43.89011...

and use the corresponding column (107) from the distance matrix for the new variable:

municipalities$dist_capital <- st_distance(municipalities)[, 107]

Exercise 4: Distance Calculation in PostGIS

Here, we repeat the distance calculation from Exercise 3 in PostGIS. First, add a new numeric field to the municipalities table in the database. Next, take a look at PostGIS’s distance calculation functions. The standard one, st_distance(), calculates distances in the coordinate system of the respective dataset. Since we have spherical coordinates (latitude, longitude), this would give us incorrect results, which is why we use the st_distanceSphere() function. Similar to its sf counterpart, the function takes as arguments two datasets and computes distances between each pair of features from them. Since we only need distances from Sarajevo, we first need to create a temporary dataset with only one entry: the “Centar” district. For this, it is convenient to use a subquery. In PostgreSQL, you can define these subqueries with the WITH keyword. Can you figure out how to do this from the PostgreSQL documentation? Your WITH query should be a simple SELECT statement that extracts the “Centar” municipality, and it should be used in an UPDATE statement that fills the new column with the distance between each municipality and the “Centar” municipality.

Solution

We first create a new column

dbExecute(db, "ALTER TABLE municipalities ADD COLUMN dist_capital double precision")

and then update this column with the distances to the “Centar” municipality:

dbExecute(db, "WITH sarajevo AS (SELECT * FROM municipalities WHERE name = 'Centar') UPDATE municipalities SET dist_capital = st_distanceSphere(sarajevo.geom, municipalities.geom) FROM sarajevo")

We can now test if the two distance calculations match:

municipalities[municipalities$name == "Neum", ]

## Simple feature collection with 1 feature and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 17.57853 ymin: 42.88908 xmax: 17.95101 ymax: 43.03375
## Geodetic CRS:  WGS 84
##   id name                           geom dist_capital
## 2 17 Neum POLYGON ((17.83138 43.01822... 103307.7 [m]

dbGetQuery(db, "SELECT dist_capital FROM municipalities WHERE name = 'Neum'")

##   dist_capital
## 1     103307.7