заливка растровым окном в R

Предположим, у меня есть растр с целочисленными значениями, описывающими время события как день года (DOY). Если в соответствующем году не было никаких событий, в ячейках устанавливается значение NA. Функция clump() пакета R 'raster' позволит обнаруживать соседние ячейки с одинаковым целочисленным значением и маркировать их уникальным идентификатором. Теперь представьте, что такие события (например, пожар) могут распространяться в пространстве с течением времени, так что ячейка (x, y) сгорела на DOY 1, а соседние ячейки (например, (x + 1, y), (x, y + 1) , ...) затем сгорели на DOY 2. Следовательно, я хотел бы идентифицировать такие события, когда соседние пиксели сгорели с разницей в DOY не более 2 дней (например, DOY (x, y) = 13 и DOY (x + 1 , y) = 15) и присвойте им уникальный идентификатор:

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 (слева) и классифицированный результат (справа)

РЕДАКТИРОВАТЬ: Ссылаясь на комментарий Лоренцо: события, где распространение, например, 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)