Search code examples
rrasterterra

How apply a function to list of SpatRaster using terra


I have a list of raster files that I can read with terra. I would like to apply simple mathematical functions that result in one SpatRaster with similar dimensions to the input. The functions (i.e. from terra) median, max, and min don't work, while mean and stdev works perfectly. Here is a minimal example:

library(terra)
f <- rast(system.file("ex/elev.tif", package="terra"))
f1<-c(f,f);f2<-c(f,f) # I'm doing this because actual rasters has multiple layers. 

rlist<-list(c(f,f),f1,f2)

# calculate different ensemble statistics
KMean <- Reduce(terra::mean, rlist)
KMedian <- Reduce(terra::median, rlist)
Ksd <- Reduce(terra::stdev, rlist)

KMax <- Reduce(terra::max, rlist)
KMin <- Reduce(terra::min, rlist)

Here are the error that I get from running median, max, and min:

> rlist
[[1]]
class       : SpatRaster 
dimensions  : 90, 95, 2  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : elev.tif  
              elev.tif  
names       : elevation, elevation 
min values  :       141,       141 
max values  :       547,       547 

[[2]]
class       : SpatRaster 
dimensions  : 90, 95, 2  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : elev.tif  
              elev.tif  
names       : elevation, elevation 
min values  :       141,       141 
max values  :       547,       547 

[[3]]
class       : SpatRaster 
dimensions  : 90, 95, 2  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : elev.tif  
              elev.tif  
names       : elevation, elevation 
min values  :       141,       141 
max values  :       547,       547 

> KMedian <- Reduce(terra::median, rlist)
Error: [median] na.rm (the second argument) must be a logical value
> KMax <- Reduce(terra::max, rlist)
Error: 'max' is not an exported object from 'namespace:terra'
> KMin <- Reduce(terra::min, rlist)
Error: 'min' is not an exported object from 'namespace:terra'

Notes:

  • All raster have the same dimensions.

Here are the session info:

> sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=de_DE.UTF-8       
 [4] LC_COLLATE=en_GB.UTF-8     LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_GB.UTF-8   
 [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] terra_1.6-17

loaded via a namespace (and not attached):
[1] compiler_4.2.1   tools_4.2.1      Rcpp_1.0.9       codetools_0.2-18

Edit:

Using the function without the namespace results in another error.

KMin <- Reduce(min, rlist)
Error: [f] unknown function argument

Solution

  • Here is the more terra-idiomatic way to do that, with tapp

    Your example data

    library(terra)
    f <- rast(system.file("ex/elev.tif", package="terra"))
    ff <- c(f,f)
    rlist <- list(ff, ff, ff)
    

    Solution:

    r <- rast(rlist)
    idx <- rep(1:2, length(rlist))
    idx
    #[1] 1 2 1 2 1 2
    
    Kmean <- tapp(r, idx, mean)
    Kmedian <- tapp(r, idx, median)
    Ksd <- tapp(r, idx, sd)
    KMin <- tapp(r, idx, min)
    KMax <- tapp(r, idx, max)
    

    Alternatively, you can create a SpatRasterDataset and use app

    s <- sds(rlist)
    Kmean <- app(s, mean)
    #etc
    

    If you want to use Reduce you can replace this

    KMin <- Reduce(min, rlist)
    #Error: [f] unknown function argument
    

    With

    KMin <- Reduce(\(i, ...) min(i), rlist)
    

    Whereas

    KMean <- Reduce(mean, rlist)
    
    

    Just works. I suppose that is related to their differences in definitions, with min being a primitive function with ellipses as first argument.

    min
    #function (..., na.rm = FALSE)  .Primitive("min")
    
    mean
    #function (x, ...) 
    #standardGeneric("mean")
    

    Finally, it is also possible to get the parallel mean, min, etc. by just doing

    RMean <- mean(ff, ff, ff)
    RMin <- min(ff, ff, ff)
    

    And, hence, (with the work-around for min and max),

    RMean <- do.call(mean, rlist)
    RMin <- do.call(\(i, ...) min, rlist)