Question: I have a dataset of about 20,000 records that represent global cities of population > 20,000. I have estimated radius which more or less describes the size of the city. It's not exactly accurate but for my purposes it will be enough. that I'm loading it into my Panda dataframe object. Below is the sample name_city,country_code,latitude,longitude,geohash,estimated_radius,population Vitry-sur-Seine,FR,48.78716,2.40332,u09tw9qjc3v3,1000,81001 Vincennes,FR,48.8486,2.43769,u09tzkx5dr13,500,45923 Villeneuve-Saint-Georges,FR,48.73219,2.44925,u09tpxrmxdth,500,30881 Villejuif,FR,48.7939,2.35992,u09ttdwmn45z,500,48048 Vigneux-sur-Seine,FR,48.70291,2.41357,u09tnfje022n,500,26692 Versailles,FR,48.80359,2.13424,u09t8s6je2sd,1000,85416 Vélizy-Villacoublay,FR,48.78198,2.19395,u09t9bmxdspt,500,21741 Vanves,FR,48.82345,2.29025,u09tu059nwwp,500,26068 Thiais,FR,48.76496,2.3961,u09tqt2u3pmt,500,29724 Sèvres,FR,48.82292,2.21757,u09tdryy15un,500,23724 Sceaux,FR,48.77644,2.29026,u09tkp7xqgmw,500,21511 Saint-Mandé,FR,48.83864,2.41579,u09tyfz1eyre,500,21261 Saint-Cloud,FR,48.84598,2.20289,u09tfhhh7n9u,500,28839 Paris,FR,48.85341,2.3488,u09tvmqrep8n,12000,2138551 Orly,FR,48.74792,2.39253,u09tq6q1jyzt,500,20528 Montrouge,FR,48.8162,2.31393,u09tswsyyrpr,500,38708 Montreuil,FR,48.86415,2.44322,u09tzx7n71ub,2000,111240 Montgeron,FR,48.70543,2.45039,u09tpf83dnpn,500,22843 Meudon,FR,48.81381,2.235,u09tdy73p38y,500,44652 Massy,FR,48.72692,2.28301,u09t5yqqvupx,500,38768 Malakoff,FR,48.81999,2.29998,u09tsr6v13tr,500,29420 Maisons-Alfort,FR,48.81171,2.43945,u09txtbkg61z,1000,53964 Longjumeau,FR,48.69307,2.29431,u09th0q9tq1s,500,20771 Le Plessis-Robinson,FR,48.78889,2.27078,u09te9txch23,500,22510 Le Kremlin-Bicêtre,FR,48.81471,2.36073,u09ttwrn2crz,500,27867 Le Chesnay,FR,48.8222,2.12213,u09t8rc3cjwz,500,29154 La Celle-Saint-Cloud,FR,48.85029,2.14523,u09tbufje6p6,500,21539 Ivry-sur-Seine,FR,48.81568,2.38487,u09twq8egqrc,1000,57897 Issy-les-Moulineaux,FR,48.82104,2.27718,u09tezd5njkr,1000,61447 Fresnes,FR,48.75568,2.32241,u09tkgenkj6r,500,24803 Fontenay-aux-Roses,FR,48.79325,2.29275,u09ts4t92cn3,500,24680 Clamart,FR,48.80299,2.26692,u09tes6dp0dn,1000,51400 Choisy-le-Roi,FR,48.76846,2.41874,u09trn12bez7,500,35590 Chevilly-Larue,FR,48.76476,2.3503,u09tmmr7mfns,500,20125 Châtillon,FR,48.8024,2.29346,u09tshnn96xx,500,32383 Châtenay-Malabry,FR,48.76507,2.26655,u09t7t6mn7yj,500,32715 Charenton-le-Pont,FR,48.82209,2.41217,u09twzu3r9hq,500,30910 Cachan,FR,48.79632,2.33661,u09tt5j7nvqd,500,26540 Bagnolet,FR,48.86667,2.41667,u09tyzzubrxb,500,33504 Bagneux,FR,48.79565,2.30796,u09tsdbx727w,500,38900 Athis-Mons,FR,48.70522,2.39147,u09tn6t2mr16,500,31225 Alfortville,FR,48.80575,2.4204,u09txhf6p7jp,500,37290 Quinze-Vingts,FR,48.84656,2.37439,u09tyh0zz6c8,500,26265 Croulebarbe,FR,48.81003,2.35403,u09tttd5hc5f,500,20062 Gare,FR,48.83337,2.37513,u09ty1cdbxcq,1000,75580 Maison Blanche,FR,48.82586,2.3508,u09tv2rz1xgx,1000,64302  Below is the visual representation of the data sample:  My goal is to find an efficient algorithm that would remove the intersecting circles and only keep the one with largest population. My initial approach was to determine which circles are intersecting using haversine formula. Problem was that to check every record for intersections with others, it needs to traverse entire dataset. The time complexity of this was too high. My second approach was to segregate the dataset by country code, and run the comparisons by chunks:  def _remove_intersecting_circles_for_country(df_country): """Helper function to remove intersections within a single country.""" indices_to_remove = set() for i in range(len(df_country)): for j in range(i + 1, len(df_country)): distance = haversine(df_country['latitude'].iloc[i], df_country['longitude'].iloc[i], df_country['latitude'].iloc[j], df_country['longitude'].iloc[j]) if distance < df_country['estimated_radius'].iloc[i] + df_country['estimated_radius'].iloc[j]: if df_country['population'].iloc[i] < df_country['population'].iloc[j]: indices_to_remove.add(df_country.index[i]) else: indices_to_remove.add(df_country.index[j]) return indices_to_remove all_indices_to_remove = set() for country_code in df['country_code'].unique(): df_country = df[df['country_code'] == country_code] indices_to_remove = _remove_intersecting_circles_for_country(df_country) all_indices_to_remove.update(indices_to_remove) new_df = df.drop(index=all_indices_to_remove) return new_df  This has significantly improved the performance because to check every record we only need to check against all the records with the same country_code. But that still makes a lot of unnecessary comparisons

Answer: Once you have the circles as polygons, determining intersections between polygons is very fast if you use a spatial index to do so. So, you can:  buffer the points to circles. In WGS84 this would be imprecise, so the buffering needs to be done via a sidestep to an equidistant projection. calculating which circles intersect can be done very fast using a spatial index. E.g. geopandas.sjoin uses a spatial index under the hood. now you can determine which circles intersect a circle representing a city with a larger population and remove those.  Result (red: original cities, blue: retained cities):  Code: from pathlib import Path import geopandas as gpd from matplotlib import pyplot as plt from pyproj import CRS, Transformer from shapely.geometry import Point from shapely.ops import transform def geodesic_point_buffer(lat, lon, distance): # Azimuthal equidistant projection aeqd_proj = CRS.from_proj4( f"+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0") tfmr = Transformer.from_proj(aeqd_proj, aeqd_proj.geodetic_crs) buf = Point(0, 0).buffer(distance) # distance in metres return transform(tfmr.transform, buf) # Read the csv file with data #csv_path = Path(__file__).resolve().with_suffix(".csv") csv_path = "https://raw.githubusercontent.com/theroggy/pysnippets/refs/heads/main/pysnippets/stackoverflow_questions/2024/Q4/2024-10-25_circles.csv" points_df = gpd.read_file(csv_path, autodetect_type=True) # Convert the points to circles by buffering them points_buffer_gdf = gpd.GeoDataFrame( points_df, geometry=points_df.apply( lambda row : geodesic_point_buffer(row.latitude, row.longitude, row.estimated_radius), axis=1 ), crs=4326, ) # Determine the intersecting city buffers (result includes self-intersections) intersecting_gdf = points_buffer_gdf.sjoin(points_buffer_gdf) # Get all city buffers that intersect a city with a larger population intersecting_larger_population_df = intersecting_gdf.loc[ intersecting_gdf.population_left < intersecting_gdf.population_right ] # Remove the city buffers that intersect with a larger population city buffer result_gdf = points_buffer_gdf[ ~points_buffer_gdf.index.isin(intersecting_larger_population_df.index) ] # Plot result ax = points_buffer_gdf.boundary.plot(color="red") result_gdf.boundary.plot(color="blue", ax=ax) plt.show()