Geospatial Analysis with R

class: center, middle, inverse, title-slide

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

---

## Today

- `terra` (raster) intro
  - agg/disagg, masking
  - assigning values
  - zonal
  - focal
  - global
  - terrain

---
- Load data

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

chirps <- terra::rast(system.file("extdata/chirps.tif", package = "geospaar"))
```

---

```r
districts <- districts %>% mutate(ID = 1:nrow(.))
distsr <- districts %>% 
  rasterize(x = ., y = raintot, field = "ID") %>% 
  print()

chirpsz <- mask(chirps, districts) 
raintot <- app(chirpsz, sum)
```
---

## Plotting issue

```r
rain_dist20 <- mask(crop(raintot, districts[20, ]), districts[20, ])

par(mar = rep(0, 4))
plot(districts$geometry)
plot(districts %>% slice(20) %>% st_geometry, add = TRUE, border = "red", 
     lwd = 4)
plot(rain_dist20, add = TRUE)

par(mar = rep(0, 4))
# plot(vect(districts), axes = FALSE)  # solution 1
plot(districts$geometry, axes = FALSE)  # solution 2
# plot(districts %>% slice(20) %>% vect(.), # solution 1
#      add = TRUE, border = "red",
#      lwd = 4)
plot(districts %>% slice(20), add = TRUE, border = "red", lwd = 4) # solution 2
plot(rain_dist20, add = TRUE, ext = districts)

par(mar = rep(0, 4))
plot(districts$geometry, axes = FALSE)  # solution 2
plot(districts %>% slice(20), add = TRUE, border = "red", lwd = 4) # solution 2
plot(rain_dist20, add = TRUE, ext = districts)
```

---
## Mosaic

```r
districts %>% 
  st_union() %>% 
  st_centroid() %>% 
  st_buffer(dist = 50000) %>% 
  st_intersects(., districts) %>% 
  unlist(.) %>% 
  slice(districts, .) -> districts_ss

plot(districts$geometry)
plot(districts_ss$geometry, add = TRUE, col = "blue")
districts %>% 
  st_union() %>% 
  st_centroid() %>% 
  st_buffer(dist = 50000) %>% plot(add = TRUE, border = "red")

rainlist <- lapply(districts_ss$geometry, function(x) {
  terra::crop(raintot, x)
})

# par(mfrow = c(2, 2))
par(mar = rep(0, 4))
plot(vect(districts))
for(i in rainlist) plot(i, add = TRUE, legend = FALSE)

plot(mosaic(sprc(rainlist)))
plot(districts_ss$geometry, add = TRUE)
```

---
## Climate data summaries and zonal stats

```r
tmax <- geodata::worldclim_country("Zambia", var = "tmax", path = tempdir())
dt <- as.Date(gsub("Y", "", names(chirpsz)), format = "%y%j")
dt <- lubridate::parse_date_time(gsub("Y", "", names(chirpsz)), orders = "yj")
wk <- lubridate::week(dt)
weekly_rainfall <- lapply(unique(wk), function(x) {
  which(wk == x)
  app(chirpsz[[which(wk == x)]], sum)
}) %>% do.call(c, .)
names(weekly_rainfall) <- paste0("wk", unique(wk))

rflim <- range(global(weekly_rainfall, range, na.rm = TRUE))
plot(weekly_rainfall, zlim = rflim)
```

---
## Distance
- plot distance from points
- points must be in same projection as model raster

```r
set.seed(1)
randsamp <- spatSample(raintot, size = 10, xy = TRUE, na.rm = TRUE) %>%
  st_as_sf(coords = c("x", "y"), crs = 4326)

ptdistr <- distance(raintot, vect(randsamp))
#> 
|---------|---------|---------|---------|
=========================================
                                          
ptdistrmsk <- mask(ptdistr, raintot)
plot(ptdistrmsk)
```