3.1 Coverage of protected areas and other effective area-based conservation measures

import datetime
now = datetime.datetime.now()
print(f"Última actualización de este notebook: {now.strftime('%d-%m-%Y %H:%M')}")

!pip install pyproj>=3.5.0
!pip install geopandas
!pip install fiona
!pip install gdown
!pip install tabulate
!pip install -U "folium>=0.12" branca matplotlib mapclassify"
!pip install geoalchemy2

import requests
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, MultiPoint, LineString, MultiLineString, Polygon, MultiPolygon, LinearRing
import gdown 
from pathlib import Path
import zipfile
from IPython.display import Markdown, display
from pyproj import CRS, Transformer, Geod
import shapely
import numpy as np
import os
from sqlalchemy import create_engine, text
from geoalchemy2 import Geometry
import geopandas as gpd

SESSION = requests.Session()
geod = Geod(ellps="WGS84")

def geodesic_area_m2(geom, crs=None, target_crs="EPSG:4326"):
    """
    Calcula área geodésica (m²) usando WGS84.
    - Si la geometría viene en CRS proyectado (metros) o distinto a lon/lat,
      la transforma a EPSG:4326 antes del cálculo.
    - `crs` puede ser: "EPSG:32719", pyproj.CRS, o el .crs de un GeoDataFrame.
    """
    if geom is None or geom.is_empty:
        return 0.0

    # Si no sabemos el CRS, asumimos que ya viene en lon/lat
    if crs is None:
        area, _ = geod.geometry_area_perimeter(geom)
        return abs(area)

    src = CRS.from_user_input(crs)
    dst = CRS.from_user_input(target_crs)

    # Si NO es geográfico (lon/lat), o no coincide con el destino, transformamos
    if (not src.is_geographic) or (src != dst):
        transformer = Transformer.from_crs(src, dst, always_xy=True)

        def _tx(x, y, z=None):
            return transformer.transform(x, y)

        # shapely >= 2
        try:
            import shapely
            geom_ll = shapely.ops.transform(_tx, geom)
        except Exception:
            # fallback por si no está el import directo
            from shapely.ops import transform
            geom_ll = transform(_tx, geom)
    else:
        geom_ll = geom

    area, _ = geod.geometry_area_perimeter(geom_ll)
    return abs(area)

def _get_json(url, params, timeout=180):
    r = SESSION.get(url, params=params, timeout=timeout)
    r.raise_for_status()
    js = r.json()
    if "error" in js:
        raise RuntimeError(js["error"])
    return js

def _get_layer_info(layer_url: str) -> dict:
    return _get_json(layer_url, {"f": "pjson"}, timeout=60)

def _get_object_ids(layer_url: str, where="1=1") -> list[int]:
    query_url = layer_url.rstrip("/") + "/query"
    js = _get_json(query_url, {
        "f": "pjson",
        "where": where,
        "returnIdsOnly": "true",
        "returnGeometry": "false",
    }, timeout=120)
    return js.get("objectIds") or []

def _ring_area(coords):
    a = 0.0
    for (x1, y1), (x2, y2) in zip(coords, coords[1:]):
        a += (x1 * y2 - x2 * y1)
    return a / 2.0

def _esri_polygon_to_shape(rings):
    if not rings:
        return None

    norm = []
    for ring in rings:
        if not ring:
            continue
        if ring[0] != ring[-1]:
            ring = ring + [ring[0]]
        norm.append(ring)

    outers, holes = [], []
    for ring in norm:
        try:
            lr = LinearRing(ring)
            if lr.is_ccw:
                holes.append(ring)
            else:
                outers.append(ring)
        except Exception:
            if _ring_area(ring) < 0:
                outers.append(ring)
            else:
                holes.append(ring)

    if not outers and norm:
        areas = [(abs(_ring_area(r)), r) for r in norm]
        areas.sort(reverse=True, key=lambda t: t[0])
        outers = [areas[0][1]]
        holes = [r for _, r in areas[1:]]

    outer_polys = [Polygon(o) for o in outers if len(o) >= 4]
    if not outer_polys:
        return None

    holes_by_outer = {i: [] for i in range(len(outer_polys))}
    for h in holes:
        hp = Polygon(h)
        if hp.is_empty:
            continue
        pt = hp.representative_point()
        for i, op in enumerate(outer_polys):
            if op.contains(pt):
                holes_by_outer[i].append(h)
                break

    polys = []
    for i, o in enumerate(outers):
        try:
            polys.append(Polygon(o, holes=holes_by_outer[i]))
        except Exception:
            polys.append(Polygon(o))

    return polys[0] if len(polys) == 1 else MultiPolygon(polys)

def _esri_geom_to_shape(geom: dict, geom_type: str):
    if geom is None:
        return None

    gt = (geom_type or "").lower()

    if "point" in gt:
        x, y = geom.get("x"), geom.get("y")
        return Point(x, y) if x is not None and y is not None else None

    if "multipoint" in gt:
        pts = geom.get("points") or []
        return MultiPoint([Point(x, y) for x, y in pts]) if pts else None

    if "polyline" in gt:
        paths = geom.get("paths") or []
        lines = [LineString(p) for p in paths if p and len(p) >= 2]
        if not lines:
            return None
        return lines[0] if len(lines) == 1 else MultiLineString(lines)

    if "polygon" in gt:
        rings = geom.get("rings") or []
        return _esri_polygon_to_shape(rings)

    return None

def _query_by_object_ids(layer_url, object_ids, out_sr=4326, out_fields="*", max_allowable_offset=None, timeout=180):
    query_url = layer_url.rstrip("/") + "/query"
    params = {
        "f": "pjson",
        "objectIds": ",".join(map(str, object_ids)),
        "outFields": out_fields,
        "returnGeometry": "true",
        "outSR": out_sr,
        "returnZ": "false",
        "returnM": "false",
    }
    if max_allowable_offset is not None:
        # Generaliza geometría en unidades del SR de salida (en 4326 es grados)
        params["maxAllowableOffset"] = max_allowable_offset

    js = _get_json(query_url, params, timeout=timeout)
    return js.get("features") or []

def _fetch_features_resilient(layer_url, ids, geom_type, out_sr=4326, timeout=180):
    """
    Intenta bajar features por IDs.
    Si falla (500 u otros), parte en 2 recursivamente.
    En caso extremo de un solo ID fallando, reintenta con generalización; si sigue fallando, lo salta.
    """
    if not ids:
        return []

    try:
        return _query_by_object_ids(layer_url, ids, out_sr=out_sr, max_allowable_offset=None, timeout=timeout)
    except Exception:
        # Si el lote es grande, bisección
        if len(ids) > 1:
            mid = len(ids) // 2
            left = _fetch_features_resilient(layer_url, ids[:mid], geom_type, out_sr=out_sr, timeout=timeout)
            right = _fetch_features_resilient(layer_url, ids[mid:], geom_type, out_sr=out_sr, timeout=timeout)
            return left + right

        # Caso: un solo OBJECTID está rompiendo el query
        bad_id = ids[0]
        # Intento 1: generalización suave (en grados, si out_sr=4326)
        for offset in (1e-5, 5e-5, 1e-4, 5e-4, 1e-3):
            try:
                feats = _query_by_object_ids(layer_url, [bad_id], out_sr=out_sr, max_allowable_offset=offset, timeout=timeout)
                return feats
            except Exception:
                continue

        print(f"[WARN] No se pudo descargar geometría para OBJECTID={bad_id} en {layer_url}. Se omite.")
        return []

def fetch_arcgis_layer_gdf(layer_url: str, where="1=1", chunk_size=100, out_sr=4326, timeout=180) -> gpd.GeoDataFrame:
    info = _get_layer_info(layer_url)
    geom_type = info.get("geometryType", "")
    oids = _get_object_ids(layer_url, where=where)

    if not oids:
        return gpd.GeoDataFrame(geometry=[], crs=f"EPSG:{out_sr}")

    rows = []
    for i in range(0, len(oids), chunk_size):
        ids_chunk = oids[i:i + chunk_size]
        feats = _fetch_features_resilient(layer_url, ids_chunk, geom_type, out_sr=out_sr, timeout=timeout)

        for ft in feats:
            attrs = ft.get("attributes") or {}
            geom = ft.get("geometry")
            attrs["geometry"] = _esri_geom_to_shape(geom, geom_type)
            rows.append(attrs)

    gdf = gpd.GeoDataFrame(rows, geometry="geometry", crs=f"EPSG:{out_sr}")
    gdf = gdf[~gdf.geometry.isna()].copy()
    return gdf

def make_engine_from_env(host_override=None):
    host = host_override or os.getenv("POSTGRES_HOST", "localhost")
    port = os.getenv("POSTGRES_PORT", "5432")
    user = os.getenv("POSTGRES_USER", "postgres")
    password = 'super_secret_pass'#os.getenv("POSTGRES_PASSWORD", "")
    db = os.getenv("POSTGRES_DB", "postgres")

    # OJO: en .env no pongas comillas dentro del valor. Debería ser sin ''.
    # Si ya tienes comillas en el valor, aquí igual funcionará en la mayoría de casos,
    # pero lo correcto es quitarlas del archivo.
    password = password.strip("'").strip('"')

    url = f"postgresql+psycopg2://{user}:{password}@{host}:{port}/{db}"
    return create_engine(url, pool_pre_ping=True)

engine = make_engine_from_env(host_override="db")  # desde tu host
schema = "public"

# Asegura CRS
def ensure_crs(gdf: gpd.GeoDataFrame, epsg=4326):
    if gdf.crs is None:
        return gdf.set_crs(epsg)
    return gdf

# -------------------------
# Descargar + unir capas
# -------------------------
gdfs = [fetch_arcgis_layer_gdf(u, chunk_size=100, out_sr=32719) for u in LAYER_AP_URLS]
areas_protegidas_capa_gpd = gpd.GeoDataFrame(gdfs[0], geometry="geometry", crs=gdfs[0].crs)
otras_designaciones_capa_gpd = gpd.GeoDataFrame(gdfs[1], geometry="geometry", crs=gdfs[1].crs)
otras_designaciones_capa_gpd = otras_designaciones_capa_gpd[otras_designaciones_capa_gpd["designacion"]!="Reserva de la Biosfera"]
conservacion_privada_capa_gpd = gpd.GeoDataFrame(gdfs[2], geometry="geometry", crs=gdfs[2].crs)

areas_protegidas_capa_gpd["area_m2"] = areas_protegidas_capa_gpd.geometry.apply(
    lambda g: geodesic_area_m2(g, crs=areas_protegidas_capa_gpd.crs)
)
areas_protegidas_capa_gpd["area_ha"] = areas_protegidas_capa_gpd["area_m2"] / 10_000
areas_protegidas_capa_gpd_area_total_ha = areas_protegidas_capa_gpd["area_ha"].sum()
areas_protegidas_capa_gpd_area_total_km2 = areas_protegidas_capa_gpd_area_total_ha / 100

otras_designaciones_capa_gpd["area_m2"] = otras_designaciones_capa_gpd.geometry.apply(
    lambda g: geodesic_area_m2(g, crs=otras_designaciones_capa_gpd.crs)
)
otras_designaciones_capa_gpd["area_ha"] = otras_designaciones_capa_gpd["area_m2"] / 10_000
otras_designaciones_capa_gpd_area_total_ha = otras_designaciones_capa_gpd["area_ha"].sum()
otras_designaciones_capa_gpd_area_total_km2 = otras_designaciones_capa_gpd_area_total_ha / 100

conservacion_privada_capa_gpd["area_m2"] = conservacion_privada_capa_gpd.geometry.apply(
    lambda g: geodesic_area_m2(g, crs=conservacion_privada_capa_gpd.crs)
)
conservacion_privada_capa_gpd["area_ha"] = conservacion_privada_capa_gpd["area_m2"] / 10_000
conservacion_privada_capa_gpd_area_total_ha = conservacion_privada_capa_gpd["area_ha"].sum()
conservacion_privada_capa_gpd_area_total_km2 = conservacion_privada_capa_gpd_area_total_ha / 100

areas_protegidas_total_gpd = gpd.GeoDataFrame(
    pd.concat([areas_protegidas_capa_gpd,otras_designaciones_capa_gpd, conservacion_privada_capa_gpd], ignore_index=True), 
    geometry="geometry", crs=gdfs[0].crs)

print("CRS (entrada):", areas_protegidas_total_gpd.crs)
print("Filas:", len(areas_protegidas_total_gpd))

areas_protegidas_total_gpd["geometry"] = areas_protegidas_total_gpd.geometry.make_valid()
areas_protegidas_total_gpd_union = areas_protegidas_total_gpd.geometry.union_all()

areas_protegidas_total_gpd["area_m2"] = areas_protegidas_total_gpd.geometry.apply(
    lambda g: geodesic_area_m2(g, crs=areas_protegidas_total_gpd.crs)
)
areas_protegidas_total_gpd["area_ha"] = areas_protegidas_total_gpd["area_m2"] / 10_000
areas_protegidas_total_gpd["area_km2"] = areas_protegidas_total_gpd["area_ha"] / 100
areas_protegidas_total_gpd_area_total_ha = areas_protegidas_total_gpd["area_ha"].sum()
areas_protegidas_total_gpd_area_total_km2 = areas_protegidas_total_gpd_area_total_ha / 100

areas_protegidas_total_gpd = gpd.read_postgis(
    sql="SELECT * FROM public.areas_protegidas_total_gpd_union_sin_validar",
    con=engine,
    geom_col="geometry",    
    crs="EPSG:32719"       
)
areas_protegidas_total_gpd_union = areas_protegidas_total_gpd.geometry.union_all()
#areas_protegidas_total_gpd["area_m2"] = efg_continental_gpd.geometry.apply(
#    lambda g: geodesic_area_m2(g, crs=efg_continental_gpd.crs)
#)
areas_protegidas_total_gpd_area_total_m2 = geodesic_area_m2(areas_protegidas_total_gpd_union, crs=areas_protegidas_total_gpd.crs)
areas_protegidas_total_gpd_area_total_ha = areas_protegidas_total_gpd_area_total_m2 / 10_000
areas_protegidas_total_gpd_area_total_km2 = areas_protegidas_total_gpd_area_total_ha / 100

ap_ha = round(areas_protegidas_total_gpd_area_total_ha, 0)
ap_ha_formateado = f"{int(ap_ha):,}".replace(",", ".")
terrestre_ha=75662600
marino_ha=340229000
totas_superficie_chile_ha=terrestre_ha+marino_ha
conservacion_porcentaje=round((ap_ha/totas_superficie_chile_ha)*100,2)

from IPython.display import Markdown as md
text = f""" 
Para el actual ciclo de reporte, la cobertura total del sistema de áreas protegidas e iniciativas de conservación reconocidas a nivel nacional alcanza las {ap_ha_formateado} hectáreas. Este volumen de conservación implica que un {conservacion_porcentaje}% de la superficie jurisdiccional total del país (combinando tierra y mar) se encuentra resguardado bajo alguna figura formal de protección.
"""
md(text)

# Cargar superficie terrestre para cruce con áreas protegidas
from pathlib import Path
import zipfile
chile_url='https://www.bcn.cl/obtienearchivo?id=repositorio/10221/10398/2/Regiones.zip'
zip_path = Path("Regiones.zip")
referer = "https://www.bcn.cl/siit/mapas_vectoriales"
headers = {
    "User-Agent": "Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/120 Safari/537.36",
    "Referer": referer,
    "Accept": "*/*",
    "Accept-Language": "es-CL,es;q=0.9,en;q=0.8",
    "Connection": "keep-alive",
}
s = requests.Session()
# 1) “Visita” la página para obtener cookies si las piden
s.get(referer, headers=headers, timeout=60)
# 2) Descarga el ZIP
with s.get(chile_url, headers=headers, stream=True, timeout=120, allow_redirects=True) as r:
    r.raise_for_status()
    with open(zip_path, "wb") as f:
        for chunk in r.iter_content(chunk_size=1024 * 1024):
            if chunk:
                f.write(chunk)

print("Descargado:", zip_path)
with zipfile.ZipFile(zip_path, "r") as z:
    names = z.namelist()
    shp_files = [n for n in names if n.lower().endswith(".shp")]
    print("SHP encontrados:", shp_files)
shp_inside = shp_files[0]  # elige el correcto si hay más de uno
regiones_gpd = gpd.read_file(f"zip://{zip_path}!{shp_inside}")
regiones_gpd = regiones_gpd.to_crs(epsg=32719)
print("CRS:", regiones_gpd.crs)
regiones_gpd["area_m2"] = regiones_gpd.geometry.apply(
    lambda g: geodesic_area_m2(g, crs=regiones_gpd.crs)
)
regiones_gpd["area_ha"] = regiones_gpd["area_m2"] / 10_000
regiones_gpd["area_km2"] = regiones_gpd["area_ha"] / 100
regiones_gpd_area_total_ha = regiones_gpd["area_ha"].sum()
regiones_gpd_area_total_km2 = regiones_gpd_area_total_ha / 100

regiones_gpd = gpd.read_postgis(
    sql="SELECT * FROM public.regiones",
    con=engine,
    geom_col="geometry",     
    crs="EPSG:32719"      
)
regiones_gpd_area_total_ha = regiones_gpd["area_ha"].sum()
regiones_gpd_area_total_km2 = regiones_gpd_area_total_ha / 100

#Cobertura por Ecosistemas Terrestres
areas_protegidas_terrestres_gpd = gpd.overlay(areas_protegidas_total_gpd, regiones_gpd, how="intersection", keep_geom_type=False)
areas_protegidas_terrestres_gpd_union = areas_protegidas_terrestres_gpd.geometry.union_all()
areas_protegidas_terrestres_gpd_area_total_m2 = geodesic_area_m2(
    areas_protegidas_terrestres_gpd_union,
    crs=areas_protegidas_terrestres_gpd.crs
)
areas_protegidas_terrestres_gpd_area_total_ha = areas_protegidas_terrestres_gpd_area_total_m2 / 10_000
areas_protegidas_terrestres_gpd_area_total_km2 = areas_protegidas_terrestres_gpd_area_total_ha / 100
porcentaje_cobertura_terrestre = round((areas_protegidas_terrestres_gpd_area_total_km2 / regiones_gpd_area_total_km2) * 100,2)

areas_protegidas_terrestres_gpd_area_total_ha = round(areas_protegidas_terrestres_gpd_area_total_ha, 0)
areas_protegidas_terrestres_gpd_area_total_ha_formateado = f"{int(areas_protegidas_terrestres_gpd_area_total_ha):,}".replace(",", ".")

text = f""" 
Ámbito Continental (Terrestre y Dulceacuícola): Se reportan {areas_protegidas_terrestres_gpd_area_total_ha_formateado} hectáreas protegidas, lo que equivale a un {porcentaje_cobertura_terrestre}% de la superficie continental oficial.
"""
#md(text)

# Cargar superficie marina para cruce con áreas protegidas
ecosistemas_marinos_gpd = gpd.read_postgis(
    sql="SELECT * FROM public.ecosistemas_marinos",
    con=engine,
    geom_col="geometry",     
    crs="EPSG:32719"      
)

#Cobertura por Ecosistemas Terrestres
import shapely
areas_protegidas_total_gpd = shapely.make_valid(areas_protegidas_total_gpd_union)

from tqdm.auto import tqdm
new_geoms = []
for g in tqdm(ecosistemas_marinos_gpd.geometry, total=len(ecosistemas_marinos_gpd), desc="Intersección"):
    new_geoms.append(g.intersection(areas_protegidas_total_gpd_union))
ecosistemas_marinos_gpd["geometry"] = new_geoms

ecosistemas_marinos_gpd_area_total_m2 = ecosistemas_marinos_gpd.geometry.apply(
    lambda g: geodesic_area_m2(g, crs=ecosistemas_marinos_gpd.crs)
).sum()

areas_protegidas_marinas_gpd_area_total_ha = ecosistemas_marinos_gpd_area_total_m2 / 10_000
porcentaje_cobertura_marina = round((areas_protegidas_marinas_gpd_area_total_ha / marino_ha) * 100,2)