Open joshcullen opened 1 year ago
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?
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.
Thanks @joshcullen, that would be very helpful to get started.
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)
@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.
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 byraster::getZ()
. In {terra}, this functionality has been kept, but adjusted (as are many of the other functions from {raster}), via theterra::time()
function instead.Would it be possible to add an
if-else
statement within theextract_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?