Geospatial Analysis with 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)

library(geospaar)
chirps <- rast(
  system.file("extdata/chirps.tif", package = "geospaar")
) 
plot_noaxes(chirps[[1]])
chirpsz <- mask(chirps, districts)
plot_noaxes(chirpsz[[1]])
plot(st_geometry(districts), add = TRUE)
# rasterVis::levelplot(chirpsz[[1:5]])

rf1 <- chirpsz[[1]] + chirpsz[[2]]
rf2 <- chirpsz[[1]] * 5

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

districts <- st_read(
  system.file("extdata/districts.geojson", package = "geospaar")
) 
districts <- districts %>% 
  mutate(ID = 1:nrow(.))
distsr <- rasterize(x = districts, y = raintot, field = "ID") 

# 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])

# 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])
}

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

set.seed(1)
randsamp <- spatSample(rast(raintotalb), size = 10, xy = TRUE) %>%
  as_tibble %>%
  st_as_sf(coords = c("x", "y")) %>% 
  drop_na()
st_crs(randsamp) <- crs(raintotalb)
# #2
ptdistr <- terra::distance(rast(raintotalb), vect(randsamp))
ptdistrmsk <- mask(ptdistr, rast(raintotalb))
plot(ptdistrmsk)