Predict suitability surfaces with observed data: binomial full interaction models

Load packages, read data and source custom scripts

rm(list = ls())
library(mgcv)

#> Loading required package: nlme

#> This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.

library(dplyr)

#> 
#> Attaching package: 'dplyr'

#> The following object is masked from 'package:nlme':
#> 
#>     collapse

#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag

#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

library(purrr)
library(stars)

#> Loading required package: abind

#> Loading required package: sf

#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

path_proj <- day2day::git_path()
path_data <- file.path(path_proj, "data")
path_cleaned <- file.path(path_data, "cleaned")
path_processed <- file.path(path_data, "processed")
path_modelled <- file.path(path_data, "modelled")
path_uk_1km_suitability <- file.path(path_modelled, "uk-1km-suitability")
dir.create(path_uk_1km_suitability, showWarnings = FALSE)

source(file.path(path_proj, "src", "52-predict-classes.R"))
source(file.path(path_proj, "src", "54-vars-to-raster.R"))

uk_bbox_1km <- read_stars(file.path(path_processed, "uk_bbox_1km.tif"))
land_cover <- fst::read_fst(file.path(path_processed, "uk_1km_dataframe_train_full.fst"))

prefix <- "binom_full_interact"
model_file <- file.path(path_modelled, paste0(sub("binom", "model_binomials", prefix), ".rds"))
data_model <- readRDS(model_file)

Predict and save surfaces

chess_label <- "chess"
file_prefix <- paste("uk_1km_suitability", prefix, chess_label, sep = "_")
file_suffix <- paste("20yr-mean-annual", "baseline", sep = "_")
classes <- sub(paste0("^", prefix, "_"), "", data_model$name)

land_cover <- predict_classes(data_model, land_cover)
predict_file <- paste(file_prefix, file_suffix, sep = "_")
predict_file <- file.path(path_uk_1km_suitability, paste0(predict_file, ".fst"))
fst::write_fst(land_cover, predict_file)

# convert each class to raster
pred_raster <- map(data_model$name, ~ vars_to_raster(land_cover, ., uk_bbox_1km))
predict_files <- paste(file_prefix, classes, file_suffix, sep = "_")
predict_files <- file.path(path_uk_1km_suitability, paste0(predict_files, ".tif"))
walk2(pred_raster, predict_files, ~ write_stars(.x, .y))

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>     user   system  elapsed 
#> 1136.090    2.747 1138.992