5

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
Georgy
  • 12,464
  • 7
  • 65
  • 73
sal
  • 1,199
  • 1
  • 13
  • 31
  • Since you're PRINTING some things to the console, could you include this output in your question please? – Hans VK Nov 24 '15 at 16:25
  • Many things in numerics are inexact. However what do you really want to do with this information? Can't you just assume that the intersection points are very close to the boundary - since they were computed exactly that way? No need to check that really, or is there? – NoDataDumpNoContribution Nov 25 '15 at 08:31

1 Answers1

7

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)
Mike T
  • 41,085
  • 18
  • 152
  • 203
  • How computationally intensive can this distance calculation be? Could a "touches" implementation including a precision model one day be much faster? What algorithms would one use for it? (Just out of curiosity.) – NoDataDumpNoContribution Nov 25 '15 at 08:33
  • Just checked the same in postgis (with GEOS 3.4.2): ```SELECT st_touches(st_geometryn(st_intersection(l,st_boundary(p)),1),st_boundary(p)) from ST_GeomFromText('LINESTRING(13.51039642756912 52.59891281441467, 13.52517380027718 52.60620240344557)',4326) l, ST_GeomFromText('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))',4326) p; ``` Same result – sal Nov 25 '15 at 08:40
  • @Trilarion I just tested if the precision model works (using JTS) for "touches", and it does not work as I thought it might. A distance-based metric for this type of predicate is actually very common under the hood of most GIS software, and it is much less computationally intensive than other operations like "buffer". – Mike T Nov 29 '15 at 19:40