Clase 12: Geoespacial con GeoPandas (Assignment)¶
Ingeniería de Datos — Universidad Católica del Uruguay
Objetivo: Implementar un pipeline geoespacial de punta a punta con GeoPandas/Shapely: carga, validación de CRS, joins espaciales, agregaciones zonales y visualización (estática e interactiva) sobre CABA.
Tiempo estimado: 120–150 minutos
📚 Lecturas mínimas (recuerdo)¶
- Brust, Cap. 6: Información geográfica y mapas
- Kaggle Geospatial (Intro, CRS, Interactive Maps, Joins, Proximity)
- GeoPandas User Guide (Introduction, CRS, Plotting)
Setup y Carga de Datos¶
Instalación rápida¶
!pip install -q geopandas shapely pyproj fiona rtree contextily folium mapclassify
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, platform
warnings.filterwarnings('ignore')
print("✅ Librerías listas")
print(f"GeoPandas: {gpd.__version__}")
print(f"Python: {platform.python_version()}")
1) Leer polígonos de radios censales (CABA)¶
Fuente GeoJSON (igual al texto de referencia):
https://bitsandbricks.github.io/data/CABA_rc.geojson
radios = gpd.read_file("https://bitsandbricks.github.io/data/CABA_rc.geojson")
print(radios.crs)
print(radios.head())
print("Filas:", len(radios))
print("Geometrías vacías:", radios.geometry.is_empty.sum())
print("Geometrías nulas:", radios.geometry.isna().sum())
2) Estándar de CRS y áreas¶
Regla práctica: para áreas/distancias, usar CRS proyectado en metros. Para CABA, una opción simple es Web Mercator (EPSG:3857) para el ejercicio. Si conoces un CRS local (p.ej., UTM 21S EPSG:32721), úsalo.
# ⚠️ COMPLETA: proyecta a CRS métrico (3857 recomendado si no sabes uno local)
radios_m = radios.to_crs(epsg=____)
radios_m["area_m2"] = radios_m.geometry.area
radios_m["densidad_hab_km2"] = radios_m["POBLACION"] / (radios_m["area_m2"] / 1e6)
radios_m[["BARRIO", "POBLACION", "AREA_KM2", "area_m2", "densidad_hab_km2"]].head()
Pista: usa epsg=3857 o un CRS local (p.ej., UTM 21S epsg=32721). Refiere: GeoPandas CRS, EPSG explorer.
Parte A — Visualización básica y normalizaciones¶
A.1 Silueta y coropleta rápida (GeoPandas.plot)¶
ax = radios.plot(figsize=(7,7), edgecolor="gray", facecolor="white")
ax.set_title("CABA — Radios Censales (WGS84)")
plt.show()
# Densidad por km² (proyectado)
ax = radios_m.plot(column="densidad_hab_km2", scheme="____", k=____,
legend=True, cmap="viridis", figsize=(8,8), linewidth=0)
ax.set_axis_off()
ax.set_title("Densidad de población (hab/km²)")
plt.tight_layout(); plt.show()
Pistas: scheme='quantiles' y k=5 (o prueba natural_breaks). Refiere: GeoPandas plot schemes, mapclassify docs.
A.2 Añadir tiles de contexto (contextily)¶
fig, ax = plt.subplots(figsize=(8,8))
radios_m.to_crs(epsg=3857).plot(ax=ax, column="densidad_hab_km2", cmap="viridis", legend=True, linewidth=0)
cx.add_basemap(ax, source=cx.providers.CartoDB.____)
ax.set_axis_off(); plt.tight_layout(); plt.show()
Pistas: Positron o PositronNoLabels (alternativa: Esri.WorldGrayCanvas). Refiere: xyzservices providers, contextily.
📝 Preguntas de reflexión — Parte A (completa los espacios)¶
1) ¿Qué esquema de clasificación fue más informativo para la coropleta (quantiles, k=5)? Por qué: _____. - Qué buscar: homogeneidad de rangos vs cortes que separen outliers; mapas legibles sin bandas vacías.
2) En la proyección elegida, la variable que normalizaste para densidad fue: POBLACION / (_____/1e6) → unidades: _____.
- Qué buscar: confirmar que el divisor está en m² y transformar a km².
3) Identifica 2 zonas con mayor densidad: Zona 1 = _; Zona 2 = . Hipótesis: __. - Qué buscar: barrios identificables por la coropleta; contrastar con conocimiento urbano.
Parte B — Attribute Join y métricas per cápita¶
Dataset de atención al ciudadano (igual a referencia):
http://bitsandbricks.github.io/data/gcba_suaci_comunas.csv
Objetivo: agregar por BARRIO, unir a geometrías y calcular per cápita.
suaci = pd.read_csv("http://bitsandbricks.github.io/data/gcba_suaci_comunas.csv", encoding="____")
suaci_barrio = suaci.groupby("BARRIO", as_index=False)["total"].sum()
# Agregado geográfico por barrio (disolver)
barrios_m = (radios_m
.dissolve(by="BARRIO", aggfunc={
"POBLACION": "____",
"VIVIENDAS": "____",
"HOGARES": "____",
"HOGARES_NBI": "____",
"area_m2": "____"
})
.reset_index())
barrios_m = barrios_m.merge(suaci_barrio, on="BARRIO", how="left")
barrios_m["contactos_pc"] = barrios_m["total"] / barrios_m["POBLACION"]
barrios_m[["BARRIO", "POBLACION", "total", "contactos_pc" ]].sort_values("contactos_pc", ascending=False).head(10)
Pistas: encoding='ISO-8859-1' (o 'latin1'). Para aggfunc usa "sum" en todas las columnas agregadas. Refiere: pandas.read_csv, GeoPandas dissolve.
# Coropleta per cápita por barrio
ax = barrios_m.plot(column="contactos_pc", cmap="____", legend=True,
figsize=(8,8), linewidth=0)
ax.set_axis_off(); ax.set_title("Contactos SUACI per cápita por barrio")
plt.tight_layout(); plt.show()
Pista: cmap='YlOrRd' o OrRd para gradientes perceptuales. Refiere: Matplotlib colormaps.
📝 Preguntas de reflexión — Parte B (completa los espacios)¶
1) La métrica per cápita calculada fue contactos_pc = total / _____. ¿Por qué es preferible a valores absolutos? _____.
- Qué buscar: sesgo por tamaño poblacional.
2) Top‑3 barrios por contactos_pc: 1) _ 2) 3) __. ¿Qué patrón común observas? _____.
- Qué buscar: centralidad, usos mixtos, actividad administrativa/comercial.
3) ¿Qué columna usarías para segmentar por tipo de solicitud si existiera? Nombre de columna: _. Métrica: ___.
- Qué buscar: equivalente a RUBRO en datasets similares; tasas por 1k habitantes.
📝 Preguntas de reflexión — Parte C (completa los espacios)¶
1) ¿Qué relación esperas entre estaciones_por_km2 y contactos_pc? Correlación esperada: _ porque ___.
- Qué buscar: relación directa/indirecta; movilidad/afluencia.
2) Barrio con menor cobertura (mayor dist_min_m): _. Posible causa: ___.
- Qué buscar: periferia, áreas portuarias/industriales, espacios verdes extensos.
3) Si tuvieras que priorizar nuevas estaciones, define un criterio: score = α·(_____) + β·(_____). Valores α, β: _____.
- Qué buscar: combinar distancia, densidad y demanda (contactos).
Parte C — Joins espaciales: SUBTE y cobertura¶
Capas:
- Líneas SUBTE: http://bitsandbricks.github.io/data/subte_lineas.geojson
- Estaciones SUBTE: http://bitsandbricks.github.io/data/subte_estaciones.geojson
lineas = gpd.read_file("http://bitsandbricks.github.io/data/subte_lineas.geojson").to_crs(barrios_m.crs)
estaciones = gpd.read_file("http://bitsandbricks.github.io/data/subte_estaciones.geojson").to_crs(barrios_m.crs)
# Conteo de estaciones por barrio (puntos dentro de polígonos)
est_x_barrio = gpd.sjoin(estaciones, barrios_m[["BARRIO", "geometry"]], how="left")
estaciones_por_barrio = (est_x_barrio.groupby("BARRIO").size()
.rename("n_estaciones").reset_index())
barrios_m = barrios_m.merge(estaciones_por_barrio, on="BARRIO", how="left")
barrios_m["n_estaciones"] = barrios_m["n_estaciones"].fillna(0).astype(int)
# Densidad estaciones por km²
barrios_m["estaciones_por_km2"] = barrios_m["n_estaciones"] / (barrios_m["area_m2"] / 1e6)
barrios_m[["BARRIO", "n_estaciones", "estaciones_por_km2"]].sort_values(["estaciones_por_km2", "n_estaciones"], ascending=False).head(12)
# Mapa: estaciones por km² + líneas del SUBTE
fig, ax = plt.subplots(figsize=(9,9))
barrios_m.plot(ax=ax, column="estaciones_por_km2", cmap="____", legend=True, linewidth=0.4, edgecolor="#999")
lineas.plot(ax=ax, color="gold", linewidth=1.2)
estaciones.plot(ax=ax, color="tomato", markersize=8)
ax.set_axis_off(); ax.set_title("Cobertura SUBTE por barrio")
plt.tight_layout(); plt.show()
Pista: cmap='PuBuGn' o prueba Blues. Refiere: Matplotlib colormaps.
C.1 Nearest: barrio → estación más cercana (opcional)¶
# ⚠️ COMPLETA: usa sjoin_nearest para distancia mínima a una estación por barrio (centroides)
barrios_centroides = barrios_m.copy()
barrios_centroides["geometry"] = barrios_centroides.geometry.centroid
nearest = gpd.sjoin_nearest(
barrios_centroides[["BARRIO", "geometry"]],
estaciones[["geometry"]],
how="____", distance_col="____"
)
barrios_m = barrios_m.merge(nearest[["BARRIO", "dist_min"]], on="BARRIO", how="left")
barrios_m["dist_min_m"] = barrios_m["dist_min"].astype(float)
barrios_m[["BARRIO", "n_estaciones", "dist_min_m"]].sort_values("dist_min_m").head(10)
Pistas: how='left' y distance_col='dist_min'. Refiere: gpd.sjoin_nearest.
Parte D — Interactivo con Folium (opcional)¶
import folium
from folium import Choropleth, Marker
barrios_ll = barrios_m.to_crs(epsg=4326)
m = folium.Map(location=[-34.61, -58.44], tiles='____', zoom_start=11)
Choropleth(
geo_data=barrios_ll.__geo_interface__,
data=barrios_ll.set_index("BARRIO")["contactos_pc"],
key_on="feature.properties.BARRIO",
fill_color='____', legend_name='Contactos per cápita'
).add_to(m)
# Puntos de estaciones
for _, row in estaciones.to_crs(epsg=4326).iterrows():
Marker([row.geometry.y, row.geometry.x], icon=folium.Icon(color='red', icon='train', prefix='fa')).add_to(m)
m
Pistas: tiles='cartodbpositron' y fill_color='YlOrRd'. Refiere: folium.Map, folium.Choropleth.
📝 Preguntas de reflexión — Parte D (completa los espacios)¶
1) ¿Qué capa usarías para permitir toggles en el mapa (encender/apagar)? Respuesta: _. ¿Por qué? ___.
- Qué buscar: folium.FeatureGroup/LayerControl.
2) ¿Qué tiles mejorarían el contraste de una coropleta intensa? Opción: _____.
- Qué buscar: CartoDB.PositronNoLabels o Esri.WorldGrayCanvas.
3) ¿Qué popup mostrarías al clickear una estación? Campos: _, , __. - Qué buscar: nombre/ID de estación, línea, barrio.
Preguntas de reflexión finales (completa los espacios)¶
1) Usar CRS proyectado (EPSG: _) cambió las áreas/distancias porque ahora están en (unidades), evitando __.
2) Normalizar por km² y per cápita evita que el mapa refleje solo _. La métrica que más cambió tu interpretación fue porque __.
3) La zona con peor cobertura (mayor dist_min_m) es _; propondría evaluar una nueva estación en (coordenadas aproximadas) considerando (criterio) __.
4) Limitaciones técnicas detectadas: (i) Formato vectorial _ para capas grandes, (ii) encoding de CSV (); mitigaciones: __, _____.
5) Checklist de calidad que aplicarías en un pipeline productivo: CRS consistente, geometrías válidas, índices espaciales, tests de conteo/áreas, documentación de fuentes/licencias. Marca lo que cumpliste hoy: [ ] CRS [ ] geometrías [ ] índices [ ] tests [ ] docs.
Tareas extra (opcional)¶
1) Hexgrid/H3 para heatmaps comparables
- Objetivo: discretizar el espacio y agregar métricas por celdas hexagonales.
- Sugerencias: usar h3/h3pandas o generar una grilla hexagonal con Shapely; calcular contactos_pc por hex y comparar contra barrios.
- Hint: normaliza por superficie de la celda (tasas por km²) y ordena por percentiles.
2) Redes y accesibilidad con OSMnx
- Objetivo: estimar acceso a estaciones por red vial (distancia/costo real) y detectar gaps de cobertura.
- Sugerencias: osmnx para descargar grafo de calles, networkx para shortest path; comparar vs dist_min_m euclidiana.
- Hint: realiza snap de centroides de barrio y estaciones a nodos del grafo y usa pesos por longitud.
3) Overlays y zonas prohibidas
- Objetivo: excluir áreas de parque/agua/industrial y recalcular cobertura y priorización de nuevas estaciones.
- Sugerencias: descargar polígonos de uso del suelo (OSM/GeoJSON), usar gpd.overlay (difference/intersection).
- Hint: valida atributos preservados y controla duplicados tras la intersección.
4) Visualización avanzada interactiva
- Objetivo: publicar un mapa con toggles de capas y estilos.
- Sugerencias: folium.LayerControl, pydeck/kepler.gl para choropleth y hexagon layers.
- Hint: exporta a HTML y adjunta al portafolio con una breve guía de uso.
5) IO y performance
- Objetivo: acelerar lectura/escritura y reducir peso de capas.
- Sugerencias: pyogrio para IO rápido, GeoParquet con pyarrow, simplificación de geometrías (simplify).
- Hint: reporta tiempos y tamaño de archivos antes/después y justifica los trade-offs.