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

# code to adapt
roads %>% mutate(km = as.numeric(st_length(.)) / 1000) %>% 
  filter(km > 200) %>% 
  st_buffer(dist = 20000) %>% 
  st_transform(crs = 4326)
st_difference(larger_object, smaller_object)
ggplot() + geom_sf

buffered_object %>% st_sample()

library(ggplot2)

p <- ggplot(districts) + 
  geom_sf() + 
  geom_sf(data = smallest_centroid, color = "red")

p1 <- p + geom_sf(data = mypoly, aes(color = as.factor(ID)),
                  fill = "transparent") + labs(color = "ID")
# small_pol2 <- pol %>% st_sfc() %>% st_sf(crs = 4326) %>% mutate(ID = 1)
p2 <- p + geom_sf(data = mypoly, aes(color = as.factor(ID)), 
                  fill = "transparent") + 
  labs(color = "ID")

gridExtra::grid.arrange(p1, p2, ncol = 2)

g1 <- cowplot::plot_grid(p1 + theme(legend.position = "none"), 
                         p2 + theme(legend.position = "none"))
cowplot::plot_grid(g1, cowplot::get_legend(p1), rel_widths = c(2, 0.3))

library(patchwork)

p1 + p2

p1 + p2 + plot_layout(guides = "collect")

p1 + p1 / p
(p1 + p2) / (p1 + p2) + plot_layout(guides = "collect")

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

plot(chirps[[3]])

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)
# values(r) <- 1:15
# plot(r)

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

# plot(vect(districts))
# plot(districts)

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

b2 <- 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(b2) <- paste0("b", 1:10)

plot(chirps[[1:10]])
# plot(b2)

b3 <- 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(b3) <- paste0("b", 1:10)
plot(b3)
plot(b3[[1]])
plot(b3[[10]])
plot(b3[[c(1, 10)]])

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

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

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

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

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

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

distsr_pol <- as.polygons(distsr, dissolve = TRUE)
distsr_pol %>% 
  st_as_sf %>% 
  st_geometry %>% 
  plot
distsr_pol %>% 
  st_as_sf %>% 
  slice(49) %>% 
  plot(add = TRUE)

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

plot_noaxes(b[[1]], legend = FALSE)
legend("right", legend = 1:10, fill = terrain.colors(10), border = FALSE, 
       bty = "n")

# plotting with larger datasets
plot_noaxes(b3, main = paste("Random", 1:10))
plot_noaxes(b3, nc = 3, nr = 4, main = paste("Random", 1:10))

rasterVis::levelplot(b3)#, #names.attr = paste("Random", 1:10))
rasterVis::levelplot(b2, names.attr = paste("Random", 1:10))
rasterVis::levelplot(b[[1]])
rasterVis::gplot(b[[1]]) + geom_tile(aes(fill = value)) + 
  scale_fill_viridis_c()

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

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

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

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

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

library(geospaar)
data(chirps)
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)