R柵格包中的交叉帶

Question

我的輸入柵格由多個圖層組成，每個圖層的圖像區域都沒有數據值。這些圖層不完全重疊，我試圖輸出一個只包含所有波段（任何圖層上沒有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) 

就如何有效地做到這一點更加任何想法/得到它與大量的圖層一起工作？

Answer 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) 

Answer 2

我第一次提出這樣的：

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

但我現在知道你想要什麼，不會讓你，因爲它會返回其中至少一層具有價值的區域（行/列）那不是NA;而不是所有的價值都不是NA的地區。

以下可能是有效的。具有至少一個NA值的所有單元格將與總和一起變爲NA（默認情況下，na.rm = FALSE）。

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

也許其次

r <- mask(r, x) 

設置所有單元格中的R以NA如果他們NA在X

Answer 3

我發現all/any是做邏輯和/或方法在一個長長的名單上。這裏展示了兩個完全相同的地塊，可以實現您在小型柵格上的目標。

#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)))