移动窗口光栅泛洪填充R

Question

假设我有一个带有整数值的光栅，将事件的时间描述为一年中的某一天(DOY)。如果在相应年份中没有事件，则将单元格设置为NA。R"光栅"包的clump()函数将允许检测具有相同整数值的相邻单元格，并用唯一ID标记它们。现在，想象这样的事件(例如火灾)可以随着时间的推移在空间中传播，因此单元格(x，y)在DOY 1和相邻单元格(例如(x+1，y)、(x，y+1)…)上燃烧然后在DOY 2上被烧伤。因此，我想识别这样的事件，其中相邻像素在最大2天的DOY差内燃烧(例如，DOY(x，y)=13和DOY(x+1，y)=15)，并为其分配一个唯一的ID:

library(raster)
m<-matrix(c(1,10,11,14,
2,2,13,NA,
20,3,25,NA,
21,25,7,NA), ncol=4, byrow = TRUE)
r<-raster(m) # raster object of matrix

应产生光栅：

res_m<-matrix(c(1,2,2,2,
1,1,2,NA,
3,1,4,NA,
3,4,5, NA), ncol=4, byrow = TRUE)
res_r<-raster(res_m)

或图形化：

par(mfrow=c(1,2))
plot(r, xlim=(c(0:1)), main="DOY")
text(r)
plot(res_r, xlim=(c(0:1)), main="classified result")
text(res_r)

图：初始DOY光栅(左)与分类结果(右)

编辑：参考Lorenzo的评论：传播为DOY1、DOY2和DOY4的事件应被视为一个事件。然而，我不知道算法会是什么样子，两个不同的事件在传播时"融化"，但仍然会被归类为两个不同事件。到目前为止，我解决的问题相当低效，如下所示：

#Round 1: find connected components
#cell indices
coli<-rep(1:ncol(r), nrow(r)) 
rowi<-rep(1:ncol(r), each= nrow(r)) 
#neighbourhood matrix (considering only NW, N, NE and W neighbours)
mat_nb <- matrix(c(1,1,1,
1,0,NA,
NA,NA,NA), nrow=3, ncol=3, byrow = T)
#create ascending class raster
cls<-1:ncell(r)
mcl<-setValues(r, cls)
#create empty raster to fill
ecl<-setValues(r, NA)
#loop through cells
for (j in 1:length(coli)){
#####get adjacent cells 
zelle<-cellFromRowCol(r, rowi[j], coli[j])
nb <- adjacent(r, zelle, directions=mat_nb, pairs=F, sorted=T)
if(is.na(r[zelle])) {next} # if cell=NA go to next cells
if(length(nb) == 0) {ecl[zelle] <- mcl[zelle]} # if no neighbours, use ascending class value
if(length(nb) > 0) {
if(all(!is.na(r[nb[]]) & r[nb[]] %in% (r[zelle]-2):(r[zelle]+2) & !(unique(ecl[nb[]])))) 
{ecl[zelle] <- ecl[nb[1]]}  # if all neighbours valid and from same class, assign to class
if(!is.na(r[nb[1]]) & r[nb[1]] %in% (r[zelle]-2):(r[zelle]+2) & is.na(ecl[zelle])) 
{ecl[zelle] <- ecl[nb[1]]} # if NW neighbour valid and zelle still unclassified, assign neighbour's class
if(!is.na(r[nb[2]]) & r[nb[2]] %in% (r[zelle]-2):(r[zelle]+2) & is.na(ecl[zelle])) 
{ecl[zelle] <- ecl[nb[2]]}  # same for N
if(!is.na(r[nb[3]]) & r[nb[3]] %in% (r[zelle]-2):(r[zelle]+2) & is.na(ecl[zelle])) 
{ecl[zelle] <- ecl[nb[3]]}  # same for NE
if(!is.na(r[nb[4]]) & r[nb[4]] %in% (r[zelle]-2):(r[zelle]+2) & is.na(ecl[zelle])) 
{ecl[zelle] <- ecl[nb[4]]} # same for W
if(all(!(r[nb[]] %in% (r[zelle]-2):(r[zelle]+2)))) {ecl[zelle] <- mcl[zelle]} # if all neighbours "invalid", assign scending class value
}
} # warnings: from pixels with less than 4 nbs
#compare result with initial raster
par(mfrow=c(1,2))
plot(r)
text(r)
plot(ecl)
text(ecl)

在第二轮中，连接的组件类被合并。

##Round 2: combine classes
ecla<-ecl #save from first recursion
# only E, SW, S and SE neighbours
mat_agg<-matrix(c(NA,NA,NA,
NA,0,1,
1,1,1), nrow=3, ncol=3, byrow = T)
for (i in 1:length(coli)){
#####get adjacent cells 
zelle<-cellFromRowCol(r, rowi[i], coli[i])
nb <- adjacent(r, zelle, directions=mat_agg, pairs=F, sorted=T)
if(is.na(r[zelle])) {next}
if(length(nb) == 0) {ecl[zelle] <- mcl[zelle]}
if(length(nb) > 0) {
if(r[nb[2]] %in% (r[zelle]-2):(r[zelle]+2) & ecla[zelle] < ecla[nb[2]]) {ecla[nb[2]] <- ecla[zelle]}
if(r[nb[2]] %in% (r[zelle]-2):(r[zelle]+2) & ecla[zelle] > ecla[nb[2]]) {ecla[zelle] <- ecla[nb[2]]}
if(r[nb[3]] %in% (r[zelle]-2):(r[zelle]+2) & ecla[zelle] < ecla[nb[3]]) {ecla[nb[3]] <- ecla[zelle]}
if(r[nb[3]] %in% (r[zelle]-2):(r[zelle]+2) & ecla[zelle] > ecla[nb[3]]) {ecla[zelle] <- ecla[nb[3]]}
if(r[nb[4]] %in% (r[zelle]-2):(r[zelle]+2) & ecla[zelle] < ecla[nb[4]]) {ecla[nb[4]] <- ecla[zelle]}
if(r[nb[4]] %in% (r[zelle]-2):(r[zelle]+2) & ecla[zelle] > ecla[nb[4]]) {ecla[zelle] <- ecla[nb[4]]}
if(r[nb[1]] %in% (r[zelle]-2):(r[zelle]+2) & ecla[zelle] < ecla[nb[1]]) {ecla[nb[1]] <- ecla[zelle]}
if(r[nb[1]] %in% (r[zelle]-2):(r[zelle]+2) & ecla[zelle] > ecla[nb[1]]) {ecla[zelle] <- ecla[nb[1]]}
} # warnings: from pixels with less than 4 nbs
}
# plot results
par(mfrow=c(1,3))
plot(ecl) # round 1 result
text(ecl)
plot(r)
text(r)
plot(ecla) # round 2 result
text(ecla)

Answer 1

这是一个棘手的老问题。我放弃了将一个复杂的函数传递给raster:：focal，所以我改为使用data.table进行处理，并使用了一系列规则。也许有更简单的方法可以做到这一点，但无论如何。

这适用于您的数据和我生成的6x6光栅。请测试一下，看看进展如何；

# packages
library(raster)
library(data.table)
# your example data
m<-matrix(c(1,10,11,14,
2,2,13,NA,
20,3,25,NA,
21,25,7,NA), ncol=4, byrow = TRUE)
r<-raster(m)
# convert raster to data.table, add cell number attribute, and zonal ID column 
# and columns for the queen's case relationships (these columns contain the cell ID 
# that holds that relationship)
df <- as.data.table(as.data.frame(r))
df[,CELL:=as.numeric(row.names(df))]
df[,ID:=0]
df[,c("TL","T","TR","R","BR","B","BL","L") := 
list(CELL-(ncol(r)+1),CELL-ncol(r),CELL-(ncol(r)-1),
CELL+1,CELL+(ncol(r)+1),CELL+ncol(r),CELL+(ncol(r)-1),CELL-1)]
# set appropriate queen's case relations to NA for all edge cells
df[c(CELL %in% seq(1,ncell(r),ncol(r))), c("TL","L","BL")] <- NA
df[c(CELL %in% seq(ncol(r),ncell(r),ncol(r))), c("TR","R","BR")] <- NA
df[c(CELL %in% 1:ncol(r)), c("TL","T","TR")] <- NA
df[c(CELL %in% (ncell(r)-nrow(r)+1):ncell(r)), c("BL","B","BR")] <- NA
# melt data table, add the cell values that correspond to each relationship,
# remove rows that wont be needed just now
dfm <- melt(df,id.vars=c("layer","CELL","ID"))
dfm[,NEIGH.VAL := r[dfm[,value]]]
dfm[,WITHIN2:=ifelse(abs(layer-NEIGH.VAL)>2,F,T)]
dfm <- dfm[!is.na(value)&!is.na(layer)&!is.na(NEIGH.VAL)][order(CELL)]
dfm <- dfm[WITHIN2==TRUE]
dfm[,NEIGH.ID := 0]
# this is a tricky loop and any errors that may occur are more than likely produced here.
# It starts at the lowest cell ID with a value and builds a lookup up vector of cell IDs
# that conform to the DOY being within 2 days, then sets the ID for them all.
# It then moves on to the next cell ID available that has not had a zone ID assigned
# and does the same thing, with a zonal ID one higher. etc.
for(n in unique(dfm[,CELL])){
while(nrow(dfm[CELL==n & NEIGH.ID==0]) > 0){
lookups <- unique(dfm[CELL==n,value])
while(all(unique(dfm[CELL %in% lookups,value]) == lookups)==F){
lookups <- unique(c(lookups,dfm[CELL %in% lookups,value]))
lookups <- unique(dfm[CELL %in% lookups,value])}
dfm[CELL %in% lookups, NEIGH.ID := max(dfm[,NEIGH.ID])+1]
}
}
# this section creates a raster of 1:ncell of the original data and assigns a zonal ID
# with reclassify. Everything that did not get a zonal ID (single 'island' cells
# and original NA cells) becomes NA
results <- unique(dfm[,list(CELL,NEIGH.ID)])
rzone <- r
rzone[] <- 1:ncell(rzone)
rc <- reclassify(rzone, results)
rc[!(rzone %in% results[[1]])] <- NA
# Now we need to determine those single 'island' cells and reinstate them, with their
# own ID, increasing incrementally from the highest extant ID based on the above analysis
final.vec <- which(!(rzone[] %in% results[[1]]) & !(rzone[] %in% which(is.na(values(r)))))
rc[final.vec] <- seq((cellStats(rc,max)+1),(cellStats(rc,max)+length(rc[final.vec])),1)
# plot to check
par(mfrow=c(1,3))
plot(r, xlim=(c(0:1)), main="DOY")
text(r)
plot(res_r, xlim=(c(0:1)), main="DOY")
text(res_r)
plot(rc, xlim=(c(0:1)), main="zones result")
text(rc)

这段代码也适用于以下内容，不过请大家玩一玩吧！(忽略警告)；

m<-matrix(c(1,3,5,14,NA,21,
2,2,17,NA,23,25,
9,15,4,5,9,NA,
14,11,14,21,25,7,
NA,NA,16,2,4,6,
NA,17,25,15,1,NA), ncol=6, byrow = TRUE)
r<-raster(m) # raster object of matrix