Como remover eficientemente círculos sobrepostos do conjunto de dados

Depois de ter os círculos como polígonos, determinar as interseções entre os polígonos será muito rápido se você usar um índice espacial para fazer isso.

Então, você pode:

• bufferize os pontos para círculos. No WGS84 isso seria impreciso, então o buffering precisa ser feito por meio de um sidestep para uma projeção equidistante.
• calcular quais círculos se cruzam pode ser feito muito rápido usando um índice espacial. Por exemplo, geopandas.sjoin usa um índice espacial por baixo dos panos.
• agora você pode determinar quais círculos cruzam um círculo que representa uma cidade com uma população maior e removê-los.

Resultado (vermelho: cidades originais, azul: cidades mantidas):

Código:

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

Da representação visual, a primeira coisa que me veio à mente foi usar uma abordagem de processamento de imagem com contornos mais externos (e mapear coordenadas geográficas para coordenadas de pixel). Não tenho certeza se é adequado para conjuntos de dados muito grandes com espaçamentos muito grandes entre cidades (mas por que alguém deveria ver se Lyon se sobrepõe a Paris?), mas nas amostras fornecidas perto de Paris, é vezes mais rápido do que a solução atual.

Estamos em O(n log(n)) + O(n) para complexidade de tempo comparado com o O(n^2 atual). No entanto, a principal compensação pode ser a criação de imagem com resolução suficiente. Seu processamento de contorno deve ser razoável considerando o OpenCV C++ otimizado sob o capô

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) 

Com o csv fornecido e o espaçamento da cidade, funcionou bem até (128, 128) na resolução da imagem. Na minha máquina, trabalhar em uma imagem (1024, 1024) ainda era mais rápido no tempo de execução do que_remove_intersecting_circles_for_country

Uma maneira comum de reduzir a complexidade é sobrepor uma grade em cima do seu mapa. Por exemplo, sobreponha uma grade de quadrados que tenham 100 milhas de lado. Crie uma estrutura de dados que lhe dirá:

1. Em que quadrado da grade uma cidade está (no centro da cidade)
2. Quais cidades (novamente, o centro do círculo da cidade) estão localizadas dentro desse quadrado da grade

Supondo que nenhuma cidade tenha mais de 100 milhas de raio, uma cidade só pode cruzar cidades cujos centros são:

• Dentro de sua própria grade quadrada
• Em um quadrado de grade adjacente
• No caso de grandes cidades, adjacente a um quadrado de grade adjacente.

Na ilustração abaixo, a grade marcada com '0' é a grade que contém a cidade que você está examinando no momento. As únicas cidades com as quais ela pode possivelmente cruzar são aquelas em posições de grade marcadas com '1' ou '2'. Novamente, assumindo que o tamanho da sua grade é maior do que o tamanho de qualquer cidade individual.

    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        

Isso reduz muito o número de cidades que você precisa olhar quando estiver tentando determinar sobreposições. Em vez de 20.000 cidades, você terá... uma dúzia? Cem?

É um pequeno gasto adicional construir essa grade para que você possa agrupar seus itens, mas o benefício no desempenho é enorme.

Pode fazer sentido, quando você estiver construindo essa estrutura de dados, incluir para cada quadrado da grade não apenas quais cidades estão centralizadas nele, mas também quais cidades o invadem. Então , se "Smallville" estiver no quadrado da grade que você rotular, por exemplo, C2, e sua circunferência invadir C3, então C2 deve conter uma entrada que diga, na verdade, "Smallville está localizada em C2." C3 tem uma entrada que diz, "Smallville invade o quadrado da grade C3."

Essa informação é bem fácil de computar em tempo de execução quando você precisa dela, no entanto. Talvez mantê-la por perto seja desnecessário.

Atualizar

Pensando bem, como você tem a lat/lon para cada cidade, você já tem a maior parte do que precisa para construir a grade. Faça sua grade com um grau de latitude por um grau de longitude. No equador, isso lhe dá quadrados de aproximadamente 69 milhas (111 km) quadrados. Isso se comprime conforme você se aproxima dos polos, mas não deve importar. Para determinar a qual grade uma cidade pertence, você apenas corta as partes decimais da sua lat/lon e indexa na grade. Você pode precisar de alguma normalização dos números (ou seja, converter de -180...180 para 0...360).

De qualquer forma, isso permite reduzir o problema a um tamanho muito mais administrável.