Shapely – Why Shapely’s Snapping (GEO Snaps) is Not Working as Expected

geopandasshapelysnapping

I'm trying to snap two lines to each other using Shapely/Geopandas but the result of the snapping is very odd.
I tried : 

import geopandas as gpd
from shapely.geometry import *
from shapely.ops import snap

lines1 = gpd.GeoDataFrame.from_file('lines1.shp')
lines1 = lines1.to_crs({'init': 'epsg:2227'})
lines2 = gpd.GeoDataFrame.from_file('lines2.shp')
lines2 = lines2.to_crs({'init': 'epsg:2227'})
res = lines1
lines2_union = lines2.geometry.unary_union
res.geometry = res.geometry.apply(lambda x: snap(x, lines2_union, 14))
res.to_file('result.shp', driver="ESRI Shapefile")

And got this result :

lines1 = red lines

lines2 = black lines

Before the snapping

After the snapping(with 14 as a tolerance) :
the blue lines are the result of the snapping

In this case the lines are correctly snapped
After the snapping

Another example where it didn't work as expected: (before the snapping)
Before the snapping

And here's the result after the snapping. Only a part is snapped to the black line (the south side). Although the original lines are pretty close and within the 14 feets
After the snapping

If I increase the tolerance I get a wrong output, something like this (after defining 20 as the tolerance of the snapping, the green line is the result):

After 20 as a tolerance

Any ideas on why the snapping is not working properly?
Any suggestions on how to solve this problem?

Best Answer

The shapely.ops.snap function snaps to the vertices of geometries only.

See the illustration below. On the left, the red vertex is within the snapping tolerance to the blue vertex, so it will snap. On the right the red vertex is outside the snapping tolerance (despite being closer to the edge!).

snapping tolerance visualisation

Shapely doesn't provide an algorithm for snapping vertices to edges. It shouldn't be too difficult to write one using shapely.ops.nearest_points though. Something like this (not tested, and not especially efficient):

from shapely.ops import nearest_points

def snap(g1, g2, threshold):
    coordinates = []
    for x, y in g1.coords:  # for each vertex in the first line
        point = Point(x, y)
        p1, p2 = nearest_points(point, g2)  # find the nearest point on the second line
        if p1.distance(p2 <= threshold):
            # it's within the snapping tolerance, use the snapped vertex
            coordinates.append(p2.coords[0])
        else:
            # it's too far, use the original vertex
            coordinates.append((x, y))
    # convert coordinates back to a LineString and return
    return LineString(coordinates)
Related Question