Geospatial Analysis with R

class: center, middle, inverse, title-slide

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

---

---
[Code for figure](https://github.com/agroimpacts/activemapper/blob/e1c238a741fd8b590e21ef1e77b89fe90a25ae4f/inst/manuscript/figures.Rmd#L454-L619)

- Optional assignment: examine the code and identify which parts of the code make each component of the figure.

---

## Today
- Homework review
- Intro to spatial, continued

---
## Homework

- Use base `R` functions to repeat the operation performed in the plot in the preceding slide, i.e. split `dat` by group, apply an `lm` to each split, and plot each split as one of three panels in a plot, with the regression line plotted over the points using `abline`.  
- If you want to go a bit further, figure out how to adorn each plot so that each series gets a different color (as with `ggplot`).

---

```r
library(sf)
library(dplyr)
sz <- 100
dat <- lapply(1:3, function(x) {
  set.seed(x)
  d <- data.frame(
    id = 1:sz,
    v1 = runif(sz, min = 2, max = 12),
    grp = paste0("g", x)
  ) %>% mutate(v2 = v1 + rnorm(sz, mean = 2, sd = 2)) %>% 
    select(id, v1, v2, grp)
}) %>% bind_rows()
```
---

```r
ggplot(dat) + 
  geom_point(aes(x = v1, y = v2)) + 
  geom_smooth(aes(x = v1, y = v2), method = "lm") + 
  facet_wrap(~grp)
```
---

## Solution

```r
# answer 1
par(mfrow = c(1, 3))
for(i in unique(dat$grp)) {
  plot(v2 ~ v1, data = dat[dat$grp == i, ])
  abline(lm(v2 ~ v1, data = dat[dat$grp == i, ]), col = "red")
}

# answer 2
layout(matrix(c(1,2,3), 1, 3, byrow = TRUE))
#code for individual plots
grouplist <- c("g1", "g2", "g3")
for(i in grouplist) {
  if(i == "g1") {
    selectedcolor <- "red"
  }
  if(i == "g2") {
    selectedcolor <- "green"
  }
  if(i == "g3") {
    selectedcolor <- "blue"
  }
  plot(x = dat$v1[dat$grp == i], y = dat$v2[dat$grp == i],
       abline(lm(dat$v2[dat$grp == i] ~ dat$v1[dat$grp == i],
                 data = dat)), col = selectedcolor)
}
```
---

Answer 2, made more compact

```r
par(mfrow = c(1, 3))
for(i in 1:3) {
  plot(v2 ~ v1, data = dat[dat$grp == unique(dat$grp)[i], ], 
       col = c("red", "green", "blue")[i], pch = 16) 
  abline(lm(v2 ~ v1, data = dat[dat$grp == unique(dat$grp)[i], ]))
}
```

---
## 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"))
farmers_s <- st_as_sf(farmers, coords = c("x", "y")) %>% 
  st_set_crs(st_crs(districts))

roads <- st_read(system.file("extdata/roads.geojson", package = "geospaar")) %>% 
  st_transform(x = roads, 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))))
```

---
## Plotting (ggplot)
- Use `geom_sf()`

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

```r
dist_areas <- districts %>% st_area()
dist_areas
dist_areas %>% units::set_units("km^2")
```

---
## Add area as column

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

```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)) %>% 
  st_transform(4326)

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