QGIS – Fixing Incorrect Geometry Check from OGR and Shapely

gdalpythonqgisshapely

I have multipolygon and trying to check for validity. Both OGR(.IsValid()) and shapely (is_valid) return True. But if I open the same geometry with QGIS and check with theCheck Geometries plugin it returns multiple errors (duplicate node, self-contact etc). Therefore, I am confused on the results. Is it possible to implement the Check Geometries plugin without using QGIS itself (through python API) or how can I get exact results as the plugin's using shapely or OGR? Or any other libraries that I can get the same results?

My geometry can be found from here

Best Answer

For me it seems that both tools give correct results. The difference is that the OGR IsValid function does not apply any tolerance but the QGIS Check Geometries tool is using a tolerance of 1E-8 by default.

For example QGIS tool reports duplicate points at place 82.762954,46.935208. Actually there are two distinct vertices with coordinates

POINT (82.76295409400032 46.935207707000075)
POINT (82.76295409399872 46.935207707)

When the tolerance of 1E-8 is applied the coordinates are the same.

enter image description here

One of the self-contact cases look like this

enter image description here

Actually the polygon ring does not touch itself but there is just a very narrow intruding spike

enter image description here

It is possible to emulate the QGIS behavior with tolerance by using ST_SnapToGrid function https://postgis.net/docs/en/ST_SnapToGrid.html that rounds the coordinates but does not care if the new geometry is valid or not. Notice the difference to ST_ReducePrecion https://postgis.net/docs/en/ST_ReducePrecision.html that returns valid geometries. I linked the PostGIS manual pages instead of the SpatiaLite reference https://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html because PostGIS pages are more detailed. Both are using the same GEOS functions and are therefore equivalent.

A simple test ogrinfo and the data that you provided

Valid geometry if tolerance is not used

ogrinfo -dialect sqlite -sql "select st_isvalid(geometry) from map" map.json
OGRFeature(SELECT):0
  st_isvalid(geometry) (Integer) = 1

Applying tolerance makes the geometry invalid

ogrinfo -dialect sqlite -sql "select st_isvalid(st_snaptogrid(geometry,0.0000001)) from map" map.json
GEOS warning: Self-intersection at or near point 82.611105102748667 47.0228712000998
OGRFeature(SELECT):0
  st_isvalid(st_snaptogrid(geometry,0.0000001)) (Integer) = 0
Related Question