Today

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

• Load data

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

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

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

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

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

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)