Superficie de humedales urbanos protegidos

import datetime
now = datetime.datetime.now()
print(f"Última actualización de este notebook: {now}")

import math
import time
import requests
import pandas as pd
import geopandas as gpd
from shapely.geometry import shape
import matplotlib.pyplot as plt
import numpy as np

import math
import time
import requests
import pandas as pd
import geopandas as gpd
from shapely.geometry import shape

# Obtener total de registros
params_count = dict(BASE_PARAMS)
params_count.update({
    "returnCountOnly": "true",
    "returnGeometry": "false",
})
r = requests.get(LAYER_URL, params=params_count, headers=HEADERS, timeout=120)
r.raise_for_status()
data = r.json()
total = int(data.get("count", 0))
total

# Función para convertir geometría ArcGIS a Shapely
# El output es una tabla que debiese quedar implementada como base de datos con metadatos para el control de verisones y cambios en la información, pero por ahora se deja como un dataframe para 
# facilitar su uso en análisis posteriores.

from shapely.geometry import Point, MultiPoint, LineString, MultiLineString, Polygon, MultiPolygon

def arcgis_geom_to_shapely(geom):
    if not geom:
        return None
    if "x" in geom and "y" in geom:
        return Point(geom["x"], geom["y"])
    if "points" in geom:
        return MultiPoint(geom["points"])
    if "paths" in geom:
        paths = geom["paths"]
        if len(paths) == 1:
            return LineString(paths[0])
        return MultiLineString(paths)
    if "rings" in geom:
        rings = geom["rings"]
        if len(rings) == 1:
            return Polygon(rings[0])
        return Polygon(rings[0], holes=rings[1:])
    return None

def fetch_page(offset: int, page_size: int = 2000) -> gpd.GeoDataFrame:
    params = dict(BASE_PARAMS)
    params.update({
        "resultOffset": offset,
        "resultRecordCount": page_size,
        "orderByFields": "OBJECTID ASC",
    })
    r = requests.get(LAYER_URL, params=params, headers=HEADERS, timeout=120)
    r.raise_for_status()
    data = r.json()

    if "error" in data:
        raise RuntimeError(f"Error ArcGIS: {data['error']}")

    feats = data.get("features", [])
    rows = []
    for f in feats:
        attrs = f.get("attributes", {})
        geom = f.get("geometry")
        attrs["geometry"] = arcgis_geom_to_shapely(geom)
        rows.append(attrs)

    return gpd.GeoDataFrame(rows, geometry="geometry", crs="EPSG:4326")

PAGE_SIZE = 2000
pages = math.ceil(total / PAGE_SIZE)

gdfs = []
for i in range(pages):
    offset = i * PAGE_SIZE
    gdf_i = fetch_page(offset, PAGE_SIZE)
    print(f"Página {i+1}/{pages}: {len(gdf_i)} filas (offset={offset})")
    gdfs.append(gdf_i)
    time.sleep(0.2)

gdf = gpd.GeoDataFrame(pd.concat(gdfs, ignore_index=True)) if gdfs else gpd.GeoDataFrame()

# Reparar geometrías si es necesario
try:
    gdf["geometry"] = gdf.geometry.make_valid()
except Exception:
    gdf["geometry"] = gdf.geometry.buffer(0)

# Calcular en CRS en metros (ajusta según zona)
gdf_m = gdf.to_crs(32719)  # UTM 19S para Chile central; usa 32718 si es más al norte

gdf["Shape_Area"] = gdf_m.area
gdf["Shape_Leng"] = gdf_m.length

gdf.head()

df=gdf
df["anio_publicacion"] = (
    df["RESOLUCION"]
    .str.extract(r"/\s*(\d{4})", expand=False)
    .astype("Int64")
)

df["anio_publicacion"].isna().any()
df["anio_publicacion"].agg(["min", "max"])

# desde df se va a calcular el total de superficie protegida para mostrarla como sin decimales

total_superficie = df["Shape_Area"].sum() / 10000  # convertir a hectáreas
total_superficie = int(total_superficie)    

df["HECTAREAS"] = pd.to_numeric(df["HECTAREAS"], errors="coerce")

tabla_anual = (
    df
    .groupby("anio_publicacion", as_index=False)["HECTAREAS"]
    .sum()
    .sort_values("anio_publicacion")
)

tabla_anual["hectareas_acumuladas"] = (
    tabla_anual["HECTAREAS"].cumsum()
)
tabla_anual

from IPython.display import Markdown as md
text = f""" 
A la fecha de reporte, el Estado de Chile ha logrado cumplir y posicionarse en el límite superior de la meta nacional:
- Superficie Protegida: Se registran {total_superficie} hectáreas de humedales urbanos reconocidos que mantienen plena vigencia legal.
- Número de Humedales: El país cuenta actualmente con 137 humedales urbanos activos (de un total de 147 declaraciones originales).
- Evolución Temporal: Desde la promulgación de la ley, la superficie acumulada ha experimentado un crecimiento sostenido: 5.471,97 ha (2021); 7.498,45 ha (2022); 9.519,14 ha (2023); 10.717,06 ha (2024); logrando el gran salto a las 16.848,28 ha en 2025."""
md(text)

#escribir tabla anual a Excel
tabla_anual.to_excel("I_3_humedales_urbanos_tabla_anual.xlsx", index=False)

data=gdf.to_crs(epsg=32719)

REGION = {
    'REGIÓN DE TARAPACÁ': "TA",
    'REGIÓN DE ANTOFAGASTA': "AN",
    'REGIÓN DE ATACAMA': "AT",
    'REGIÓN DE COQUIMBO': "CO",
    'REGIÓN DE VALPARAÍSO': "VA",
    'REGIÓN METROPOLITANA DE SANTIAGO': "RM",
    "REGIÓN DEL LIBERTADOR GENERAL BERNARDO O'HIGGINS": "LI",
    'REGIÓN DEL MAULE': "ML",
    'REGIÓN DEL ÑUBLE': "NB",
    'REGIÓN DEL BIOBÍO': "BI",
    'REGIÓN DE LA ARAUCANÍA': "AR",
    'REGIÓN DE LOS RÍOS': "LR",
    'REGIÓN DE LOS LAGOS': "LL",
    'REGIÓN DE AYSÉN DEL GENERAL CARLOS IBÁÑEZ DEL CAMPO': "AI",
    'REGIÓN DE MAGALLANES Y DE LA ANTÁRTICA CHILENA': "MA",
}

data["COD_REG"] = data["REGION"].map(REGION)

data_region = pd.DataFrame(
    list(REGION.items()),
    columns=["Región", "Código Región"]
)
data_region["geometry"] = None
data_region["Hectáreas"] = 0.0
data_region["Hectáreas (geometría)"] = 0.0
data_region["Número de humedales"] = 0

for region in REGION.values():
    tmp = data[data["COD_REG"] == region]
    tmp_geometry = tmp.geometry.union_all()
    cur_index = data_region[data_region["Código Región"] == region].index[0]
    data_region.loc[cur_index, "geometry"] = tmp_geometry
    data_region.loc[cur_index, "Hectáreas"] = tmp["HECTAREAS"].sum()
    data_region.loc[cur_index, "Hectáreas (geometría)"] = tmp_geometry.area / 10000
    data_region.loc[cur_index, "Número de humedales"] = len(tmp)
    del tmp, tmp_geometry

data_region

fig, ax = plt.subplots(figsize=(10, 8))

position = np.arange(len(data_region))
width = 0.35

area_bar = ax.bar(position - width / 2, data_region["Hectáreas (geometría)"], width, label="Área [Ha]", color="steelblue")
ax.set_xlabel("Regiones")
ax.set_ylabel("Área [Ha]")

ax2 = ax.twinx()
number_bar = ax2.bar(position + width / 2, data_region["Número de humedales"], width, label="N° humedales", color="coral")
ax2.set_ylabel("N° humedales")

ax.set_xticks(position)
ax.set_xticklabels(data_region["Región"], rotation=45, ha='right')

lines1, labels1 = ax.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax.legend(lines1 + lines2, labels1 + labels2, loc='upper right')
fig.savefig("humedales_por_region.png", dpi=300, bbox_inches="tight")
plt.tight_layout()

from IPython.display import Markdown as md
total = 17000
text = f""" 
El logro de las **{total_superficie} hectáreas** demuestra una política pública altamente efectiva en su fase de declaración. La **Ley N° 21.202** ha funcionado como un catalizador institucional que trasciende la mera conservación de la biodiversidad, obligando a que cualquier proyecto que pueda afectar a estos ecosistemas ingrese al **Sistema de Evaluación de Impacto Ambiental (SEIA)** y forzando a carteras enfocadas en obras civiles (como el **MOP**) a generar nuevos estándares de infraestructura pública.
"""
md(text)