遍历光栅砖
Looping through a raster brick
我知道这里有大量关于遍历光栅块的问题,但其中 none 提供了我正在寻找的 answer/advice。
我有一个大的(17.2GB,7901 层)netcdf
文件,我已将其作为 RasterBrick
导入到 R
。
> KK10Brick
class : RasterBrick
dimensions : 2160, 4320, 9331200, 7901 (nrow, ncol, ncell, nlayers)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : D:\LandUse\KK10.nc
names : X8000, X7999, X7998, X7997, X7996, X7995, X7994, X7993, X7992, X7991, X7990, X7989, X7988, X7987, X7986, ...
z-value : 100, 8000 (min, max)
varname : land_use
文件中的每一层代表 1 年,我需要为砖块中的每个像素创建时间移动平均值。尽管变量看起来是绝对的 (land_use
),但它实际上是一个百分比覆盖率。
我想创建一个 30 年移动平均线,10 年滑动 window。例如第一个 window 将生成图层 1:30
的平均值栅格,下一个 window 将生成另一个图层 11:40
...[=19 的平均值栅格=].
我当时认为 for 循环可能是完成此操作的最佳方式,但我不确定我是否走上了正确的道路,例如
for (i in 1:7901){
subsetLayers <- code to subset relevant layers
out <- stackApply(KK10Brick, indices = subsetLayers, fun = "mean", na.rm = TRUE, filename = paste("./Output/", "meanLU_window_", i, ".tif", sep = ""))
rm(out)}
我卡住的地方是编写代码来为 subsetLayers
生成序列。任何帮助将不胜感激。
编辑。
library(raster)
exBrick <- brick(nrow = 180, ncol = 360, nl = 100)
values(exBrick) <- runif(ncell(exBrick) * nlayers(exBrick))
crs(exBrick) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
exBrick
这应该适用于您的示例数据。我不确定它在速度和 RAM 使用方面如何扩展到您的非常大的 netcdf 数据 - 如果它适用于大数据,请告诉我。
starts = seq(1, nlayers(exBrick)-30, 10)
nout = length(starts)
out = brick(nrow = 180, ncol = 360, nl = nout)
values(out) = NA
crs(out) = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
for (i in 1:nout) {
start = starts[i]
out[[i]] = mean(exBrick[[start:(start+30)]])
}
如果 RAM 使用是分配大砖块来存储结果的限制因素,我们可以通过将每个结果图层保存到磁盘,一次一个栅格,以牺牲一些速度为代价来节省 RAM:
for (i in starts) {
out = mean(exBrick[[i:(i+30)]])
writeRaster(out, filename=paste0("out",i,".grd"), overwrite=TRUE)
}
我知道这里有大量关于遍历光栅块的问题,但其中 none 提供了我正在寻找的 answer/advice。
我有一个大的(17.2GB,7901 层)netcdf
文件,我已将其作为 RasterBrick
导入到 R
。
> KK10Brick
class : RasterBrick
dimensions : 2160, 4320, 9331200, 7901 (nrow, ncol, ncell, nlayers)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : D:\LandUse\KK10.nc
names : X8000, X7999, X7998, X7997, X7996, X7995, X7994, X7993, X7992, X7991, X7990, X7989, X7988, X7987, X7986, ...
z-value : 100, 8000 (min, max)
varname : land_use
文件中的每一层代表 1 年,我需要为砖块中的每个像素创建时间移动平均值。尽管变量看起来是绝对的 (land_use
),但它实际上是一个百分比覆盖率。
我想创建一个 30 年移动平均线,10 年滑动 window。例如第一个 window 将生成图层 1:30
的平均值栅格,下一个 window 将生成另一个图层 11:40
...[=19 的平均值栅格=].
我当时认为 for 循环可能是完成此操作的最佳方式,但我不确定我是否走上了正确的道路,例如
for (i in 1:7901){
subsetLayers <- code to subset relevant layers
out <- stackApply(KK10Brick, indices = subsetLayers, fun = "mean", na.rm = TRUE, filename = paste("./Output/", "meanLU_window_", i, ".tif", sep = ""))
rm(out)}
我卡住的地方是编写代码来为 subsetLayers
生成序列。任何帮助将不胜感激。
编辑。
library(raster)
exBrick <- brick(nrow = 180, ncol = 360, nl = 100)
values(exBrick) <- runif(ncell(exBrick) * nlayers(exBrick))
crs(exBrick) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
exBrick
这应该适用于您的示例数据。我不确定它在速度和 RAM 使用方面如何扩展到您的非常大的 netcdf 数据 - 如果它适用于大数据,请告诉我。
starts = seq(1, nlayers(exBrick)-30, 10)
nout = length(starts)
out = brick(nrow = 180, ncol = 360, nl = nout)
values(out) = NA
crs(out) = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
for (i in 1:nout) {
start = starts[i]
out[[i]] = mean(exBrick[[start:(start+30)]])
}
如果 RAM 使用是分配大砖块来存储结果的限制因素,我们可以通过将每个结果图层保存到磁盘,一次一个栅格,以牺牲一些速度为代价来节省 RAM:
for (i in starts) {
out = mean(exBrick[[i:(i+30)]])
writeRaster(out, filename=paste0("out",i,".grd"), overwrite=TRUE)
}