extract_covariates_var_time() not compatible w/ {terra} SpatRaster

[joshcullen]

I'm working to extract a set of values for my tracks from 4 raster covariates, where 3 of these must be time-matched. Due to the eventual ending of support for {raster}, I've transitioned much of my analyses to {terra}. While trying to use the extract_covariates_var_time() function, I've run into an issue caused by raster::getZ(). In {terra}, this functionality has been kept, but adjusted (as are many of the other functions from {raster}), via the terra::time() function instead.

Would it be possible to add an if-else statement within the extract_covariates_var_time() function based on whether the class of the object is of type 'RasterStack' or 'RasterBrick' for {raster} and 'SpatRaster' for {terra}? This would provide an automated way to extract covariate values without needing to change the class from a SpatRaster to a RasterStack or RasterBrick.

Also, since I am working with monthly raster data, would it be possible to add an argument that extracts covariate values based on the name of a specified column corresponding to the name/time of each raster layer (e.g., year-month) rather than the amount of time before or after the time of the raster layer?

[jmsigner]

Thanks for your feedback. I will transition amt to terra and sf within the next weeks, so hopefully this will solve issue 1. I have to give issue 2 some thought. E.g., you would like to extract from a raster called March if a column called month has the value march, right?

[joshcullen]

Thanks! And yes regarding the second topic.

I've written a custom function to do this for my own purposes, so I can add that here if you think it would be helpful as a place to start.

[jmsigner]

Thanks @joshcullen, that would be very helpful to get started.

[joshcullen]

Here are my internal and wrapper functions for extracting covariates in parallel across IDs. I've also created a reprex below to demonstrate. Hope that helps!

extract.covars.internal = function(data, layers, dyn_names, ind, p) {
  #data = data frame containing at least the id, coordinates (x_,y_), and indicator column (ind) to match to the dynamic raster layer(s) 
  #layers = a {terra} SpatRaster object containing environ covars
  #dyn_names = vector of names dynamic raster layers (in same order as layers); NULL by default
  #ind = character/integer. The name or column position of the indicator column of data to be
  #matched w/ names of a dynamic raster
  #p = a stored 'progressr' object for creating progress bar
  
  
    # Extract environ covariate values by points (CRS assumed to be the same across 'data' and 'layers')
    pts <- data %>%
      sf::st_as_sf(., coords = c('x_','y_'), crs = terra::crs(layers[[1]]))
    
    #Extract values from each line segment
    for (j in 1:n_distinct(pts[[ind]])) {
      
      #create subsetted data.frame of original given selected month.year
      data.sub <- data[data[[ind]] == unique(pts[[ind]])[j],]
      
      
      # Create time-matched raster stack
      time.ind <- purrr::map(layers[dyn_names], ~which(names(.x) == unique(pts[[ind]])[j]))
      layers.tmp <- layers
      
      layers.tmp[dyn_names] <- purrr::map2(layers.tmp[dyn_names], time.ind, ~{.x[[.y]]})
      layers.tmp <- rast(layers.tmp)
      
      
      tmp1 <- terra::extract(layers.tmp,
                            terra::vect(pts[pts[[ind]] == unique(pts[[ind]])[j],]),
                            along = FALSE, cells = FALSE)
      
      
      
      
      # add extracted covariates to original dataset
      tmp2 <- tmp1 %>%
        dplyr::select(-ID) %>%
        cbind(data.sub, .)
      
      
      
      if (j == 1) {
        extr.covar <- tmp2
      } else {
        extr.covar <- rbind(extr.covar, tmp2)
      }
    }
  
  p()  #plot progress bar
  return(extr.covar)
}



#--------------------------------------



extract.covars = function(data, layers, dyn_names = NULL, ind) {
  ## data must be a data frame with "id" column, coords labeled "x_" and "y_" ; optionally can have column that specifies dynamic covar names
  
  message("Prepping data for extraction...")
  tictoc::tic()
  
  dat.list <- split(data, data$id)
  
  ## Make raster data (stored in `layers`) usable in parallel
  .layers <- purrr::map(layers, terra::wrap)
  
  
  
  message("Extracting environmental values for IDs...")
  
  progressr::with_progress({
    #set up progress bar
    p <- progressr::progressor(steps = length(dat.list))
    
    pts <- furrr::future_map(dat.list,
                              ~extract.covars.internal(data = .x, layers = map(.layers, terra::rast),
                                                       dyn_names = dyn_names, ind = ind,
                                                       p = p),
                              .options = furrr_options(seed = TRUE))
  })
  
  pts <- dplyr::bind_rows(pts, .id = "id")
  
  tictoc::toc()
  
  
  return(pts)
}


#--------------------------------------


library(tidyverse)
library(lubridate)
library(amt)
library(terra)
library(furrr)
library(future)

data("amt_fisher")
data("amt_fisher_covar")

# Add indicator column for month
amt_fisher2 <- amt_fisher %>% 
  mutate(month = month.abb[month(t_)])
unique(amt_fisher2$month)  #December through March



# Convert from Raster to SpatRaster
amt_fisher_covar2 <- sapply(amt_fisher_covar, rast)

# Make sure that all covars have same extent; probably want to crop by specified extent in real analysis
for (var in c("landuse", "elevation")) {
  amt_fisher_covar2[[var]] <- terra::resample(amt_fisher_covar2[[var]],
                                              amt_fisher_covar2$popden,
                                              method = "average")
}

# Create object w/ 4 layers of landuse representing 4 different months
amt_fisher_covar2$landuse <- terra::app(amt_fisher_covar2$landuse, function(x) rep(x, 4))
names(amt_fisher_covar2$landuse) <- c('Feb','Mar','Dec','Jan')



### Extract time-matched covars by month ###

#for running in parallel; define number of cores to use
plan(multisession, workers = 4)
amt_fisher3 <- extract.covars(data = amt_fisher2, layers = amt_fisher_covar2, dyn_names = "landuse",
                       ind = "month")
#takes 4 sec to run
plan(sequential)


head(amt_fisher3)

[jmsigner]

@joshcullen the first question should now be addressed in the development version of amt, which does not rely on raster anymore. I will now work on your second question.