Geospatial Analysis with 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"))
chirpsz <- mask(chirps, districts)

raintot <- app(chirpsz, sum)

plot_noaxes(chirpsz[[1]])
plot(st_geometry(districts), add = TRUE)

library(dplyr)
library(sf)
library(terra)
wcmeantemp <- geodata::worldclim_country(var = "tavg", res = 2.5,
                                         country = "Zambia", path = tempdir())
districts <- st_read(
  system.file("extdata/districts.geojson", package = "geospaar"),
)
zamtmean <- mask(app(wcmeantemp, mean), districts)
trng <- global(zamtmean, range, na.rm = TRUE)
reclmat <- cbind(c(floor(trng[1]), 20, 24), c(20, 24, ceiling(trng[2])), 1:3)

zamtclass <- classify(x = zamtmean, rcl = as.integer(reclmat), 
                      include.lowest = TRUE)

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

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)

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)
districts %>% 
  slice(20) %>% 
  st_geometry() %>% 
  plot(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
districts %>% 
  slice(20) %>% 
  st_geometry() %>% 
  plot(add = TRUE, border = "red", lwd = 4) # solution 2
plot(rain_dist20, add = TRUE, ext = districts)

library(parallel)

# terra
r <- rast(ext(30, 31, -14, -13), res = 0.01, 
          crs = "+proj=longlat +datum=WGS84")
values(r) <- 1:ncell(r)

rp <- terra::wrap(r)
system.time(b <- mclapply(1:2000, function(x, r = rp) {
  r <- unwrap(r)
  values(r) <- sample(1:10, size = ncell(r), replace = TRUE)
  set.seed(x)
  r3 <- r2 <- r
  values(r2) <- values(r) * runif(n = ncell(r), 0.9, 1.2)
  values(r3) <- values(r) * runif(n = ncell(r), 0.8, 1.3)
  wrap(c(r3, r2, r))
}) %>% lapply(., unwrap))