多级栅格列表 - 在 R 中执行代数

Question

我有一个多级栅格列表，如何有效地对此类结构进行操作总是让我头疼。我想根据索引从列表第三层的栅格中取平均值。例如：

# make some dummy rasters
a <- raster(xmn=0,xmx=5,ymn=0,ymx=5,res=1)
a[] <- sample(1:5,25,replace=T)
b <- raster(xmn=0,xmx=5,ymn=0,ymx=5,res=1)
b[] <- sample(1:5,25,replace=T)
c <- raster(xmn=0,xmx=5,ymn=0,ymx=5,res=1)
c[] <- sample(1:5,25,replace=T)
d <- raster(xmn=0,xmx=5,ymn=0,ymx=5,res=1)
d[] <- sample(1:5,25,replace=T)
e <- raster(xmn=0,xmx=5,ymn=0,ymx=5,res=1)
e[] <- sample(1:5,25,replace=T)
f <- raster(xmn=0,xmx=5,ymn=0,ymx=5,res=1)
f[] <- sample(1:5,25,replace=T)
g <- raster(xmn=0,xmx=5,ymn=0,ymx=5,res=1)
g[] <- sample(1:5,25,replace=T)
h <- raster(xmn=0,xmx=5,ymn=0,ymx=5,res=1)
h[] <- sample(1:5,25,replace=T)
i <- raster(xmn=0,xmx=5,ymn=0,ymx=5,res=1)
i[] <- sample(1:5,25,replace=T)
j <- raster(xmn=0,xmx=5,ymn=0,ymx=5,res=1)
j[] <- sample(1:5,25,replace=T)
k <- raster(xmn=0,xmx=5,ymn=0,ymx=5,res=1)
k[] <- sample(1:5,25,replace=T)
l <- raster(xmn=0,xmx=5,ymn=0,ymx=5,res=1)
l[] <- sample(1:5,25,replace=T)

)

list_abcd <- list(a,b,c,d)
list_efgh <- list(e,f,g,h)
list_ijkl <- list(i,j,k,l)

list_all <- list(list_abcd,list_efgh,list_ijkl)

我想对栅格执行计算以获得列表中各个第一个栅格的平均值（a,e and i 的平均值），第二 (b,f and j)、第三 (c,g,k) 和第四 (d, h and l)。

结果将是一个包含 4 个栅格的列表。

实际上，我有一个这样结构的列表，包含 10 个组，每个组有 100 多个栅格，因此解决方案需要可扩展。

其他信息 我使用 NetCDF 文件，它们的结构如下所示（如下）。我创建了多级列表来处理它，但后来很头疼如何在没有太多费力的 for 循环的情况下很好地处理它。

3 variables:
    double var1[lon,lat,years,irr]
    double var2[lon,lat,years,irr]
    double var3[lon,lat,years,irr]
 4 dimensions:
    lon  Size:720
        units: degrees_east
        long_name: lon
    lat  Size:360
        units: degrees_north
        long_name: lat
    years  Size:100
        units: mapping
        long_name: years
    irr  Size:2
        units: mapping
        long_name: rainfed/irrigated

Answer 1

您的示例数据

library(raster)
r <- raster(xmn=0,xmx=5,ymn=0,ymx=5,res=1)
set.seed(123)
list_abcd <- replicate(4, setValues(r, sample(1:5,25,replace=T)))
list_efgh <- replicate(4, setValues(r, sample(1:5,25,replace=T)))
list_ijkl <- replicate(4, setValues(r, sample(1:5,25,replace=T)))   
list_all <- list(list_abcd,list_efgh,list_ijkl)
names(list_all) <- c("abcd", "efgh", "ijkl")

解决方案，创建单个 RasterStack 并使用 stackApply

x <- stack(unlist(list_all))
i <- rep(1:4, 3)
s <- stackApply(x, i, mean)
s
#class      : RasterBrick 
#dimensions : 5, 5, 25, 4  (nrow, ncol, ncell, nlayers)
#resolution : 1, 1  (x, y)
#extent     : 0, 5, 0, 5  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : memory
#names      :  index_1,  index_2,  index_3,  index_4 
#min values : 1.333333, 2.000000, 1.000000, 1.666667 
#max values : 4.666667, 5.000000, 4.666667, 4.333333

现在让我们用 terra（raster 的替换）来做这件事。我们从三个 SpatRaster 开始（类似于 RasterStack）

library(terra)
rr <- rast(xmin=0,xmax=5,ymin=0,ymax=5,res=1)
set.seed(123)
abcd <- replicate(4, setValues(rr, sample(1:5,25,replace=T))) |> rast()
efgh <- replicate(4, setValues(rr, sample(1:5,25,replace=T))) |> rast()
ijkl <- replicate(4, setValues(rr, sample(1:5,25,replace=T))) |> rast()
# the above is equivalent to: abcd <- c(a,b,c,d)

解决方案 1：合并为一个 SpatRaster 并使用 tapp

rall <- c(abcd, efgh, ijkl) 
z1 <- tapp(rall, rep(1:4, 3), "mean")

解决方案 2：组合成 SpatRasterDataset 并使用 app

rsd <- sds(abcd, efgh, ijkl) 
z2 <- app(rsd, "mean")

方案三：直接使用mean

z3 <- mean(abcd, efgh, ijkl)

稍后：从你的评论中我了解到你犯的错误确实是在创建这些列表时。如果我理解你的话，你有 n 个代表模型的文件，每个文件有 2*100 层；并且您想要跨模型平均每一层的值（年份和水状况）；这样您最终会得到一个包含 200 个图层的栅格数据集（或 2 个包含 100 个图层的数据集）。你可以通过这样的方式实现：

 library(terra)
 ff <- list.files(pattern="nc$")  
 sd <- sds(ff)
 x <- app(sd, mean)

或者，如果 irrigated/rainfed 是分开的 sub-datasets，比如

 library(terra)
 ff <- list.files(pattern="nc$")  
 sd <- lapply(ff, \(i) rast(i, "irrigated")) |> sds()
 # or sd <- lapply(ff, \(i) rast(i, 1)) |> sds()
 x <- app(sd, mean)

一些示例数据

f <- system.file("ex/logo.tif", package="terra")
sd <- sds(c(f,f,f))
x <- app(sd, mean)

多级栅格列表 - 在 R 中执行代数

Multi-level list of rasters - performing algebra in R

r

list

stack

lapply

raster