#323 closed defect (fixed)
A point interpolated from a line does not always intersect the same line
Reported by: | starsareblue | Owned by: | sgillies |
---|---|---|---|
Priority: | major | Milestone: | 3.4.0 |
Component: | SWIG Python | Version: | 3.2.0 |
Severity: | Unassigned | Keywords: | history |
Cc: | sgillies |
Description
Using the shapely Python bindings to GEOS, the following assertion fails.
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
Change History (22)
comment:1 by , 14 years ago
Resolution: | → invalid |
---|---|
Status: | new → closed |
comment:2 by , 14 years ago
Resolution: | invalid |
---|---|
Status: | closed → reopened |
That is true, but in the following case, the projection point (0, -1) is exactly representable. Is it perhaps an issue with the MultiLineString?
I would expect the intersection of the following MultiLineString with the Point to produce either (0, -1) or (1, 0) but instead it produces (-1, 1). The distance() function works correctly though. So interpolation/projection onto a MultiLineString must be done manually for now for each individual LineString.
import shapely.geometry as g point = g.Point(1, 0) l1 = g.LineString([(0, 2), (0,0)]) l2 = g.LineString([(-1, 1), (1, 1)]) m = g.MultiLineString([l1.coords, l2.coords]) assert m.interpolate(m.project(point)).coords[0] == (-1, 1) # True, but wrong assert m.interpolate(m.project(point)).coords[0] == (0, -1) # Fails assert m.interpolate(m.project(point)).coords[0] == (1, 0) # Fails
comment:4 by , 14 years ago
The interpolate and project methods correspond to the linear referencing methods in the GEOS 3.2 API. You can see similar strange behavior in the GeoDjango package ( http://code.djangoproject.com/attachment/ticket/11948/contrib-gis-linearref.patch ) so I assumed the problem had to do with GEOS itself.
Here is another example. The first interpolation is correct. However, the second interpolation is wrong. It looks like some kind of wraparound is happening incompletely for the edge case when the interpolation result for a MultiLineString happens on the endpoint of one of its LineStrings.
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)'
Roy Hyunjin Han http://invisibleroads.com Empowering people to create since 2008
comment:5 by , 14 years ago
Roy, are you certain you don't have a coordinate mixup in your code?
The situation is a multilinestring consisting of 2 perpendicular lines that cross at (0,1), right? Interpolating (via Shapely) gives me a result that I expect: as the normalized length approaches 0.5, the interpolated point approaches (0, 0), the end of the first linestring. Go past length 0.5 and there's a jump to (-1, 1), the starting point of the second linestring.
I'm not seeing any wrap-around when the normalized length > 1: I always get the end of the second linestring. But with normalized length < 0, the interpolation appears to go in reverse, starting from the end of the second linestring and working its way back to the start of the first linestring.
Paul, that's an intended feature, right?
comment:6 by , 14 years ago
My problem is that I am trying to find the coordinates of the projection of a Point onto a MultiLineString. For example, the coordinates of the point (2,1) projected onto the line segment (0,0)-(5,0) should be (2,0).
I understand that the interpolate() method for a MultiLineString composed of two line segments wraps into the second line segment if the distance argument is longer than the first line segment. That seems to be working properly.
Then what exactly is MultiLineString.project(Point) returning?
I use LineString.interpolate(LineString.project(Point)) to find the coordinates of the point projected onto the LineString. So I assumed that MultiLineString.interpolate(MultiLineString.project(Point)) would do the same. However, sometimes it works and sometimes it doesn't, giving strange results as in the two examples above.
Sean, if the interpolate() method is working as intended, is project() for MultiLineStrings working as it is intended? What do 8.0 and 4.0 mean in the following example?
In [1]: import shapely.geometry as g In [2]: m = g.MultiLineString([[(0, -2), (0, 2)], [(-2, 0), (2, 0)]]) In [3]: m.project(g.Point(2,1.9)) Out[3]: 8.0 In [4]: m.project(g.Point(2,2.1)) Out[4]: 4.0
Currently, my workaround is to first compute the distance of the Point to the MultiLineString and then compute the ProjectedPoint using ProjectedPoint = MultiLineString.interpolate(MultiLineString.project(Point)). If the distance from the ProjectedPoint to the Point does not match the expected distance, then I redo the projection for each LineString inside MultiLineString until I find the ProjectedPoint whose distance matches the first distance.
comment:7 by , 14 years ago
Sorry, I didn't understand that you were trying to project a point not on the geometry. I have no idea what project is supposed to do in that case.
comment:8 by , 13 years ago
Component: | Default → SWIG Python |
---|---|
Milestone: | → 3.2.3 |
Owner: | changed from | to
Status: | reopened → new |
Sean: are you willing to take a look at this ? I've no 'project' interface in the C++ Geometry API so dunno what is this all about.
comment:9 by , 13 years ago
See use of geos::linearref::LengthIndexedLine::project in http://trac.osgeo.org/geos/browser/trunk/capi/geos_ts_c.cpp#L5665.
Is linear interpolation supported for all multilinestrings or just ones that form single curves?
comment:10 by , 13 years ago
That's just for single curves, as the return is a distance from first point along the line.
comment:11 by , 13 years ago
The problem really lies in Project returning a single number (a distance from start point of first component), which isn't appropriate to address multi-component geometries (multilinestring). We can't change the Project C-API interface, but eventually GEOSInterpolate could be changed to use other parts of the C++ API where a componentNumber/distance object is used for this kind of operations.
Anyone feels like working on a patch for it ?
comment:12 by , 13 years ago
Oops, GEOSInterpolate also take a single value. I'm afraid we can then only forbid using multi-component geometries, or deprecate the whole projection API :/
comment:13 by , 13 years ago
Milestone: | 3.2.3 → GEOS Future |
---|
comment:14 by , 13 years ago
Milestone: | GEOS Future → 3.3.1 |
---|
Fix committed in JTS: http://jts-topo-suite.svn.sourceforge.net/viewvc/jts-topo-suite?view=revision&revision=463 http://jts-topo-suite.svn.sourceforge.net/viewvc/jts-topo-suite?view=revision&revision=465
plus 464 and 467 as side minor tweaks.
The testcase is interestingly expecting the result of not interpolating.
comment:15 by , 13 years ago
Milestone: | 3.3.1 → 3.4.0 |
---|
The fix makes use of a yet-to-be-ported LocationIndexedLine
comment:16 by , 13 years ago
r3484 is a port of the JTS fix in trunk (to be 3.4.0).
Not sure we want this back-ported to 3.3.1 as it doesn't seem that the fix results in what the original bug reported was looking for.
Comments welcome.
comment:17 by , 12 years ago
Resolution: | → fixed |
---|---|
Status: | new → closed |
since nobody answered I'll consider this fine as is (fixed in 3.4.0)
comment:18 by , 12 years ago
Thank you so much for this fix!
Looking back at the comments, I retried the following lines and got the answer that I expected if I draw the same situation with pen and paper.
In [1]: import shapely.geometry as g In [2]: m = g.MultiLineString([[(0, -2), (0, 2)], [(-2, 0), (2, 0)]]) In [3]: print m.interpolate(m.project(g.Point(2,1.9))).wkt In [4]: print m.interpolate(m.project(g.Point(2,2.1))).wkt ## -- End pasted text -- POINT (2.0000000000000000 0.0000000000000000) POINT (0.0000000000000000 2.0000000000000000)
Judging from this single test case, I will also assume that this problem is fixed. Thanks again.
comment:19 by , 11 years ago
Keywords: | history added |
---|
The point isn't on the line, it's just very very very close. But in double precision, the coordinate space is chunky, not infinitely dividable. Think about diagonal line in graph paper. If you can only represent points where the graph lines intersect, how many places on the diagonal line are exactly representable?