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

---
- Creating multi-layer brick
- what are the dimensions of this raster? resolution?
- how are the values filled?

```r
b <- lapply(1:10, function(x) {
  r <- rast(ext(30, 31, -14, -13), resolution = 0.1, 
            crs = "EPSG:4326")
  set.seed(x)
  values(r) <- sample(1:10, size = ncell(r), replace = TRUE)
  r
}) %>% do.call(c, .)
names(b) <- paste0("b", 1:10)

plot(b)
```

---
- stars package

```r
library(stars)
stars::st_as_stars(b) %>% plot(col = viridis::viridis(10))
stars::st_as_stars(b) %>% plot()

ggplot() + 
  geom_stars(data = stars::st_as_stars(b[[1]]))
```

---
## Raster to other types

```r
plot(b)
b[1:10]
b[1:ncell(b)]
as.matrix(b)
as.data.frame(b) %>% as_tibble()
```

---

## Classification

```r
plot(chirps)
chirpsz <- mask(chirps, districts)
plot(chirpsz[[1:4]])
raintot <- app(chirpsz, sum)
plot(raintot)

qtiles <- unlist(global(raintot, quantile, na.rm = TRUE))
plot(terra::classify(raintot, rcl = qtiles))

rclmat <- rbind(c(10, 37, 1), 
                c(37, 48, 2), 
                c(48, 60, 3), 
                c(60, 160, 4))
rcldat <- data.frame(rclmat)
rcldat$X3 <- letters[1:4]

rclmat <- cbind(qtiles[-length(qtiles)], 
                qtiles[-1], 1:(length(qtiles) - 1))

plot(terra::classify(raintot, rcl = rclmat, include.lowest = TRUE))
# plot(terra::classify(raintot, rcl = rcldat))
```

---

```r
r1 <- raintot < 50
r2 <- (raintot > 50 & raintot < 100) * 2
r3 <- (raintot > 100) * 3
plot(app(c(r1, r2, r3), sum))
```

---
## Pre-processing
- Aggregating/disaggregating
- what is default method for aggregating?

```r
plot(b, main = paste0("Brick ", 1:10))
aggregate(b, fact = 2) %>% plot(., main = paste0("Brick agg", 1:10))
# aggregate(b, fun = min, fact = 2) %>% plot
# aggregate(b, fun = sd, fact = 2) %>% plot
```

---

## Pre-processing
- what is default method for disaggregating?

```r
plot(b, main = paste0("Brick ", 1:10))
disagg(b, fact = 2) %>% 
  plot(., main = paste0("Brick disagg", 1:10))
disagg(b, fact = 2, method = 'bilinear') %>% 
  plot(., main = paste0("Brick disagg bilinear", 1:10))
```

---
- Masking
- `mask` makes values `NA` outside of mask area.

```r
plot_noaxes(chirps[[1]])
chirpsz <- mask(chirps, districts)
plot_noaxes(chirpsz[[1]])
plot(st_geometry(districts), add = TRUE)
# rasterVis::levelplot(chirpsz[[1:5]])
```

---

## Calculations

- raster algebra
- statistics
- z dimension stats

---
## Calculations
- math operations work cell-wise

```r
rf1 <- chirpsz[[1]] + chirpsz[[2]]
rf2 <- chirpsz[[1]] * 5
plot(c(rf1, rf2), main = paste0("r", 1:2))
```

```r
# dev.off()
plot(rf1 / 1000)
plot(((log10(rf1 + 1) * 10) / 5)^2)
plot(((log10(chirps[[1:5]] + 1) * 10) / 5)^2)
```

---
## Terrain analysis with zonal statistics

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

# Download DEM
# dem <- getData(name = "alt", country = "ZMB", path = tempdir())
dem <- geodata::elevation_30s(country = "Zambia",  path = tempdir())

# calculate slope
slope <- terrain(x = dem, v = 'slope', unit = 'degrees')
plot(slope)

# calculate mean by district
distsr_rs <- resample(x = distsr, y = slope, method = "near")
plot(distsr_rs)
zoneslope <- zonal(x = slope, z = distsr_rs, fun = "mean", na.rm = TRUE)
hist(zoneslope[, 2])
zoneelevation <- zonal(x = dem, z = distsr_rs, fun = "mean", na.rm = TRUE)
hist(zoneelevation[, 2])
```
---

## Subst

```r
# map zonal statistics
distr_slopes <- subst(x = distsr_rs, from = zoneslope$ID, 
                      to = zoneslope$slope)
distr_elevzone <- subst(x = distsr_rs, from = zoneelevation$ID, 
                        to = zoneelevation$ZMB_elv_msk)

# plot
plot(s[[2]])
s <- c(distr_elevzone, distr_slopes)
titles <- c("Elevation", "Slope")
par(mfrow = c(1, 2))
for(i in 1:length(s)) {
  plot_noaxes(s[[i]], main = titles[i])
}
```

---
## Climate data summaries and zonal stats

```r
# rainfall <- getData(name = "worldclim", res = 10, 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)
#   calc(chirpsz[[which(wk == x)]], sum)
# }) %>% stack(.)
# names(weekly_rainfall) <- paste0("wk", unique(wk))
# 
# rflim <- range(cellStats(weekly_rainfall, range))
# plot(stack(weekly_rainfall), zlim = rflim)
```

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

```r
# set.seed(1)
# randsamp <- sampleRandom(raintotalb, size = 10, xy = TRUE) %>% 
#   as_tibble %>%  
#   st_as_sf(coords = c("x", "y")) 
# st_crs(randsamp) <- crs(raintotalb)
# # #2
# ptdistr <- distanceFromPoints(object = raintotalb, xy = as_Spatial(randsamp))
# ptdistrmsk <- mask(ptdistr, raintotalb)
# plot(ptdistrmsk)
```