从平均值创建新的光栅图像

Creating a new raster image from averages

我有一长串光栅,我需要取一定数量的平均值并创建一个新图像。例如,如果我有栅格 r1 r2 r3 r4 r5 r6 r7 r8 我想取 r1 和 r2 的平均值来给我一个图像,比方说 new1。然后我想要 r3 和 r4 的平均值给我图像 new2。我尝试在 caTools 中使用 runmean,但无法获得所需的输出。如果我有 8 个光栅图像,那么使用两个 window 应该给我留下四个光栅图像。我知道栅格通常属于 GIS 站点,但我需要代码方面的帮助,所以我希望它在这里没问题。

假设您将所有栅格都放在一个文件夹中:rasdir(此文件夹中除了要循环的栅格之外没有其他任何东西),设置环境变量:

rasdir="myrasters/"
raspaths <- list.files(path=rasdir, full.names=T)

假设所有栅格具有相同的范围和分辨率,它们可以堆叠:

rascube <- stack(raspaths)

创建执行某些功能的函数,例如跨频带均值

rascube 是要循环的图像堆栈,win window 大小,outdir 是输出目录

rasfun <- function(x=rascube, win=2, outdir=getwd()){

#Sanity check
if(!(length(raspaths)/win)%%1==0){stop("Number of rasters must be divisible by window size")}

#Create ```mat``` , an index that determines how rasters in ```rascube``` are aggregated: 
#indices in the same row refer to rasters to be averaged into the ith output.

mat <- matrix(data=1:length(raspaths), ncol=win, byrow=T)

#Loop over ```rascube```, calculating the moving average as controlled by ```mat```

for (i in 1:nrow(mat)){

#Compute ith moving mean, You can alter this to compute a moving "whatever you like"
#Note the usage of ```[[ ]]``` to subset raster bands: see ```raster``` docu.
#Also Note the usage of ```na.rm=T```, just in case your images have NA's you dont care about

res_i <- sum(x[[ mat[i,1]:mat[i,win] ]], na.rm=T)/win #

#Write output to file: note how output filename is derived from the respective input rasters
#makes it possible to trace the outputs back to their source rasters.

writeRaster(x=res_i, filename=paste(mat[i,1:win], collapse=""),
format="GTiff", overwrite=T)

}
}

#Run newly created function on your stack of inputs with whatever args:
rasfun(x=rascube, win=2, outdir="moving_mean_rasters/")

注意:栅格的数量必须能被 window 大小整除,例如试图在 window 大小为 2 的 7 个栅格上进行 运行 移动 window 计算将被健全性检查失败。当然,您可以更改函数以按照您认为最适合您的用例的方式运行。干杯!