Line Intersection produces a point which results in invalid value in line_locate_point

[sharon92]

Expected behavior and actual behavior.

Creating an intersection point of a river cross-section [LineString with Z values] with the river [LineString 2D]
and then projecting the intersection point on the same river cross-section produces invalid value in line_locate_point.

Coordinates of river:

river2d = LineString([[535947.0545920385, 5384870.337515112],
 [535948.4767302983, 5384868.485428078],
 [535949.2374089044, 5384866.368757173],
 [535949.4314629119, 5384864.317398101],
 [535949.9447410172, 5384856.532624327],
 [535950.4611092675, 5384851.452841926],
 [535951.2879338376, 5384844.474442553]])

Coordinates of cross-section:

segment = LineString([(535947.7556526819, 5384855.732786809, 368.795007324),
 (535948.0260501119, 5384855.83494383, 368.795007324),
 (535948.0260501119, 5384855.83494383, 368.244989014),
 (535948.1377846444, 5384855.877157496, 366.774987793),
 (535948.1564934427, 5384855.884225741, 366.665002441),
 (535949.9360188602, 5384856.556536182, 366.665),
 (535951.3422283962, 5384857.087806678, 366.655),
 (535951.8095666496, 5384857.264368578, 366.675012207),
 (535951.8183376329, 5384857.267682284, 366.765008545),
 (535951.8511331676, 5384857.280072543, 368.165002441),
 (535951.8839287009, 5384857.292462801, 368.554986572)])

intepolated_point = river2d.intersection(segment)
print(intepolated_point.coords[0])
>> (535949.9429907533, 5384856.559170186, 366.66495042066657)

segment.project(interpolated_point)

  File "C:\Python311\Lib\site-packages\shapely\linear.py", line 89, in line_locate_point
    return lib.line_locate_point(line, other)

RuntimeWarning: invalid value encountered in line_locate_point

Steps to reproduce the problem.

(For example, a script with required data)

Operating system

Windows 10

Shapely version and provenance

2.0.1 from PyPI using pip

[mwtoews]

Issue confirmed with main shapely branch with latest GEOS. Here is a short version:

from shapely import geos, wkt

print(geos.geos_version_string)
A = wkt.loads("LINESTRING (0 0, 1 1, 1 1, 2 2)")
B = wkt.loads("POINT (0 1)")
print(A.project(B))

3.11.2-CAPI-1.17.2
/path/to/src/shapely/shapely/linear.py:90: RuntimeWarning: invalid value encountered in line_locate_point
return lib.line_locate_point(line, other)
0.7071067811865476

The issue is the coordinate that has repeated XY coords giving a zero segment distance in 2d (Z dim is ignored). The issue is likely deep within C++ of the GEOS library.

[jorisvandenbossche]

produces invalid value in line_locate_point.

Just to be clear, the "invalid value" is only a warning (that needs to be avoided/suppressed in GEOS). But is the actual result correct?

[sharon92]

@jorisvandenbossche yes the value is correct. I meant to say that it produces the warning invalid value in line_locate_point.

@mwtoews is right, it ignores the Z dimension. Perhaps duplicate 2d coordinates are not allowed in GEOS

[chau-intl]

Would it be possible from within shapely to expose a bit more of the stack trace?

I see the same warning but it only shows the part of the stack trace which is related to shapely and not the rest of the stack trace which would be helpful to direct me to the place in MY code where the warning is related to. If this is possible it would be easier to check if the result is still valid despite the warning.

[dr-jts]

@mwtoews are you able to reproduce this in GEOS directly?

[mwtoews]

Hmm, I'm not seeing much with:

./bin/geosop -a "LINESTRING (0 0, 1 1, 1 1, 2 2)" -b "POINT (0 1)" project

it could potentially be within https://github.com/shapely/shapely/blob/2.0.2/src/ufuncs.c (look for line_locate_point); further investigation is needed to pinpoint it.

[mwtoews]

This was fixed with GEOS 3.12.0.

While I'm not sure which commit was the fix, I've added a test for it libgeos/geos@34b29f8