Faster extract

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: pastclim
 Type: Package
 Title: Manipulate Time Series of Climate Reconstructions 
-Version: 2.2.0.9000
+Version: 2.2.0.9001
 Authors@R: c(
        person("Michela", "Leonardi", role = "aut"),
        person(c("Emily","Y."), "Hallet", role = "ctb"),
@@ -25,7 +25,7 @@ BugReports: https://github.com/EvolEcolGroup/pastclim/issues
 Encoding: UTF-8
 LazyData: true
 LazyDataCompression: xz
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Roxygen: list(markdown = TRUE)
 Depends: 
   R (>= 4.0.0),

R/bathy_to_spatraster.R

@@ -1,9 +1,9 @@
 #' Cast `bathy` to `SpatRaster`
 #'
-#' This function converts a [`marmap::bathy`][`marmap::as.bathy()`] object to
+#' This function converts a `marmap::bathy` object to
 #' a [`terra::SpatRaster`].
 #'
-#' @param bathy a [`marmap::bathy`][`marmap::as.bathy()`] to convert
+#' @param bathy a `marmap::bathy` to convert
 #' @returns a [`terra::SpatRaster`] with the relief for the chosen region
 #'
 #' @keywords internal

R/download_etopo_subset.R

@@ -8,7 +8,7 @@
 #' fly from the NOAA server. If you plan to use the ETOPO2022 dataset
 #' extensively, it is worthwhile downloading it permanently to your computer
 #' with [download_etopo()], but beware that it is a large file (>1Gb). This
-#' function uses [marmap::getNOAA.bathy()] to download the data, and then
+#' function uses `marmap::getNOAA.bathy()` to download the data, and then
 #' converts them into a [`terra::SpatRaster`] formatted to be compatible with
 #' `pastclim`. NOTE: this function does not save the relief, it returns a
 #' [`terra::SpatRaster`]. If you plan to reuse this relief multiple times, it
@@ -17,7 +17,7 @@
 #' @param rast_template a [`terra::SpatRaster`] providing the extent and
 #'   resolution to be downloaded. This raster needs to have identical vertical
 #'   and horizontal resolution, and standard lat/long projection.
-#' @param ... additional parameters to be passed to [marmap::getNOAA.bathy()] to
+#' @param ... additional parameters to be passed to `marmap::getNOAA.bathy()` to
 #'   customise how files are stored. See the manpage for that function for
 #'   details
 #' @returns a [`terra::SpatRaster`] with the relief for the chosen region

R/download_paleoclim.R

@@ -47,7 +47,7 @@ download_paleoclim <- function(dataset, bio_var, filename = NULL) {
   paleoclim_path[1] <- file.path(paleoclim_path[1], resolution)
   paleoclim_path[8] <- file.path(paleoclim_path[8], resolution)
   # create a vrt for each variable
-  for (i in seq_len(length(band_vector))) {
+  for (i in seq_along(band_vector)) {
     # build the vsizip paths
     paleoclim_vsizip <- paste0("/vsizip/", file.path(
       paleoclim_path,

R/download_worldclim_future.R

@@ -78,7 +78,7 @@ download_worldclim_future <- function(dataset, bio_var, filename = NULL) {
 
 
   # create a vrt for each variable
-  for (i in seq_len(length(band_vector))) {
+  for (i in seq_along(band_vector)) {
     vrt_path <- file.path(get_data_path(), paste0(
       dataset, "_",
       band_vector[i], "_v",

R/download_worldclim_present.R

@@ -64,7 +64,7 @@ download_worldclim_present <- function(dataset, bio_var, filename) {
   }
 
   # create a vrt for each variable
-  for (i in seq_len(length(band_vector))) {
+  for (i in seq_along(band_vector)) {
     # build the vsizip paths
     if (!grepl("altitude", bio_var)) {
       worldclim_vsizip <- paste0("/vsizip/", file.path(

R/get_biome_classes.R

@@ -34,7 +34,7 @@ get_biome_classes <- function(dataset) {
     biomes_string <- biomes_string[-length(biomes_string)]
     biomes_string <- substr(biomes_string, 4, nchar(biomes_string))
     biome_categories <- data.frame(
-      id = seq_len(length(biomes_string)),
+      id = seq_along(biomes_string),
       category = biomes_string
     )
   }

R/get_land_mask.R

@@ -68,7 +68,6 @@ get_land_mask <- function(time_bp = NULL, time_ce = NULL, dataset) {
   ), terra::nlyr(land_mask))
 
 
-
   if (is.null(time_ce)) {
     names(land_mask) <- paste("land_mask", time_bp(land_mask), sep = "_")
   } else {

R/koeppen_geiger.R

@@ -251,7 +251,6 @@ methods::setMethod(
 )
 
 
-
 #' @param filename filename to save the raster (optional).
 #' @rdname koeppen_geiger-methods
 #' @export

R/location_series.R

@@ -65,15 +65,21 @@ location_series <-
 
     check_dataset_path(dataset = dataset, path_to_nc = path_to_nc)
 
-    # if we are using standard datasets, check whether a variables exists
-    # and get the times
-    if (dataset != "custom") {
-      check_var_downloaded(bio_variables, dataset)
-      times <- get_time_bp_steps(dataset = dataset, path_to_nc = path_to_nc)
-    } else { # else check that the variables exist in the custom nc
-      check_var_in_nc(bio_variables, path_to_nc)
-      times <- get_time_bp_steps(dataset = "custom", path_to_nc = path_to_nc)
+    # get the region series for this dataset
+    climate_brick <- region_series(
+      bio_variables = bio_variables,
+      dataset = dataset,
+      path_to_nc = path_to_nc
+    )
+
+    # get all available times
+    times <- time_bp(climate_brick)
+
+    # if time_bp is NULL, get all times from the region series
+    if (is.null(time_bp)) {
+      time_bp <- time_bp(climate_brick)
     }
+
     time_bp_i <- time_bp_to_i_series(
       time_bp = time_bp,
       time_steps = times
@@ -105,21 +111,30 @@ location_series <-
 
     # now copy over the times to match the coordinates
     time_bp <- rep(time_bp, each = n_loc)
-    # and now feed the info to location_slice
-    location_ts <- location_slice(
-      x = x, time_bp = time_bp, coords = coords, bio_variables = bio_variables,
-      dataset = dataset, path_to_nc = path_to_nc,
-      nn_interpol = nn_interpol, buffer = buffer,
+
+    # now simply wrap around location_slice_from_region_series
+    location_ts <- location_slice_from_region_series(
+      x = x,
+      time_bp = time_bp,
+      time_ce = NULL,
+      coords = coords,
+      region_series = climate_brick,
+      nn_interpol = nn_interpol,
+      buffer = buffer,
       directions = directions
     )
 
     # TODO if we had time_ce, we should convert back from time_bp
+    if (!is.null(time_ce)) {
+      location_ts$time_ce <- location_ts$time_bp + 1950
+      # remove the time_bp column
+      location_ts <- location_ts[, !names(location_ts) %in% "time_bp"]
+    }
 
     return(location_ts[, !names(location_ts) %in% "time_bp_slice"])
   }
 
 
-
 #' Extract a time series of bioclimatic variables for one or more locations.
 #'
 #' Deprecated version of [location_series()]

R/location_slice_from_region_series.R → R/location_slice_from_region_series_fast.R

@@ -55,10 +55,10 @@ location_slice_from_region_series <- # nolint
     bio_variables <- names(region_series)
 
     time_bp <- check_time_vars(time_bp = time_bp, time_ce = time_ce)
-    # boolean whether we will need to readd time_ce instead of time_bp
-    readd_ce <- FALSE
+    # boolean whether we will need to read time_ce instead of time_bp
+    read_ce <- FALSE
     if (any(!is.null(time_ce), "time_ce" %in% names(x))) {
-      readd_ce <- TRUE
+      read_ce <- TRUE
     }
 
     # if we have a data.frame
@@ -123,78 +123,97 @@ location_slice_from_region_series <- # nolint
       time_bp = locations_data$time_bp, time_steps = times
     )
     locations_data$time_bp_slice <- times[time_indeces]
-    unique_times <- unique(locations_data$time_bp_slice)
 
-    for (i_time in unique_times) {
-      this_slice <- slice_region_series(climate_brick,
-        time_bp = i_time
+    # if we are not using a buffer, we try to get all the values from
+    # specific locations directly
+    if (!buffer) {
+      locations_climate <- terra::extract(climate_brick,
+        y = locations_data[coords],
+        layer = time_indeces
       )
+      # for each bio_variable (the list element), extract the value  column from
+      # the data.frame and assign to locations_data
+      for (var in bio_variables) {
+        locations_data[[var]] <- locations_climate[[var]]$value
+      }
+    } else {
+      # if we are using a buffer, set to NA as we will compute them later
+      for (var in bio_variables) {
+        locations_data[[var]] <- NA
+      }
+    }
 
-      this_slice_indeces <- which(locations_data$time_bp_slice == i_time)
-      if (!buffer) { # get the specific values for those locations
-        this_climate <- terra::extract(
-          x = this_slice,
-          y = locations_data[locations_data$time_bp_slice == i_time, coords]
+    # if we interpolate or use a buffer, we have to find the incomplete cases
+    if (nn_interpol || buffer) {
+      loc_id_to_move <- which(!stats::complete.cases(locations_data))
+      # only do something if we have some locations for which we have no data
+      if (length(loc_id_to_move) != 0) {
+        # for each id, get the appropriate list of neighbours
+        rast_id_to_move <- terra::cellFromXY(
+          climate_brick[[1]],
+          as.matrix(coords_df[loc_id_to_move, ])
         )
-        # factors don't behave nicely when adding new elements, cast to
-        # character
-        if ("biome" %in% names(this_climate)) {
-          this_climate$biome <- as.character(this_climate$biome)
-        }
-        # sort out the indexing here
-        locations_data[locations_data$time_bp_slice == i_time, bio_variables] <-
-          this_climate[
-            ,
-            bio_variables
-          ]
-      } else { # set to NA as we will compute them with a buffer
-        locations_data[this_slice_indeces, ] <- NA
-      }
+        neighbours_list <- terra::adjacent(
+          climate_brick[[1]],
+          rast_id_to_move,
+          directions = directions,
+          pairs = FALSE
+        )
+        # convert from matrix (one row per location) to data.frame where first
+        # column is focal location, second column is id of each neighbour
+        neighbours_list <- data.frame(
+          focal_id = rep(as.numeric(rownames(neighbours_list)),
+            each = ncol(neighbours_list)
+          ),
+          neighbour_id = as.vector(t(neighbours_list)),
+          layer = rep(time_indeces[loc_id_to_move],
+            each = ncol(neighbours_list)
+          ),
+          unique_id = rep(seq_len(nrow(neighbours_list)),
+            each = ncol(neighbours_list)
+          )
+        )
+        # there is a BUG terra does not seem to cope with using layer when
+        # extracting from sds if y is cellID, y has to be a data.frame
+        neighbours_coords <- as.data.frame(terra::xyFromCell(
+          climate_brick[[1]],
+          neighbours_list$neighbour_id
+        ))
+        # extract the climate for all neighbours
+        neighbours_values <-
+          terra::extract(
+            x = climate_brick,
+            y = neighbours_coords,
+            layer = neighbours_list$layer
+          )
 
-      if (nn_interpol || buffer) {
-        locations_to_move <- this_slice_indeces[
-          this_slice_indeces %in%
-            which(!stats::complete.cases(locations_data))
-        ]
-        if (length(locations_to_move) == 0) {
-          next
-        }
-        for (i in locations_to_move) {
-          if (inherits(x, "data.frame")) {
-            cell_id <-
-              terra::cellFromXY(this_slice, as.matrix(coords_df[
-                i,
-              ]))
+        # for each variable, compute the mean across neighbours
+        # and replace the vales in locations_data
+        for (i_var in bio_variables) {
+          if (i_var == "biome") {
+            # for factors, compute the mode across neighbours
+            neighbours_mode <- tapply(
+              neighbours_values[[i_var]]$value,
+              neighbours_list$unique_id,
+              mode
+            )
+            locations_data[loc_id_to_move, i_var] <-
+              neighbours_mode
           } else {
-            cell_id <- coords_df[i]
-          }
-          neighbours_ids <-
-            terra::adjacent(this_slice, cell_id,
-              directions = directions, pairs = FALSE
+            neighbours_mean <- tapply(
+              neighbours_values[[i_var]]$value,
+              neighbours_list$unique_id,
+              mean,
+              na.rm = TRUE
             )
-
-          neighbours_values <-
-            terra::extract(
-              x = this_slice,
-              y = neighbours_ids[1, ]
-            ) # [, bio_variables]
-
-          neighbours_values_mean <- colMeans(
-            neighbours_values[, !names(neighbours_values) %in% "biome",
-              drop = FALSE
-            ],
-            na.rm = TRUE
-          )
-          if ("biome" %in% bio_variables) {
-            neighbours_values_mean["biome"] <-
-              mode(as.character(neighbours_values[, "biome"]))
+            locations_data[loc_id_to_move, i_var] <-
+              neighbours_mean
           }
-          locations_data[i, bio_variables] <-
-            neighbours_values_mean[bio_variables]
         }
       }
     }
-    # is.nan has not method for a data.frame
+
+    # is.nan has no method for a data.frame
     # nolint start
     is.nan.data.frame <- function(x) {
       do.call(cbind, lapply(x, is.nan))
@@ -205,7 +224,7 @@ location_slice_from_region_series <- # nolint
 
     locations_data <- locations_data[order(orig_id), ]
 
-    if (readd_ce) {
+    if (read_ce) {
       locations_data$time_ce <- locations_data$time_bp + 1950
       locations_data$time_ce_slice <- locations_data$time_bp_slice + 1950
       locations_data <- locations_data[

R/region_series.R

@@ -104,22 +104,26 @@ region_series <-
         this_var_longname <- NULL
         this_var_units <- NULL
       }
+
+      # 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
+      )
+
       # 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)
+        # get times from the var_brick
+        times <- time_bp(var_brick)
+        # convert it to indeces
         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)

R/sample_region_series.R

@@ -130,7 +130,7 @@ sample_rs_variable <- function(x, size, method = "random", replace = FALSE,
   # create list to store samples for each time step
   sample_list <- list()
   t_steps <- time_bp(x[1])
-  for (i in seq_len(length(size))) {
+  for (i in seq_along(size)) {
     if (size[i] > 0) {
       x_step <- slice_region_series(x, t_steps[i])
       # spatSample samples additional points to make sure it has enough points