Fastest way to extract the neighbours of a GPS point among millions of other points?

Question

I have millions of GPS point located in the same city. I would like to make spatial bins that contain adjacent points together (GPS points that are within a specific range of each other). I have a shapefile of that city and the GPS points located in it. Given the enormous amount of data, I'm looking for the least-cost way to do so.

I have already considered making a distances matrix but its size would grow quadratically with the size of the dataset.

I have also considered making a bounding box around that point and compare its top-left and bottom-right coordinates to all of the points in the dataset but it proved high temporal complexity.

I'm quite new to GeoPandas and I was wondering if there's a way to extract neighbours of a given point from the plot itself? Like an interactive map or something like this.

PS: The plot below contains a sample of data (not all of it)

Answer 1

There are many clustering algorithms you can use, for example DBSCAN:

Very large n_samples, medium n_clusters

import geopandas as gpd
import pandas as pd
from sklearn.cluster import DBSCAN
import datetime
start = datetime.datetime.now()


in_shapefile = r'/home/bera/Desktop/gistemp/1000000_random_points.shp'
out_shapefile = r'/home/bera/Desktop/gistemp/1000000_random_points_clustered.shp'

df = gpd.read_file(in_shapefile)

df['xcoord'] = df.geometry.apply(lambda x: x.x) #Add a x coordinate column
df['ycoord'] = df.geometry.apply(lambda x: x.y) #and y

coords = df[['xcoord','ycoord']].values #Create a numpy array of x and y coordinate
db = DBSCAN(eps=30, min_samples=1).fit(coords) #cluster them. eps is distance between points, adjust it
cluster_labels = pd.Series(db.labels_).rename('cluster')  #.labels_ creates a numpy array. Convert to a series and rename it

df = pd.concat([df, cluster_labels], axis=1) #Add the cluster series to your dataframe

df2 = df.loc[(df.xcoord<666148) & (df.ycoord>6587010)] #I create a small sample to visualize in QGIS by selecting points in the upper left corner
df2[['geometry','cluster']].to_file(out_shapefile.replace('.shp','_sample_.shp')) #Export

stop = datetime.datetime.now()
print((stop-start).seconds, ' seconds')

With 1 million points it takes 50 s.

My points are completely random so alot of them ends up in one cluster: