Пересечение полос в растровом пакете 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

Я впервые предложил это:

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 2

Спасибо за ответы, они дали мне хорошие идеи. Я придумал это, что гораздо быстрее:

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