Land suitability surfaces: visualization and diagnostics

Load packages, read data and source custom scripts

rm(list = ls())
library(dplyr)

#> 
#> Attaching package: 'dplyr'

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

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

library(stars)

#> Loading required package: abind

#> Loading required package: sf

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

library(purrr)

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

source(file.path(path_proj, "src", "54-vars-to-raster.R"))
source(file.path(path_proj, "src", "61-vis-stars.R"))
source(file.path(path_proj, "src", "63-tidy-names.R"))
source(file.path(path_proj, "src", "64-vis-suitability.R"))

uk_bbox_1km <- read_stars(file.path(path_processed, "uk_bbox_1km.tif"))
uk <- st_read(file.path(path_cleaned, "uk_simple.gpkg"), "union")

#> Reading layer `union' from data source `/home/rstudio/Cloud/lancs/projects/land-suitability/data/cleaned/uk_simple.gpkg' using driver `GPKG'
#> Simple feature collection with 1 feature and 0 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -116.1923 ymin: 7054.099 xmax: 655644.8 ymax: 1218625
#> projected CRS:  OSGB 1936 / British National Grid

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

suffixes <- c("binom_full_simple", "binom_full_interact",
              "binom_maskout_simple", "binom_maskout_interact")
predict_files <- file.path(
    path_uk_1km_suitability,
    paste("uk_1km_suitability", suffixes, "chess", "20yr-mean-annual_baseline.fst",
          sep = "_")
)
land_cover <- left_join(land_cover, read_sce(predict_files), by = "id")

Custom function, prepare data and visualize

# prepare data
land_names <- get_land_classes(names(land_cover), suffixes[1])
data_plot <- suitability_prepare(land_cover)
data_plot <- rename_with(data_plot, ~ sub("binom", "binomial", .x))

# Custom function
custom_plot <- function (land_name, nbreaks = 11, ...) {
    cat("\n\n")
    cat(paste0("## ", tidy_make_classes(land_name, 2), "\n\n"))
    cat("**Proportional breaks:**\n")
    suitability_check(data_plot, uk_bbox_1km, land_name, nbreaks, uk, shortname = TRUE, ...)
    cat("\n")
    cat("**Quantile breaks:**", "\n\n")
    suitability_check(data_plot, uk_bbox_1km, land_name, floor(1.3 * nbreaks), uk,
                      shortname = TRUE, breaks = "quantile", ...)
}

# visualize
walk(land_names, ~ custom_plot(.))

(1) Arable

Proportional breaks: Quantile breaks:

(2) Wetland

Proportional breaks: Quantile breaks:

(3) Improved grassland

Proportional breaks: Quantile breaks:

(4) Forest

Proportional breaks: Quantile breaks:

(5) Semi natural grassland

Proportional breaks: Quantile breaks:

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#> 211.620   5.668 217.942