How to convert self-intersecting polygons into distinctive, well behaving polygons?

My polygon inputs consist of self-intersecting polygons, which fails the st_isvalid test, ie: the below statement returns some self-intersection rings:

ogrinfo -sql "select geom, st_isvalidreason(geom) from out" intersectingpolygon.shp -q

How can I sanitize the problematic polygons in intersectingpolygon.shp by breaking them into a few distinctive, well-behaving polygons, using GDAL?


There is a similar question on GIS, but I don't find it useful because the answers seem to deal with how to avoid the bad inputs in the first place, and how to convert MultiPolygon into a few Polygons. Nothing about how to convert self-intersecting polygons into well-behaving ones.

  • Could you explain from what source you receive selfintersecting polygons. The problem is that such polygons can be just the result of an incorrect input and appear to be kind of mess. The problem of accuracy of input can be involved. Maybe you need first to eliminate some extra vertices and line that lie within some threshold of accuracy. – Vadym Dec 25 '18 at 11:07
  • @Vadym , self intersecting polygons come from another library, which I can't fix. The library output is buggy and I can't do anything about ( let's assume it this way). Let's just concentrate on the original question instead of asking where those incorrect inputs come from – Graviton Dec 25 '18 at 11:17
  • Use ogr2ogr with sqlite dialect and -sql "select st_makevalid(geometry), attribute1, attribute2 from layer". That will turn self-intersecting polygons into multipolygons. Use -explodecollections if you want to continue to single-part polygons. – user30184 Dec 25 '18 at 12:58
  • @user30184, your idea doesn't really work because a lot of st_makevalid function is simply not available in ogr2ogr, see this post. – Graviton Dec 27 '18 at 4:09
  • You are right that GDAL with SpatiaLite with LWGEOM features is needed but it is a compile time decision to take them or not so "simply not available" is too strict. By your comment to a linked question you seem to use Windows and then you can use the OSGeo4W version of GDAL. See my answer that includes a simple test as a proof. – user30184 Dec 27 '18 at 9:37

You must have a GDAL version that is compiled with a SpatiaLite version that is compiled with LWGEOM features. At least OSGeo4W delivers such GDAL.

A simple test with a self-intersecting geometry and GDAL version

GDAL 2.3.2, released 2018/09/21

POLYGON (( 8 21, 15 21, 8 13, 15 13, 8 21 ))

I stored the polygon into OpenJUMP JML format.

First test with ogrinfo

ogrinfo -dialect sqlite -sql "select ST_MakeValid(geometry) from self_intersection" self_intersection.jml
Had to open data source read-only.
INFO: Open of `self_intersection.jml'
      using driver `JML' successful.

Layer name: SELECT
Geometry: Unknown (any)
Feature Count: 1
Extent: (8.000000, 13.000000) - (15.000000, 21.000000)
Layer SRS WKT:
Geometry Column = ST_MakeValid(geometry)
  MULTIPOLYGON (((11.5 17.0,15 13,8 13,11.5 17.0)),((8 21,15 21,11.5 17.0,8 21)))

Self-intersecting polygon was fixed into a multipolygon. Next test with ogr2ogr and -explodecollections.

ogr2ogr -f jml  -dialect sqlite -sql "select ST_MakeValid(geometry) from self_intersection" split_to_polygons.jml self_intersection.jml -explodecollections

Check the result with ogrinfo

ogrinfo split_to_polygons.jml -al
INFO: Open of `split_to_polygons.jml'
      using driver `JML' successful.

Layer name: split_to_polygons
Geometry: Unknown (any)
Feature Count: 2
Extent: (8.000000, 13.000000) - (15.000000, 21.000000)
Layer SRS WKT:
R_G_B: String (0.0)
  R_G_B (String) = (null)
  POLYGON ((11.5 17.0,15 13,8 13,11.5 17.0))

  R_G_B (String) = (null)
  POLYGON ((8 21,15 21,11.5 17.0,8 21))

Self-intersecting polygon has been converted into two simple polygons.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.