Merge pull request #46 from EvolEcolGroup/biome4_classes

Treat biome classes as factors

R/data.R

@@ -46,3 +46,10 @@
 #'
 #' @format A data.frame with multiple columns to describe.
 "koeppen_classes"
+
+#' BIOME4 classes.
+#'
+#' A data.frame defining the details of each class
+#'
+#' @format A data.frame with multiple columns to describe.
+"biome4_classes"

R/get_ice_mask.R

@@ -33,6 +33,18 @@ get_ice_mask <- function(time_bp = NULL, dataset) {
     ice_mask <- climate_series["biome"]
     ice_mask[ice_mask != 28] <- NA
     ice_mask[ice_mask == 28] <- 1
+    # sort out categories
+      # we pass a list so that each level if turned into a categorical variable
+      levels(ice_mask) <-
+        rep(list(data.frame(id=1, category= c("ice"))),
+            terra::nlyr(ice_mask))
+      terra::coltab(ice_mask) <- rep (list(
+        data.frame(
+          values = c(1),
+          cols = c("#B5D1E0")
+        )
+      ), terra::nlyr(ice_mask))
+
     names(ice_mask) <- paste("ice_mask", time_bp(ice_mask), sep = "_")
     varnames(ice_mask) <- "ice_mask"
     return(ice_mask)

R/get_land_mask.R

@@ -50,7 +50,20 @@ get_land_mask <- function(time_bp = NULL, time_ce = NULL, dataset) {
   } else {
     stop("no method yet for this dataset")
   }
-
+  # sort out categories
+  # we pass a list so that each level if turned into a categorical variable
+  levels(land_mask) <-
+    rep(list(data.frame(id=1, category= c("land"))),
+        terra::nlyr(land_mask))
+  terra::coltab(land_mask) <- rep (list(
+    data.frame(
+      values = c(1),
+      cols = c("#1C540F")
+    )
+  ), terra::nlyr(land_mask))
+
+
+
   if (is.null(time_ce)) {
     names(land_mask) <- paste("land_mask", time_bp(land_mask), sep = "_")
   } else {

R/location_slice_from_region_series.R

@@ -109,7 +109,6 @@ location_slice_from_region_series <-
       coords_df <- locations_data[, c("cell_number")]
     }
 
-
     # now sort out the time slices corresponding to each location
     times <- time_bp(climate_brick)
     time_indeces <- time_bp_to_index(
@@ -129,7 +128,7 @@ location_slice_from_region_series <-
           x = this_slice,
           y = locations_data[locations_data$time_bp_slice == i_time, coords]
         )
-        # sor tout the indexing here
+        # sort out the indexing here
         locations_data[locations_data$time_bp_slice == i_time, bio_variables] <-
           this_climate[
             ,
@@ -139,6 +138,10 @@ location_slice_from_region_series <-
         locations_data[this_slice_indeces, ] <- NA
       }
 
+      # factors don't behave nicely when adding new elements, cast to character
+      if ("biome" %in% names(locations_data)){
+        locations_data$biome <- as.character(locations_data$biome)
+      }
       if (nn_interpol | buffer) {
         locations_to_move <- this_slice_indeces[this_slice_indeces %in%
           which(!stats::complete.cases(locations_data))]
@@ -165,12 +168,12 @@ location_slice_from_region_series <-
               y = neighbours_ids[1, ]
             ) # [, bio_variables]
 
-          neighbours_values_mean <- apply(neighbours_values, 2,
+          neighbours_values_mean <- apply(neighbours_values[,!names(neighbours_values) %in% "biome"], 2,
             mean,
             na.rm = T
           )
           if ("biome" %in% bio_variables) {
-            neighbours_values_mean["biome"] <- mode(neighbours_values[, "biome"])
+            neighbours_values_mean["biome"] <- mode(as.numeric(neighbours_values[, "biome"]))
           }
           locations_data[i, bio_variables] <-
             neighbours_values_mean[bio_variables]
@@ -191,6 +194,12 @@ location_slice_from_region_series <-
       locations_data$time_ce_slice <- locations_data$time_bp_slice + 1950
       locations_data <- locations_data[, !names(locations_data) %in% c("time_bp", "time_bp_slice")]
     }
+
+    # reintroduce the factor
+    if ("biome" %in% bio_variables){
+      locations_data$biome <- factor(levels(region_series$biome)[[1]]$category[as.numeric(locations_data$biome)],
+                                     levels = levels(region_series$biome)[[1]]$category)
+    }
     return(locations_data)
   }
 

R/region_series.R

@@ -121,6 +121,21 @@ region_series <-
         var_brick <- terra::crop(var_brick, crop)
       }
 
+      # special treatment for biome variables
+      if (this_var_nc == "biome") {
+        # we pass a list so that each level if turned into a categorical variable
+        levels(var_brick) <-
+          rep(list(get_biome_classes(dataset = dataset)),
+              terra::nlyr(var_brick))
+        terra::coltab(var_brick) <- rep (list(
+          data.frame(
+            values = get_biome_classes(dataset = dataset)$id,
+            cols = pastclim::biome4_classes$colour[match(get_biome_classes(dataset = dataset)$id,
+                                                         pastclim::biome4_classes$id)]
+          )
+        ), terra::nlyr(var_brick))
+      }
+
       climate_spatrasters[[this_var]] <- var_brick
 
       varnames(climate_spatrasters[[this_var]]) <- this_var

data-raw/data_files/biome4_classes.csv

@@ -0,0 +1,30 @@
+id,category,colour
+0,Water bodies,47 61 73
+1,Tropical evergreen forest,11 33 6
+2,Tropical semi-deciduous forest,40 57 3
+3,Tropical deciduous forest/woodland,68 49 12
+4,Temperate deciduous forest,33 92 28
+5,Temperate conifer forest,10 51 44
+6,Warm mixed forest,0 0 40
+7,Cool mixed forest,79 100 79
+8,Cool conifer forest,0 60 9
+9,Cold mixed forest,44 14 20
+10,Evegreen taiga/montane forest,0 12 79
+11,Deciduous taiga/montane forest,40 70 100
+12,Tropical savanna,73 100 21
+13,Tropical xerophytic shrubland,55 63 15
+14,Temperate xerophytic shrubland,10 87 79
+15,Temperate sclerophyll woodland,55 63 15
+16,Temperate broadleaved savanna,45 100 20
+17,Open conifer woodland,100 60 87
+18,Boreal parkland,40 49 100
+19,Tropical grassland,100 73 21
+20,Temperate grassland,100 88 60
+21,Desert,97 100 79
+22,Steppe tundra,90 90 9
+23,Shrub tundra,40 100 60
+24,Dwarf shrub tundra,48 52 29
+25,Prostrate shrub tundra,82 62 60
+26,Cushion forb lichen moss tundra,60 40 100
+27,Barren,73 71 67
+28,Land ice,71 82 88

data-raw/make_data/biome4_classes.R

@@ -0,0 +1,7 @@
+## code to prepare `koeppen_classes` dataset
+## source this script in RStudio
+biome4_classes <- read.csv("./data-raw/data_files/biome4_classes.csv")
+biome4_classes$colour <- sapply(strsplit(biome4_classes$colour, " "), function(x) {
+  rgb(x[1], x[2], x[3], maxColorValue = 100)
+})
+usethis::use_data(biome4_classes, overwrite = TRUE)

tests/testthat/test_location_slice.R

@@ -58,7 +58,7 @@ test_that("location_slice", {
     dataset = "Example", nn_interpol = TRUE
   )
   # checked value by hand
-  expect_true(this_climate_biome$biome[4] == 17)
+  expect_true(as.numeric(this_climate_biome$biome[4]) == 18)
 
   # now use the full dataframe for pretty labelling
   this_climate_df <- location_slice(

vignettes/a0_pastclim_overview.Rmd

@@ -564,20 +564,25 @@ detailing the extension of biomes.
 get_biome_classes("Example")
 ```
 
-If we need to plot the extent of a specific biome, for example the
-desert, we first extract the variable and then subset it to just the
-class we are interested in using its ID (21, in this case):
-
-```{r, fig.width=6, fig.height=2.5}
+We can get the biome for 20k years ago and plot it with:
+```{r, fig.width=8, fig.height=6}
 biome_20k <- region_slice(
   time_bp = -20000,
   bio_variables = c("biome"),
   dataset = "Example"
 )
+plot(biome_20k)
+```
+
+Note that the legend is massive. When plotting multiple time slices, it is best
+to use `legned=FALSE` in the plotting statement to avoid having the legend.
+If we need to plot the extent of a specific biome, for example the
+desert, we can simply set the other levels to NA:
+
+```{r, fig.width=6, fig.height=2.5}
 biome_20k$desert <- biome_20k$biome
-biome_20k$desert[biome_20k$desert != 21] <- FALSE
-biome_20k$desert[biome_20k$desert == 21] <- TRUE
-terra::plot(biome_20k)
+biome_20k$desert[biome_20k$desert != 21] <- NA
+terra::plot(biome_20k$desert)
 ```
 
 The climate reconstructions do not show areas under permanent ice. Ice