fix self intersecting polygons using ogr2ogr or other open source librar

Question

I have a shapefile that includes some self intersecting polygons. What would be the easiest way to fix them using open source Python/command line? Is it possible to use ogr2ogr to do this (answers using fiona, shapely, geopandas etc. are also welcome) ? I want to do this because Google Earth Engine cannot handle self-intersecting polygons.

The answer using spatialite seems like a good option, however my GDAL doesn't currently support spatialite functions. I built GDAL using Docker.

ogrinfo --formats sqlite | grep 'spatialite' -i

Answer 1

Spatiallite has a MakeValid method which should fix the issues with invalid polygons. You can run it against your dataset using the sql parameter in ogr2ogr like so:

ogr2ogr -dialect SQLITE \
-sql "select PFAF_ID, ST_MakeValid(Geometry) as Geometry from hybas_lev06_v1c_merged_fiona_upstream_downstream_FAO_V01" \
valid_corrected.shp \
hybas_lev06_v1c_merged_fiona_upstream_downstream_FAO_V01.shp

Answer 2

Since you mentioned you had access to docker, you can utilize this postgis image to clean/export your geometries in a single docker run command:

docker run -it --rm -e POSTGRES_PASSWORD=mysecretpassword \
 --mount type=bind,source=<host_shapefile_location>,target=/tmp/shapedir \
 mdillon/postgis bash -c "
apt-get update && apt-get install gdal-bin && tee /docker-entrypoint-initdb.d/z.sh <<EOF
psql -c $'CREATE DATABASE ogr WITH TEMPLATE = 'template_postgis''
ogr2ogr -f 'PostgreSQL' -a_srs 'EPSG:<YOUR_SRID>' PG:'host=localhost user=postgres dbname=ogr password=mysecretpassword' /tmp/shapedir/myshp.shp -nlt MULTIPOLYGON -lco GEOMETRY_NAME=geom -nln myshp
ogr2ogr -overwrite -f 'ESRI Shapefile' /tmp/shapedir/output.shp PG:'host=localhost user=postgres dbname=ogr password=mysecretpassword' -sql 'SELECT myfield, st_collectionextract(st_makevalid(st_buffer(st_makevalid(geom),0)),3) FROM myshp'
EOF
/docker-entrypoint.sh postgres"

This command does the following:

1. Bind mounts the host directory containing the shapefile to /tmp/shapedir in the container.
2. Installs "vanilla" GDAL into the container.
3. Creates a Postgis database suitable for processing the shapefile.
4. Uses ogr2ogr to import the shapefile into the newly created database (and coerces the geometry to MULTIPOLYGON for the sake of simplicity and sanity)
5. Cleans geometries using st_makevalid(geom) and st_buffer(geom,0) and then applies st_collectionextract(geom,3) to ensure that only multipolygons are returned.
6. Exports the cleaned geometry to the specified host directory.

Once your cleaned shapefile has successfully exported, just ctrl +c in your terminal to stop and remove the temporary container.

Note: Make sure to explicitly provide the correct EPSG code via a_srs -OGR is not great at autodetecting it! You can look up the correct projection by pasting the contents of the .prj file here