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,36 @@
 #' 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")
-  }
-  if (any(time_calkaBP > 798)) {
-    stop("the dataset of sea level reconstructions stops at 798ky BP")
+get_sea_level <- function(time_bp, dataset = "spratt2016") {
+  dataset <- match.arg(dataset, c("spratt2016", "clark2025"))
+  if (dataset == "spratt2016") {
+    # check that time is not too old for the dataset
+    if (any(time_bp < -798000)) {
+      stop("spratt2016 only reached -798,000 years BP")
+    }
+    sea_level_info <- spratt2016
+  } else if (dataset == "clark2025") {
+    if (any(time_bp < -4882000)) {
+      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]
+  sea_level <- sea_level - sea_level_info$sea_level[sea_level_info$time_bp == 0]
   return(sea_level)
 }

R/make_land_mask.R

@@ -16,20 +16,30 @@
 #'
 #' @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.
 #' @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) {
-  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")
+make_land_mask <- function(relief_rast, time_bp, 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)) {