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)