rterra

How to create maximum value composite using terra r package?


I want to create a maximum value composite raster by taking maximum values from each band of a raster stack. I am using the following code

library(terra)
library(RStoolbox)

rast1 <- rast(lsat)

#Create two raster stacks
rast2 <- rast(lsat)
random_nums <- runif(length(rast2), min = 0.3, max = 0.5) # a set of random numbers the size of the image
rast2[] <- rast2[] + random_nums 

rast3 <- rast(lsat)
random_nums <- runif(length(rast3), min = 0.1, max = 0.4) # a set of random numbers the size of the image
rast3[] <- rast3[] - random_nums 

#Create a stack of all the rasters
s <- c(rast1, rast2, rast3)

#Create a maximum value composite 
b1 <- tapp(s, index=c(1,8,15,2,9,16,3,10,17,4,11,18,5,12,19,6,13,20,7,14,21), 
           fun=max)
b1

Now it is giving me 21 layers. But the output should have 7 bands with each band made of by taking the maximum of rast1, rast2 and rast3 i.e. the band1 of b1 should take the maximum of band1 of rast1, 2 and 3. Likewise for b2 to b7.


Solution

  • Example data

    library(terra)
    r1 <- rast(ncol=10, nrow=10, nlyr=7)
    set.seed(1)
    r1 <- init(r1, runif)
    r2 <- init(r1, runif)
    r3 <- init(r1, runif)
    

    The most straightforward way to get the parallel maximum cell values (the maximum value for each cell and each layer, across SpatRasters):

    m1 <- max(r1, r2, r3)
    

    Or create a SpatRasterDataset and use app

    d <- sds(r1, r2, r3)
    m2 <- app(d, max)
    

    I would not recommend it as it is unnecessarily complex, but you can also do that with tapp like this:

    i <- rep(1:7, 3)
    i
    # [1] 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7
    s <- c(r1, r2, r3)
    m3 <- tapp(s, i, fun=max)
    
    m3
    #class       : SpatRaster 
    #dimensions  : 10, 10, 7  (nrow, ncol, nlyr)
    #resolution  : 36, 18  (x, y)
    #extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
    #coord. ref. : lon/lat WGS 84 
    #source      : memory 
    #names       :             X1,             X2,             X3,             X4,             X5,             X6, ... 
    #min values  :   1.808664e-01,   0.000000e+00,  2.803072e-309, -6.884876e-311,   0.000000e+00,   0.000000e+00, ... 
    #max values  :   9.926841e-01,  3.131513e-294,  1.536768e+277,  3.202245e+307,  3.046950e+294,  1.047578e+296, ... 
    

    Note that the layers with the same index are grouped. So for what you want you should only have seven numbers, each occurring three times.