我有一個多層次的柵格串列,它總是讓我頭疼如何有效地操作這些結構。我想根據索引從串列的第 3 級的柵格中獲取平均值。例如:
# 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
uj5u.com熱心網友回復:
您的示例資料
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)來做這件事。我們從三個SpatRasters 開始(類似于 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")
解決方案3: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 是單獨的子資料集,例如
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)
轉載請註明出處,本文鏈接:https://www.uj5u.com/shujuku/481197.html
上一篇:從包含重復串列值的字典中獲取鍵,然后將其對應的鍵存盤在嵌套串列中
下一篇:無法洗掉字串串列的特定值
