trying to designate dynamic LU covariates as categorical and plot with LU colors

00_TUTORIAL_script.R

@@ -18,10 +18,12 @@
 
 
 
-# Load required packages --------------------------------------------------
+# Empty memory and workspace; Load required packages ----------------------
+gc()
+rm(list=ls())
 
 pkgs <- c("tidyverse", "foreach", "raster", "terra", "sf", "mapview", "ranger",
-          "caret", "viridis")
+          "caret", "viridis", "hexbin")
 # make sure 'mapview' pkg installed from github to avoid pandoc error:
 # remotes::install_github("r-spatial/mapview")
 lapply(pkgs, library, character.only = TRUE)
@@ -362,45 +364,113 @@ r_stack_LU_1910_present <- stack(
 names(r_stack_LU_1910_present) <- paste0("LU_", seq(1910, format(Sys.Date(), "%Y"), 1))
 
 # calculating mode using raster pkg takes 100 times longer than using terra,
-# so convert to spatRaster using terra and then perform terra::modal()
+# so convert to SpatRaster using terra and then perform terra::modal()
 sr_LU_1910_present <- rast(r_stack_LU_1910_present)
 
 # Obtain land use for every location, with coordinates x and y for every year t
 # between 1953 and 2022, by assigning the same class as in the temporally nearest
 # year for which a map is available.
-# LU_xyt_1km for year 1953 & 2022
-sr_LU_xyt_1km <- imap(c(1953, 2022), ~ {
+# Compute LU_xyt_1km, here only for 2 of the 70 years, e.g. 1990 & 2022
+sr_LU_xyt_1km <- imap(c(1990, 2022), ~ {
   i <- .x
-  sr_LU_1910_present[[paste0("LU_", i)]]
-}) %>%  rast(.)
+  output <- sr_LU_1910_present[[paste0("LU_", i)]]
+  setNames(output, paste0("LU_", i))
+}) %>% rast(.)
 # Sidenote: since output of foreach is list, rast() combines into multilayer SpatRaster
 
 # In the same manner, assign the land use class that occurred most frequently in
 # the 5, 10, 20 and 40 years prior to and including year t...
-# LU_xyt_delta5_1km for year 1953 & 2022
-sr_LU_xyt_delta5_1km <- imap(c(1953, 2022), ~ {
+# Compute LU_xyt_delta5_1km, here only for 2 of the 70 years, e.g. 1990 & 2022
+sr_LU_xyt_delta5_1km <- imap(c(1990, 2022), ~ {
   i <- .x
-  modal(sr_LU_1910_present[[paste0("LU_", seq(i - 4, i))]], na.rm = TRUE)
+  output <- modal(sr_LU_1910_present[[paste0("LU_", seq(i - 4, i))]], na.rm = TRUE)
+  setNames(output, paste0("LU_", i, "_delta5"))
 }) %>% rast(.)
 
-# LU_xyt_delta10_1km for year 1953 & 2022
-sr_LU_xyt_delta10_1km <- imap(c(1953, 2022), ~ {
+# Compute LU_xyt_delta10_1k, here only for 2 of the 70 years, e.g. 1990 & 2022
+sr_LU_xyt_delta10_1km <- imap(c(1990, 2022), ~ {
   i <- .x
-  modal(sr_LU_1910_present[[paste0("LU_", seq(i - 9, i))]], na.rm = TRUE)
+  output <- modal(sr_LU_1910_present[[paste0("LU_", seq(i - 9, i))]], na.rm = TRUE)
+  setNames(output, paste0("LU_", i, "_delta10"))
 }) %>% rast(.)
 
-# LU_xyt_delta20_1km for year 1953 & 2022
-sr_LU_xyt_delta20_1km <- imap(c(1953, 2022), ~ {
+# Compute LU_xyt_delta20_1km, here only for 2 of the 70 years, e.g. 1990 & 2022
+sr_LU_xyt_delta20_1km <- imap(c(1990, 2022), ~ {
   i <- .x
-  modal(sr_LU_1910_present[[paste0("LU_", seq(i - 19, i))]], na.rm = TRUE)
+  output <- modal(sr_LU_1910_present[[paste0("LU_", seq(i - 19, i))]], na.rm = TRUE)
+  setNames(output, paste0("LU_", i, "_delta20"))
 }) %>% rast(.)
 
-# LU_xyt_delta40_1km for year 1953 & 2022
-sr_LU_xyt_delta40_1km <- imap(c(1953, 2022), ~ {
+# Compute LU_xyt_delta40_1km, here only for 2 of the 70 years, e.g. 1990 & 2022
+sr_LU_xyt_delta40_1km <- imap(c(1990, 2022), ~ {
   i <- .x
-  modal(sr_LU_1910_present[[paste0("LU_", seq(i - 39, i))]], na.rm = TRUE)
+  output <- modal(sr_LU_1910_present[[paste0("LU_", seq(i - 39, i))]], na.rm = TRUE)
+  setNames(output, paste0("LU_", i, "_delta40"))
 }) %>% rast(.)
 
+# combine dynamic LU covariates for 1990 and designate as categorical covariates
+sr_LU_xyt_1990_1km <- as.factor(c(sr_LU_xyt_1km$LU_1990,
+                                  sr_LU_xyt_delta5_1km$LU_1990_delta5,
+                                  sr_LU_xyt_delta10_1km$LU_1990_delta10,
+                                  sr_LU_xyt_delta20_1km$LU_1990_delta20,
+                                  sr_LU_xyt_delta40_1km$LU_1990_delta40))
+
+# combine dynamic LU covariates for 2022
+sr_LU_xyt_2022_1km <- c(sr_LU_xyt_1km$LU_2022,
+                        sr_LU_xyt_delta5_1km$LU_2022_delta5,
+                        sr_LU_xyt_delta10_1km$LU_2022_delta10,
+                        sr_LU_xyt_delta20_1km$LU_2022_delta20,
+                        sr_LU_xyt_delta40_1km$LU_2022_delta40)
+
+# In the following lines of code, we designate proper names to the LU categories/classes...
+# read in tables with reclassified values for dynamic LU covariates
+ls_tbl_recl <- foreach(tbl = 1:length(v_cov_dyn_names)) %do%
+  readr::read_csv(paste0("data/covariates/organism/",
+                         gsub(".grd", "", v_cov_dyn_names)[tbl], "_reclassify_dyn.csv"))
+
+# tibble of all LU classes used in dynamic LU covariates
+tbl_LU_xyt_classes <- distinct(ls_tbl_recl[[length(ls_tbl_recl)]][,3:4]) %>% 
+  drop_na() %>% 
+  rename(id = 1, landuse = 2)
+
+# list of tibbles of LU classes specific to each dynamic LU covariate
+# (not all contain all classes)
+ls_tbl_LU_xyt_classes <- map2(
+  rep(list(tbl_LU_xyt_classes), length(nlyr(sr_LU_xyt_1990_2022_1km))),
+  map(as.list(sr_LU_xyt_1990_2022_1km), ~unique(.x)[,1]),
+  ~filter(.x, id %in% .y)
+)
+
+# designate as categorical variables and assign proper classes
+foreach(i = 1:length(ls_tbl_LU_xyt_classes)) %do% {
+  levels(sr_LU_xyt_1990_2022_1km[[i]]) <- as.data.frame(ls_tbl_LU_xyt_classes[[i]])
+}
+
+# list of cols for plotting
+cols <- c("#57c4ad", "#e6e1bc", "#f6d3e8", "#999b99", "#035f05", "#e68ae2",
+          "#eae205", "#e18c0c", "#7809c7", "#f8050f", "#0206fa")
+
+cols[pull(ls_tbl_LU_xyt_classes[[1]], id)]
+
+ls_cols <- map2(rep(list(cols), length(nlyr(sr_LU_xyt_decades_1km))),
+                ls_tbl_LU_xyt_classes, ~.x[pull(.y, id)])
+
+# plot and save plots to disk
+foreach(i = 1:nlyr(sr_LU_xyt_decades_1km)) %do% {
+  pdf(paste0("out/maps/explorative/covariates/dynamic/m_",
+             names(sr_LU_xyt_decades_1km)[i], ".pdf"), width = 9, height = 9)
+  plot(sr_LU_xyt_decades_1km[[i]],
+       col = ls_cols[[i]])
+  dev.off()
+}
+
+# visualize dynamic LU covariate
+mapView(sr_LU_xyt_delta40_1km[[1]],
+        #att = "description",
+        col.regions = c("", "", "", "#006164", "#db4325",
+                        "white", "#386cb0", "#ffff99"),
+        layer.name = "Land use around 1960")
+
 
 
 # Dynamic covariates: peat occurrence -------------------------------------