6

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.

Improve this question
3
  • 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
    Commented Oct 13, 2017 at 11:22
  • It is a bug of the version. My old GDAL 1.10.1 for your code returns True
    – nadya
    Commented Oct 15, 2017 at 0:03
  • Thanks nadya. As user30184 suggested, I have raised this issue in the mailing list.
    – Roberto Ribeiro
    Commented Oct 16, 2017 at 15:41

1 Answer 1

Reset to default
1

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

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.