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')])