Measuring distance between two lat long points in 1000's in Python

Sort both dataframes by latitude. This will allow not to calculate distance between pairs of points if their latitudes differ significantly. In the best case you can get speedup of 580 times.

The idea is that you loop through rows of df1 and for each row of that array find left and right indices of the second array, which latitudes are not far from that row

df1.sort_values(by='latitude')
df2.sort_values(by='latitude')
n1 = df1.shape[0]
n2 = df2.shape[0]
left = 0
right = 0
threshold = 400
lat_threshold = threshold / radius # latitude difference that corresponds to 400 m
for i in range(n1):
    row1 = df1.iloc[[i]]
    lat1 = row1['latitude']
    lon1 = row1['longitude']
    while left < n2 and df2.iloc[[left]]['latitude'] < lat1 - lat_threshold:
        left += 1
    while right < n2 and df2.iloc[[right]]['latitude'] < lat1 + lat_threshold:
        right += 1
    for j in range(left, right):
        row2 = df2.iloc[[j]]
        lat2 = row2['latitude']
        lon2 = row2['longitude']
        k = distance((lat1, lon1), (lat2, lon2))
        if (k <= threshold):
            ll.append(row2)
        

You may try this using geopandas :

import geopandas as gpd
import pandas as pd
import pyproj

df1 = pd.read_csv("/home/ubuntu/maid80.csv")
df2 = pd.read_csv("/home/ubuntu/iodr.csv")

gdf1 = gpd.GeoDataFrame(df1, geometry=gpd.points_from_xy(df1['longitude'], df1['latitude']), crs=pyproj.CRS.from_epsg(4326))
gdf2 = gpd.GeoDataFrame(df2, geometry=gpd.points_from_xy(df2['longitude'], df2['latitude']), crs=pyproj.CRS.from_epsg(4326))

radius = 400
for gdf in [gdf1, gdf2]:
  gdf.to_crs(pyproj.CRS.from_epsg(3857), inplace=True)

gdf1['geometry'] = gdf1['geometry'].buffer(radius)
gdf2['IS_WITHIN_400M'] = 1

gdf = gpd.sjoin(gdf1, gdf2['geometry'], how='left')
print(gdf[gdf.IS_WITHIN_400M_right==1].head())

Some explanations :

Geopandas will allow you to use a GeoDataFrame, on which you can "buffer" your points using a radius (really fast). The points_from_xy function is also pretty fast and will allow you to construct those objects efficiently.

The sjoin method (stands for spatial join) is also fast. I suspect this as something to do with algorithm including bounding boxes and/or sorting coordinates... I've had some good results using this method.

Warning :

I projected the datasets into EPSG 3857, which is global AND has cartesian coordinates (in meters). Regarding to any "real" project, you have to chose the projection carefully (ie chose the best "metric-system-friendly" projection in your area) to avoid any distorsion of the buffer...

You can only use numpy functions for your distance function and vectorize it. That should be a lot faster:

from __future__ import print_function

import pandas as pd
import math

import numpy as np


def distance(origin: pd.DataFrame, lat2, lon2):
'''Measure distance not for a pair but for the whole dataframa against one point'''
    lat1 = origin['latitude']
    lon1 = origin['longitude']
    radius = 6371 * 1000  # km
    dlat = np.radians(lat2 - lat1)
    dlon = np.radians(lon2 - lon1)
    a = np.sin(dlat / 2) * np.sin(dlat / 2) + np.cos(np.radians(lat1)) \
        * np.cos(np.radians(lat2)) * np.sin(dlon / 2) * np.sin(dlon / 2)
    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
    d = radius * c
    return d


def main():
    df1 = pd.read_csv("/home/ubuntu/maid80.csv")
    df2 = pd.read_csv("/home/ubuntu/iodr.csv")
    ll = []
    for index, row in df2.iterrows():
        #because you can test the whole dataframe gainst one point you can remove    one loop.
        mask= distance(df1,row['latitude'],row['longitude'])<400.0
        ll.extend(df1.index[mask].to_list()) #only add points to the list where the distance is <400

    
    myset = set(ll)
    print(myset)

Maybe you have to switch the dataframes. I dont know which one is the one you want to collect the ids from.

You can use BallTree with HaversineDistance metric. At first construct the tree with coordinates from the first table and then query coordinates from the second table from that tree

from sklearn.neighbors import BallTree, DistanceMetric

radius = 6371 * 1000
max_distance = 400 / radius

# ensure that format of array is [latitude, longitude]
rows1 = np.deg2rad(df1[['latitude', 'longitude']].to_numpy())
rows2 = np.deg2rad(df2[['latitude', 'longitude']].to_numpy())

# haversine metric accepts latitude and longitude only in radians and returns distance
# on unit sphere
tree = BallTree(rows1, metric=DistanceMetric.get_metric('haversine'))

count = tree.query_radius(rows2, r=max_distance, count_only=True)
print(df2['id'].iloc[np.nonzero(count)[0]])