5

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.

  • 2
    This feels pretty much like a bug. With POINT (10, 0) poly.Touches(point) returns true, as well as poly.Intersects(point) if the point is (5, 1). Write a mail to gdal-devel mailing list and report this. – user30184 Oct 13 '17 at 11:22
  • It is a bug of the version. My old GDAL 1.10.1 for your code returns True – nadya Oct 15 '17 at 0:03
  • Thanks nadya. As user30184 suggested, I have raised this issue in the mailing list. – Roberto Ribeiro Oct 16 '17 at 15:41
1

This has been fixed per GDAL ticket https://trac.osgeo.org/gdal/ticket/7091

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.