Movimientos en masa en Los Rosales — cuaderno integrado - SIDRo

Resumen¶

Los movimientos en masa son la amenaza característica de las laderas de los Cerros Orientales. Este cuaderno superpone la zonificación oficial del POT Alcaldía Mayor de Bogotá, 2021 con las edificaciones y los equipamientos del barrio, con la lógica de exposición de CCDR-tools Amadio, 2024. Luego analiza la lluvia que detona los deslizamientos con índices ETCCDI Zhang et al., 2011, sus niveles de retorno y su cambio bajo escenarios climáticos.

Resultados clave (ejecución del 2026-10-27 con datos reales):

El área con búfer tiene 24 ha en amenaza alta (2,9 %) y 187 ha en amenaza media (22,4 %).
992 edificaciones están en amenaza alta y 5.386 en amenaza media. Tres colegios y 2,8 km de vías quedan en amenaza alta. Dentro del barrio, unos 3.500 de sus 13.700 residentes estimados viven en amenaza alta y unos 10.000 en amenaza media.
La lluvia de un día con periodo de retorno de 100 años es de unos 102 mm (IC 95 %: 77–123 mm).
La mediana de los modelos CMIP6 proyecta menos lluvia media (−7 % a −14 %). La dispersión es grande: el percentil 90 entre modelos llega a +20 %.

Las páginas de detalle a, b, c y d explican cada paso con más extensión.

1. Marco¶

1.1 La amenaza y el barrio¶

Los movimientos en masa se activan cuando se combinan pendiente, material (geología y geomorfología), agua e intervención humana. Los Rosales está en el contacto entre la ciudad consolidada y la ladera. La zonificación del POT, derivada de los estudios básicos de IDIGER a escalas de 1:2.000 a 1:5.000, es la fuente principal. SIDRo no construye un modelo propio de susceptibilidad.

1.2 Cadena de análisis¶

Como la zonificación es categórica y no da probabilidades, no se calcula un impacto anual esperado. Se reporta la exposición por categoría.

2. Configuración del entorno¶

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from IPython.display import Markdown, display

from sidro import process, viz
from sidro.analisis import clima, espacial, masa
from sidro.aoi import aoi_bufer
from sidro.config import parametros
from sidro.fetch import obtener_capa

amenaza = intentar(obtener_capa, "amenaza_mov_masa",
                   mensaje="La zonificación de amenaza por movimientos en masa no está disponible")
if amenaza is not None and len(amenaza):
    amenaza = masa.preparar(amenaza)   # traduce los códigos numéricos del POT a Baja / Media / Alta
    campo = "categoria"
else:
    campo = None

3. Amenaza¶

3.1 Zonificación del POT¶

if amenaza is not None and len(amenaza) and campo:
    m = viz.mapa_base()
    viz.capa(m, amenaza, campo, nombre="Amenaza por movimientos en masa", categorica=True, cmap="RdYlGn_r")
    viz.cerrar(m)
else:
    if amenaza is not None:
        display(Markdown(":::{warning}\nLa capa no tiene polígonos con categoría en el área de interés: el barrio puede "
                         "estar fuera de la zonificación urbana consultada (revise el suelo rural).\n:::"))
    m = None
m

if amenaza is not None and len(amenaza) and campo:
    area = aoi_bufer().geometry.iloc[0]
    z = amenaza.to_crs(aoi_bufer().crs)[[campo, "geometry"]].copy()
    z = z.assign(categoria=z[campo].astype(str).str.title())
    z["geometry"] = z.geometry.intersection(area)
    tabla = (z.dissolve(by="categoria").area / 1e4).rename("área (ha)").to_frame()
    tabla["% del área con búfer"] = 100 * tabla["área (ha)"] / (area.area / 1e4)
    display(tabla.round(1))

3.2 Estudios básicos: geología y geomorfología¶

geo = intentar(obtener_capa, "geologia_ingenieria", mensaje="Geología para ingeniería no disponible")
geom = intentar(obtener_capa, "geomorfologia", mensaje="Geomorfología no disponible")
for nombre, g in (("Geología para ingeniería", geo), ("Geomorfología", geom)):
    if g is not None and len(g):
        texto = [c for c in g.columns if c != "geometry" and g[c].dtype == "string" and g[c].nunique() < 40]
        print(f"{nombre}: {len(g)} polígonos; campos descriptivos: {', '.join(texto[:6]) or 'ninguno detectado'}")

Geología para ingeniería: 29 polígonos; campos descriptivos: LEY_GEO, DESCRIPCIO
Geomorfología: 25 polígonos; campos descriptivos: UNIDAD, GEOMORFOES, PROVINCIA, REGION, SUBUNIDAD, COMPONENTE

if geo is not None and len(geo):
    texto = [c for c in geo.columns if c != "geometry" and geo[c].dtype == "string" and geo[c].nunique() < 40]
    m = viz.mapa_base()
    viz.capa(m, geo, texto[0] if texto else None, nombre="Geología", categorica=bool(texto), cmap="Set3")
    viz.cerrar(m)
else:
    m = None
m

3.3 Pendiente e índice topográfico de humedad¶

ruta_dem = intentar(obtener_capa, "dem", mensaje="No se pudo obtener el modelo digital de elevación")
if ruta_dem is not None:
    pend = process.pendiente_grados(process.abrir_raster(ruta_dem))
    v = pend.values[np.isfinite(pend.values)]
    display(pd.Series({"pendiente media (°)": v.mean(), "pendiente máxima (°)": v.max(),
                       "área con pendiente > 25° (%)": 100 * (v > 25).mean(),
                       "área con pendiente > 35° (%)": 100 * (v > 35).mean()}).round(1).to_frame("área con búfer"))

El TWI, $\ln(a/\tan\beta)$ , es alto donde el agua converge (vaguadas, pie de ladera) Beven & Kirkby, 1979. Se calcula con la secuencia de morfometría de InundationRisk Krifai, 2025 y solo da contexto.

if ruta_dem is not None:
    from sidro.aoi import construir_aoi

    humedad = espacial.twi(process.abrir_raster(ruta_dem))
    fig, ax = plt.subplots(1, 2, figsize=(10, 4.2))
    for a_, da, titulo, cmap in ((ax[0], pend, "Pendiente (°)", "magma_r"), (ax[1], humedad, "TWI", "YlGnBu")):
        da.plot(ax=a_, cmap=cmap, cbar_kwargs={"shrink": 0.75})
        construir_aoi().boundary.plot(ax=a_, color="#111", lw=1.5)
        a_.set_title(titulo), a_.set_aspect("equal"), a_.set_xlabel(""), a_.set_ylabel("")
    fig.tight_layout()
    fig

Lectura. La amenaza alta ocupa poca área (24 ha), pero el terreno es empinado: un 27,5 % del área con búfer tiene pendiente mayor de 25° y un 8,9 % mayor de 35°. Con un DEM de 30 m, el TWI solo resuelve vaguadas mayores que una manzana y no ve el drenaje urbano.

3.4 Lluvia extrema histórica y niveles de retorno¶

serie = intentar(obtener_capa, "clima_observado_diario", mensaje="No se pudo obtener la serie climática diaria")
if serie is not None and "precipitation_sum" in serie.columns:
    idx_lluvia = clima.indices_precipitacion(serie["precipitation_sum"])
    display(idx_lluvia[parametros()["movimientos_masa"]["indices_lluvia"]].describe().loc[["mean", "min", "max"]].round(1))
    fig, ax = plt.subplots(figsize=(7, 3.2))
    for c in parametros()["movimientos_masa"]["indices_lluvia"][:2]:
        ax.plot(idx_lluvia.index, idx_lluvia[c], label=c)
    ax.set_xlabel("Año"), ax.set_ylabel("mm"), ax.legend(), ax.grid(alpha=0.3)
    fig
elif serie is not None:
    display(Markdown(":::{warning}\nLa serie no tiene la columna `precipitation_sum`.\n:::"))

retorno = None
if serie is not None and "precipitation_sum" in serie.columns:
    mx = idx_lluvia["rx1day"]
    retorno = {d: clima.niveles_retorno_gev(mx, distribucion=d) for d in ("gev", "gumbel")}
    display(retorno["gumbel"].round(1).rename(columns={"periodo": "periodo de retorno (años)", "nivel": "Rx1day (mm)"})
            .set_index("periodo de retorno (años)"))
    fig, ax = plt.subplots(figsize=(6.5, 3.4))
    for d, color in (("gev", "#3a6ea5"), ("gumbel", "#d1495b")):
        r = retorno[d]
        ax.plot(r["periodo"], r["nivel"], color=color, label=d.upper())
        ax.fill_between(r["periodo"], r["ic_inf"], r["ic_sup"], color=color, alpha=0.15)
    orden = np.sort(mx.dropna().values)[::-1]
    ax.plot((len(orden) + 1) / np.arange(1, len(orden) + 1), orden, "ko", ms=3, label="máximos observados")
    ax.set_xscale("log"), ax.set_xlabel("Periodo de retorno (años)"), ax.set_ylabel("Rx1day (mm)")
    ax.legend(), ax.grid(alpha=0.3)
    fig

Lectura. En promedio, la lluvia máxima de un día es de 43 mm y la de cinco días de 93 mm. El ajuste Gumbel da 67 mm para 10 años y 102 mm para 100 años. El nivel de 100 años es una extrapolación con unos 30 años de datos. ERA5 (~25 km) suaviza los aguaceros convectivos, así que estos valores son un piso.

4. Exposición y riesgo¶

import geopandas as gpd
import numpy as np
import pandas as pd
from IPython.display import Markdown, display

from sidro import exposicion, process, viz
from sidro.analisis import indices, masa
from sidro.aoi import aoi_bufer
from sidro.config import rutas
from sidro.fetch import obtener_capa

amenaza = intentar(obtener_capa, "amenaza_mov_masa", mensaje="La zonificación de amenaza no está disponible")
edif = intentar(exposicion.edificaciones, mensaje="No se pudo construir la exposición de edificaciones")
if amenaza is not None and len(amenaza):
    amenaza = masa.preparar(amenaza)   # traduce los códigos numéricos del POT a Baja / Media / Alta
    campo = "categoria"
else:
    campo = None
listo = edif is not None and campo is not None
if amenaza is not None and campo is None:
    display(Markdown(":::{warning}\nLa capa de amenaza no tiene polígonos o un campo de categoría reconocible.\n:::"))

4.1 Edificaciones por categoría¶

if listo:
    edif["categoria"] = masa.clasificar(edif, amenaza, campo)
    tabla = masa.resumen_exposicion(edif, edif["categoria"], {
        "edificaciones": ("id_edif", "count"), "personas": ("ocupantes_est", "sum"),
        "valor_millones_cop": ("valor_cop", lambda s: s.sum() / 1e6)})
    tabla.loc["Total"] = tabla.sum()
    display(tabla.round(0).astype(int))

if listo:
    m = viz.mapa_base()
    viz.capa(m, amenaza, campo, nombre="Amenaza", categorica=True, cmap="RdYlGn_r")
    viz.capa(m, edif, "categoria", nombre="Edificaciones", categorica=True, cmap="RdYlGn_r",
             tooltip=["categoria", "pisos", "anio", "uso"])
    viz.cerrar(m)
else:
    m = None
m

4.2 Vías y equipamientos críticos¶

vias = intentar(obtener_capa, "red_vial", mensaje="No se pudo obtener la red vial de OpenStreetMap")
equip = intentar(obtener_capa, "equipamientos_criticos", mensaje="No se pudo obtener los equipamientos críticos")

if listo and vias is not None and len(vias):
    vias_aoi = gpd.clip(vias, aoi_bufer().geometry.iloc[0])
    km = masa.longitud_por_categoria(vias_aoi, amenaza, campo)
    display(km.round(2).rename("km de vía").to_frame())
if listo and equip is not None and len(equip):
    equip = equip.copy()
    equip["categoria"] = masa.clasificar(equip, amenaza, campo)
    cols = [c for c in ("name", "amenity", "categoria") if c in equip.columns]
    display(equip[equip.categoria.isin(masa.orden_tabla()[:2])][cols].reset_index(drop=True))

Lectura. El 5,9 % de las edificaciones está en amenaza alta y casi un tercio en amenaza media. Los colegios Nueva Granada, Jordán de Sajonia y Rosario de Santo Domingo están en amenaza alta. De la red vial del área con búfer, 2,8 km están en amenaza alta y 35 km en amenaza media. El inventario de equipamientos de OpenStreetMap es incompleto y hay que contrastarlo con IDIGER. «Sin categoría» no significa libre de amenaza.

La medida de riesgo por unidad son las personas estimadas en amenaza alta.

if listo:
    orden = masa.orden_tabla()
    edif["personas_alta"] = edif["ocupantes_est"].where(edif["categoria"] == orden[0], 0.0)
    edif["personas_media"] = edif["ocupantes_est"].where(edif["categoria"] == orden[1], 0.0)
    edif["edif_alta"] = (edif["categoria"] == orden[0]).astype(int)
    unidades = process.unidades_analisis()
    agg = process.agregar(edif, unidades, {
        "n_edificaciones": ("id_edif", "count"), "edificaciones_alta": ("edif_alta", "sum"),
        "personas_alta": ("personas_alta", "sum"), "personas_media": ("personas_media", "sum")})
    agg["vulnerabilidad_social"] = indices.vulnerabilidad_social(agg)
    print(f"Unidades: {len(agg)} (origen: {unidades.attrs.get('origen')})")
    display(agg.drop(columns="geometry").sort_values("personas_alta", ascending=False).head(10).round(1))

WARNING sidro.fetch: manzanas_censo: falló nivel nacional/usuario -> Archivo de usuario no encontrado: /Users/mauricio/Documents/Data-science-local/SIDRo/data/user/dane_mgn_manzanas_cnpv2018.gpkg

WARNING sidro.process: Unidades: manzanas_censo no disponible (No se pudo obtener 'manzanas_censo'. Intentos: [{'nivel': 'nacional', 'tipo': 'usuario', 'error': 'Archivo de usuario no encontrado: /Users/mauricio/Documents/Data-science-local/SIDRo/data/user/dane_mgn_manzanas_cnpv2018.gpkg'}])

Unidades: 57 (origen: estratificacion)

if listo:
    biv = indices.bivariado(agg["personas_alta"].fillna(0), agg["vulnerabilidad_social"])
    agg = agg.join(biv)
    m = viz.mapa_base()
    viz.capa_bivariada(m, agg)
    viz.cerrar(m)
else:
    m = None
m

Leyenda: exposición a movimientos en masa en el eje vertical, vulnerabilidad social en el horizontal.

if listo:
    process.guardar(agg[["id_unidad", "geometry", "personas_alta", "edificaciones_alta", "personas_media",
                         "vulnerabilidad_social"]].rename(columns={"personas_alta": "riesgo"}), "riesgo_mov_masa")
    print("Guardado:", (rutas().processed / "riesgo_mov_masa.parquet").relative_to(rutas().raiz))

Guardado: data/processed/riesgo_mov_masa.parquet

5. Escenarios climáticos¶

La zonificación no cambia con el clima, pero la lluvia detonante sí. Se aplica el cambio de la precipitación de NEX-GDDP-CMIP6 Thrasher et al., 2022 a la serie local bajo SSP2-4.5 y SSP5-8.5 Riahi et al., 2017O'Neill et al., 2017.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from IPython.display import display

from sidro.analisis import clima
from sidro.config import aoi_cfg, parametros
from sidro.fetch import obtener_capa

MODELOS = None  # p. ej. ["ACCESS-CM2", "EC-Earth3"] para una prueba rápida
cfg = aoi_cfg()["clima"]
serie = intentar(obtener_capa, "clima_observado_diario", mensaje="No se pudo obtener la serie climática diaria")
pr = serie["precipitation_sum"] if serie is not None and "precipitation_sum" in serie.columns else None
deltas = None
if pr is not None:
    deltas = intentar(clima.deltas_tabla, ["pr"], cfg["escenarios"], cfg["horizontes"], MODELOS,
                      cache="deltas_mov_masa",
                      mensaje="No se pudieron obtener las proyecciones CMIP6 (requiere acceso a Planetary Computer)")
if deltas is not None:
    display(deltas.groupby(["escenario", "horizonte"])["delta"].describe()[["count", "min", "50%", "max"]].round(3))

cambios = None
if deltas is not None:
    indices_lluvia = parametros()["movimientos_masa"]["indices_lluvia"]
    q = cfg["percentiles_ensamble"]
    filas = []
    for (esc, hor), g in deltas.groupby(["escenario", "horizonte"]):
        razon = pd.DataFrame([clima.cambio_lluvia_extrema(pr, d) for d in g["delta"]])[indices_lluvia]
        for ind in indices_lluvia:
            filas.append({"escenario": esc, "horizonte": hor, "índice": ind,
                          **{f"p{x}": 100 * (np.percentile(razon[ind], x) - 1) for x in q}})
    cambios = pd.DataFrame(filas)
    display(cambios.round(1))

if cambios is not None:
    fig, ax = plt.subplots(figsize=(7, 3.4))
    r1 = cambios[cambios["índice"] == "rx1day"]
    x = np.arange(len(r1))
    med = r1[f"p{q[1]}"].values
    ax.bar(x, med, yerr=[med - r1[f"p{q[0]}"].values, r1[f"p{q[2]}"].values - med], color="#3a6ea5", capsize=4)
    ax.set_xticks(x, [f"{e}\n{h}" for e, h in zip(r1.escenario, r1.horizonte)], fontsize=8)
    ax.set_ylabel("Cambio en Rx1day (%)"), ax.axhline(0, color="#555", lw=0.8)
    fig

niveles_futuros = None
if cambios is not None and "precipitation_sum" in serie.columns:
    hist = clima.niveles_retorno_gev(clima.indices_precipitacion(serie["precipitation_sum"])["rx1day"], distribucion="gumbel")
    r1 = cambios[cambios["índice"] == "rx1day"].set_index(["escenario", "horizonte"])
    cols = {"histórico": hist.set_index("periodo")["nivel"]}
    for (esc, hor), fila in r1.iterrows():
        cols[f"{esc} {hor}"] = hist.set_index("periodo")["nivel"] * (1 + fila[f"p{q[1]}"] / 100)
    niveles_futuros = pd.DataFrame(cols).rename_axis("periodo de retorno (años)")
    display(niveles_futuros.round(0))
    fig, ax = plt.subplots(figsize=(6.5, 3.4))
    for col in niveles_futuros:
        ax.plot(niveles_futuros.index, niveles_futuros[col], label=col, lw=2 if col == "histórico" else 1.2,
                color="#333" if col == "histórico" else None)
    ax.set_xscale("log"), ax.set_xlabel("Periodo de retorno (años)"), ax.set_ylabel("Rx1day (mm)")
    ax.legend(fontsize=7), ax.grid(alpha=0.3)
    fig

Lectura. La mediana de los cinco modelos reduce la lluvia, y por construcción también Rx1day y Rx5day (−7 % a −14 %). El rango p10–p90 cruza el cero en tres de los cuatro casos: no hay una señal robusta de cambio en la lluvia detonante. Solo SSP5-8.5 a final de siglo es negativo en todo el rango. Como el método no capta la intensificación de los extremos, no debe leerse como una disminución de la amenaza.

if cambios is not None:
    from sidro import export

    export.guardar_resumen("mov_masa_escenarios", {"cambio_lluvia_pct": cambios.to_dict("records")})

6. Tabla resumen para el informe ejecutivo¶

Indicador	Valor (ejecución 2026-10-25)
Área en amenaza alta / media (área con búfer)	24,3 ha (2,9 %) / 187,2 ha (22,4 %)
Edificaciones en amenaza alta / media	992 / 5.386
Residentes en amenaza alta (barrio / búfer)	≈ 3.500 / ≈ 7.500
Valor expuesto en amenaza alta	≈ 2,6 billones COP
Vías en amenaza alta / media	2,8 km / 35,3 km
Equipamientos en amenaza alta (OSM)	3 colegios
Rx1day, periodo de retorno 10 / 100 años	67 mm / 102 mm
Cambio de Rx1day, SSP5-8.5 2081–2100 (p10 – p50 – p90)	−24 % – −14 % – −7 %

7. Conclusiones¶

La exposición a movimientos en masa se concentra en las manzanas del borde oriental. Las edificaciones y los colegios en amenaza alta son la prioridad para vigilancia y estudios de detalle.
La lluvia detonante del barrio es del orden de 70 a 100 mm en un día para eventos de 10 a 100 años.
Las proyecciones no muestran un aumento de la lluvia media. Con el método actual eso no descarta aguaceros más intensos.

7.1 Limitaciones¶

La zonificación es categórica y no da probabilidad: solo se reporta exposición.
Hay dos capas del POT con el mismo nombre (ids 69 y 80). Hay que confirmar cuál es la vigente.
Un elemento toma la categoría más alta que toca, lo que sobrestima un poco la exposición de los bordes.
Los residentes se estiman desde el área habitacional (110 m² brutos por residente), sin censo por manzana.
No se modela la humedad antecedente ni los umbrales de lluvia de cada ladera. El resultado es indicativo.

7.2 Próximos pasos¶

Usar una estación pluviométrica de IDEAM o de la red SAB de IDIGER como nivel local de la serie.
Estimar la intensificación de los extremos con un método que no escale todos los días por igual.

References¶

Alcaldía Mayor de Bogotá. (2021). Decreto 555 de 2021, por el cual se adopta la revisión general del Plan de Ordenamiento Territorial de Bogotá D.C. https://services7.arcgis.com/lsxbLWF2l19Rmhqj/ArcGIS/rest/services/POT_Bogota_Decreto_555_2021/FeatureServer
Amadio, M. (2024). CCDR tools data and methodology. Global Facility for Disaster Reduction and Recovery (GFDRR), World Bank. https://gfdrr.github.io/CCDR-tools
Zhang, X., Alexander, L., Hegerl, G. C., & others. (2011). Indices for monitoring changes in extremes based on daily temperature and precipitation data. WIREs Climate Change, 2, 851–870. 10.1002/wcc.147
Beven, K. J., & Kirkby, M. J. (1979). A physically based, variable contributing area model of basin hydrology. Hydrological Sciences Bulletin, 24(1), 43–69. 10.1080/02626667909491834
Krifai, M. (2025). InundationRisk: modular framework for regional flood risk assessment. https://github.com/mkrifai/inundation-risk-model
Thrasher, B., Wang, W., Michaelis, A., Melton, F., Lee, T., & Nemani, R. (2022). NASA Global Daily Downscaled Projections, CMIP6. Scientific Data, 9, 262. 10.1038/s41597-022-01393-4
Riahi, K., van Vuuren, D. P., Kriegler, E., & others. (2017). The Shared Socioeconomic Pathways and their energy, land use, and greenhouse gas emissions implications: An overview. Global Environmental Change, 42, 153–168. 10.1016/j.gloenvcha.2016.05.009
O’Neill, B. C., Kriegler, E., Ebi, K. L., & others. (2017). The roads ahead: Narratives for shared socioeconomic pathways describing world futures in the 21st century. Global Environmental Change, 42, 169–180. 10.1016/j.gloenvcha.2015.01.004