R栅格数据包中的带交集 [英] Intersection of bands in R raster package

问题描述

我的输入栅格包含多个图层，每个图层的图像区域均不包含任何数据值.这些层没有完全重叠，我正在尝试输出仅由所有波段的交点(在任何层上没有NoData值的区域)组成的文件.

My input raster consists of multiple layers, each with an image area surrounded by no data values. These layers don't overlap completely, and I am trying to output a file which consists of only the intersection of all bands (the zone which has no NoData values on any layer).

以下内容适用于几层，但不适用于我的真实文件中的50层以上(至少3000x3000像素):

The following works for a few layers, but not for the 50+ that I have in my real files (at least 3000x3000 pixels):

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)

关于如何更有效地执行此操作/使其在大量图层上工作的任何想法?

Any ideas on how to do this more efficiently / get it working with a large number of layers?

推荐答案

感谢您的回答，他们给了我很好的主意.我想出了这个，速度更快:

Thanks for the answers, they gave me good ideas. I came up with this, which is much faster:

myraster = stack(fn, bands) # You get the idea
NAvalue(myraster) = 0

# Tranform to 1 where there is data    
logical_raster = as.logical(myraster)

# Make a raster with 1 in the zone of intersection
a = subset(myraster, 1)
values(a) = TRUE
for(i in 1:nlayers(myraster)) {
  a = a & logical_raster[[i]]
}

# Apply the "mask" and trim to intersection extent
myraster = myraster * a
intersect_only = trim(myraster)

这篇关于R栅格数据包中的带交集的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！