Home > Back-end >  How apply a function to list of SpatRaster using terra
How apply a function to list of SpatRaster using terra

Time:09-20

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

CodePudding user response:

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)
  • Related