Search code examples
rrasterterra

Add layers of two raster stacks element wise and store in new stack


I have two raster stacks each with ca. 80000 layers. The stacks have the same extent, the same CRS, the same resolution and the same number of layers. I want to add the layers of the stacks element wise to each other (e.g. stack3[[n]] <- stack1[[n]]+stack2[[n]]). I tried the following approaches so far:

  1. Using a loop, which should work but takes way too long with over 80000 layers.
library(terra)
nlayers <- nlyr(stack1)
stack3 <- rast()
for (layer_index in 1:nlayers){
    res <- stack1[[layer_index]]+stack2[[layer_index]]
    stack3[[layer_index]] <- res
}
  1. Using simply '+'. This operation doesn't produce the desired result.
stack3 <- stack1+stack2
identical(stack3[[1]],stack1[[1]]+stack2[[1]])
[1] FALSE
  1. Using terra::app() throws an error message:
stack3 <- rast()
sum_layers <- function(layer1, layer2) {
   layer1 + layer2
}

stack3 <- app(stack1,stack2, sum_layers)

[1] Error in get(as.character(FUN), mode = "function", envir = envir) : 
  object 'fun' of mode 'function' was not found

Any help, tricks and suggestions appreciated!


Solution

  • The simplest approach is to use +.

    library(terra)
    r1 <- rast(system.file("ex/logo.tif", package="terra"))   
    r2 <- r1 / 2
    x <- r1 + r2
    

    This does work:

    y <- r1[[1]] + r2[[1]]
    all.equal(x[[1]], y)
    #[1] TRUE
    
    # no differences in the values
    plot(x[[1]] - y)
    

    An alternative approach would be

    s <- sds(r1, r2)
    a <- app(s, sum)
    

    I think the base::identical method sees a difference in the encapsulated C++ object (it is not the same pointer). You could use isTRUE(all.equal(x, y)) instead. To avoid the confusion, I have now also added a specialized identical<SpatRaster,SpatRaster> method to terra 1.7-67.