¿Por qué el ajuste de Shapely (instantáneas GEO) no funciona como se esperaba?

Estoy tratando de ajustar dos líneas entre sí usando Shapely / Geopandas, pero el resultado del ajuste es muy extraño. Lo intenté :

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")

Y obtuve este resultado:

líneas1 = líneas rojas

líneas2 = líneas negras

Después del ajuste (con 14 como tolerancia): las líneas azules son el resultado del ajuste

En este caso, las líneas se ajustan correctamente

Otro ejemplo donde no funcionó como se esperaba: (antes del ajuste)

Y aquí está el resultado después del ajuste. Solo una parte se ajusta a la línea negra (el lado sur). Aunque las líneas originales están bastante cerca y dentro de los 14 pies

Si aumento la tolerancia obtengo un resultado incorrecto, algo como esto (después de definir 20 como la tolerancia del ajuste, el resultado es la línea verde):

¿Alguna idea de por qué el ajuste no funciona correctamente? ¿Alguna sugerencia sobre cómo resolver este problema?

La shapely.ops.snapfunción se ajusta solo a los vértices de las geometrías.

Vea la siguiente ilustración. A la izquierda, el vértice rojo está dentro de la tolerancia de ajuste al vértice azul, por lo que se ajustará. A la derecha, el vértice rojo está fuera de la tolerancia de ajuste (¡a pesar de estar más cerca del borde!).

Shapely no proporciona un algoritmo para ajustar vértices a los bordes. Sin shapely.ops.nearest_pointsembargo , no debería ser demasiado difícil escribir uno usando . Algo como esto (no probado, y no especialmente eficiente):

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)