Geospatial Analysis with R

class: center, middle, inverse, title-slide

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

---

---

## Today

- `terra` (raster) intro
  - classes
  - dummy rasters
  - reading and writing
  - rasterize, raster to vector
  - agg/disagg, masking
  - assigning values
  - zonal
  - focal
  - global

---
- 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"))
mypoly <- st_polygon(list(cbind(x = c(26.5, 27.5, 27, 26, 26.5), 
                                y = c(-15.5, -16.5, -17, -16, -15.5)))) %>%
  st_sfc(crs = 4326) %>% st_intersection(., districts) %>% st_as_sf() %>% 
  mutate(ID = seq(1:nrow(.)))
#st_crs(mypoly) <- 4326
smallest_centroid <- districts %>% 
  filter(distName == "Lusaka") %>% 
  st_centroid()
```

## terra objects
- `SpatRaster` object - 1 layer to many layers, from multiple or single file sources
- `SpatVector` object - Vector data

---
## accessing bands
- use list `[[ ... ]]` notation

```r
chirps <- rast(system.file("extdata/chirps.tif", package = "geospaar"))
chirps[[3]] ## access by index
#chirps[["Y16301"]] ## access by band name
```

```r
plot(chirps[[3]])
```

---

## Creating raster from scratch

- how does `extent` work?

```r
library(terra)
r <- rast(ext(30, 31, -14, -13), resolution = 0.1, crs = "EPSG:4326")
values(r) <- sample(1:10, size = ncell(r), replace = TRUE)
r2 <- r
# values(r2) <- rnorm(n = 10)
values(r2) <- rnorm(n = ncell(r2))
# values(r) <- 1:15
# plot(r)

par(mar = c(0, 0, 0, 0))
plot(districts %>% st_geometry)
# plot(r, add = TRUE)
# plot(r2, add = TRUE)
# plot(r - r2, add = TRUE)

plot(r - r2)
```

---
- Creating multi-layer stack

```r
s <- c(r, log10(r))
names(s) <- c("dummy", "log10dummy")
plot(s)
print(class(s))
```

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

---
### A large brick
- how does this differ from above?

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

plot(b2)
plot(b2[[1]])
plot(b2[[10]])
plot(b2[[c(1, 10)]])
plot(b2[["b3"]])
plot(b2["b3"])

```

---
## read and write

```r
## write raster
writeRaster(r, filename = file.path(tempdir(), "mydummy.tif"), overwrite = TRUE)
writeRaster(r, filename = file.path(tempdir(), "mydummy2.tif"), 
            overwrite = TRUE, datatype = "INT1U")
dir(tempdir())

## read in raster
r <- rast(file.path(tempdir(), "mydummy.tif"))
plot(r)
```

---
## read and write

- `brick` function can be used both to read and write

```r
## write file
writeRaster(s, filename = file.path(tempdir(), "mybrick.tif"), overwrite = TRUE)

## read file
newbrick <- rast(file.path(tempdir(), "mybrick.tif"))
```

---

## Projecting
- `project` 
- similar to `st_transform()`, use crs of existing layer to project
- need to define resolution with `res` argument

```r
roads <- read_sf(system.file("extdata/roads.geojson", package = "geospaar"))
chirpsz_alb <- project(x = chirps, y = crs(roads), res = 5000) 
chirpsz_alb %>% print()
plot(chirpsz_alb[[1:10]])
```

---
## Vector <--> Raster
- `rasterize` converts to raster. Need to pass argument that gives raster dimensions.
- In below example, we use `chirpsz` for dimensions

```r
distsr <- rasterize(districts, chirps, field = "distName")
distsr %>% plot_noaxes
str(distsr)
values(distsr)
ext(distsr)
# stack(distsr, chirpsz[[1:2]]) %>% plotRGB(stretch = "lin")

```

---
## Vector <--> Raster
- `as.polygons`, use dissolve to merge common raster values

```r
distsr_pol <- as.polygons(distsr)
# plot(as.polygons(distsr))
# plot(distsr_pol)
distsr_pol %>% 
  st_as_sf %>% 
  st_geometry %>% 
  plot
distsr_pol %>% 
  st_as_sf %>% 
  slice(49) %>% 
  plot(add = TRUE)
```

---
## Plotting
- run examples below individually

```r
plot(b)
# plot(b2)
plot_noaxes(b)
plot_noaxes(b, main = paste("Random", 1:10))
plot(b2, main = paste("Random", 1:10))
plot_noaxes(b, legend = FALSE)
```

---

```r
# png("external/mydem.png", width = 6, height = 4, res = 300, units = "in")
# plot_noaxes(b[[1]], legend = FALSE, mar = c(0, 0, 0, 6), oma = c(0, 0, 0, 4))
plot_noaxes(b[[1]], legend = "bottom", mar = c(1, 1, 1, 6), oma = c(0, 0, 0, 4))
# add_legend(x = 30.5, y = -13.5, legend = 1:10, fill = terrain.colors(10), 
#            border = FALSE, bty = "n")
# dev.off()
```

---

```r
# plotting with larger datasets
plot_noaxes(b2, main = paste("Random", 1:10))
plot_noaxes(b2, nc = 2, nr = 5, main = paste("Random", 1:10))
```

---
- ggplot() works with `geom_raster()`

```r
ggplot2::ggplot(as.data.frame(b[[1]], xy = TRUE)) + 
  geom_raster(aes(x = x, y = y, fill = b1)) + scale_fill_viridis_c()
```

---
- stars package

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

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

---

```r
quantile(values(chirps[[1]]))
#>        0%       25%       50%       75%      100% 
#>  0.000000  0.000000  0.000000  3.119399 21.333221
plot(terra::classify(chirps[[1]], rcl = c(0, 5, 10, 25)))
```

---
## Raster to other types

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

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

```r
plot(b, main = paste0("Brick ", 1:2) )
aggregate(b, fact = 2) %>% plot(., main = paste0("Brick agg", 1:2) )
# 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:2) )
disagg(b, fact = 2) %>% 
  plot(., main = paste0("Brick disagg", 1:2) )
disagg(b, fact = 2, method = 'bilinear') %>% 
  plot(., main = paste0("Brick disagg bilinear", 1:2) )
```

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

```r

# dev.off()
plot(rf1 / 1000)
plot(((log10(rf1 + 1) * 10) / 5)^2)

```