R/region_series.R

In pastclim: Manipulate Time Series of Climate Reconstructions

#' Extract a time series of climate variables for a region
#'
#' This function extracts a time series of one or more climate variables for
#' a given
#' dataset covering a region (or the whole world). The function returns a
#' [`terra::SpatRasterDataset`] object, with
#' each variable as a sub-dataset.
#'
#' @param time_bp time slices in years before present (negative values represent
#' time before present, positive values time in the future). This parameter can
#' be a vector of times (the slices need
#' to exist in the dataset), a list with a min and max element setting the
#' range of values, or left to NULL to retrieve all time steps.
#' To check which slices are available, you can use
#' [get_time_bp_steps()].
#' @param time_ce time slices in years CE (see `time_bp` for options).
#' For available time slices in years CE, use [get_time_ce_steps()].
#' Only one of `time_bp` or `time_ce` should be used.
#' @param bio_variables vector of names of variables to be extracted
#' @param dataset string defining the dataset to use. If set to "custom",
#' then a single nc file is used from "path_to_nc"
#' @param path_to_nc the path to the custom nc file containing the palaeoclimate
#' reconstructions. All the variables of interest need to be included in
#' this file.
#' @param ext an extent, coded as numeric vector (length=4;
#' order= xmin, xmax, ymin, ymax) or a [terra::SpatExtent] object. If NULL,
#' the full extent of the reconstruction is given.
#' @param crop a polygon used to crop the reconstructions (e.g. the outline
#' of a continental mass). A [`sf::sfg`][sf::st] or a [terra::SpatVector] object
#' is used to define the polygon.
#' @returns a
#' [`terra::SpatRasterDataset`] object, with
#' each variable as a sub-dataset.
#'
#' @import terra
#' @export

region_series <-
  function(time_bp = NULL,
           time_ce = NULL,
           bio_variables,
           dataset,
           path_to_nc = NULL,
           ext = NULL,
           crop = NULL) {
    time_bp <- check_time_vars(time_bp = time_bp, time_ce = time_ce)

    check_dataset_path(dataset = dataset, path_to_nc = path_to_nc)
    if (dataset == "Barreto2023") {
      message("This is a large dataset, it might take a while...")
    }

    if (!is.null(ext)) {
      if (!any(
        inherits(ext, "SpatExtent"),
        all(inherits(ext, "numeric"), length(ext) == 4)
      )) {
        stop(
          "ext should be a numeric vector of length 4 or a ",
          "terra::SpatExtent object created terra::ext"
        )
      }
      if (inherits(ext, "numeric")) {
        ext <- terra::ext(ext)
      }
    }
    if (!is.null(crop)) {
      if (!any(
        inherits(crop, "SpatVector"),
        inherits(crop, "sfg")
      )) {
        stop(
          "crop should be a sf::sfg or a terra::SpatVector ",
          "object created with terra::vect"
        )
      }
      if (inherits(crop, "sfg")) {
        crop <- terra::vect(crop)
      }
    }

    # check whether the variables exist for this dataset
    if (dataset != "custom") { # if we are using standard datasets
      check_var_downloaded(bio_variables, dataset)
    } else { # else check that the variables exist in the custom nc
      check_var_in_nc(bio_variables, path_to_nc)
    }

    time_index <- NULL
    climate_spatrasters <- list()

    for (this_var in bio_variables) {
      # get name of file that contains this variable
      if (dataset != "custom") {
        var_meta <- get_var_meta(this_var, dataset)
        this_file <- var_meta$file_name
        this_file <- file.path(get_data_path(), this_file)
        this_var_orig <- get_varname(variable = this_var, dataset = dataset)
        this_var_longname <- var_meta$long_name
        this_var_units <- var_meta$units
      } else {
        this_file <- file.path(path_to_nc)
        this_var_orig <- this_var
        this_var_longname <- NULL
        this_var_units <- NULL
      }
      # figure out the time indeces the first time we run this
      if (is.null(time_index)) {
        # as we have the file name, we can us the same code for custom and
        # standard datasets.
        times <- get_time_bp_steps(dataset = "custom", path_to_nc = this_file)
        time_index <- time_bp_to_i_series(
          time_bp = time_bp,
          time_steps = times
        )
      }
      # retrieve time axis for virtual file
      var_brick <- pastclim_rast(
        x = this_file, bio_var_orig = this_var_orig,
        bio_var_pastclim = this_var, var_longname = this_var_longname,
        var_units = this_var_units
      )
      # subset to time steps
      if (!is.null(time_bp)) {
        var_brick <- terra::subset(var_brick, subset = time_index)
      }

      # subset extent
      if (!is.null(ext)) {
        var_brick <- terra::crop(var_brick, ext)
      }
      # subset to crop
      if (!is.null(crop)) {
        terra::crs(crop) <- terra::crs(var_brick)
        var_brick <- terra::mask(var_brick, crop)
        var_brick <- terra::crop(var_brick, crop)
      }

      # special treatment for biome variables
      if (this_var == "biome") {
        if (dataset != "custom") {
          # 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))
        } else {
          # for custom datasets, we just turn them into a factor
          var_brick <- terra::as.factor(var_brick)
        }
      }

      climate_spatrasters[[this_var]] <- var_brick
    }
    climate_sds <- terra::sds(climate_spatrasters)
    return(climate_sds)
  }