Geospatial Analysis with R

class: center, middle, inverse, title-slide

.title[
# Geospatial Analysis with R
]
.subtitle[
## Class 17
]

---

## Today
### Spatial continued
#### Data

```r
districts <- st_read(system.file("extdata/districts.geojson", 
                                 package = "geospaar"))
farmers_sf <- read_csv(
  system.file("extdata/farmer_spatial.csv", package = "geospaar")
) %>% st_as_sf(coords = c("x", "y")) %>% 
  st_set_crs(st_crs(districts))

roads <- st_read(system.file("extdata/roads.geojson", package = "geospaar")) %>%   
  st_transform(crs = st_crs(districts))

pt <- st_point(x = c(28, -14))
pts <- st_multipoint(x = cbind(x = c(27.5, 28, 28.5), y = c(-14.5, -15, -15.5)))
sline <- st_linestring(cbind(x = c(27, 27.5, 28), y = c(-15, -15.5, -16)))
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))))
```

---

## Follow-up from last class

### 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()`

---

- Mixing `%>%` and `+` in `ggplot`

- Previously

```r
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()
```

- Single pipeline

```r
ggplot(districts) + geom_sf() +
  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_minimal()
```

---

## Spatial Operations
Data setup

```r
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

```r
par(mar = rep(0, 4))
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

```r
st_intersects(x = roads_rpj, y = districts)
```

---
## `st_intersects`

```r
highway <- roads_rpj[401, ]
dists_int <- districts %>% 
  slice(st_intersects(x = highway, y = districts)[[1]])
#> 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)

```r
pol_int_dists <- st_intersection(x = pol, y = districts)

par(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)
```

---
## `st_difference`
- Similar to "erase" in ArcGIS
- Returns first object without second object.

```r
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

```r
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)

```r
## 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`

```r
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)
```

---
## Exercise
Find all districts within 50 km of sampled points?

---
## Exercise
How many farmers are within 100km of one of the "long roads"?

---
## 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

```r
# 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

```r
buffered_object %>% st_sample()
```