Vectorizing function which calls a dataframe for use in raster::overlay()

Question

I have a function which alters a RasterLayer given conditions specified in a dataframe. I would like to call this function using the overlay() function to alter a RasterStack but am having trouble vectorizing my function.

Here is an example of what I am doing with dummy data:

library(raster)
library(dplyr)

#generate dummy dataframe
df <- data.frame(matrix(nrow=3,ncol=43))
colnames(df) <- c("date_start", "date_end", "pixel_id")
df$date_start <- seq(as.Date("2024-01-01"), as.Date("2024-01-03"), by="+1 day")
df$date_end <- seq(as.Date("2024-02-01"), as.Date("2024-02-03"), by="+1 day")
df$pixel_id <- c(1,2,3)

#generate raster with pixel ids
id_raster <- raster(nrow=3, ncol=3, xmn=-123, xmx=-121, ymn=43, ymx=45) %>% setValues(seq(1,9))

#make a raster stack for editing
stackable <- id_raster / id_raster # make one blank raster
daylist <- seq(as.Date("2024-01-01"), as.Date("2024-02-03"), by="+1 day")
stack <- stack(replicate(length(daylist), stackable)) %>% setZ(daylist) # make a stack of rasters with one layer for every day of interest

#set function to apply to each layer
fun <- function(layer, id_raster,...) {
    day <- getZ(layer)
    day_data <- df %>% filter(date_start <= as.Date(day)) %>% filter(date_end >= as.Date(day))
    new_layer <- Which((id_raster %in% day_data$pixel_id), na.rm=F, cells=F)
    new_layer <- new_layer + layer
    return(new_layer)}

#use overlay to apply my function to each layer of the stack
new_stack <- overlay(stack, id_raster, fun=Vectorize(fun))

When I run this, the final line returns the error: "cannot use this formula, probably because it is not vectorized." So I assume I am using the Vectorize() function incorrectly somehow.

This question gets close to mine but the function calls on a list while mine calls on a dataframe. Somehow this must make a difference because simply using the Vectorize() call isn't working for me like it did for them.

I have tried adding a vectorize.args argument to the Vectorize() call using various syntaxes but I either still get the same error or my code seems to run endlessly (which shouldn't happen with such a small dataset).

Answer 1

I do not think this can be easily done with raster::overlay. I show a more direct way using "terra", as that package has replaced "raster".

Your example data

df <- data.frame(matrix(nrow=3,ncol=3))
colnames(df) <- c("date_start", "date_end", "pixel_id")
df$date_start <- seq(as.Date("2024-01-01"), as.Date("2024-01-03"), by="+1 day")
df$date_end <- seq(as.Date("2024-02-01"), as.Date("2024-02-03"), by="+1 day")
df$pixel_id <- c(1,2,3)

library(terra)
id_raster <- rast(nrow=3, ncol=3, xmin=-123, xmax=-121, ymin=43, ymax=45, vals=1:9)

s <- rast(id_raster, nlyr=34, vals=1)
time(s) <- seq(as.Date("2024-01-01"), as.Date("2024-02-03"), by="+1 day")

Solution

days <- time(s)
out <- lapply(1:nlyr(s), \(i) {
            d <- days[i]
            ids <- df$pixel_id[df$date_start <= d & df$date_end >= d]
            if (length(ids) > 0) {
                m <- subst(id_raster, ids, 1, NA)
                mask(s[[i]], m)
            } else { 
                NULL 
            }
        })

x <- sum(rast(out), na.rm=TRUE)