多层射手列表 - 在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和i的平均值）， 2nd（b，f和j），第三（c，g，k）和第4（d，h和h和l）。

结果将是包含4个射手的列表。

实际上，我有一个针对10个组的列表，每个列表都有100多个raster，因此解决方案需要可扩展。

---

附加信息我使用NetCDF文件及其结构的关注（如下）。我创建多层列表来处理它，但是如果没有太多费力的循环，那么如何处理它如何处理它。

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 s（类似于栅格挡板）

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")

解决方案3：使用平均直接

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)

或者，如果灌溉/雨水是单独的子数据集，则类似于

 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)