Doing a topological analysis in GDAL/OGR returned an unexpected result. When I started investigating, it seems possible that the Intersects() method is badly implemented. Specifically, points along a polygon's boundary returns False, when it should return True:
from osgeo import ogr poly = ogr.Geometry(ogr.wkbPolygon) ring = ogr.Geometry(ogr.wkbLinearRing) ring.AddPoint(0, 0) ring.AddPoint(5, 5) ring.AddPoint(10, 0) ring.CloseRings() poly.AddGeometry(ring) point = ogr.Geometry(ogr.wkbPoint) point.AddPoint(10, 0) poly.Intersects(point) >>> False poly.Disjoint(point) >>> False
Not only does the point correctly intersect the polygon, Intersects and Disjoint should be opposite subsets, and never return the same value. Furthermore, the Intersection() method, which should only work when geometries intersect, works normally:
inter = poly.Intersection(point) inter.ExportToWkt() >>> 'POINT (10 0 0)'
I'd like to know if perhaps I'm not understanding this correctly, or if my GDAL compilation is problematic, or if it's GDAL itself that has a bug.
Info: GDAL 2.2.1 64-bits, Tamas Szekeres' VS2008 version.