Find pairs of near features with geopandas, fiona, GDAL etc

Question

I have a feature class of points, a few thousands I have now, but I expect to deal with hundreds of thousands of points. I would like find pair of points that is within a specific distance. The points I have is typically spread apart spatially so I don't expect to find many pairs. How I can find those pairs.

Right now the points are originally from ESRI shapefile, and read into geopandas geodataframe. I tried to tbl = df['geometry'].distance(df['geometry']) see if I can get distance of all pairs.

But I got the error "python: GeometryComponentFilter.cpp:34: virtual void geos::geom::GeometryComponentFilter::filter_ro(const geos::geom::Geometry*): Assertion `0' failed.". Beside

But, I don't want to get full set of distance as I am only interested in points that are close to each other.

I have a version of code which uses ArcGIS, and GenerateNearTable_analysis is what I used for the purpose.

I want to run code with open source code now.

My solution has to be python and should run on Linux.

Answer 1

You can use Geopandas with spatial index nearest:

import geopandas as gpd
import numpy as np

#Create a point dataframe with 1000 points
np.random.seed(42)
coords = np.random.uniform(low=0, high=1000, size=(1000, 2))
df = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=coords[:,0], y=coords[:,1]))
df["id"] = range(df.shape[0])

#Create the spatial index and query nearest
si = df.sindex

(in_geometry_indices, tree_geom_indices), distance =  si.nearest(df.geometry, 
                                                           exclusive=True,
                                                           return_distance=True)
#Point with index 0 nearest point is index 69, at a distance of 23.7:
print(in_geometry_indices[0])
#[0]
print(tree_geom_indices[0])
#[69]
print(distance[0])
#[23.72531493]

#Create a dataframe from the three arrays
df2 = gpd.pd.DataFrame(data=[in_geometry_indices, tree_geom_indices, distance]).T
df2.columns = ["in_index", "si_index", "distance"]

print(df2.head())
#    in_index  si_index   distance
# 0       0.0      69.0  23.725315
# 1       1.0     231.0  20.765461
# 2       2.0     700.0  16.020174
# 3       3.0     264.0   2.243012
# 4       4.0     699.0  16.082228

#Find all points that are within 50 distance from eachother
df3 = df2.loc[df2.distance<=50].copy()

#Merge with the start df by id-in_index or id-si_index
point_pairs = gpd.pd.concat([gpd.pd.merge(df, df3, left_on='id', right_on='in_index'),
           gpd.pd.merge(df, df3, left_on='id', right_on='si_index')])