CODEWITHSUNDEEP
rrasterterra
Jaken
Jaken

Reputation: 521

Raster app functions with multiple stacks and shared index

I have two raster stacks, both of which represent a sequence of months. Stack1 is a repeating sequence across several years (e.g., Jan, Feb, March, ..., Jan, Feb, March, ..., Jan, Feb, March, ...) and Stack2 is a single sequence of reference values for those same months (e.g., Jan, Feb, March, ...). I want to perform an operation combining the two stacks using their common index - e.g., Stack1/Stack2, but specific to the index (so Stack1Jan1/Stack2Jan, Stack1Feb1/Stack2Feb, ... Stack1Jan2/Stack2Jan).

I'm quite sure this is possible using one of the apply-based methods (app, tapp, rapp, etc.), but despite spending a lot of time on the documentation, I can't figure out which one. I know that using tapp, I can apply a function to different layers of the raster using an index, as in the example below, but I'm not sure how to make that same index apply when involving a second stack.

#sample data
r <- rast(ncols=9, nrows=9)
values(r) <- 1:ncell(r)
s <- c(r, r, r, r, r, r)
s <- s * 1:6
Stack1 <- c(s, s, s)
Stack2 <- rast(ncols=9, nrows=9, nlyrs=6)
values(Stack2) <- rep(1, ncell(Stack2)*6)
Stack2 <- Stack2 * 1:6

#tapp lets you apply an index for layers within a single stack
idx <- rep(1:6, 3)
s1_mean <- tapp(Stack1, idx, fun=mean)

# what I'd like would be something like:
myfun <- function(x) {x/Stack2[[idx]]}
s1_out <- tapp(Stack1, index=idx, fun=myfun)

I know I could do this by subsetting each raster into their months and then doing the calculation that way, but I'm sure there's a more efficient way!

Upvotes: 1

Views: 33

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47481

It seems that what you are looking for is

s1_out <- Stack1/Stack2

The above recycles the shorter object. That is expected in R, for example:

(1:12) / (1:3)
# [1]  1.0  1.0  1.0  4.0  2.5  2.0  7.0  4.0  3.0 10.0  5.5  4.0

It is possible to do this with an *app function like this

out <- lapp(sds(Stack1, Stack2), \(x, y) x/y, recycle=TRUE)

In more complex cases you could match with the indices. In this case that would be

out <- Stack1[[1:18]]/Stack2[[1:6]]
# or even 
out <- Stack1[[1:18]]/Stack2[[rep(1:6, 3)]]

But I wouldn't use these less efficient alternatives if you don't need them.

Upvotes: 1

Related Questions