A point interpolated from a line does not always intersect the same line

[starsareblue]

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

[pramsey]

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?

[starsareblue]

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

[pramsey]

Not sure what these interpolate and project methods are...

[starsareblue]

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

[sgillies]

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?

[starsareblue]

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.

[sgillies]

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.

[strk]

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.

[sgillies]

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?

[strk]

That's just for single curves, as the return is a distance from first point along the line.

[strk]

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 ?

[strk]

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 :/

[strk]

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.

[strk]

The fix makes use of a yet-to-be-ported LocationIndexedLine?

[strk]

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.

[strk]

since nobody answered I'll consider this fine as is (fixed in 3.4.0)

[starsareblue]

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.