shapely / shapely

Manipulation and analysis of geometric objects
https://shapely.readthedocs.io/en/stable/
BSD 3-Clause "New" or "Revised" License
3.85k stars 565 forks source link

Floating point issue with splitter points #1058

Open AugustBrenner opened 3 years ago

AugustBrenner commented 3 years ago

Expected behavior and actual behavior.

When using ops.nearest_points to get a splitter Point to split a LineString, the returned nearest_point on the line will not split the LineString because the nearest point has a distance of ~1.0E-16 away from the line.

Steps to reproduce the problem.

nearest_point, _ = nearest_points(edge.line, event_point) edge_segments = list(split(edge.line, nearest_point))

Operating system

(For example, Mac OS X 10.15.7)

Shapely version and provenance

1.7.1 installed from PyPI using poetry

tomplex commented 3 years ago

I wonder if this is related to #841? @AugustBrenner, can you share the geometries which are specifically causing this problem so we can test it?

AugustBrenner commented 3 years ago

I don't think this is related.

POINT (-111.9737684950623 40.38074884410073)

LINESTRING (-111.9641265714514 40.38234418201984, -111.9655417114647 40.38210678889413, -111.9679690859342 40.38169955832259, -111.9679690859342 40.38169955832259, -111.9685463679903 40.38160274258401, -111.9704676571729 40.38128086854964, -111.9704676571729 40.38128086854964, -111.9707731021402 40.381229746174, -111.9726739229372 40.38091143077548, -111.9726739229372 40.38091143077548, -111.9728509113205 40.3808817437128, -111.9742219007016 40.38068317441761)

Jeremiah-England commented 3 years ago

I've encountered similar issues splitting linestrings by points as well. I'm not sure why the nearest_points method sometimes returns points for which this is true:

https://github.com/Toblerity/Shapely/blob/4fcf83523dbe96fcc51fe12b6826a6ecdc4b2a66/shapely/ops.py#L450

That is, the interior of the splitter (a Point in this case) does not touch the interior of the line. It really seems like this should work after nearest_points. May be worth looking into further.

One solution is to snap the LineString to Point which results from nearest_points (a tip I got from here). This inserts the Point as a vertex of the LineString. It should split in that situation.

There is a very similar issue in #934 and another solution is offered in #630.

At least we could update our documentation to note that you can snap the LineString to the Point. Currently it says to snap the Point to the LineString, but if the vertices of the LineString are sparse where you want to split this doesn't work well.

It may be convenient to snap the splitter with low tolerance to the geometry. For example in the case of splitting a line by a point, the point must be exactly on the line, for the line to be correctly split. (link)

Edit: Another related issue: #561

AugustBrenner commented 3 years ago

Ah, yea. I also tried this using snap, as well as interpolate(project()). All methods on the shapely contract boundary that return points on the line don't seem to be able to split that line for complex floating point coordinates.

I imagine some sort of tolerance mechanism needs to be added for this to work.

My current work around is using a buffered point as the splitter, which returns 2 to 3 LineString segments depending on the geometry and the method I use (nearest point, snap, interpolate-project), then adding a conditional for 2 and 3 segments respectively, and reconstructing the valid point split segments using the original splitter point, and the first/last LineStrings of the GeometryCollection.

This could actually be a useful and robust mechanism to add a tolerance to a split. In several thousand of the geometries I have yet to come across a distance that was more than 1E-15 when snapping/nearest_point/interpolate-project-ing. I currently set my buffer to 1E-14 when doing my workaround hack.

edge_segments = list(split(edge.line, nearest_point.buffer(1.0E-14)))
entry_segment = LineString(list(edge_segments[0].coords) + list(nearest_point.coords))
etc...
timovanasten commented 3 years ago

I also have this issue.

Here is an example case that does not behave as expected in def _split_line_with_point(line, splitter) with: line: LINESTRING (2.1847328 41.401595, 2.1848937 41.4014736, 2.184928 41.4014477, 2.1849834 41.4014026) splitter: POINT (2.1848937 41.4014736) (second coordinate of the line)

Returns: LINESTRING (2.1847328 41.401595, 2.1848937 41.4014736, 2.1848937 41.4014736) <- last coordinate duplicated LINESTRING (2.1848937 41.4014736, 2.184928 41.4014477, 2.1849834 41.4014026)