包裹:
数据:
- 具有 10 个波段的 rasterStack。
- 每个波段都包含一个由 NA 包围的图像区域
- 波段是合乎逻辑的,即图像数据为“1”,周围区域为“0”/NA
- 每个波段的“图像区域”彼此不完全对齐,尽管大多数有部分重叠
客观的:
- 编写一个快速函数,该函数可以返回每个“区域”的 rasterLayer 或像元编号,例如,仅包含带 1 和 2 数据的像素位于区域 1,仅包含带 3 和 4 数据的像素位于区域 2等。如果返回 rasterLayer,我需要稍后能够将区域值与波段编号匹配。
第一次尝试:
# Possible band combinations
values = integer(0)
for(i in 1:nlayers(myraster)){
combs = combn(1:nlayers(myraster), i)
for(j in 1:ncol(combs)){
values = c(values, list(combs[,j]))
}
}
# Define the zone finding function
find_zones = function(bands){
# The intersection of the bands of interest
a = subset(myraster, 1)
values(a) = TRUE
for(i in bands){
a = a & myraster[[i]]
}
# Union of the remaining bands
b = subset(myraster, 1)
values(b) = FALSE
for(i in seq(1:nlayers(myraster))[-bands]){
b = b | myraster[[i]]
}
#plot(a & !b)
cells = Which(a & !b, cells=TRUE)
return(cells)
}
# Applying the function
results = lapply(values, find_zones)
我当前的功能需要很长时间才能执行。你能想出更好的方法吗?请注意,我不只是想知道每个像素有多少波段有数据,我还需要知道哪些波段。这样做的目的是在之后以不同的方式处理不同的区域。
另请注意,实际场景是 3000 x 3000 或更大的栅格,可能有超过 10 个波段。
编辑
由 10 个偏移图像区域组成的一些示例数据:
# Sample data
library(raster)
for(i in 1:10) {
start_line = i*10*1000
end_line = 1000000 - 800*1000 - start_line
offset = i * 10
data = c(rep(0,start_line), rep(c(rep(0,offset), rep(1,800), rep(0,200-offset)), 800), rep(0, end_line))
current_layer = raster(nrows=1000, ncols=1000)
values(current_layer) = data
if(i == 1) {
myraster = stack(current_layer)
} else {
myraster = addLayer(myraster, current_layer)
}
}
NAvalue(myraster) = 0 # You may not want to do this depending on your solution...