Consider the following binary raster, with patches in various shapes:
Using terra
to get separate polygons for each patch, I converted to vector via as.polygons()
and then disaggregated via disagg()
. After dropping polygons with 0 values, the result is this:
I would like patches that connect diagonally, via a single shared vertex, to be the same geometry. So the end result would be:
My working solution is to create a data frame that stores the polygon ID as well as a new field that is the value to ultimately merge/aggregate with. Then iterate over the original disaggregated vector data:
I found that cropping the vector data down - then checking for touching only among these "close" polygons - ran faster than checking all polygons for if they touch the current isolated polygon.
Is there a way to speed up this process? Or a function I'm missing that would do all of this but more quickly? I'm sure there are many solutions, but the end use is on 3k+ vector files, each with a few thousand polygons to potentially merge.
Code to reproduce, and my solution, here:
### setup
library(terra)
m <- matrix(
c(1,1,0,1,1,0,1,
1,0,0,1,1,0,0,
0,0,0,0,0,1,0,
0,1,0,0,0,0,0,
0,1,1,0,1,1,1,
0,0,0,1,0,0,0,
0,0,0,0,0,1,1),
nrow=7,
ncol=7,
byrow=TRUE
)
r <- rast(m)
plot(r)
v <- as.polygons(r)
v <- v[which(v$lyr.1 == 1),]
plot(v, "lyr.1")
v <- disagg(v)
v$id <- 1:nrow(v)
plot(v)
text(v, labels="id")
### solution
nv <- nrow(v)
mv_df <- data.frame(id=1:nv, mval=1:nv)
for (i in 1:nv) {
i_poly <- v[i,]
i_ext <- extend(ext(i_poly), 1)
i_crop <- crop(v, i_ext)
i_touch <- is.related(i_crop, i_poly, relation="touches")
i_update <- append(i_crop$id[i_touch], i_poly$id)
mv_df$mval[i_update] <- min(mv_df$mval[i_update])
}
v <- merge(v, mv_df, by.x="id", by.y="id")
v <- aggregate(v, "mval")
plot(v, "mval")
text(v, labels="mval")
With your example data
library(terra)
r <- rast(matrix(nrow=7, ncol=7, c(1,1,0,0,0,0,0,1,0,0,1,1,0,0,0,0,0,0,1,0,0,1,1,0,0,0,1,0,1,1,0,0,1,0,0,0,0,1,0,1,0,1,1,0,0,0,1,0,1)))
I would do
x <- patches(r, directions=8, zeroAsNA=TRUE)
p <- as.polygons(x)
See:
plot(p, col=rainbow(5))
text(p, halo=TRUE, inside=TRUE)
The second polygon is a multi-polygon. The two parts should not share the same vertex in a single polygon as that would make the polygon invalid.