Predict land suitability UKCP18 SPEED scenarios: binomial mask-out interaction models


Predicting land suitability surfaces for each land class using generalised additive 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"))
source(file.path(path_proj, "src", "57-predict-scenario.R"))

uk_bbox_1km <- read_stars(file.path(path_processed, "uk_bbox_1km.tif"))

speed_tag <- "01"
land_sce <- fst::read_fst(
    file.path(path_processed, paste0("uk_1km_dataframe_ukcp18-speed_", speed_tag, ".fst"))
)

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

Select period

# data_model
vars_no_speed <- c("elev", "slope_nb8")
vars_all <- unique(unlist(map(data_model$form, ~ all.vars(.[[3]]))))
vars_speed <- setdiff(vars_all, vars_no_speed)
vars_sce <- paste0("rcp\\d+_(", paste(vars_speed, collapse = "|"), ")_") %>%
    grep(names(land_sce), value = TRUE)
periods <- sub("^rcp\\d+_[[:alpha:]_]+_([0-9]{6}).*_([0-9]{6}).*", "\\1-\\2", vars_sce) %>%
    unique()
rcps <- unique(sub("^(rcp\\d+)_.+$", "\\1", vars_sce))

Predict and save surfaces

for (rcp in rcps) {
    walk(periods,
        ~ predict_scenario(land_sce, uk_bbox_1km, rcp, ., vars_speed,
                           data_model, path_uk_1km_suitability, speed_tag)
    )
}

Time to execute the task

Only useful when executed with Rscript.

proc.time()
#>      user    system   elapsed 
#> 10747.733     7.747 10763.379