Sea level

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: pastclim
 Type: Package
 Title: Manipulate Time Series of Climate Reconstructions 
-Version: 2.2.0.9003
+Version: 2.2.0.9004
 Authors@R: c(
        person("Michela", "Leonardi", role = "aut"),
        person(c("Emily","Y."), "Hallet", role = "ctb"),

NEWS.md

@@ -2,6 +2,7 @@
 * improve the speed of `location_series()`.
 * fix for missing time step in CHELSA-TraCE21k (-1300 BP is missing for min
   and max temperature
+* Add Clark2025 as a possible sea level dataset.
 
 # pastclim 2.2.0
 * Add functions to perform delta downscale of climate data.

R/get_sea_level.R

@@ -6,34 +6,43 @@
 #' year ago).
 #'
 #' @param time_bp the time of interest
+#' @param dataset the dataset to use, either "Spratt2016" or "Clark2025"
 #' @returns a vector of sea levels in meters from present level
 #'
 #' @keywords internal
 
 
-get_sea_level <- function(time_bp) {
-  # get sea level from Spratt 2016
-  sea_level_info <- utils::read.table(
-    system.file("extdata/sea_level_spratt2016.txt",
-      package = "pastclim"
-    ),
-    header = TRUE
-  )
-  time_calkaBP <- -time_bp / 1000 # nolint
-  if (any(time_calkaBP < 0)) {
-    stop("this function only supports times in the past")
+get_sea_level <- function(time_bp, dataset = "Spratt2016") {
+  dataset <- match.arg(dataset, c("Spratt2016", "Clark2025"))
+  if (any(time_bp > 0, na.rm = TRUE)) {
+    stop("time_bp should be in the past")
   }
-  if (any(time_calkaBP > 798)) {
-    stop("the dataset of sea level reconstructions stops at 798ky BP")
+  if (dataset == "Spratt2016") {
+    # check that time is not too old for the dataset
+    if (any(time_bp < -798000, na.rm = TRUE)) {
+      stop("Spratt2016 only reached -798,000 years BP")
+    }
+    sea_level_info <- spratt2016
+  } else if (dataset == "Clark2025") {
+    if (any(time_bp < -4882000, na.rm = TRUE)) {
+      stop("Clark2025 only reached -4,882,000 years BP")
+    }
+    sea_level_info <- clark2025
   }
+
+
   ## TODO this is not safe, we should be getting the closest values
   ## or even better interpolate
   sea_level <- stats::approx(
-    x = sea_level_info$age_calkaBP,
-    y = sea_level_info$SeaLev_longPC1,
-    xout = time_calkaBP
+    x = sea_level_info$time_bp,
+    y = sea_level_info$sea_level,
+    xout = time_bp
   )$y
   # rescale to have 0 for 0kBP
-  sea_level <- sea_level - sea_level_info$SeaLev_longPC1[1]
+  baseline_sea_level <- sea_level_info$sea_level[sea_level_info$time_bp == 0]
+  if (length(baseline_sea_level) != 1) {
+    stop("sea level dataset should include a single time_bp == 0 entry")
+  }
+  sea_level <- sea_level - baseline_sea_level
   return(sea_level)
 }

R/location_slice_from_region_series_fast.R

@@ -239,8 +239,10 @@ location_slice_from_region_series <- # nolint
     if ("biome" %in% bio_variables) {
       biome_levels <- levels(region_series$biome)[[1]]$category
       # get the levels from the numeric values
-      biome_numeric <- match(locations_data$biome,
-                             levels(region_series$biome)[[1]]$id)
+      biome_numeric <- match(
+        locations_data$biome,
+        levels(region_series$biome)[[1]]$id
+      )
       locations_data$biome <- factor(biome_levels[biome_numeric],
         levels = biome_levels
       )

R/make_land_mask.R

@@ -16,20 +16,36 @@
 #'
 #' @param relief_rast a [`terra::SpatRaster`] with relief
 #' @param time_bp the time of interest
-#' @param sea_level sea level at the time of interest (if left to NULL, this is
-#'   computed using Spratt 2016)
+#' @param sea_level sea level at the time of interest. It can be set to
+#' "Spratt2016" (the default) or "Clark2025" to automatically compute the
+#' level from one of those two datasets, or to a numeric vector of sea levels
+#' with the same length as `time_bp`. `NULL` is treated as "Spratt2016" for
+#' backwards compatibility.
 #' @returns a [`terra::SpatRaster`] of the land masks (with land as 1's and sea
 #'   as NAs), where the layers are different times
 #'
 #' @export
 
-make_land_mask <- function(relief_rast, time_bp, sea_level = NULL) {
+make_land_mask <- function(relief_rast, time_bp, sea_level = "Spratt2016") {
   if (is.null(sea_level)) {
-    sea_level <- get_sea_level(time_bp = time_bp)
-  } else { # check that we have as many sea level estimates as times
-    if (length(time_bp) != length(sea_level)) {
-      stop("time_bp and sea_level should have the same number of elements")
+    sea_level <- "Spratt2016"
+  }
+
+  # if sea_level is a character, check that it is either spratt or clark
+  if (is.character(sea_level)) {
+    if (!sea_level %in% c("Spratt2016", "Clark2025")) {
+      stop("sea_level should be either 'Spratt2016' or 'Clark2025'")
     }
+    sea_level <- get_sea_level(time_bp = time_bp, dataset = sea_level)
+  }
+
+  # now sea level should be numeric and the same length as time_bp
+  if (!is.numeric(sea_level)) {
+    stop("sea_level should be numeric")
+  }
+
+  if (length(time_bp) != length(sea_level)) {
+    stop("time_bp and sea_level should have the same number of elements")
   }
   land_mask <- NULL
   for (i in seq_along(time_bp)) {