2

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

  • I am aware of the PostGIS solution but this takes some additional steps to setup. – RutgerH Sep 15 '17 at 21:20
2

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
  • With GDAL 2.1.1, released 2016/07/07 I get the error: ERROR 1: In ExecuteSQL(): sqlite3_prepare(select PFAF_ID, ST_MakeValid(Geometry) as Geometry from hybas_lev06_v1c_merged_fiona_upstream_downstream_FAO_V01): no such function: ST_MakeValid .....Looking for a solution now – RutgerH Sep 16 '17 at 3:53
  • ogrinfo --formats does not show SQLITE as a driver. Ubuntu 16.10 GDAL via apt-get – RutgerH Sep 16 '17 at 4:03
  • You may need to install liblwgeom (from PostGIS), then configure/install libspatialite (--enable-lwgeom), then install GDAL with libspatialite support, to get access to the newest functions. – James Mar 28 at 19:37
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

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.