1

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

g1.distance(g2)
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.

8
  • 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