I have two polylines whose endpoints are near to each other, but do not overlap. I've put them into this shapefile. I want to use shapely's ops.snap() to snap their endpoints, but it's not working as expected.

from shapely import ops
import geopandas as gpd

wontsnap = gpd.read_file(path_to_wontsnap)
g1 = wontsnap.geometry.values[0]
g2 = wontsnap.geometry.values[1]

If I check the distance between the geometries, I get

Out[661]: 0.0006421029901831172

This agrees (roughly) with the measure tool in QGIS:

enter image description here

Now if I try to snap the two geometries with a tolerance that is greater than the distance between their endpoints:

snapped = ops.snap(g1, g2, 0.001)

The snapped object contains only the g1 geometry. I.e., they don't snap.

If I increase the tolerance by 10x:

snapped = ops.snap(g1, g2, 0.01)

The snapped object includes both geometries, snapped correctly.

If I trim the g1 and g2 geometries to the last and first 10 points, respectively, snapping works as expected.

I have tried various tolerances between 0.001 and 0.01--the geometries will snap at 0.003 but not 0.002. At this point I can only assume this is a bug in shapely--any ideas?

I have posted this as a bug report in the shapely repository. It appears to be a bug.

closed as off-topic by PolyGeo Jul 3 at 23:35

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "Questions that relate to the business of vendors and product teams, such as licensing, pricing, release dates, submission of bug reports and enhancement requests, etc, are off-topic here and should be directed instead to the official software support teams." – PolyGeo
If this question can be reworded to fit the rules in the help center, please edit the question.

  • please provide a short WKT that demonstrates the issue – Ian Turton Oct 6 '18 at 8:54
  • The projection of your shapefile is WGS84 and the coordinates are in longitude, latitude (angles in degrees). Shapely uses a cartesian plane system for computing geometries (distance = euclidean distance). That means that all your distance calculations are wrong without reprojecting your layer (What units are used by GeoPandas/Shapely area and distance functions?) – gene Oct 6 '18 at 10:03
  • @Gene, then why does shapely's distance() function return the correct distance? Shapely doesn't care about your projection, and euclidean distance between lat/lon is fine with me. I'm not trying to compute a meaningful distance, I'm trying to snap points that are within a threshold distance. Cartesian is fine. – Jon Oct 6 '18 at 14:45
  • @Ian, I can't do it, that's why I uploaded the shapefile. Like I said, when I try to snip the polylines down to 10 points, the snapping works fine. So I can't provide a shortened version of the polylines, which are hundreds of points. I can paste in a wkt of the entire thing if you want, but they're already provided as a shapefile so I'm not inclined to do that. – Jon Oct 6 '18 at 14:47
  • "and euclidean distance between lat/lon is fine with me": Did you realize what you're writing ? – gene Oct 6 '18 at 16:29