1

我的输入栅格由多个图层组成，每个图层都有一个没有数据值包围的图像区域。这些图层不完全重叠，我正在尝试输出一个文件，该文件仅包含所有波段的交集（在任何图层上都没有 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)

关于如何更有效地做到这一点/让它与大量层一起工作的任何想法？

4

3 回答 3

2

我首先提出这个：

x <- trim(inband, filename='fout')

但是我现在意识到这不会让你得到你想要的，因为它会返回至少一层的值不是 NA 的区域（行/列）；而不是所有值都不是 NA 的区域。

以下可能是有效的。具有至少一个 NA 值的所有单元格将变为 NA 和 sum（默认情况下，na.rm=FALSE）。

x <- sum(inband)
x <- trim(x)
r <- crop(inband, x)

也许紧随其后

r <- mask(r, x)

如果它们在 x 中为 NA，则将 r 中的所有单元格设置为 NA

于 2011-04-09T03:58:26.620 回答
1

感谢您的回答，他们给了我很好的想法。我想出了这个，这要快得多：

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)
于 2011-04-15T16:55:37.970 回答
0

我发现all/any是对一长串进行逻辑 AND/OR 的方法。这里展示了两个相同的图，它们可以在小栅格上实现您的目标。

#dummy data
m1 = cbind(matrix(NA, nrow = 4, ncol = 2), matrix(1, nrow = 4, ncol = 2))
m2 = t(m1)
m3 = matrix(rep(c(1, NA), 8), nrow = 4)
inbands = stack(lapply(list(m1, m2, m3), raster))

# first method, using & operator 
plot(inbands[[1]] & inbands[[2]] & inbands[[3]])

# second method, using `all`
plot(all(inbands))

在我的系统上，有将数字强制为逻辑的警告。以下两种方法可以避免警告但可能会更慢？它们在逻辑上是等价的，但您可以将它们相互对立，并与上面的第二种方法进行对比。

plot(all(!is.na(inbands)))
plot(!any(is.na(inbands)))
于 2011-04-09T09:24:02.503 回答