R/read_soilmap.R

In inbo/n2khab: Providing Preprocessed Reference Data for Flemish Natura 2000 Habitat Analyses

#' Return the \code{soilmap} or \code{soilmap_simple} data source as an
#' \code{sf} multipolygon layer
#'
#' Returns either the raw data source \code{soilmap}
#' or (by default) the
#' processed data source \code{soilmap_simple}
#' as a standardized \code{sf} multipolygon layer (tidyverse-styled,
#' internationalized) in the Belgian Lambert 72 CRS (EPSG-code
#' \href{https://epsg.io/31370}{31370}).
#' Given the size of these data sources (especially the raw one), this function
#' takes a bit longer than usual to run.
#'
#' The raw data source is published
#' [at DOV](https://www.dov.vlaanderen.be/geonetwork/srv/dut/catalog.search#/metadata/5c129f2d-4498-4bc3-8860-01cb2d513f8f)
#' (Databank Ondergrond Vlaanderen) and
#' is discussed by Van Ranst & Sys (2000) and Dudal et al. (2005).
#' A 'pure' (single) dataformat of the raw data source (no metadatafiles etc.)
#' has also been stored (with versioning) at
#' Zenodo (\doi{10.5281/zenodo.3387007}) - which we refer to as the
#' `soilmap` data source - in order
#' to support the `read_soilmap()` function and to sustain long-term workflow
#' reproducibility.
#'
#' The processed data source `soilmap_simple` is a GeoPackage, available at
#' \href{https://doi.org/10.5281/zenodo.3732903}{Zenodo}.
#'
#' Note that factors are generated with implicit \code{NA} values (i.e. there is
#' no factor level to represent the missing values).
#' If you want this category to appear in certain results, you can convert
#' such variables with
#' [forcats::fct_explicit_na()].
#'
#' In case the raw data source \code{soilmap} is used
#' (\code{use_processed = FALSE}), it is possible to
#' manually perform the standardization for coastal plain features and/or the
#' simplification,
#' both of which were applied in the \code{soilmap_simple} data source.
#' See \emph{Arguments} for more information.
#'
#' See R-code in the \href{https://github.com/inbo/n2khab-preprocessing}{
#' n2khab-preprocessing} repository for the creation of the
#' \code{soilmap_simple} data source from
#' the \code{soilmap} data source.
#'
#' @md
#'
#' @param file The absolute or relative file path of the _processed_ data
#' source `soilmap_simple`.
#' Used only if \code{use_processed = TRUE} (= default).
#' The default value follows the data management advice in the
#' vignette on data storage (run \code{vignette("v020_datastorage")}).
#' It uses the first \code{n2khab_data} folder that is found when
#' sequentially climbing up 0 to 10 levels in the file system hierarchy,
#' starting from the working directory.
#' @param file_raw Same as `file`, to define the filepath of the _raw_
#' datasource `soilmap`.
#' Used only if \code{use_processed = FALSE}.
#' @param use_processed Logical.
#' If \code{TRUE} (the default), load and return the processed data source
#' \code{soilmap_simple}, instead of the raw data source \code{soilmap}.
#' The central layer of
#' `soilmap_simple` can be manually generated by reading the raw `soilmap` data
#' source with `standardize_coastalplain=TRUE`, `simplify=TRUE` and
#' `explan=FALSE`, but this takes some time.
#' @param version_processed Version ID of the `soilmap_simple` data source.
#' Only used with \code{use_processed = TRUE} (the default).
#' Defaults to the latest available version defined by the package.
#' @param standardize_coastalplain Logical.
#' Only applied with \code{use_processed = FALSE}.
#' If \code{TRUE}, fill the values of the morphogenetic substrate,
#' texture and drainage variables (`bsm_mo_substr`, `bsm_mo_tex` and
#' `bsm_mo_drain` + their `_explan` counterparts) where possible,
#' _for features with a geomorphological soil type code_.
#' This largely applies to features in the 'coastal plain' area.
#' - To derive morphogenetic texture and drainage levels from the geomorphological
#' soil types, a conversion table by Bruno De Vos &
#' Carole Ampe is applied (for earlier work on this, see Ampe 2013).
#' - Substrate classes are copied over from `bsm_ge_substr` into `bsm_mo_substr`
#' (`bsm_ge_substr` already follows the categories of `bsm_mo_substr`).
#' - A logical variable is added to the output to mark conversions
#' (see section _Value_).
#'
#' These steps coincide with the approach that was taken to construct
#' `bsm_mo_soilunitype` in the raw data source.
#' @param simplify Logical.
#' Only applied with \code{use_processed = FALSE}.
#' If \code{TRUE}, only a limited number of variables that are most
#' useful for analytical work are returned.
#' @param explan Logical, defaults to `FALSE`.
#' Should the `_explan` variables accompanying `bsm_mo_xxx` variables be
#' returned in the _simplified_ result?
#' If \code{use_processed = FALSE}: only has effect if `simplify=TRUE`.
#' (With `simplify=FALSE`, the `_explan` variables are always returned.)
#' If \code{use_processed = TRUE}: is applied in returning `soilmap_simple`.
#'
#' @return
#' A Simple feature collection of geometry type \code{MULTIPOLYGON},
#' representing either the processed data source `soilmap_simple` (default) or
#' the raw data source `soilmap`.
#'
#' Besides the standardization for the coastal plain areas, `soilmap_simple`
#' contains
#' only a subset of the `soilmap` variables (marked with an asterisk below).
#'
#' The `soilmap` attribute variables all start with prefix `bsm_` (referring
#' to the 'Belgian soil map'), in order to distinguish from similar attributes
#' derived from other maps or field observations.
#'
#' Most attributes represent categories and are returned as
#' factors.
#' When a variable is a one-to-one translation of another (e.g. code vs.
#' explanation), the order of factor levels is aligned.
#'
#' Three types of data frame variables are returned when reading `soilmap`:
#' - **variables with `mo_` in their name**: their categories follow the
#' Belgian Morphogenetic System.
#'   - With `standardize_coastalplain = FALSE`, these are only available _outside
#' the coastal plain areas_ except for `bsm_mo_soilunitype` (which is
#' standardized already in the raw data source).
#' - **variables with `ge_` in their name**: their categories follow the
#' Belgian Geomorphological System.
#' (Note however, that `bsm_ge_substr` does follow the Belgian Morphogenetic
#' System as well.)
#'   - Values are typically available _within the coastal plain areas_,
#'   but some geomorphological soil types (starting with letter `O`)
#'   have a wider distribution across Flanders.
#'   - They are _not_ included in `soilmap_simple`.
#'   - A special variable is `bsm_ge_typology`, which is `TRUE` if
#'   `bsm_soiltype` follows the geomorphological typology,
#'   and `FALSE` otherwise.
#' - **variables without `mo_` or `ge_` in their name** are:
#'   - either _system-agnostic_ metadata (first two + last four variables:
#'   `bsm_poly_id`, `bsm_map_id`, `bsm_map_url`, `bsm_book_url`,
#'   `bsm_detailmap_url`, `bsm_profloc_url`),
#'   - or _mixed_ (representing `mo_` categories within and `ge_` categories
#'   outside coastal plains): the other ones, like `bsm_region`, `bsm_legend`,
#'   `bsm_soiltype` and `bsm_soilseries`.
#'   - A special variable is `bsm_converted`, returned only if
#'   `standardize_coastalplain = TRUE`.
#'
#' Many variables have a 'counterpart variable' with suffix `_explan`:
#' they provide a more elaborate textual explanation.
#' They are not listed below.
#'
#' Short explanation of attributes is given below.
#' More elaborate explanations can be found in the references and in metadata
#' [at DOV](https://www.dov.vlaanderen.be/geonetwork/srv/dut/catalog.search#/metadata/5c129f2d-4498-4bc3-8860-01cb2d513f8f).
#' 1. Meaning of the main non-metadata variables:
#'     - `bsm_region` (*): name of the region
#'     - `bsm_ge_region`: code of the region within the coastal plain area
#'     - `bsm_legend`: generalised (simplified) legend key (37 levels)
#'     - `bsm_legend_title` and `bsm_legend_explan`:
#'     the legend keys and text of Van Ranst & Sys (2000) (833 and 622 levels,
#'     respectively)
#'     - `bsm_soiltype`: the soil type of the Belgian soil map (mixed nature:
#'     morphogenetic & geomorphological codes).
#'     `bsm_soiltype_id` represents a numeric code for each level.
#'     - `bsm_ge_typology`: Logical.
#'     Does the soiltype code follow the geomorphological typology?
#'     - `bsm_soiltype_region`: `bsm_soiltype`, followed by a code representing
#'     `bsm_region`
#'     - `bsm_soilseries`: either the morphogenetic soil series (outside the
#'     coastal plain areas),
#'     which is the three core characters of `bsm_soiltype`,
#'     or just `bsm_soiltype` if the latter has a geomorphological code.
#'     - `bsm_converted` (*): Logical.
#'     Were morphogenetic texture and drainage variables (`bsm_mo_tex` and
#'     `bsm_mo_drain`) derived from a conversion table?
#'     This is equivalent with the question: does `bsm_mo_soilunitype` differ
#'     from `bsm_soiltype`?
#'     Value `TRUE` is largely confined to the 'coastal plain' areas.
#'     Only returned if `standardize_coastalplain = TRUE`.
#'     (Note: the variable is not included in version `soilmap_simple_v1`.)
#'     - `bsm_mo_soilunitype` (*): as `bsm_soiltype`, but applying morphogenetic
#'     codes within the coastal plain areas in most cases
#'     (see the `standardize_coastalplain`
#'     argument for more information about this conversion)
#'     - `bsm_mo_substr` (*), `bsm_ge_substr`: code of the soil substrate
#'     - `bsm_mo_tex` (*): code of the soil texture category
#'     - `bsm_mo_drain` (*): code of the soil drainage category
#'     - `bsm_mo_prof` (*): code of the soil profile category
#'     - `bsm_mo_parentmat` (*): code of a variant regarding the parent material
#'     - `bsm_mo_profvar` (*): code of a variant regarding the soil profile
#'     - `bsm_mo_phase`: code of the soil phase (i.e. additional soil
#'     properties).
#'     They are explained in the book that accompanies the specific analog map
#'     identified by `bsm_map_id`.
#'     - `bsm_ge_series`: the geomorphological soil series
#'     - `bsm_ge_subseries`: the geomorphological soil subseries
#' 1. Meaning of the metadata variables:
#'     - `bsm_poly_id` (*): unique polygon ID (numeric)
#'     - `bsm_map_id`: code of the analog map covering this area
#'     - `bsm_map_url`: hyperlink to the scanned analog map scale 1:20000 (pdf),
#'     identified by `bsm_map_id`
#'     - `bsm_bookurl`: hyperlink to the scanned book (pdf), accompanying the
#'     analog map identified by `bsm_map_id`
#'     - `bsm_detailmap_url`: hyperlink to the scanned maps at scale 1:5000
#'     (zip-file with jpg files) belonging to the map identified by
#'     `bsm_map_id`
#'     - `bsm_profloc_url`: hyperlink to the scanned maps with the profile
#'     locations
#'     (zip-file with jpg files) belonging to the map identified by
#'     `bsm_map_id`
#'
#' (*) Included in the `soilmap_simple` data source.
#'
#' @family functions returning environmental data sets
#'
#' @references
#' - Ampe C. (2013). Databank aardewerk Vlaanderen 2010.
#' Omzetten (zeer) oude legende bodemkartering naar legende bodemkaart
#' Kuststreek. Vlaamse Landmaatschappij Regio West, Bruges, 45 p.
#' - Dudal R., Deckers J., Van Orshoven J. & Van Ranst E. (2005). Soil survey
#' in Belgium and its applications. In: Bullock P., Jones R.J.A., Montanarella
#' L. (editors). Soil Resources of Europe. Office for Official Publications of
#' the European Communities, Luxembourg, p. 63–71.
#' URL: <http://hdl.handle.net/1854/LU-368514>.
#' - Van Ranst E. & Sys C. (2000). Eenduidige legende van de digitale
#' bodemkaart van Vlaanderen (schaal 1: 20000). Universiteit Gent, Laboratorium
#' voor Bodemkunde, Ghent, 361 p. URL: <http://hdl.handle.net/1854/LU-125899>.
#'
#' @examples
#' \dontrun{
#' # This example supposes that your working directory or a directory up to 10
#' # levels above has the 'n2khab_data' folder AND that the latest version of
#' # the 'soilmap_simple'
#' # data source is present in the default subdirectory.
#' # In all other cases, this example won't work but at least you can
#' # consider what to do.
#'
#' library(dplyr)
#' soilmap_simple <- read_soilmap()
#' soilmap_simple
#' soilmap_simple %>%
#'   filter(!is.na(bsm_mo_substr)) %>%
#'   glimpse()
#' soilmap_simple %>%
#'   filter(bsm_converted) %>%
#'   glimpse()
#' }
#'
#' @importFrom assertthat
#' assert_that
#' is.flag
#' noNA
#' is.string
#' @importFrom sf
#' read_sf
#' @importFrom git2rdata
#' read_vc
#' @importFrom dplyr
#' %>%
#' select
#' mutate
#' mutate_at
#' mutate_if
#' filter
#' filter_at
#' left_join
#' vars
#' contains
#' recode
#' matches
#' distinct
#' pull
#' everything
#' @importFrom stats
#' setNames
#' @importFrom rlang .data
#' @export
read_soilmap <-
  function(file = file.path(
             locate_n2khab_data(),
             "20_processed/soilmap_simple/soilmap_simple.gpkg"
           ),
           file_raw = file.path(locate_n2khab_data(), "10_raw/soilmap"),
           use_processed = TRUE,
           version_processed = "soilmap_simple_v2",
           standardize_coastalplain = FALSE,
           simplify = FALSE,
           explan = FALSE) {
    assert_that(is.flag(simplify), noNA(simplify))
    assert_that(
      is.flag(standardize_coastalplain),
      noNA(standardize_coastalplain)
    )
    assert_that(is.flag(use_processed), noNA(use_processed))
    assert_that(is.flag(explan), noNA(explan))
    assert_that(is.string(version_processed))

    ####### 1. Reading soilmap_simple ####
    ######################################

    if (use_processed) {
      soilmap_simple_path <- file
      assert_that(file.exists(soilmap_simple_path))

      soilmap_simple <-
        read_sf(
          soilmap_simple_path,
          "soilmap_simple"
        ) %>%
        mutate_if(is.character, factor)

      if (explan & version_processed == "soilmap_simple_v1") {
        warning(
          "Version soilmap_simple_v1 is not supported for ",
          "adding explanatory variables."
        )
      }

      if (explan & version_processed != "soilmap_simple_v1") {
        suppressWarnings(
          explanations <-
            read_sf(
              soilmap_simple_path,
              "explanations"
            ) %>%
            split(.$subject) %>%
            lapply(function(df) {
              select(df, -.data$subject)
            })
        )

        soilmap_simple <-
          soilmap_simple %>%
          mutate(
            bsm_mo_substr_explan =
              namelist_factor(.data$bsm_mo_substr,
                codelist = explanations[["bsm_mo_substr"]]
              ),
            bsm_mo_tex_explan =
              namelist_factor(.data$bsm_mo_tex,
                codelist = explanations[["bsm_mo_tex"]]
              ),
            bsm_mo_drain_explan =
              namelist_factor(.data$bsm_mo_drain,
                codelist = explanations[["bsm_mo_drain"]]
              ),
            bsm_mo_prof_explan =
              namelist_factor(.data$bsm_mo_prof,
                codelist = explanations[["bsm_mo_prof"]]
              ),
            bsm_mo_parentmat_explan =
              namelist_factor(.data$bsm_mo_parentmat,
                codelist = explanations[["bsm_mo_parentmat"]]
              ),
            bsm_mo_profvar_explan =
              namelist_factor(.data$bsm_mo_profvar,
                codelist = explanations[["bsm_mo_profvar"]]
              )
          ) %>%
          select(
            .data$bsm_poly_id:.data$bsm_mo_soilunitype,
            !!(c(5, 12) + rep(0:5, each = 2))
          )
      }

      if (version_processed == "soilmap_simple_v1") {
        soilmap_simple <-
          select(
            soilmap_simple,
            -.data$bsm_ge_coastalplain
          )
        suppressWarnings(st_crs(soilmap_simple) <- 31370)
      }

      return(soilmap_simple)
    } else {
      ####### 2. Reading soilmap        ####
      ######################################

      soilmap_path <- file_raw
      assert_that(file.exists(soilmap_path))

      suppressWarnings(
        soilmap <- read_sf(soilmap_path,
          crs = 31370
        )
      )

      soilmap <-
        soilmap %>%
        convertdf_enc(from = "latin1", to = "UTF-8") %>%
        select(
          bsm_poly_id = .data$gid,
          bsm_map_id = .data$Kaartbldnr,
          bsm_region = .data$Streek,
          bsm_ge_region = .data$Streek_c,
          bsm_legend = .data$Grove_leg,
          bsm_legend_title = .data$Uitleg_tit,
          bsm_legend_explan = .data$Uitleg,
          bsm_soiltype_id = .data$codeid,
          bsm_soiltype = .data$Bodemtype,
          bsm_ge_typology = .data$Type_class,
          bsm_soiltype_region = .data$Bodtypstr,
          bsm_soilseries = .data$Bodemser_c,
          bsm_soilseries_explan = .data$Bodemserie,
          bsm_mo_soilunitype = .data$Unitype,
          bsm_mo_substr = .data$Substr_V_c,
          bsm_mo_substr_explan = .data$SubstraatV,
          bsm_mo_tex = .data$Textuur_c,
          bsm_mo_tex_explan = .data$Textuur,
          bsm_mo_drain = .data$Drainage_c,
          bsm_mo_drain_explan = .data$Drainage,
          bsm_mo_prof = .data$Profontw_c,
          bsm_mo_prof_explan = .data$Profontw,
          bsm_mo_parentmat = .data$Varimoma_c,
          bsm_mo_parentmat_explan = .data$Varimoma,
          bsm_mo_profvar = .data$Variprof_c,
          bsm_mo_profvar_explan = .data$Variprof,
          bsm_mo_phase = .data$Fase_c,
          bsm_ge_substr = .data$Substr_p_c,
          bsm_ge_substr_explan = .data$Substr_pol,
          bsm_ge_series = .data$Serie_c,
          bsm_ge_series_explan = .data$Serie,
          bsm_ge_subseries = .data$Subserie_c,
          bsm_ge_subseries_explan = .data$Subserie,
          bsm_map_url = .data$Scan_kbl,
          bsm_book_url = .data$Scan_boek,
          bsm_detailmap_url = .data$Scan_5000,
          bsm_profloc_url = .data$Scan_stip
        ) %>%
        mutate(bsm_ge_typology = .data$bsm_ge_typology == "Zeepolders") %>%
        mutate_at(
          .vars = vars(
            -.data$bsm_poly_id,
            -.data$bsm_soiltype_id,
            -.data$bsm_ge_typology,
            -.data$bsm_soiltype_id,
            -.data$geometry,
            -matches(".+_mo_.+_explan")
          ),
          .funs = factor
        )

      # setting factor levels of mo_*_explan variables in the same
      # order as mo_*

      keyvars <- c(
        "bsm_mo_substr",
        "bsm_mo_tex",
        "bsm_mo_drain",
        "bsm_mo_prof",
        "bsm_mo_parentmat",
        "bsm_mo_profvar"
      )
      keys <- list()
      soilmap_df <-
        soilmap %>%
        st_drop_geometry()
      for (i in keyvars) {
        temp_df <-
          soilmap_df %>%
          select(matches(str_c(i, ".*"))) %>%
          select(1:2) %>%
          filter_at(1, function(x) !is.na(x)) %>%
          distinct()
        keys[[i]] <-
          setNames(
            temp_df %>% pull(2),
            temp_df %>% pull(1)
          )
      }

      soilmap <-
        soilmap %>%
        mutate(
          bsm_mo_substr_explan = recode(
            .data$bsm_mo_substr,
            !!!keys[["bsm_mo_substr"]]
          ),
          bsm_mo_tex_explan = recode(
            .data$bsm_mo_tex,
            !!!keys[["bsm_mo_tex"]]
          ),
          bsm_mo_drain_explan = recode(
            .data$bsm_mo_drain,
            !!!keys[["bsm_mo_drain"]]
          ),
          bsm_mo_prof_explan = recode(
            .data$bsm_mo_prof,
            !!!keys[["bsm_mo_prof"]]
          ),
          bsm_mo_parentmat_explan = recode(
            .data$bsm_mo_parentmat,
            !!!keys[["bsm_mo_parentmat"]]
          ),
          bsm_mo_profvar_explan = recode(
            .data$bsm_mo_profvar,
            !!!keys[["bsm_mo_profvar"]]
          ),
        )

      ####### STANDARDIZATION FOR COASTAL PLAIN FEATURES ################

      if (standardize_coastalplain) {
        transl <- read_vc(
          file = "soil_translation_coastalplain",
          root = pkgdatasource_path(
            "textdata/soil_translation_coastalplain", ".yml"
          )
        ) %>%
          mutate(soiltype_orig = factor(.data$soiltype_orig,
            levels =
              levels(soilmap$bsm_soiltype)
          )) %>%
          filter(!is.na(.data$texture_transl)) %>%
          mutate(
            tex_explan_transl = recode(
              .data$texture_transl,
              !!!keys[["bsm_mo_tex"]]
            ),
            drain_explan_transl = recode(
              .data$drainage_transl,
              !!!keys[["bsm_mo_drain"]]
            ),
            bsm_converted = TRUE
          )

        stand_vars <- c(
          "bsm_mo_substr",
          "bsm_mo_tex",
          "bsm_mo_drain"
        )
        soilmap <-
          soilmap %>%
          left_join(transl, by = c("bsm_soiltype" = "soiltype_orig")) %>%
          mutate_at(
            c(stand_vars, paste0(stand_vars, "_explan")),
            as.character
          ) %>%
          mutate(
            bsm_mo_substr = ifelse(is.na(.data$bsm_mo_substr) &
              !is.na(.data$bsm_ge_substr),
            as.character(.data$bsm_ge_substr),
            .data$bsm_mo_substr
            ) %>%
              factor(levels = levels(soilmap$bsm_mo_substr)),
            bsm_mo_substr_explan = ifelse(is.na(.data$bsm_mo_substr_explan) &
              !is.na(.data$bsm_ge_substr_explan),
            as.character(.data$bsm_ge_substr_explan),
            .data$bsm_mo_substr_explan
            ) %>%
              factor(levels = levels(soilmap$bsm_mo_substr_explan)),
            bsm_mo_tex = ifelse(is.na(.data$bsm_mo_tex) &
              !is.na(.data$texture_transl),
            .data$texture_transl,
            .data$bsm_mo_tex
            ) %>%
              factor(levels = levels(soilmap$bsm_mo_tex)),
            bsm_mo_tex_explan = ifelse(is.na(.data$bsm_mo_tex_explan) &
              !is.na(.data$tex_explan_transl),
            .data$tex_explan_transl,
            .data$bsm_mo_tex_explan
            ) %>%
              factor(levels = levels(soilmap$bsm_mo_tex_explan)),
            bsm_mo_drain = ifelse(is.na(.data$bsm_mo_drain) &
              !is.na(.data$drainage_transl),
            .data$drainage_transl,
            .data$bsm_mo_drain
            ) %>%
              factor(levels = levels(soilmap$bsm_mo_drain)),
            bsm_mo_drain_explan = ifelse(is.na(.data$bsm_mo_drain_explan) &
              !is.na(.data$drain_explan_transl),
            .data$drain_explan_transl,
            .data$bsm_mo_drain_explan
            ) %>%
              factor(levels = levels(soilmap$bsm_mo_drain_explan)),
            bsm_converted = ifelse(is.na(.data$bsm_converted),
              FALSE,
              .data$bsm_converted
            )
          ) %>%
          select(-contains("transl")) %>%
          select(
            .data$bsm_poly_id:.data$bsm_soilseries_explan,
            .data$bsm_converted,
            everything()
          )
      }

      ########## SIMPLIFICATION ############

      if (simplify) {
        soilmap <-
          soilmap %>%
          {
            if (standardize_coastalplain) {
              .
            } else {
              mutate(.,
                bsm_converted = NA
              )
            }
          } %>%
          select(
            .data$bsm_poly_id,
            .data$bsm_region,
            .data$bsm_converted,
            .data$bsm_mo_soilunitype,
            .data$bsm_mo_substr,
            .data$bsm_mo_substr_explan,
            .data$bsm_mo_tex,
            .data$bsm_mo_tex_explan,
            .data$bsm_mo_drain,
            .data$bsm_mo_drain_explan,
            .data$bsm_mo_prof,
            .data$bsm_mo_prof_explan,
            .data$bsm_mo_parentmat,
            .data$bsm_mo_parentmat_explan,
            .data$bsm_mo_profvar,
            .data$bsm_mo_profvar_explan
          ) %>%
          {
            if (explan) . else select(., -matches("_explan"))
          } %>%
          {
            if (standardize_coastalplain) {
              .
            } else {
              select(., -.data$bsm_converted)
            }
          }
      }

      return(soilmap)
    }
  }