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.

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Jon
  • 2,894
  • 2
  • 19
  • 35
  • please provide a short WKT that demonstrates the issue – Ian Turton Oct 06 '18 at 08: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 06 '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 06 '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 06 '18 at 14:47
  • "and euclidean distance between lat/lon is fine with me": Did you realize what you're writing ? – gene Oct 06 '18 at 16:29
  • Indeed, I do. It is not a distance that is meaningful in terms of length. If it makes it easier for people to understand my question, pretend that the points aren't lon/lat, but just x,y coordinates on a plane. I intentionally left out that they were coordinates because it's irrelevant to the problem I'm asking about. – Jon Oct 06 '18 at 17:09
  • If the issue goes away when you reduce the number of points it is probably a validity issue not a bug. – Ian Turton Oct 06 '18 at 17:36
  • @Ian, well, remember that the snapping does work, but only when the tolerance is well beyond what should be required. Additionally, is_valid returns True for both g1 and g2. – Jon Oct 06 '18 at 19:47

0 Answers0