On 01/29/2010 03:51 AM, Sean Gillies wrote:
>
> On Jan 28, 2010, at 11:55 PM, Roy Hyunjin Han wrote:
>
>> from shapely.wkt import loads
>>
>> m = loads("""MULTILINESTRING ((336700.9845364579814486
>> 1738886.2085458301007748, 343584.4261802079854533
>> 1738555.8341833299491554), (343584.4261802079854533
>> 1738555.8341833299491554, 350539.0000000000000000
>> 1736135.0000000000000000), (350539.0000000000000000
>> 1736135.0000000000000000, 353532.0000000000000000
>> 1733674.0000000000000000))""")
>>
>> p = loads("""POINT (344705.7704083333374001
>> 1737219.7509562498889863)""")
>>
>> assert m.intersects(m.interpolate(m.project(p))) == True
>
> I just got a pair of email about GEOS #323. Point and line
> intersection is fraught with numerical precision problems. Using a
> geometry's distance method, you can implement something like an
> "almostIntersects":
>
> bool(a.distance(b)<= tolerance)
>
> Would be nice if GEOS projection and interpolation were exactly
> inverse, but that's not the case.
>
> --
> Sean
>
> _______________________________________________
> Community mailing list
> [email protected]
> http://lists.gispython.org/mailman/listinfo/community
Thanks, Sean. Since this is a GEOS issue, the conversation has been
moved there: http://trac.osgeo.org/geos/ticket/323
In [1]: import shapely.geometry as g
In [2]: m = g.MultiLineString([[(0, -2), (0, 2)], [(-2, 0), (2, 0)]])
In [3]: m.interpolate(m.project(g.Point(2,1.9))).wkt # Correct
Out[3]: 'POINT (2.0000000000000000 0.0000000000000000)'
In [4]: m.interpolate(m.project(g.Point(2,2.1))).wkt # Wrong
Out[4]: 'POINT (-2.0000000000000000 0.0000000000000000)' # != (2,0)
_______________________________________________
Community mailing list
[email protected]
http://lists.gispython.org/mailman/listinfo/community
