1

I use GeoPandas to link a point feature to the nearest line feature with the code below. I thought I might share it for anyone to comment on it or offer a better solution. Pay attention to the comment "sindex i is NOT the df_line index". It makes sense now, but caused me a lot of headaches (and nasty bugs) before I figured that one out:

from math import hypot
import geopandas as gpd
import pandas as pd

# link_nearest_line
def link_nearest_line(point_geom, df_line):
    x, y = point_geom.x, point_geom.y
    sindex = df_line.sindex
    i = next(sindex.nearest((x, y, x, y)))
    line_geom = df_line.iloc[i].geometry
    snap_geom = line_geom.interpolate(line_geom.project(point_geom))
    snap_dist = hypot(x - snap_geom.x, y - snap_geom.y)
    nearest = []
    for i in sindex.intersection((x - snap_dist, y - snap_dist, x + snap_dist, y + snap_dist)):
        line_geom = df_line.iloc[i].geometry
        snap_geom = line_geom.interpolate(line_geom.project(point_geom))
        snap_dist = hypot(x - snap_geom.x, y - snap_geom.y)
        # sindex i is *NOT* the df_line index
        # map it using df_line.index[i] before returning the result !!!
        nearest.append((df_line.index[i], snap_geom, snap_dist))
    return min(nearest, key=lambda t: t[2])

# main
df_point = gpd.read_file('point.shp')
df_line = gpd.read_file('line.shp')

df_result = df_point.apply(
    lambda row: pd.Series(link_nearest_line(row.geometry, df_line),
                          index=('line_index', 'projected_point_geom', 'distance')),
    axis='columns', result_type='expand')
df_point = pd.concat([df_point, df_result], axis='columns')

# geodataframe df_point now has fields:
# line_index: the index of the nearest line
# projected_point_geom: geometry of the point projected onto the nearest line
# distance: distance from point to projected_point
  • 4
    I'm voting to close this question as off-topic because it's an attempt to start a code review discussion. – Vince Dec 18 '19 at 1:18

Browse other questions tagged or ask your own question.