# Ticket #323 (closed defect: fixed)

Opened 3 years ago

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

Reported by: Owned by: starsareblue sgillies major 3.4.0 SWIG Python 3.2.0 Unassigned sgillies

### Description

Using the shapely Python bindings to GEOS, the following assertion fails.

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))""")

assert m.intersects(m.interpolate(m.project(p))) == True

## Change History

### Changed 3 years ago by pramsey

• status changed from new to closed
• resolution set to invalid

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?

### Changed 3 years ago by starsareblue

• status changed from closed to reopened
• resolution invalid deleted

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
```

### Changed 3 years ago by pramsey

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

### Changed 3 years ago by 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

### Changed 3 years ago by 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?

### Changed 3 years ago by 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.

### Changed 3 years ago by 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.

### Changed 23 months ago by strk

• owner changed from geos-devel@… to sgillies
• status changed from reopened to new
• component changed from Default to SWIG Python
• milestone set to 3.2.3

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.

### Changed 23 months ago by 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?

### Changed 23 months ago by strk

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

### Changed 20 months ago by 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 ?

### Changed 20 months ago by 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 :/

### Changed 20 months ago by strk

• milestone changed from 3.2.3 to GEOS Future

### Changed 20 months ago by strk

• milestone changed from GEOS Future to 3.3.1

plus 464 and 467 as side minor tweaks.

The testcase is interestingly expecting the result of not interpolating.

### Changed 20 months ago by strk

• milestone changed from 3.3.1 to 3.4.0

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

### Changed 20 months ago by 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.

### Changed 12 months ago by strk

• status changed from new to closed
• resolution set to fixed

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

### Changed 12 months ago by 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.

Note: See TracTickets for help on using tickets.