Menu
I have an X by Y grid with cells containing 1 if a certain criteria is met or 0 if it is not. Now I want to identify features in the grid where there are at least N contiguous cells containing a 1. Contiguous cells can be adjacent side by side, or adjacent diagonally. I made a picture to illustrate the problem (see link), with N = 5. For clarity I omitted marking the 0s, and they are in the unmarked cells. Red 1s belong to features I want to identify, and black 1s do not. The desired result would be as shown in the picture, but with all the black 1s changed to 0s. I use R, so solutions using that language would be thoroughly appreciated, but I'll happily settle for others. I couldn't find anything in the R libraries (such as rgeos) specifically, but maybe I'm missing something. Any help appreciated, thanks!
Here is a small reproducible example created
input.mat <- structure(c(1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L,
0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L,
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L,
1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L,
0L, 1L, 1L, 1L), .Dim = c(15L, 15L), .Dimnames = list(NULL, NULL))
input.mat
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 1 1 0 0 0 0 0 0 0 0 0 0 1 0 0
[2,] 1 1 0 0 1 1 1 0 0 1 0 0 0 1 0
[3,] 0 0 1 0 0 0 0 0 0 1 1 0 1 0 1
[4,] 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0
[5,] 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
[6,] 1 0 0 0 0 0 0 0 0 0 1 0 1 1 0
[7,] 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0
[8,] 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
[9,] 1 0 0 0 0 1 0 1 0 0 0 1 1 1 0
[10,] 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0
[11,] 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1
[12,] 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0
[13,] 0 0 1 0 1 0 0 0 1 0 0 0 0 0 1
[14,] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1
[15,] 1 1 1 1 1 0 0 0 1 1 0 0 0 0 1
output.mat <- structure(c(1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L,
0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L,
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L,
1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L,
1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L,
0L, 0L, 0L, 0L), .Dim = c(15L, 15L), .Dimnames = list(NULL, NULL))
output.mat
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 1 1 0 0 0 0 0 0 0 0 0 0 1 0 0
[2,] 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0
[3,] 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1
[4,] 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0
[5,] 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
[6,] 1 0 0 0 0 0 0 0 0 0 1 0 1 1 0
[7,] 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0
[8,] 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
[9,] 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0
[10,] 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0
[11,] 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1
[12,] 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0
[13,] 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0
[14,] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
[15,] 1 1 1 1 1 0 0 0 1 1 0 0 0 0 0
Created on 2021-05-27 by the reprex package (v2.0.0)
227
asked Sep 11 '25 07:09
Using terra functions:
Convert matrix to raster (rast). Identify patches of 1s, surrounded by zeros (zeroAsNA = TRUE). Consider also diagonal neighbors when defining contiguity (directions = 8). Count number of cells in each patch (freq). Check which patches have a count of < 5. At these indices, set cells to NA. Coerce raster to matrix and check which values are NA. At these indices, set original matrix values to 0.
library(terra)
m = input.mat
p = patches(rast(input.mat), directions = 8, zeroAsNA = TRUE)
p[p %in% which(freq(p)[ , "count"] < 5)] = NA
m[is.na(as.matrix(p, wide = TRUE))] = 0
all.equal(m, output.mat)
# [1] TRUE
Patches in original input.mat (plot(p)):
After removal of patches with < 5 cells:
Related posts: Combining polygons and calculating their area (i.e. number of cells) in R; Obtaining connected components in R
140
With data.table non equi-join to find neighbouring points and igraph:
library(igraph)
library(data.table)
# index of pixels fulfilling criteria
idx <- which(input.mat==1)
# Coordinates of pixels
coord <- data.table(arrayInd(idx,dim(input.mat)))
setnames(coord,c("x","y"))
coord[,c('xmin','xmax','ymin','ymax'):=.(x-1,x+1,y-1,y+1)]
# Find neighbours indices
neighbours <- coord[coord,.(x.x,x.y,i.x,i.y),on=.(x>=xmin,x<=xmax,y>=ymin,y<=ymax)][!(i.x==x.x&i.y==x.y)][
,.(start = nrow(input.mat)*(x.y-1)+x.x,
end = nrow(input.mat)*(i.y-1)+i.x)]
g <- graph_from_data_frame(neighbours)
g
#> IGRAPH 503ba64 DN-- 53 120 --
#> + attr: name (v/c)
#> + edges from 503ba64 (vertex names):
#> [1] 2 ->1 16 ->1 17 ->1 1 ->2 16 ->2 17 ->2 7 ->6 22 ->6
#> [9] 6 ->7 8 ->7 22 ->7 23 ->7 7 ->8 9 ->8 22 ->8 23 ->8
#> [17] 8 ->9 23 ->9 30 ->15 1 ->16 2 ->16 17 ->16 1 ->17 2 ->17
#> [25] 16 ->17 33 ->17 6 ->22 7 ->22 8 ->22 23 ->22 7 ->23 8 ->23
#> [33] 9 ->23 22 ->23 15 ->30 45 ->30 17 ->33 49 ->33 57 ->41 57 ->43
#> [41] 30 ->45 60 ->45 33 ->49 41 ->57 43 ->57 71 ->57 73 ->57 45 ->60
#> [49] 75 ->60 77 ->62 57 ->71 57 ->73 60 ->75 62 ->77 92 ->77 77 ->92
#> [57] 134->133 147->133 133->134 135->134 150->134 134->135 150->135 138->137
#> + ... omitted several edges
# Find clusters
clust <- clusters(g)
# Minimum size
kept <- clust$membership[clust$membership %in% which(clust$csize >= 5)]
idx_kept <- as.numeric(names(kept))
M <- input.mat*0
M[idx_kept]<-1
M
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
#> [1,] 1 1 0 0 0 0 0 0 0 0 0 0 1
#> [2,] 1 1 0 0 0 0 0 0 0 0 0 0 0
#> [3,] 0 0 1 0 0 0 0 0 0 0 0 0 1
#> [4,] 0 0 0 1 0 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0 0 0 0 1 0
#> [6,] 1 0 0 0 0 0 0 0 0 0 1 0 1
#> [7,] 1 1 0 0 0 0 0 0 0 0 0 1 0
#> [8,] 1 1 0 0 0 0 0 0 0 0 0 0 0
#> [9,] 1 0 0 0 0 0 0 0 0 0 0 1 1
#> [10,] 0 0 0 0 0 0 0 0 0 0 0 1 1
#> [11,] 0 0 1 0 1 0 0 0 0 0 0 0 0
#> [12,] 0 0 0 1 0 0 0 0 0 1 0 0 0
#> [13,] 0 0 1 0 1 0 0 0 1 0 0 0 0
#> [14,] 0 0 0 0 0 0 0 0 1 0 0 0 0
#> [15,] 1 1 1 1 1 0 0 0 1 1 0 0 0
#> [,14] [,15]
#> [1,] 0 0
#> [2,] 1 0
#> [3,] 0 1
#> [4,] 1 0
#> [5,] 0 0
#> [6,] 1 0
#> [7,] 0 0
#> [8,] 0 0
#> [9,] 1 0
#> [10,] 1 0
#> [11,] 0 1
#> [12,] 0 0
#> [13,] 0 0
#> [14,] 0 0
#> [15,] 0 0
all.equal(output.mat,M)
#[1] TRUE
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
