如何有效地从数据集中删除重叠的圆

一旦将圆变成多边形，如果使用空间索引来执行此操作，则确定多边形之间的交点会非常快。

因此，您可以：

• 将点缓冲到圆。在 WGS84 中，这会不精确，因此需要通过回避等距投影来完成缓冲。
• 使用空间索引可以非常快速地计算哪些圆相交。例如geopandas.sjoin在后台使用空间索引。
• 现在您可以确定哪些圆圈与代表人口较多城市的圆圈相交，然后将其删除。

结果（红色：原有城市，蓝色：保留城市）：

代码：

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()

从视觉表现来看，首先想到的是使用具有最外轮廓的图像处理方法（并将地理坐标映射到像素坐标）。不确定是否适合城市之间间距很大的非常大的数据集（但为什么要查看里昂是否与巴黎重叠？），但在提供的近巴黎样本上，它的时间比当前解决方案更快。

与当前的 O(n^2) 相比，我们的时间复杂度为 O(n log(n)) + O(n)。然而，主要的权衡可能是创建具有足够分辨率的图像。考虑到底层优化的 C++ OpenCV，它的轮廓处理应该还不错

import numpy as np
import cv2
import pandas as pd
from typing import Tuple
from io import StringIO

def detect_non_overlapping_circles(df: pd.DataFrame, image_size:Tuple[int, int] = (128, 128)):
    img = np.zeros(image_size, dtype=np.uint8) # Black image

    lat_min, lat_max = df['latitude'].min(), df['latitude'].max()
    lon_min, lon_max = df['longitude'].min(), df['longitude'].max()

    def geo_to_pixel(lat: float, lon: float):
        x = int((lon - lon_min) / (lon_max - lon_min) * (image_size[0] - 100) + 50)
        y = int((lat - lat_min) / (lat_max - lat_min) * (image_size[1] - 100) + 50)
        return (x, y)

    # Convert radius to pixels based on the image size and coordinate range
    def radius_to_pixel(radius_meters: float):
        # Using an approximate conversion factor
        earth_circumference = 40075000  # meters
        meters_per_degree = earth_circumference / 360
        lat_range = lat_max - lat_min
        pixels_per_meter = (image_size[1] - 100) / (lat_range * meters_per_degree)
        return int(radius_meters * pixels_per_meter)

    df_sorted = df.sort_values('population', ascending=False)

    kept_indices = []

    for idx, row in df_sorted.iterrows():
        test_img = img.copy()
    
        # Current circle
        center = geo_to_pixel(row['latitude'], row['longitude'])
        radius = radius_to_pixel(row['estimated_radius'])
    
        # Find outermost contour at each time
        cv2.circle(test_img, center, radius, 255, -1)    
        contours, _ = cv2.findContours(test_img, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    
        # If exactly one contour, this circle dont overlap with any one
        if len(contours) == len(kept_indices) + 1:
            kept_indices.append(idx)
            img = test_img  # Update the main image

    return df.loc[kept_indices]

def main(df):
    result_df = detect_non_overlapping_circles(df)
    
    print(f"Original number of cities: {len(df)}")
    print(f"Cities after removing overlaps: {len(result_df)}")
    print(result_df[['name_city', 'population', 'estimated_radius']].to_string())

    return result_df

# Not sure if you're working with csv but from what provided
data = 
"""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
""" # ...

df = pd.read_csv(StringIO(data))
df.dropna()

proc_df = main(df) 

使用提供的 csv 和城市间距，它可以很好地处理 (128, 128) 图像分辨率。在我的计算机上，处理 (1024, 1024) 图像的执行时间仍然比_remove_intersecting_circles_for_country

降低复杂性的常用方法是在地图上叠加网格。例如，叠加一个边长为 100 英里的方格网格。创建一个数据结构来告诉您：

1. 城市位于哪个网格区域（城市中心）
2. 哪些城市（同样是城市圆的中心）位于该网格范围内

假设没有一个城市的半径大于 100 英里，那么一个城市只能与以下中心的城市相交：

• 在自己的网格方格内
• 在相邻的网格方块中
• 对于大城市来说，相邻的网格方块相邻。

在下图中，标记为“0”的网格是包含您当前正在检查的城市的网格。它唯一可能与之相交的城市是网格位置标记为“1”或“2”的城市。同样，假设您的网格大小大于任何单个城市的规模。

    X X X X X X X
    X 2 2 2 2 2 X
    X 2 1 1 1 2 X
    X 2 1 O 1 2 X
    X 2 1 1 1 2 X
    X 2 2 2 2 2 X
    X X X X X X X        

这大大减少了您在尝试确定重叠时需要查看的城市数量。您将有...十几个？一百个？而不是 20,000 个城市？

构建该网格以便您可以分类存放物品的开销很小，但性能优势却是巨大的。

在构建该数据结构时，为每个网格方格不仅包含哪些城市位于其中心，还包含哪些城市侵占了它，这可能是有意义的。因此，如果“超人前传”位于您标记为 C2 的网格方格中，并且其周长侵占了 C3，则 C2 应该包含一个条目，实际上显示“超人前传位于 C2”。C3 有一个条目显示“超人前传侵占了网格方格 C3。”

不过，在需要时，这些信息在运行时很容易计算。也许保留它是不必要的。

更新

想想看，因为您有每个城市的纬度/经度，所以您已经拥有了构建网格所需的大部分信息。让您的网格纬度为一度，经度为一度。在赤道处，您的方格面积约为 69 英里（111 公里）见方。当您靠近两极时，面积会压缩，但这并不重要。要确定城市属于哪个网格，您只需将纬度/经度和索引的小数部分切掉即可。您可能需要对数字进行一些标准化（即从 -180...180 转换为 0...360）。

无论如何，这可以让你将问题缩小到更易于管理的规模。