我的输入栅格由多个图层组成,每个图层都有一个没有数据值包围的图像区域。这些图层不完全重叠,我正在尝试输出一个文件,该文件仅包含所有波段的交集(在任何图层上都没有 NoData 值的区域)。
以下适用于几层,但不适用于我在真实文件中拥有的 50 层以上(至少 3000x3000 像素):
library(raster)
fin = "D:\\temp\\all_modes.pix"
fout = "D:\\temp\\test.pix"
inbands = stack(fin, bands = c(3:20))
NAvalue(inbands) = 0
# Not great:
#out = all(is.na(inbands) == FALSE) * inbands
#writeRaster(out, filename=fout, format="PCIDSK", dtype="INT2U", overwrite=TRUE, NAflag=0)
# A little better:
#mymask = all(as.logical(inbands))
#mask(inbands, mymask, filename=fout, format="PCIDSK", dtype="INT2U", overwrite=TRUE, NAflag=0)
# Even better, don't need to keep everything (but still not efficient):
#trim(all(as.logical(inbands)) * inbands, filename=fout, format="PCIDSK", dtype="INT2U", overwrite=TRUE, NAflag=0)
# Even better, calculations get smaller as we progress (is it possible to do even better?)
for(i in 1:nlayers(inbands)){
band_i = subset(inbands, i)
inbands = trim(as.logical(band_i) * inbands)
}
writeRaster(inbands, filename=fout, format="PCIDSK", dtype="INT2U", overwrite=TRUE, NAflag=0)
关于如何更有效地做到这一点/让它与大量层一起工作的任何想法?