3.2.b Evaluación de la amenaza — Incendios forestales

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 bbox_wgs84
from sidro.config import parametros
from sidro.fetch import obtener_capa

pot = intentar(obtener_capa, "amenaza_incendio_forestal",
               mensaje="La capa de amenaza por incendios forestales del POT no está disponible")

if pot is not None and len(pot):
    pot = masa.preparar(pot)   # traduce los códigos numéricos del POT (1 Baja, 2 Media, 3 Alta ...)
    m = viz.mapa_base()
    viz.capa(m, pot, "categoria", nombre="Amenaza POT", categorica=True, cmap="OrRd")
    viz.cerrar(m)
elif pot is not None:
    display(Markdown("La capa del POT no tiene polígonos en el área de interés: la interfaz urbano-forestal del barrio "
                     "se evalúa con la distancia a la vegetación de [3.2.c](c_riesgo.md)."))
    m = None
else:
    m = None
m

ruta_cob = intentar(obtener_capa, "cobertura_suelo", mensaje="No se pudo obtener la cobertura del suelo")
if ruta_cob is not None:
    cob = process.abrir_raster(ruta_cob)
    clases = parametros()["incendio_forestal"]["clases_bosque_worldcover"]
    bosque = espacial.poligonos_de_clase(cob, clases)
    area_ha = bosque.area.sum() / 1e4
    print(f"Vegetación combustible en el área con búfer: {area_ha:,.1f} ha "
          f"({100 * np.isin(cob.values, clases).sum() / np.isfinite(cob.values).sum():.0f} % de las celdas)")
    m = viz.mapa_base()
    viz.capa(m, bosque, nombre="Vegetación combustible", color="#2d6a4f")
    viz.cerrar(m)
else:
    bosque, m = None, None
m

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

focos = intentar(obtener_capa, "focos_calor",
                 mensaje="Sin focos de calor: defina la variable de entorno FIRMS_MAP_KEY para consultarlos (opcional)")
if focos is not None:
    print(f"{len(focos)} focos de calor detectados por VIIRS en el área con búfer en los últimos días consultados")

import matplotlib.pyplot as plt

clima_df = intentar(obtener_capa, "clima_observado_diario", mensaje="No se pudo obtener la serie climática diaria")
if clima_df is not None:
    lat = (bbox_wgs84()[1] + bbox_wgs84()[3]) / 2
    fwi = intentar(clima.fwi_diario, clima_df, lat, mensaje="No se pudo calcular el FWI (faltan variables en la serie)")
else:
    fwi = None

if fwi is not None:
    umbral = parametros()["incendio_forestal"]["fwi_umbral_alto"]
    dias = clima.dias_fwi_alto(fwi)
    resumen = pd.DataFrame({"FWI medio": fwi.groupby(fwi.index.year).mean(),
                            "FWI máximo": fwi.groupby(fwi.index.year).max(), "días FWI alto": dias})
    display(resumen.describe().loc[["mean", "min", "max"]].round(1))

    fig, ax = plt.subplots(1, 2, figsize=(10, 3.4))
    mensual = fwi.groupby(fwi.index.month).mean()
    ax[0].bar(mensual.index, mensual.values, color="#bc6c25")
    ax[0].set_xlabel("Mes"), ax[0].set_ylabel("FWI medio"), ax[0].set_title("Ciclo anual")
    ax[1].plot(resumen.index, resumen["FWI máximo"], color="#bc6c25")
    ax[1].axhline(umbral, ls="--", color="#555")
    ax[1].set_xlabel("Año"), ax[1].set_ylabel("FWI máximo anual"), ax[1].set_title(f"Umbral alto = {umbral}")
    fig.tight_layout()
    fig