moving-window raster flood fill in R 中
moving-window raster flood fill in R
假设我有一个带有整数值的栅格，将事件的时间描述为一年中的某一天 (DOY)。如果相应年份没有事件，则单元格设置为 NA。 R 'raster' 包的 clump() 函数将允许检测具有相同整数值的相邻单元格并使用唯一 ID 标记它们。现在，想象一下这样的事件（例如火灾）可以随着时间的推移在 space 中传播，因此单元格 (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)
plot: initial DOY raster (left) vs. classified result (right)
编辑：
参考 Lorenzo 的评论：事件，其中传播例如DOY1、DOY2 和 DOY4 应视为一个事件。但是，我无法弄清楚算法会是什么样子，其中两个不同的事件 "melt" 在传播时仍然会被归类为两个不同的事件。
到目前为止，我解决的问题效率很低，如下所示：
#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)
这是一个棘手的老问题。我放弃了将复杂函数传递给 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
假设我有一个带有整数值的栅格，将事件的时间描述为一年中的某一天 (DOY)。如果相应年份没有事件，则单元格设置为 NA。 R 'raster' 包的 clump() 函数将允许检测具有相同整数值的相邻单元格并使用唯一 ID 标记它们。现在，想象一下这样的事件（例如火灾）可以随着时间的推移在 space 中传播，因此单元格 (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)
plot: initial DOY raster (left) vs. classified result (right)
编辑： 参考 Lorenzo 的评论：事件，其中传播例如DOY1、DOY2 和 DOY4 应视为一个事件。但是，我无法弄清楚算法会是什么样子，其中两个不同的事件 "melt" 在传播时仍然会被归类为两个不同的事件。 到目前为止，我解决的问题效率很低，如下所示：
#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)
这是一个棘手的老问题。我放弃了将复杂函数传递给 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