pypro.mx

Optimización de consultas espaciales

Con unos cientos de polígonos todo es rápido. Con cientos de miles de puntos, un código que funcionaba en segundos puede tardar horas. La buena noticia es que casi toda la lentitud viene de cuatro hábitos: leer más de lo necesario, recorrer filas con bucles, ignorar el índice espacial y trabajar con geometrías más detalladas de lo que hace falta. Esta lección corrige cada uno y mide la diferencia.

Medir antes de optimizar

Optimiza lo que es lento de verdad, no lo que parece lento. time.perf_counter() basta para comparar dos versiones del mismo paso. Para los ejemplos generamos 20,000 puntos aleatorios sobre México con una semilla fija:

import time
import numpy as np
import geopandas as gpd

NE = "https://naciscdn.org/naturalearth/10m/cultural/"
estados = gpd.read_file(NE + "ne_10m_admin_1_states_provinces.zip",
                        where="iso_a2 = 'MX'", columns=["name", "iso_a2"])
estados = estados.dropna(subset=["name"])

rng = np.random.default_rng(42)
xmin, ymin, xmax, ymax = estados.total_bounds
puntos = gpd.GeoDataFrame(
    geometry=gpd.points_from_xy(rng.uniform(xmin, xmax, 20_000),
                                rng.uniform(ymin, ymax, 20_000)), crs=4326)

1. Leer solo lo que necesitas

Con el motor pyogrio (el predeterminado en GeoPandas 1.x), read_file filtra durante la lectura: where descarta filas con SQL, columns descarta columnas y bbox descarta todo lo que cae fuera de un rectángulo. Con use_arrow=True la lectura usa Apache Arrow (requiere pyarrow) y suele ser más rápida en archivos grandes. Leer un archivo completo y filtrar después en pandas mueve a memoria datos que vas a tirar.

2. Vectorizar en lugar de recorrer filas

Desde Shapely 2, las operaciones de GeoPandas se ejecutan en C sobre todo el arreglo de geometrías. Un bucle con iterrows o apply vuelve a Python en cada fila y pierde esa ventaja:

puntos_m = puntos.to_crs(6372)
cdmx = gpd.GeoSeries(gpd.points_from_xy([-99.13], [19.43]), crs=4326).to_crs(6372).iloc[0]

t = time.perf_counter()
lento = puntos_m.geometry.apply(lambda g: g.distance(cdmx))
t_bucle = time.perf_counter() - t

t = time.perf_counter()
rapido = puntos_m.distance(cdmx)
t_vector = time.perf_counter() - t
print(f"apply: {t_bucle:.3f} s   vectorizado: {t_vector:.4f} s")
print(f"{t_bucle / t_vector:.0f} veces más rápido;", np.allclose(lento, rapido))

Las dos versiones calculan la distancia de cada punto a la Ciudad de México y dan el mismo resultado, pero la vectorizada suele ser decenas de veces más rápida. En operaciones costosas por sí mismas, como un buffer, la diferencia es menor porque domina el cálculo geométrico.

3. Usar el índice espacial

Para saber en qué estado cae cada punto, la forma ingenua compara cada punto con cada polígono. El índice espacial (un árbol STR que guarda el rectángulo envolvente de cada geometría) descarta casi todas las comparaciones antes de hacer la prueba exacta. gpd.sjoin lo usa automáticamente:

t = time.perf_counter()
asignados = []
for p in puntos.geometry[:2000]:
    dentro = estados[estados.contains(p)]
    asignados.append(dentro["name"].iloc[0] if len(dentro) else None)
t_ingenuo = (time.perf_counter() - t) * 10   # estimado para 20,000

t = time.perf_counter()
union = gpd.sjoin(puntos, estados, predicate="within")
t_sjoin = time.perf_counter() - t
print(f"bucle (estimado): {t_ingenuo:.1f} s   sjoin: {t_sjoin:.3f} s")
print(len(union), "puntos dentro de algún estado")

El índice también se consulta directamente: estados.sindex.query(geom, predicate="intersects") devuelve las posiciones de las geometrías que cumplen el predicado. Es útil cuando necesitas candidatos rápidos dentro de tu propio algoritmo.

4. Simplificar geometrías cuando el detalle no importa

Una costa a escala 1:10 millones tiene miles de vértices. Para un mapa del país o para una prueba de pertenencia aproximada, simplify reduce los vértices manteniendo la forma general. La tolerancia está en las unidades del CRS, por eso conviene proyectar primero:

import shapely

estados_m = estados.to_crs(6372)
simple = estados_m.simplify(2000)        # tolerancia de 2 km
print("vértices originales:", shapely.get_num_coordinates(estados_m.geometry).sum())
print("vértices simplificados:", shapely.get_num_coordinates(simple.geometry).sum())
cambio = abs(simple.area.sum() / estados_m.area.sum() - 1)
print(f"cambio en el área total: {cambio:.2%}")

Usa simplify(..., preserve_topology=True) (el valor por omisión) para no crear polígonos inválidos. No simplifiques antes de medir distancias o áreas que deban ser exactas.

5. Elegir un formato rápido

GeoParquet (to_parquet y read_parquet) es columnar y comprimido: se lee mucho más rápido que un Shapefile o un GeoJSON y permite leer solo algunas columnas. GeoJSON es cómodo para compartir datos pequeños en la web, pero es texto y crece mucho.

Trampas comunes

iterrows con geometrías. Casi siempre existe una versión vectorizada o un sjoin.
Uniones con CRS distintos. sjoin avisa, pero el resultado será incorrecto si sigues adelante.
Buffers en grados. Además de incorrectos, generan polígonos deformes que hacen más lentas las operaciones siguientes.
Simplificar sin pensar en la escala. Una tolerancia de 2 km es razonable para un mapa nacional y desastrosa para un mapa de colonias.
Optimizar a ciegas. Mide cada paso; a veces el cuello de botella es la lectura, no el análisis.

Cierre

El orden de ataque es: leer menos, vectorizar, dejar que sjoin use el índice y simplificar solo para dibujar. Con esos cuatro cambios, la mayoría de los análisis pasan de minutos a fracciones de segundo.

Visualización Geospacial Avanzada con GeoPandas

Optimización de Consultas Espaciales