Visualización Geospacial Avanzada con GeoPandas

Implementación del proyecto final

En este notebook construyes, paso a paso, el proyecto final: un atlas urbano de México. Respondes tres preguntas con datos: ¿qué estados concentran más población urbana por kilómetro cuadrado?, ¿qué parte del territorio queda a menos de 100 km de las cinco ciudades más grandes?, y ¿cómo se reparte la población urbana entre tres franjas de latitud? Los datos son los estados y las ciudades de Natural Earth.

Una advertencia sobre los datos

La población de cada ciudad en Natural Earth (pop_max) es una estimación de su área metropolitana, y la capa solo incluye las ciudades principales (189 en México). Sirve para practicar las técnicas y comparar estados, pero no sustituye al censo del INEGI. Dilo en tus conclusiones: reconocer los límites de los datos es parte de un buen análisis.

Paso 1: datos limpios y proyectados

import geopandas as gpd
import pandas as pd

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"]).to_crs(6372)
estados["name"] = estados["name"].replace({"Distrito Federal": "Ciudad de México"})
estados["area_km2"] = estados.area / 1e6
ciudades = gpd.read_file(NE + "ne_10m_populated_places_simple.zip",
                         where="iso_a2 = 'MX'",
                         columns=["name", "pop_max", "iso_a2"]).to_crs(6372)
ciudades["name"] = ciudades["name"].replace({"Mexico City": "Ciudad de México"})
print(len(estados), "estados,", len(ciudades), "ciudades")

Paso 2: unión espacial y densidad

sjoin asigna a cada ciudad el estado que la contiene; después agregas por estado y unes el resultado a la capa de polígonos. Si algún estado no tuviera ciudades en la capa, fillna(0) le asigna cero en lugar de dejarlo vacío (con estos datos los 32 tienen al menos una, pero el paso protege el análisis si cambias el filtro):

union = gpd.sjoin(ciudades, estados[["name", "geometry"]],
                  predicate="within", lsuffix="ciudad", rsuffix="estado")
por_estado = union.groupby("name_estado").agg(
    ciudades=("name_ciudad", "size"), poblacion=("pop_max", "sum"))
atlas = estados.merge(por_estado, left_on="name", right_index=True, how="left")
atlas[["ciudades", "poblacion"]] = atlas[["ciudades", "poblacion"]].fillna(0)
atlas["densidad"] = atlas["poblacion"] / atlas["area_km2"]
print(atlas.nlargest(5, "densidad")[["name", "densidad"]].round(1))

Paso 3: la coropleta

La Ciudad de México tiene una densidad tan alta que, con una escala continua, todos los demás estados quedarían del mismo color. Por eso usas scheme="NaturalBreaks" (o "FisherJenks", su versión exacta y reproducible): la capital ocupa su propia clase y el resto se distingue. Encima dibujas las diez ciudades más pobladas.

Paso 4: área de influencia de las grandes ciudades

Un buffer de 100 km alrededor de cada una de las cinco ciudades más pobladas, unido con union_all() para no contar dos veces las zonas que se traslapan, y recortado al país. La proporción del territorio cubierto es un cociente de áreas:

pais = estados.dissolve()
zonas = ciudades.nlargest(5, "pop_max").buffer(100_000).union_all()
zonas = gpd.GeoSeries([zonas], crs=6372).clip(pais)
cobertura = zonas.area.sum() / pais.area.sum()
print(f"Territorio a menos de 100 km de las 5 mayores ciudades: {cobertura:.1%}")

Paso 5: franjas de latitud con dissolve

Divides el país en tres franjas según la latitud del centroide de cada estado (sur, por debajo de 18.5° N; centro, hasta 23° N; norte, el resto) y las unes con dissolve, sumando población y área. Son franjas geométricas, no regiones oficiales: la península de Yucatán, por ejemplo, cae en la franja central. La tabla resultante permite comparar la densidad urbana de cada franja:

latitud = atlas.centroid.to_crs(4326).y
atlas["region"] = pd.cut(latitud, [0, 18.5, 23, 33], labels=["sur", "centro", "norte"])
regiones = atlas.dissolve(by="region", aggfunc={"poblacion": "sum", "area_km2": "sum"},
                          observed=True)
regiones["densidad"] = regiones["poblacion"] / regiones["area_km2"]
print(regiones[["poblacion", "area_km2", "densidad"]].round(1))

Paso 6: mapa interactivo y conclusiones

Exportas el atlas con explore() a un archivo HTML, con la densidad como color y el nombre, las ciudades y la población en la ventana emergente. Cierras con tres a cinco frases que respondan las preguntas iniciales, cada una con su número.

Trampas comunes

• Hacer el sjoin con capas en CRS distintos.
• Sumar áreas de buffers sin unirlos: las zonas traslapadas se cuentan dos veces.
• Olvidar fillna(0) tras la unión: con estos datos los 32 estados tienen al menos una ciudad, pero si filtras ciudades (por ejemplo, solo las de más de un millón de habitantes) los estados sin coincidencias quedan con NaN y desaparecen de la coropleta.
• Presentar pop_max como población oficial.