1

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?

PS:

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.

Improve this question
5
  • 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
    Commented Dec 25, 2018 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
    Commented Dec 25, 2018 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
    Commented Dec 25, 2018 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
    Commented Dec 27, 2018 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
    Commented Dec 27, 2018 at 9:37

1 Answer 1

Reset to default
2

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:
(unknown)
Geometry Column = ST_MakeValid(geometry)
OGRFeature(SELECT):0
  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:
(unknown)
R_G_B: String (0.0)
OGRFeature(split_to_polygons):0
  R_G_B (String) = (null)
  POLYGON ((11.5 17.0,15 13,8 13,11.5 17.0))

OGRFeature(split_to_polygons):1
  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.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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