rgispolygonrasterterra

Merging diagonal-touching polygons from raster using terra


Consider the following binary raster, with patches in various shapes:

Binary raster showing patches, where some patches join by a single shared diagonal vertex

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:

Polygons are shown where diagonally touching polygons are separate geometries.

I would like patches that connect diagonally, via a single shared vertex, to be the same geometry. So the end result would be:

Polygons are shown but those touching diagonally are the same geometry.

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:

  1. Select one polygon
  2. Get its extent and extend it very minimally
  3. Crop original vector data based on this extended extent
  4. See which geometries touch the original geometry
  5. Update the merge values for these geometries, and the original, to be the same

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")

Solution

  • 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)
    

    enter image description here

    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.