3.3.b Evaluación de la amenaza — Movimientos en masa - SIDRo

Esta página muestra la zonificación oficial de amenaza por movimientos en masa Alcaldía Mayor de Bogotá, 2021, los estudios básicos que la sustentan, la pendiente del terreno y el régimen histórico de lluvias extremas.

Zonificación de amenaza¶

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

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

Área por categoría de amenaza¶

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))

Las áreas se calculan sobre el AOI con su búfer; si las zonas de distintas categorías se traslapan, los porcentajes pueden sumar más que la superficie cubierta.

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

Estos estudios (escala 1:5.000) son el insumo técnico de la zonificación; aquí se muestran como contexto y no entran al cálculo del riesgo.

Pendiente¶

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"))

Índice topográfico de humedad (TWI)¶

El TWI, $\ln(a/\tan\beta)$ , combina el área que drena hacia cada punto ( $a$ ) con la pendiente local ( $\beta$ ): es alto donde el agua converge y se acumula (vaguadas, pie de ladera) y bajo en laderas empinadas y cumbres. Es un indicador de dónde se concentra la humedad del suelo que acompaña a la lluvia detonante Beven & Kirkby, 1979. Se calcula del DEM de 30 m con un encaminamiento D8, siguiendo la secuencia de morfometría de InundationRisk Krifai, 2025; no entra al cálculo del riesgo, 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

Con un DEM de 30 m el TWI solo resuelve vaguadas mayores que una manzana. El drenaje urbano (sumideros, alcantarillado) no está en el DEM, de modo que el índice describe el terreno natural, no el agua pluvial real.

Lluvia extrema histórica¶

Los movimientos en masa se activan con lluvia intensa o prolongada. Se calculan los índices ETCCDI Zhang et al., 2011: Rx1day (máximo diario), Rx5day (máximo de cinco días) y R95pTOT (lluvia en días muy húmedos).

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:::"))

Niveles de retorno del Rx1day¶

¿Qué lluvia de un día se excede, en promedio, una vez cada 10, 50 o 100 años? Se ajusta una distribución de valores extremos (GEV, y Gumbel como alternativa más estable en series cortas) a los máximos anuales, con intervalo de confianza del 95 % por remuestreo Coles, 2001.

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

La tabla es la del ajuste Gumbel. Con ~30 años de datos, el intervalo para 50 y 100 años es amplio; el nivel de 100 años es una extrapolación, no una observación. La lluvia es la de ERA5 (~25 km), que suaviza los picos convectivos: tome estos valores como piso, no como techo.

Una serie de reanálisis de ~25 km suaviza los aguaceros convectivos locales: los máximos diarios reales en la ladera pueden ser mayores. Si se dispone de una estación pluviométrica (IDEAM o la red SAB de IDIGER), úsela como nivel local de clima_observado_diario.

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
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
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
Coles, S. (2001). An Introduction to Statistical Modeling of Extreme Values. Springer. 10.1007/978-1-4471-3675-0