Geospatial Analysis with R

class: center, middle, inverse, title-slide

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

---

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

---
## Add area as column

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

---
## Plot districts (exercise)
- Use `geom_sf()` to plot area
- Fill districts based on area column.
- Use `scale_fill_continuous` , with `type = 'viridis'`

```r
ggplot(districts) + geom_sf(aes(fill = area)) + 
  scale_fill_continuous(type = "viridis")
```

---
## Compare areas

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

ggplot(districts) + 
  geom_histogram(aes(x = area - area2))

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

# calculate (area - area2) as a percent of area, and assignn it to a new variable areadiff. Plot areadiff against area in a scatterplot using ggplot
```
---

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

---

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

---

```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_viridis_c() + 
  theme_minimal()
```

---
## ggplot steps
- Start with `ggplot()`. You can include data if you're only using one data set.
- Add plot type (`geom_point`, `geom_line`, `geom_sf`, `geom_histogram`)
- Add aesthetics in `aes()`. These are columns used for plotting. 
- For `geom_point`, `geom_line` you need 'x' and 'y' aesthetics. 
- For `geom_sf`, spatial context used for 'x' and 'y'. But you might want 'fill' or 'color'