[GIS] How to convert self intersecting polygon to multipolygon

algorithmgdalpolygon

Working with the union and polygonize methods in GDAL/OGR, I often end up with self intersecting polygons. The case is one polygon which essentially consists of two polygons:

Self intersecting polygon example

In the above example, the polygon has two vertices at the intersection point, but whether the lines meet and goes seperate ways at the point or crosses each other, is unknown.

Are there methods for splitting a self intersecting polygon like above, into two polygons or a multipolygon?
Alternatively are there any existing algorithms I could implement my self (I would rather not invent another wheel)?

Best Answer

It's strange to get such self intersecting polygons because it violates the polygon definition in OGC Simple feature access specification and may not be viewed as legal polygons.

I don't use GDAL that often, but am quite familiar with GEOS and if I remember clearly GDAL union/polygonize are built based on GEOS library. In GEOS, if you union two polygons that share one vertex coordinates, like the case in your given picture, it will return you a multipolygon. In case you might want to look into the code, here is how GEOS do it:

http://svn.osgeo.org/geos/branches/3.3/src/geom/

The relevant functions may be the Geometry::Union in file Geometry.cpp

And also the GeometryFactory::buildGeometry in file GeometryFactory.cpp

Related Solutions

Polygon Validity – Detecting Invalid Polygons (Self-Intersecting or Touching)

If you don't want to write something from scratch, you could possibly call Check Geometry from ArcToolbox. I've used it a few times for certain geometries. Not sure if all of the above cases will be returned as bad geometry. You may want to test the tool on a test dataset. Please keep us posted.

R – Dissolving/Unifying Ill-Behaved Polygons in R

Without your original data, I can't be sure this will work, but I thought it might help you out. I didn't bring it all the way there, this solution still likely needs some level of automation, but might give you a general way forward

First, I create some spatial polygons

polypoints1 <- matrix(c(1,2,2,1,1,2,2,1,1,2),ncol=2)
polypoints2 <- matrix(c(1,3,3,1,1,3,3,1,1,3),ncol=2)
polypoints3 <- matrix(c(1,2,2,1,1,2,2,1,1,2)+1.1,ncol=2)
polypoints4 <- matrix(c(1,2,2,1,1,2,2,1,1,2)+0.5,ncol=2)

p1 <- Polygon(polypoints1)
ps1 <- Polygons(list(p1),1)
sps1 <- SpatialPolygons(list(ps1))

p2 <- Polygon(polypoints2)
ps2 <- Polygons(list(p2),2)
sps2 <- SpatialPolygons(list(ps2))

p3 <- Polygon(polypoints3)
ps3 <- Polygons(list(p3),3)
sps3 <- SpatialPolygons(list(ps3))

p4 <- Polygon(polypoints4)
ps4 <- Polygons(list(p4),4)
sps4 <- SpatialPolygons(list(ps4))

I plotted them just to see

plot(sps2,col='green')
plot(sps1,add=T,col='blue')
plot(sps3,add=T,col='yellow')
plot(sps4,add=T,col='purple')

I merged them into an spdf

data=data.frame(c(x=rep(1,4)),row.names=c(1:4))
sps <- SpatialPolygons(list(ps1,ps2,ps3,ps4))
spdf <- SpatialPolygonsDataFrame(sps,data)

You can identify which polygon overlaps which like so:

gIntersects(spdf,spdf,byid =T)

From the above command you could create some kind of loops to do the overlapping combinations below (I'm just ignoring sps4 for brevity at this point)

poly2a <- gIntersection(spdf[2,],spdf[1,],drop_lower_td=T)
poly2a <- SpatialPolygonsDataFrame(poly2a,data.frame(c(x=1),row.names=c(1)))
plot(poly2a,add=T,col='red')

This time we need to change the ID since we're going to rbind these later

poly2b <- gIntersection(spdf[2,],spdf[3,],drop_lower_td=T)
poly2b <- spChFIDs(poly2b,"2")
poly2b <- SpatialPolygonsDataFrame(poly2b,data.frame(c(x=1),row.names=c(2)))
plot(poly2b,add=T,col='red')

Merge the overlapping polygons into another spdf

spdf_overlaps <- rbind(poly2a,poly2b)
poly2 <- unionSpatialPolygons(spdf_overlaps,rep(1,2))
plot(poly2,add=T,col='blue')

Now we have poly2 which is where we have 2 layers overlapping (except combinations with sps4) then to figure out 3 layers, we just have to check out where poly2 and spdf overlap (if you make a more automated version of this, you'll need to make sure that 'poly2' does not include sps4 as in this example)

gIntersects(poly2,spdf[4,],byid =T)

poly3 <- gIntersection(poly2,spdf[4,],drop_lower_td=T)
plot(poly3,add=T,col="red")

Check it out

gIsValid(poly2)
gIsValid(poly3)

Alternatively, you could always do a pseudo rasterization, much easier, but you loose some detail depending on your cell size:

First make the grid:

bb <- bbox(spdf)
cs <- c(0.1,0.1)  # cell size
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)


sp_grd <- SpatialGridDataFrame(grd,
                           data=data.frame(id=1:prod(cd)))

Then, make grid into a polygon which used for overlap

library(Grid2Polygons)
grid <- Grid2Polygons(sp_grd)
plot(grid)

Then count the number of polygons that overlap each grid cell

count <- apply(gContains(spdf,grid,byid=T),1,sum)

Finally, plot it!

plot(grid)
for(i in 1:length(grid)){
    plot(grid[i,],col=rev(heat.colors(3))[count[i]],add=T)
}

Best Answer

Related Solutions

Polygon Validity – Detecting Invalid Polygons (Self-Intersecting or Touching)

R – Dissolving/Unifying Ill-Behaved Polygons in R

Related Question