Snap points to polygon/lines in python shapely?

Question

I have a shapefile of point features

And a shapefile of polygons features

I would like to snap these points to the polygons, as they are often just outside the polygons and cannot be spatially joined.

How do I snap the points to the polygon edge if they are within a range of around 100m or so? One thing I tried was to convert the polygon boundary into a line feature and then snap(), but it didn't modify the geometry no mater what the tolerance was

res = point
lines2_union = parcel_as_lines.geometry.unary_union
res.geometry = point.geometry.apply(lambda x: snap(x, lines2_union, 999999999999999999999999999))

Answer 1

sjoin_nearest to find any polygons within your tolerance. Measure distance from point to polygon. If the distance is below your tolerance and more than 0, snap. Else leave the point unchanged:

import geopandas as gpd
from shapely.ops import snap, nearest_points

poly = gpd.read_file(r"/home/bera/Desktop/GIStest/1400_buildings.geojson")
point = gpd.read_file(r"/home/bera/Desktop/GIStest/500_points.geojson")
snapdist = 100 #The distance from point to polygon to search and snap
output = r"/home/bera/Desktop/GIStest/500_points_snapped.geojson"

#To measure distances in meters, the data needs to be in a projected crs with meters as units.
#  My data is in Sweden so I choose EPSG:3006
originalcrs = point.crs
poly = poly.to_crs(3006)
point = point.to_crs(3006)

#Create unique ids
poly["polyid"] = range(poly.shape[0]) 
point["pointid"] = range(point.shape[0])

#For each point find the closest polygon(s)
#  If there are more than one polygon within snapdistance, the points will be duplicated
poly["polygeom"] = poly.geometry #Save the geometry, or it is lost in sjoin_nearest
sj = gpd.sjoin_nearest(left_df=point, right_df=poly, how="left", max_distance=snapdist)

#Measure distances. If there is no polygon within snapdistance, return None
sj["distance"] = sj.apply(lambda x: x.geometry.distance(x.polygeom) if x.polygeom is not None else None, axis=1)

#Sort by distance and drop any duplicated points
sj = sj.sort_values(by=["pointid","distance"], ascending=True, na_position="last")
sj = sj.drop_duplicates(subset="pointid", keep="first")

print(sj[["pointid","polyid","distance"]].head(2))
#    pointid  polyid  distance
# 0        0  1422.0  47.02335
# 1        1     NaN       NaN

#So point 0 should snap to polygon 1422 if the distance < snapdist. 
#Point 1 should not snap. Any other points with 0 distance (they are already intersecting a polygon) should
#  stay unchanged

#Create/find the point on the polygon to snap to
sj["nearestpoint"] = sj.apply(lambda x: nearest_points(x.geometry, x.polygeom)[1] if (x.polygeom is not None or x.distance==0) else None, axis=1)
print(sj.head(2))
#    id                        geometry  ...  distance                    nearestpoint
# 0   0  POINT (680677.291 6577536.845)  ...  47.02335  POINT (680703.749 6577497.971)
# 1   1  POINT (682069.722 6590459.711)  ...       NaN                            None

sj["geometry"] = sj.apply(lambda x: snap(x.geometry, x.nearestpoint, snapdist) 
                          if x.nearestpoint is not None else x.geometry, axis=1) #Snap
#Export the result
sj = sj.to_crs(originalcrs)
sj = sj[[c for c in sj.columns if c in point.columns]] #Drop processing columns
sj.to_file(output)