4

GeoPandas has this really convenient tool sjoin_nearest() to find the nearest features to, say, each point in a dataframe. It can also return the distance:

distances = gpd.sjoin_nearest(points_gdf, lines_gdf, distance_col="distances")

But what if I want to find out the exact coordinates of the nearest node? I guess I could do a quite expensive operation exploding the lines to points, but as the sjoin_nearest operation found the distance to the nearest line (and did it pretty fast, by Python and GeoPandas standards) somehow, I'm guessing there must be a more efficient way?

1 Answer 1

4

You can use shapely.ops nearest points:

Returns a tuple of the nearest points in the input geometries. The points are returned in the same order as the input geometries.

import geopandas as gpd
from sqlalchemy import create_engine
from shapely.ops import nearest_points

#Read some local PostGIS data into two data frames    
con = create_engine(url="postgresql://postgres@localhost/data")
    road = gpd.read_postgis(sql="select geom from lmv.ok_vl_riks", con=con) #Line layer with roads
station = gpd.read_postgis(sql="select geom from lmv.ok_js_riks", con=con) #Point layer with railway stations
station["stationid"] = range(0, station.shape[0]) #Create a unique id, to use for dropping duplicates

#Join the points to the lines
station["pointgeom"] = station.geometry #Save the point geometry or it is lost in spatial join
sj = gpd.sjoin_nearest(left_df=road, right_df=station, how="inner", max_distance=500, distance_col="distance")

#Duplicate points/rows are created when there are multiple join features within max_distance. Drop the duplicates, keep the nearest
sj = sj.sort_values(by=["stationid", "distance"], ascending=True).drop_duplicates(subset="stationid", keep="first")
sj["nearest_point"] = sj.apply(lambda x: nearest_points(g1=x["geom"], g2=x["pointgeom"])[0], axis=1) #Calculate nearest points on line and the joined point. Keep index [0], which is the point on the line nearest to the joined point

#Plot a subset
ymin, xmin = 6945207.5,645918.3
ymax, xmax = 6951569.9,652789.3
subset = sj.cx[xmin:xmax, ymin:ymax]
ax = subset.plot(figsize=(20,20), color="black")
subset["pointgeom"].plot(ax=ax, marker = "s", color="red")
subset["nearest_point"].plot(ax=ax, marker="*", color="blue", markersize=200)

Black lines = roads, red squares = railway stations, blue stars = nearest point on line to red square.

enter image description here

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.