python - Calculate Distance to Nearest Feature with Geopandas

Question

I'm looking to do the equivalent of the ArcPy Generate Near Table using Geopandas / Shapely. I'm very new to Geopandas and Shapely and have developed a methodology that works, but I'm wondering if there is a more efficient way of doing it.

I have two point file datasets - Census Block Centroids and restaurants. I'm looking to find, for each Census Block centroid, the distance to it's closest restaurant. There are no restrictions in terms of same restaurant being the closest restaurant for multiple blocks.

The reason this becomes a bit more complicated for me is because the Geopandas Distance function calculates elementwise, matching based on index. Therefore, my general methodology is to turn the Restaurants file into a multipoint file and then set the index of the blocks file to all be the same value. Then all of the block centroids and the restaurants have the same index value.

import pandas as pd
import geopandas as gpd
from shapely.geometry import Polygon, Point, MultiPoint

Now read in the Block Centroid and Restaurant Shapefiles:

Blocks=gpd.read_file(BlockShp)
Restaurants=gpd.read_file(RestaurantShp)

Since the Geopandas distance function calculates distance elementwise, I convert the Restaurant GeoSeries to a MultiPoint GeoSeries:

RestMulti=gpd.GeoSeries(Restaurants.unary_union)
RestMulti.crs=Restaurants.crs
RestMulti.reset_index(drop=True)

Then I set the index for the Blocks equal to 0 (the same value as the Restaurants multipoint) as a work around for the elementwise calculation.

Blocks.index=[0]*len(Blocks)

Lastly, I use the Geopandas distance function to calculate the distance to the nearest restaurant for each Block centroid.

Blocks['Distance']=Blocks.distance(RestMulti)

Please offer any suggestions on how any aspect of this could be improved. I'm not tied to using Geopandas or Shapely, but I am looking to learn an alternative to ArcPy.

Thanks for the help!

Answer 1

If I understand correctly your issue, Blocks and Restaurants can have very different dimensions. For this reason, it's probably a bad approach to try to force into a table format by reindexing.

I would just loop over blocks and get the minimum distance to restaurants (just as @shongololo was suggesting).

I'm going to be slightly more general (because I already have this code written down) and do a distance from points to lines, but the same code should work from points to points or from polygons to polygons. I'll start with a GeoDataFrame for the points and I'll create a new column which has the minimum distance to lines.

%matplotlib inline
import matplotlib.pyplot as plt
import shapely.geometry as geom
import numpy as np
import pandas as pd
import geopandas as gpd

lines = gpd.GeoSeries(
    [geom.LineString(((1.4, 3), (0, 0))),
        geom.LineString(((1.1, 2.), (0.1, 0.4))),
        geom.LineString(((-0.1, 3.), (1, 2.)))])

# 10 points
n  = 10
points = gpd.GeoSeries([geom.Point(x, y) for x, y in np.random.uniform(0, 3, (n, 2))])

# Put the points in a dataframe, with some other random column
df_points = gpd.GeoDataFrame(np.array([points, np.random.randn(n)]).T)
df_points.columns = ['Geometry', 'Property1']

points.plot()
lines.plot()

Now get the distance from points to lines and only save the minimum distance for each point (see below for a version with apply)

min_dist = np.empty(n)
for i, point in enumerate(points):
    min_dist[i] = np.min([point.distance(line) for line in lines])
df_points['min_dist_to_lines'] = min_dist
df_points.head(3)

which gives

    Geometry                                       Property1    min_dist_to_lines
0   POINT (0.2479424516236574 2.944916965334865)    2.621823    0.193293
1   POINT (1.465768457667432 2.605673714922998)     0.6074484   0.226353
2   POINT (2.831645235202689 1.125073838462032)     0.657191    1.940127

---- EDIT ----

(taken from a github issue) Using apply is nicer and more consistent with how you'd do it in pandas:

def min_distance(point, lines):
    return lines.distance(point).min()

df_points['min_dist_to_lines'] = df_points.geometry.apply(min_distance, df_lines)

EDIT: As of at least 2019-10-04 it seems that a change in pandas requires a different input in the last code block, making use of the args parameters in .apply():

df_points['min_dist_to_lines'] = df_points.geometry.apply(min_distance, args=(df_lines,))