[SOLVED] Splitting a multipolygon into polygons with edges that don't intersect

Splitting a multipolygon into polygons with edges that don't intersect

I am working with spatial data and I need to split a multipolygon into polygons with no intersecting edges. For instance, consider this multipolygon mp:

p <- st_polygon(list(rbind(
     c(0, 0),
     c(2, 2),
     c(0, 2),
     c(2, 0),
     c(0, 0))))
q <- st_polygon(list(rbind(
     c(-4, -4),
     c(-2, -2),
     c(-4, -2),
     c(-2, -4),
     c(-4, -4))))
mp <- st_multipolygon(list(p,q))

How can I turn mp into these four polygons?

a <- st_polygon(list(rbind(
     c(1, 1),
     c(2, 2),
     c(0, 2),
     c(1, 1))))

b <- st_polygon(list(rbind(
     c(1, 1),
     c(0, 0),
     c(0, 2),
     c(1, 1))))

c <- st_polygon(list(rbind(
     c(-3, -3),
     c(-2, -4),
     c(-4, -4),
     c(-3, -3))))

d <- st_polygon(list(rbind(
     c(-3, -3),
     c(-4, -2),
     c(-2, -2),
     c(-3, -3))))

I already tried to use st_cast, but there are edges that cross

Solution

Valid simple features should be non-self intersecting - https://r-spatial.github.io/sf/articles/sf1.html#simple-feature-geometry-types - so not sure if it helps with your real data or not, but passing that mp through st_make_valid() turns it into a mulitpolygon of 4 polygons and that can be casted:

library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
p <- st_polygon(list(rbind(
  c(0, 0),
  c(2, 2),
  c(0, 2),
  c(2, 0),
  c(0, 0))))
q <- st_polygon(list(rbind(
  c(-4, -4),
  c(-2, -2),
  c(-4, -2),
  c(-2, -4),
  c(-4, -4))))
mp <- st_multipolygon(list(p,q))

st_is_valid(mp)
#> [1] FALSE

(mp2 <- st_make_valid(mp))
#> MULTIPOLYGON (((-4 -4, -3 -3, -2 -4, -4 -4)), ((-3 -3, -4 -2, -2 -2, -3 -3)), ((0 0, 1 1, 2 0, 0 0)), ((1 1, 0 2, 2 2, 1 1)))

st_cast(st_sfc(mp2), "POLYGON")
#> Geometry set for 4 features 
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -4 ymin: -4 xmax: 2 ymax: 2
#> CRS:           NA
#> POLYGON ((-4 -4, -3 -3, -2 -4, -4 -4))
#> POLYGON ((-3 -3, -4 -2, -2 -2, -3 -3))
#> POLYGON ((0 0, 1 1, 2 0, 0 0))
#> POLYGON ((1 1, 0 2, 2 2, 1 1))

^{Created on 2023-02-20 with reprex v2.0.2}