Spatial continued

Data

districts <- st_read(system.file("extdata/districts.geojson", 
                                 package = "geospaar"))
Reading layer `districts' from data source 
  `/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/geospaar/extdata/districts.geojson' 
  using driver `GeoJSON'
Simple feature collection with 72 features and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 21.99978 ymin: -18.0751 xmax: 33.6875 ymax: -8.226213
Geodetic CRS:  WGS 84
farmers <- read_csv(system.file("extdata/farmer_spatial.csv", 
                                package = "geospaar"))

Projections, cont.

Setting

farmers_sf <- st_as_sf(farmers, coords = c("x", "y")) %>% 
  st_set_crs(st_crs(districts))

Transform

roads <- st_read(system.file("extdata/roads.geojson", package = "geospaar")) %>% 
  st_transform(crs = st_crs(districts))
Reading layer `roads' from data source 
  `/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/geospaar/extdata/roads.geojson' 
  using driver `GeoJSON'
Simple feature collection with 473 features and 1 field
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: -292996.9 ymin: -2108848 xmax: 873238.8 ymax: -1008785
Projected CRS: Africa_Albers_Equal_Area_Conic

Features from scratch

pt <- st_point(x = c(28, -14))
pt
#> POINT (28 -14)
#
# #2
pts <- st_multipoint(x = cbind(x = c(27.5, 28, 28.5), y = c(-14.5, -15, -15.5)))
pts
#> MULTIPOINT ((27.5 -14.5), (28 -15), (28.5 -15.5))
#
# #3
sline <- st_linestring(cbind(x = c(27, 27.5, 28), y = c(-15, -15.5, -16)))
sline
#> LINESTRING (27 -15, 27.5 -15.5, 28 -16)
#
# #4
pol <- st_polygon(list(cbind(x = c(26.5, 27.5, 27, 26, 26.5), 
                             y = c(-15.5, -16.5, -17, -16, -15.5))))
pol

Plotting

  • sf’s generic plot() (use add = TRUE)
par(mar = c(0, 0, 0, 0))
plot(districts %>% st_geometry(), col = "grey")
plot(pt, add = TRUE, col = "red", pch = 20)
plot(pts, add = TRUE, col = "blue", pch = 20)
plot(sline, add = TRUE, col = "cyan", lwd = 2)
plot(pol, add = TRUE, col = "orange", lwd = 2)

Plotting (ggplot)

  • Use geom_sf()
ggplot(districts) + geom_sf() +
  geom_sf(data = roads, col = "red") + 
  geom_sf(data = farmers_sf, col = "blue")

Area and Length

  • st_area(), st_length. Default units are in meters
  • How many square m in a square km?
dist_areas <- districts %>% st_area()
dist_areas
dist_areas %>% units::set_units("km^2")

Add area as column

districts <- districts %>% 
  mutate(area = st_area(.) %>% units::set_units("km^2") %>% as.numeric()) 
districts

Plot districts (exercise)

  • Use geom_sf() to plot area
  • Fill districts based on area column.
  • Use scale_fill_continuous , with type = 'viridis'
ggplot(districts) + geom_sf(aes(fill = area)) + 
  scale_fill_continuous(type = "viridis")

Compare areas

districts <- districts %>% 
  st_transform("ESRI:102022") %>% 
  mutate(area2 = as.numeric(st_area(.) / 10000)) %>% 
  st_transform(4326)

ggplot(districts) + 
  geom_histogram(aes(x = area - area2))
districts <- districts %>%
  mutate(areadiff = (area - area2) / area * 100)

ggplot(districts) +
  geom_point(aes(x = area2, y = areadiff))

districts %>%
  mutate(areadiff = ((area - area2) / area) * 100) %>%
  ggplot() + geom_point(aes(x = area2, y = areadiff))

districts %>%
  mutate(areadiff = ((area - area2) / area) * 100) %>% 
  summarize(mu = mean(areadiff))
districts %>%
  mutate(areadiff = ((area - area2) / area) * 100) %>%
  ggplot() + geom_sf(aes(fill = areadiff))

districts %>%
  mutate(areadiff = ((area - area2) / area) * 100) %>% 
  summarize(mu = mean(areadiff)) %>% 
  ggplot() + geom_sf(aes(fill = mu))

districts %>%
  mutate(areadiff = ((area - area2) / area) * 100) %>% 
  st_drop_geometry() %>% 
  summarize(mu = mean(areadiff)) #%>% 
  # ggplot() + geom_sf(aes(fill = mu))

districts %>%
  mutate(areadiff = ((area - area2) / area) * 100) %>% 
  as_tibble() %>% 
  summarize(mu = mean(areadiff)) #%>%

Calculate length (exercise)

  • Use st_length() to add column ‘road_length’ to roads
  • ‘road_length’ should be in km
  • Use ggplot to plot districts with black outline (color), and empty fill (fill = NA)
  • Then plot roads based on length (use scale_color_continuous() )
  • Plot roads with line width 1 (lwd = 1)
  • Use theme_bw()
roads <- st_read(
  system.file("extdata/roads.geojson", package = "geospaar")
) 
roads <- roads %>% 
  mutate(road_length = as.numeric(units::set_units(st_length(.), "km"))) 
ggplot() +
  geom_sf(data = districts, color = "black", fill = NA) +
  geom_sf(data = roads, aes(color = road_length), lwd = 1) + 
  scale_color_continuous(type = "viridis") + 
  theme_bw()
  • Mixing %>% and + in ggplot
  • Single pipeline
ggplot(districts) + 
  geom_sf(color = "black", fill = NA) +
  st_read(system.file("extdata/roads.geojson", package = "geospaar")) %>% 
  mutate(road_length = as.numeric(units::set_units(st_length(.), "km"))) %>% 
  geom_sf(data = ., aes(color = road_length), lwd = 1) + 
  scale_color_continuous(type = "viridis") +
  theme_bw()
  • Calculate length of roads in km
  • select roads between 50 and 500 km length
  • plot the selected roads over districts using ggplot. Make the roads red
p <- ggplot(districts) + geom_sf()
p + roads %>% 
  mutate(road_length = as.numeric(units::set_units(st_length(.), "km"))) %>% 
  filter(., road_length <= 500 & road_length > 49) %>% 
  geom_sf(data = ., color = 'red') + theme_minimal()

roads_ss <- roads %>% 
  mutate(road_length = as.numeric(units::set_units(st_length(.), "km"))) %>% 
  filter(., road_length <= 500 & road_length > 49)

p + geom_sf(data = roads_ss, color = 'red') + 
  theme_minimal()

Spatial Operations

Data setup

coords <- cbind("x" = c(25, 25.1, 26.4, 26.4, 26.4, 25), 
                "y" = c(-15.5, -14.5, -14.6, -14.8, -15.6, -15.5))
pol <- st_polygon(x = list(coords)) %>% 
  st_sfc %>% 
  st_sf(ID = 1, crs = 4326)
roads_rpj <- st_transform(roads, crs = st_crs(districts))

Plotting

par(mar = rep(0, 4))
# plot(districts)
# plot(districts$geometry, col = "grey")
plot(st_geometry(districts), col = "grey")
plot(pol, col = "purple", add = TRUE)
plot(roads_rpj, col = 'red', add = TRUE)

st_intersects

  • Similar to “Select by Spatial Location”
  • Finds indices of second object that intersect each element of first object
st_intersects(x = roads_rpj, y = districts)

st_intersects

highway <- roads_rpj[401, ]
dists_int <- districts %>% 
  slice(st_intersects(x = highway, y = .)[[1]])
# unlist(st_intersects(x = highway, y = districts))
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
dists_int
plot(st_geometry(dists_int), col = 'grey')
plot(highway, col = 'red', add = TRUE)

st_intersection()

  • finds intersecting pairs between two objects
  • Each row in result is one “intersection object” (i.e. an object from x and y that intersect)
pol_int_dists <- st_intersection(x = pol, y = districts)

par(mfrow = c(1, 2), mar = rep(0, 4))
plot(st_geometry(districts), col = "grey")
plot(st_geometry(pol_int_dists), col = rainbow(n = nrow(pol_int_dists)), 
     add = TRUE)
plot(st_geometry(districts), col = "grey")
districts %>% 
  slice(unlist(st_intersects(x = pol, y = .))) %>% 
  plot(add = TRUE)

st_difference

  • Similar to “erase” in ArcGIS
  • Returns first object without second object.
d1 <- st_difference(x = districts, y = pol)
d2 <- st_difference(x = pol, y = districts)

par(mfrow = c(1, 2), mar = c(0, 0, 1, 0))
d1 %>% 
  st_geometry %>% 
  plot(col = "grey")
d2 %>% 
  st_geometry %>%
  plot(col = "grey")

st_union

  • Combines objects without resolving borders.
  • Really takes all possible combinations of first and 2nd object
  • Switching order of inputs changes column order
uni1 <- st_union(x = pol, y = districts)
uni1

uni2 <- st_union(x = districts, y = pol)
uni2

par(mfrow = c(1, 2), mar = c(0, 0, 1, 0))
plot(uni1[2], reset = FALSE, main = "Union of pol with districts")
plot(uni2[1], reset = FALSE, main = "Union of districts with pol")

st_buffer()

  • Best to use projected coordinates (e.g. Albers)
## project everything to Albers
farmers_alb <- farmers_sf %>% 
  st_transform(st_crs(roads))
long_roads <- roads %>% 
  filter(as.numeric(st_length(.)) > 500000)
districts_alb <- districts %>% 
  st_transform(st_crs(roads))

roads_buff_20km <- st_buffer(long_roads, 20000)
plot(st_geometry(districts_alb), col = "grey")
plot(st_geometry(roads_buff_20km), col = "blue", add = TRUE)
plot(st_geometry(long_roads), col = "red", add = TRUE)

st_sample

par(mar = rep(0, 4))
districts %>% 
  st_union %>% 
  plot(col = "grey")
set.seed(1)
districts %>% 
  st_union %>% 
  st_sample(size = 25, exact = TRUE) %>% 
  plot(pch = 20, add = TRUE)
par(mar = rep(0, 4))
districts %>% 
  st_geometry %>% 
  plot(col = "grey")
set.seed(1)
districts %>% 
  st_sample(size = rep(5, nrow(.)), by_polygon = TRUE) %>% 
  plot(pch = 20, add = TRUE)

Exercise

Place 3 random points in Zambia (union districts to make single shape of Zambia). Find all districts within 50 km of sampled points?

Adapt this

set.seed(1)
pts <- st_union() %>% 
  st_sample() %>% 
  st_transform(st_crs())
dists_in_pts <- districts_alb %>% 
  slice(unlist(st_intersects(st_buffer(, dist = 50000), .)))
ggplot() + 
  geom_sf(data = ) + 
  geom_sf(data = , fill = "blue") + 
  geom_sf(data = )

Exercise

How many farmers are within 100km of one of the “long roads”?

Adapt

farmers_alb %>% 
  slice(unlist(st_intersects(st_buffer(), .))) %>% 
  ?

Further Exercises

Round 1

  • Calculate the length of roads in kilometers
  • Buffer roads > 100 km by 30 kilometers, save as new object
  • Buffer roads > 100 km 10 kilometers, save as new object
  • Take the difference between the 30 km and 10 km buffers, i.e. erase/difference the two
  • Plot the resulting difference using ggplot with the difference in red over Zambia’s districts
# code to adapt
roads %>% 
  mutate(km = as.numeric(st_length(.)) / 1000) %>% 
  filter(km > 200) %>% 
  st_buffer(dist = 20000) %>% 
  st_transform(crs = 4326)
st_difference(larger_object, smaller_object)
ggplot() + geom_sf

Round 2

  • Randomly sample 100 points within the 30 km road buffer (use a seed of 1)
  • Plot the 30 km buffer (red) the points in it (blue) over the districts
buffered_object %>% st_sample()