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)
