Análisis Geoespacial con GeoPandas: Estudio de Radios Censales en CABA

Objetivos de Aprendizaje

• Trabajar con datos geoespaciales utilizando GeoPandas y Shapely
• Comprender y gestionar CRS (Coordinate Reference Systems) para cálculos correctos
• Implementar joins espaciales para combinar datos por ubicación
• Realizar agregaciones zonales para resumir información por áreas
• Crear visualizaciones estáticas e interactivas con mapas
• Aplicar análisis de proximidad y operaciones geométricas
• Discretizar el espacio usando grillas hexagonales (H3)

Contexto de Negocio

El dataset contiene información de radios censales de la Ciudad Autónoma de Buenos Aires (CABA). Los radios censales son unidades geográficas utilizadas por los censos para dividir el territorio y recopilar información demográfica.

Como analista de datos urbanos, necesitamos:

• Analizar la distribución poblacional por áreas de la ciudad
• Identificar zonas con alta densidad habitacional
• Estudiar relaciones espaciales entre diferentes variables
• Crear visualizaciones informativas para toma de decisiones

Descripción del Dataset

• Fuente: GeoJSON de radios censales de CABA
• Registros: ~3,554 radios censales
• CRS Original: EPSG:4326 (WGS 84 - coordenadas geográficas)
• Geometrías: MultiPolygon

Estructura:

• RADIO_ID: Identificador único del radio censal
• BARRIO: Nombre del barrio
• COMUNA: Comuna a la que pertenece
• POBLACION: Cantidad de habitantes
• VIVIENDAS: Número de viviendas
• HOGARES: Cantidad de hogares
• HOGARES_NBI: Hogares con Necesidades Básicas Insatisfechas
• AREA_KM2: Área en kilómetros cuadrados
• geometry: Geometría del polígono

Proceso de Análisis

1. Preparación del Entorno

import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import contextily as cx
from shapely.ops import unary_union
import warnings
warnings.filterwarnings('ignore')

# Configurar visualización
plt.style.use('seaborn-v0_8-darkgrid')
pd.set_option('display.max_columns', None)

print("✅ Librerías listas")

Librerías clave:

• geopandas: Extensión de pandas para datos geoespaciales
• shapely: Manipulación y análisis de geometrías
• contextily: Agregar mapas base (tiles) a visualizaciones
• pyproj: Conversión entre sistemas de coordenadas

2. Carga y Exploración de Datos

# Cargar datos geoespaciales desde URL
url = "https://bitsandbricks.github.io/data/CABA_rc.geojson"
radios = gpd.read_file(url)

# Inspeccionar datos
print(f"Shape: {radios.shape}")
print(f"CRS: {radios.crs}")
print(f"Geometry type: {radios.geometry.type.unique()}")

# Verificar integridad
print(f"Geometrías vacías: {radios.geometry.is_empty.sum()}")
print(f"Geometrías nulas: {radios.geometry.isna().sum()}")

Resultados:

• 3,554 radios censales en total
• CRS: EPSG:4326 (coordenadas latitud/longitud)
• Sin geometrías vacías o nulas (dataset limpio)

2.1 Análisis Exploratorio

# Estadísticas básicas
print(radios[['POBLACION', 'AREA_KM2', 'HOGARES']].describe())

# Ver primeras filas
print(radios.head())

Decisión tomada: El CRS original (EPSG:4326) usa grados, lo cual no es adecuado para cálculos de área y distancia. Necesitamos proyectar a un CRS métrico.

3. Gestión de CRS (Coordinate Reference Systems)

¿Por qué es importante el CRS?

El Coordinate Reference System define cómo las coordenadas se mapean a ubicaciones en la Tierra:

• Geográficas (EPSG:4326): Latitud/Longitud en grados
• Proyectadas (EPSG:3857, UTM): Coordenadas X/Y en metros

Para cálculos de área, distancia y buffer, necesitamos un CRS proyectado.

# Proyectar a Web Mercator (EPSG:3857) - coordenadas en metros
radios_m = radios.to_crs(epsg=3857)

# Calcular área real en metros cuadrados
radios_m['area_m2'] = radios_m.geometry.area

# Calcular densidad poblacional (habitantes por km²)
radios_m['densidad_hab_km2'] = radios_m['POBLACION'] / (radios_m['area_m2'] / 1e6)

print(radios_m[['BARRIO', 'area_m2', 'densidad_hab_km2']].head())

Resultados:

• Áreas correctamente calculadas en m²
• Densidad poblacional muestra alta variabilidad entre radios
• Algunos radios tienen densidad > 30,000 hab/km²

4. Visualización Básica

4.1 Silueta de CABA

fig, ax = plt.subplots(figsize=(10, 10))
radios.plot(ax=ax, color='lightblue', edgecolor='black', linewidth=0.3)
ax.set_title('Radios Censales de CABA', fontsize=14)
ax.axis('off')
plt.tight_layout()
plt.show()

4.2 Mapa Coroplético (Densidad Poblacional)

fig, ax = plt.subplots(figsize=(12, 10))

# Mapa coroplético
radios_m.plot(
    column='densidad_hab_km2',
    cmap='YlOrRd',
    legend=True,
    ax=ax,
    edgecolor='black',
    linewidth=0.1,
    legend_kwds={'label': 'Densidad (hab/km²)', 'shrink': 0.8}
)

ax.set_title('Densidad Poblacional por Radio Censal - CABA', fontsize=14)
ax.axis('off')
plt.tight_layout()
plt.show()

Insight: Las zonas céntricas muestran mayor densidad poblacional, mientras que áreas periféricas tienen densidades más bajas.

5. Agregación por Barrio

# Disolver geometrías por barrio (unificar radios del mismo barrio)
barrios = radios_m.dissolve(by='BARRIO', aggfunc={
    'POBLACION': 'sum',
    'HOGARES': 'sum',
    'VIVIENDAS': 'sum',
    'area_m2': 'sum'
}).reset_index()

# Recalcular densidad a nivel barrio
barrios['densidad_hab_km2'] = barrios['POBLACION'] / (barrios['area_m2'] / 1e6)

# Ordenar por densidad
barrios_top = barrios.nlargest(10, 'densidad_hab_km2')
print(barrios_top[['BARRIO', 'POBLACION', 'densidad_hab_km2']])

Resultados:

• Barrios más densos: San Nicolás, Montserrat, Constitución
• Menos densos: Villa Soldati, Villa Riachuelo, Parque Avellaneda

6. Joins Espaciales

Los spatial joins permiten combinar datos basándose en relaciones geométricas (intersección, contención, proximidad).

6.1 Ejemplo: Puntos de Interés

# Crear puntos de interés ficticios
pois = gpd.GeoDataFrame({
    'nombre': ['Obelisco', 'Casa Rosada', 'Congreso'],
    'tipo': ['Monumento', 'Político', 'Político']
}, geometry=gpd.points_from_xy(
    [-58.3816, -58.3703, -58.3927],  # longitud
    [-34.6037, -34.6082, -34.6098]   # latitud
), crs='EPSG:4326')

# Proyectar a mismo CRS
pois_m = pois.to_crs(epsg=3857)

# Spatial join: asignar cada punto al radio que lo contiene
pois_con_radio = gpd.sjoin(pois_m, radios_m, how='left', predicate='within')
print(pois_con_radio[['nombre', 'BARRIO', 'COMUNA']])

Resultado: Cada punto queda asociado a su radio censal y barrio correspondiente.

7. Buffers y Análisis de Proximidad

# Crear buffer de 500m alrededor del Obelisco
obelisco = pois_m[pois_m['nombre'] == 'Obelisco']
buffer_500m = obelisco.geometry.buffer(500)

# Encontrar radios que intersectan el buffer
radios_cercanos = radios_m[radios_m.geometry.intersects(buffer_500m.iloc[0])]

print(f"Radios dentro de 500m del Obelisco: {len(radios_cercanos)}")

Aplicación: Útil para análisis de cobertura de servicios, zonas de influencia o estudios de accesibilidad.

8. Visualización con Mapa Base

fig, ax = plt.subplots(figsize=(12, 10))

# Plotear choropleth
radios_m.plot(
    column='densidad_hab_km2',
    cmap='YlOrRd',
    alpha=0.7,
    ax=ax,
    edgecolor='black',
    linewidth=0.1
)

# Agregar mapa base (tiles de OpenStreetMap)
cx.add_basemap(ax, crs=radios_m.crs, source=cx.providers.OpenStreetMap.Mapnik, alpha=0.5)

ax.set_title('Densidad Poblacional - CABA (con mapa base)', fontsize=14)
ax.axis('off')
plt.tight_layout()
plt.show()

9. Grillas Hexagonales (H3)

La discretización con grillas hexagonales permite agregar datos en celdas uniformes, útil para análisis espaciales y visualizaciones.

import h3

# Convertir geometrías a hexágonos H3 (resolución 9)
def geometria_a_h3(geometry, resolucion=9):
    """Convierte un polígono a conjunto de hexágonos H3"""
    coords = list(geometry.exterior.coords)
    # Invertir lat/lon para H3
    geojson_coords = [[coord[1], coord[0]] for coord in coords]
    hexagons = h3.polyfill({'type': 'Polygon', 'coordinates': [geojson_coords]}, resolucion)
    return list(hexagons)

# Aplicar a algunos radios
sample_radios = radios.head(10)
sample_radios['h3_hexagons'] = sample_radios.geometry.apply(geometria_a_h3)

print("Ejemplo de hexágonos H3:")
print(sample_radios[['RADIO_ID', 'h3_hexagons']].head())

Ventajas de H3:

• Tamaño uniforme de celdas
• Vecindad bien definida
• Eficiente para agregaciones espaciales
• Compatible con análisis a múltiples escalas (resoluciones)

Análisis Crítico y Decisiones

Decisiones de CRS

¿Por qué EPSG:3857 (Web Mercator)?

• Pro: Compatible con la mayoría de servicios de tiles (OpenStreetMap, Google Maps)
• Pro: Coordenadas en metros facilitan cálculos de área/distancia
• Contra: Distorsión en latitudes altas (no es problema para CABA)

Alternativa: EPSG:22185 (POSGAR 2007 / Argentina zona 5) sería más preciso para Argentina, pero menos compatible con mapas base.

Decisión: Usar EPSG:3857 para balance entre precisión y compatibilidad.

Elección de Visualizaciones

1. Mapas coropléticos: Efectivos para mostrar distribución de variables continuas
2. Clasificación de colores: Esquema YlOrRd funciona bien para densidad (amarillo = bajo, rojo = alto)
3. Transparencia al agregar tiles: Permite ver tanto datos como contexto geográfico

Agregación por Barrio vs Radio

• Radio censal: Máximo detalle, pero 3,554 unidades pueden ser abrumadoras
• Barrio: Balance entre detalle y legibilidad (48 barrios)
• Comuna: Menor granularidad pero útil para decisiones administrativas (15 comunas)

Decisión: Mostrar múltiples niveles de agregación según el análisis específico.

Extra: Mapas Interactivos con Folium

Aunque GeoPandas es excelente para visualizaciones estáticas, a veces necesitamos mapas interactivos para explorar datos con mayor detalle. Folium es una librería de Python que facilita la creación de mapas interactivos basados en Leaflet.js.

¿Por qué usar Folium?

• Interactividad: Zoom, paneo y popups con información detallada.
• Mapas base: Fácil acceso a OpenStreetMap, Stamen, CartoDB, etc.
• Capas: Posibilidad de superponer múltiples capas de datos (Heatmaps, Choropleths, Marcadores).

Ejemplo: Mapa de Calor y Marcadores

A continuación, mostramos cómo crear un mapa interactivo que visualice la densidad de puntos (simulados) y permita inspeccionar detalles individuales.

import folium
from folium import plugins

# 1. Crear el mapa base centrado en CABA
m = folium.Map(location=[-34.6037, -58.3816], zoom_start=12, tiles='CartoDB positron')

# 2. Convertir GeoDataFrame a formato compatible (lat, lon)
# Asumiendo que tenemos un GeoDataFrame 'gdf_puntos' con geometría de puntos
heat_data = [[point.xy[1][0], point.xy[0][0]] for point in gdf_puntos.geometry]

# 3. Añadir Mapa de Calor
plugins.HeatMap(heat_data, radius=15).add_to(m)

# 4. Añadir Marcadores con Popups
for idx, row in gdf_puntos.iterrows():
    folium.Marker(
        location=[row.geometry.y, row.geometry.x],
        popup=f"Categoría: {row['category']}",
        icon=folium.Icon(color='blue', icon='info-sign')
    ).add_to(m)

# 5. Guardar o mostrar
m.save('mapa_interactivo.html')
m

Pie de figura: Esta imagen es una representación estática de lo que sería el mapa interactivo. En la versión HTML real, podrías hacer zoom para ver calles específicas y hacer clic en cada marcador para ver la categoría y otros datos del punto. Los mapas de calor (heatmap) son especialmente útiles para identificar "zonas calientes" de actividad o densidad poblacional sin saturar el mapa con miles de puntos individuales.

Conclusión

El análisis geoespacial con GeoPandas permite trabajar con datos geográficos con la misma facilidad que pandas, realizando cálculos espaciales precisos mediante el manejo correcto de CRS. Los joins espaciales facilitan la combinación de datos de diferentes fuentes basándose en relaciones geométricas, mientras que las visualizaciones coropléticas y mapas interactivos comunican efectivamente patrones espaciales complejos. Las grillas hexagonales (H3) ofrecen una solución escalable para análisis de grandes volúmenes de datos geoespaciales, manteniendo uniformidad y eficiencia computacional.

Próximos pasos: Implementar análisis de autocorrelación espacial (Moran's I) para detectar clustering de variables demográficas. Explorar técnicas de interpolación espacial (Kriging, IDW) para estimar valores en áreas sin datos. Integrar datos temporales para análisis espacio-temporales de evolución urbana. Aplicar machine learning espacial para predicción de variables demográficas en áreas no censadas. Desarrollar dashboards interactivos con Plotly Dash o Streamlit que permitan exploración dinámica de múltiples variables geoespaciales simultáneamente.