Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Shapely intersections vs shapely relationships - inexact?

I wonder if I am thinking the wrong way or if this is a bug:

I have a linestring and a polygon, I create the intersection points of the line and the polygon's boundary

enter image description here

These intersection points should intersect (at least touch) the polygon's boundary, right?

from shapely import geometry,wkt
line = geometry.LineString([(13.51039642756912, 52.598912814414675), (13.525173800277184, 52.60620240344557)])
poly = geometry.Polygon ([(13.52072838433517, 52.61735554606274), (13.52233276805985, 52.59511541819082), (13.51312087418833, 52.59394589806786),( 13.51526963068252, 52.60338701649216),( 13.51836560008325 ,52.6009395669487), (13.52072838433517, 52.61735554606274)])

ips = line.intersection(poly.boundary)
for i in ips:
    print i.touches(poly.boundary)  # should touch but it doesnt!!!!

>>>False
like image 520
sal Avatar asked Nov 24 '15 01:11

sal


1 Answers

It's not a bug, but this is a frequent question.

Without a precision model, all floating point calculations are limited by the machine epsilon. The intersected points are interpolated from each geometry, and are seldom exact (unless you have right angles). All of the DE-9IM predicates like 'touches' currently require exact noding (unless we have a precision model, which might happen one day UPDATE: testing with JTS Topology Suite, the DE-9IM predicates don't use the precision model, therefore it is unlikely that the GEOS clone will work any different).

A more robust strategy is to test the distance between the two, which should be less than the machine epsilon if they intersect. E.g.:

EPS = 1e-15
for i in ips:
    print(i.distance(poly) < EPS)
like image 154
Mike T Avatar answered Oct 20 '22 01:10

Mike T