3  RENABAP

Análisis de la exposición poblacional a peligros de inundación en el Partido de La Plata

3.1 Resumen ejecutivo

Los datos procesados, incluyendo las estimaciones de exposición por barrio popular para los tres métodos de estimación, desglosados por nivel de peligro, están disponibles para descarga en formato CSV. Haga clic aquí para descargar los datos por barrio. Los datos agregados por cuenca también están disponibles, incluyendo una columna “eje” que permite resumir los datos al nivel de eje de cuenca. Haga clic aquí para descargar los datos por cuenca.

3.2 Objetivos

El objetivo principal de este análisis es calcular con mayor precisión la exposición poblacional en la región del Partido de La Plata, comparando el enfoque actual de interpolación por área de datos del RENABAP con enfoques de mapeo dasimétrico de mayor resolución utilizando datos de la Capa Global de Asentamientos Humanos (GHSL) y el conjunto de datos de edificios abiertos de Google-Microsoft OpenStreetMap.

Para lograr este objetivo, necesitamos:

1. Explicar y comparar metodologías: Desarrollar una comprensión clara de las diferencias entre interpolación por área y mapeo dasimétrico
2. Evaluar fuentes de datos: Analizar las ventajas y limitaciones de cada conjunto de datos utilizado
3. Crear estimaciones robustas: Generar un rango de posibles exposiciones poblacionales para informar la toma de decisiones
4. Priorizar recursos: Identificar áreas donde se requiera recopilación de datos más precisa

3.3 Contexto

import matplotlib.pyplot as plt
import contextily as ctx


from io import BytesIO, StringIO
from owslib.wfs import WebFeatureService
from shapely.geometry import box
import geopandas as gpd
import requests
import pandas as pd
import os

import boto3
import duckdb


import numpy as np
import s2sphere
from botocore.config import Config
from rasterstats import zonal_stats
import rasterstats
from rasterio.features import rasterize
import numpy.ma as ma
import itables
from itables import show
import json
import requests
import matplotlib.cm as cm
from IPython.display import HTML, display

import seaborn as sns
import rioxarray

# =============================================================================
# ITABLES SPANISH CONFIGURATION
# =============================================================================

# Configure Argentine Spanish for itables
try:
    spanish_url = "https://cdn.datatables.net/plug-ins/2.3.3/i18n/es-AR.json"
    response = requests.get(spanish_url)
    response.raise_for_status()
    spanish_config = response.json()
    itables.options.language = spanish_config
except Exception:
    pass  # Fall back to English if configuration fails

# Configure smaller font size for all itables
css = """
.dt-container {
  font-size: small;
}
"""
display(HTML(f"<style>{css}</style>"))

# Helper function to round numeric columns for display
def round_numeric_columns(df, decimals=0):
    """Round all numeric columns in a DataFrame to specified decimal places."""
    df_display = df.copy()
    numeric_columns = df_display.select_dtypes(include=[np.number]).columns
    df_display[numeric_columns] = df_display[numeric_columns].round(decimals)
    return df_display

# =============================================================================
# CONSTANTS AND CONFIGURATION
# =============================================================================

# Coordinate Reference Systems
USE_CRS = "EPSG:5349"  # POSGAR 2007 / Argentina 4
WEB_MERCATOR_CRS = "EPSG:3857"  # Web Mercator for visualization
WGS84_CRS = "EPSG:4326"  # WGS84 for API calls

# File paths
BASE_PATH = "/home/nissim/Documents/dev/fulbright/ciut-riesgo"
DATA_PATH = f"{BASE_PATH}/notebooks/data"
PELIGRO_PATH = f"{DATA_PATH}/la_plata_pelig_2023_datos_originales.geojson"
PARTIDOS_PATH = f"{DATA_PATH}/pba_partidos.geojson"
CUENCAS_PATH = f"{BASE_PATH}/notebooks/cuencas_buenos_aires.geojson"
BUILDINGS_PATH = f"{BASE_PATH}/notebooks/buildings_filtered.parquet"
GHSL_PATH = "/home/nissim/Downloads/spatial/GHS_POP_E2025_GLOBE_R2023A_54009_100_V1_0_R14_C13/GHS_POP_E2025_GLOBE_R2023A_54009_100_V1_0_R14_C13.tif"

# Data URLs
RENABAP_URL = (
    "https://www.argentina.gob.ar/sites/default/files/renabap-2023-12-06.geojson"
)
PARTIDOS_WFS_URL = "https://geo.arba.gov.ar/geoserver/idera/wfs"
CUENCAS_API_URL = "https://services1.arcgis.com/atxllciEI8CHWvwW/ArcGIS/rest/services/Cuencas_BuenosAires_2023/FeatureServer/0/query"

# Data processing constants
HAZARD_LEVELS = ["baja", "media", "alta"]
METHOD_NAMES = ["edificios", "ghsl", "areal"]
EXPOSURE_COLUMNS = [
    "fam_exp_edificios",
    "fam_exp_ghsl",
    "fam_exp_areal",
]
NON_HAZARD_VALUE = "none"
NODATA_VALUE = -200

# Column mappings and renaming
COLUMN_MAPPINGS = {
    "buildings_to_edificios": {"fam_expuestas_buildings": "fam_expuestas_edificios"},
    "method_cleanup_prefix": "fam_expuestas_",
}

# Basic visualization settings (only for repeated values)
DEFAULT_FIGSIZE = (12, 10)
MAP_PADDING = 500
PLASMA_CMAP = plt.cm.plasma

# Color schemes for visualization
PELIGROSIDAD_COLORS = {
    "alta": PLASMA_CMAP(0.8),
    "media": PLASMA_CMAP(0.5),
    "baja": PLASMA_CMAP(0.2),
}

METHOD_COLORS = {
    "fam_exp_areal": PLASMA_CMAP(0.8),
    "fam_exp_ghsl": PLASMA_CMAP(0.5),
    "fam_exp_edificios": PLASMA_CMAP(0.2),
}

# Eje mapping for watershed analysis
EJE_MAPPING = {
    "noreste": ["Area de Bañados", "Cuenca Arroyo Rodriguez-Don Carlos"],
    "noroeste": ["Cuenca Arroyo Martín-Carnaval", "Cuenca Arroyo Pereyra"],
    "central": ["Cuenca Arroyo del Gato"],
    "sudoeste": ["Cuenca A° Maldonado", "Cuenca Río Samborombón"],
    "sudeste": ["Cuenca Arroyo El Pescado"],
}


def setup_base_map(figsize=None, bounds=None, padding_x=None, padding_y=None):
    """Create figure and set up basic map boundaries with padding."""
    if figsize is None:
        figsize = DEFAULT_FIGSIZE
    if padding_x is None:
        padding_x = MAP_PADDING
    if padding_y is None:
        padding_y = MAP_PADDING

    if bounds is None:
        bounds = renabap_pba_intersect.total_bounds

    # Convert bounds to Web Mercator for basemap compatibility
    if bounds is not None:
        # Create a temporary GeoDataFrame with the bounds to reproject
        temp_bounds = gpd.GeoDataFrame(
            geometry=[box(bounds[0], bounds[1], bounds[2], bounds[3])], crs=USE_CRS
        )
        bounds_3857 = temp_bounds.to_crs(WEB_MERCATOR_CRS).total_bounds
    else:
        bounds_3857 = bounds

    fig, ax = plt.subplots(figsize=figsize)
    ax.set_xlim(bounds_3857[0] - padding_x, bounds_3857[2] + padding_x)
    ax.set_ylim(bounds_3857[1] - padding_y, bounds_3857[3] + padding_y)
    return fig, ax


def add_basemap(ax, zoom=13):
    """Add CartoDB basemap to the axes."""

    ctx.add_basemap(
        ax,
        source=ctx.providers.CartoDB.PositronNoLabels,
        zorder=0,
        zoom=zoom,
    )

    return ax


def add_north_arrow(ax, x=0.95, y=0.05, arrow_length=0.04):
    """Add a north arrow to the map."""
    # Add north arrow, https://stackoverflow.com/a/58110049/604456
    ax.annotate(
        "N",
        xy=(x, y),
        xytext=(x, y - arrow_length),
        arrowprops=dict(facecolor="black", width=3, headwidth=10),
        ha="center",
        va="center",
        fontsize=14,
        xycoords=ax.transAxes,
    )


def add_la_plata_outline(ax):
    """Add the outline of Partido de La Plata to a map."""
    la_plata_3857 = la_plata.to_crs(WEB_MERCATOR_CRS)
    la_plata_3857.plot(
        ax=ax,
        facecolor="none",
        edgecolor="black",
        linewidth=0.5,
        linestyle="--",
        legend=False,
        zorder=5,
    )


def create_consistent_map(title, bounds=None):
    """Create a map with consistent styling and basemap."""
    fig, ax = setup_base_map(bounds=bounds)

    add_basemap(ax)

    add_north_arrow(ax)

    add_la_plata_outline(ax)

    ax.set_title(title, fontsize=16, fontweight="bold", pad=20)

    ax.set_axis_off()

    return fig, ax


def wfs_to_gdf(
    wfs_url: str, layer_name: str, srs: str = "EPSG:4326"
) -> gpd.GeoDataFrame:
    """
    Descarga una capa WFS y la devuelve como GeoDataFrame.

    Args:
        wfs_url (str): URL del servicio WFS.
        layer_name (str): Nombre de la capa (typename).
        srs (str): Código EPSG del sistema de referencia de coordenadas.

    Returns:
        gpd.GeoDataFrame: Capa descargada como GeoDataFrame.
    """
    wfs = WebFeatureService(url=wfs_url, version="2.0.0")
    response = wfs.getfeature(typename=layer_name, srsname=srs)
    gdf = gpd.read_file(BytesIO(response.read()))
    return gdf


def create_exposure_tidy_data(
    data,
    id_column,
    peligrosidad_column,
    method_suffix,
    exposure_values,
    exclude_zero=True,
):
    """
    Create tidy exposure dataset in a standardized format.

    Args:
        data: DataFrame containing the base data
        id_column: Column name for the identifier (e.g., 'id_renabap', 'Cuenca', 'eje')
        peligrosidad_column: Column name for hazard level
        method_suffix: Suffix for the exposure column (e.g., 'areal', 'ghsl', 'edificios')
        exposure_values: Series or array of exposure values matching data rows
        exclude_zero: Whether to exclude zero exposure values

    Returns:
        pd.DataFrame: Tidy format dataframe with id, peligrosidad, and exposure columns
    """
    tidy_data = []
    for idx, (_, row) in enumerate(data.iterrows()):
        exposure_value = (
            exposure_values.iloc[idx]
            if hasattr(exposure_values, "iloc")
            else exposure_values[idx]
        )

        if exclude_zero and exposure_value <= 0:
            continue

        tidy_data.append(
            {
                id_column: row[id_column],
                "peligrosidad": row[peligrosidad_column],
                f"fam_expuestas_{method_suffix}": exposure_value,
            }
        )

    return pd.DataFrame(tidy_data)


def create_wide_exposure_dataframe(
    areal_data, ghsl_data, buildings_data, id_columns, exclude_hazard_value="none"
):
    """
    Create wide format exposure dataframe by merging tidy datasets.

    Args:
        areal_data: Tidy dataframe with areal interpolation results
        ghsl_data: Tidy dataframe with GHSL dasymetric results
        buildings_data: Tidy dataframe with buildings dasymetric results
        id_columns: List of columns to merge on (e.g., ['id_renabap', 'peligrosidad'])
        exclude_hazard_value: Hazard value to exclude from results

    Returns:
        pd.DataFrame: Wide format dataframe with all exposure methods
    """
    # Filter out non-hazard values
    areal_filtered = areal_data[areal_data["peligrosidad"] != exclude_hazard_value]
    ghsl_filtered = ghsl_data[ghsl_data["peligrosidad"] != exclude_hazard_value]
    buildings_filtered = buildings_data[
        buildings_data["peligrosidad"] != exclude_hazard_value
    ]

    # Apply column mapping for buildings if needed
    if "fam_expuestas_buildings" in buildings_filtered.columns:
        buildings_filtered = buildings_filtered.rename(
            columns={"fam_expuestas_buildings": "fam_expuestas_edificios"}
        )

    # Merge all datasets
    wide_data = areal_filtered.merge(ghsl_filtered, on=id_columns, how="outer").merge(
        buildings_filtered, on=id_columns, how="outer"
    )

    # Fill NaN values with 0
    wide_data = wide_data.fillna(0)
    
    # Rename columns to shorter format
    column_mapping = {
        "fam_expuestas_areal": "fam_exp_areal",
        "fam_expuestas_ghsl": "fam_exp_ghsl", 
        "fam_expuestas_edificios": "fam_exp_edificios"
    }
    wide_data = wide_data.rename(columns=column_mapping)

    return wide_data


response = requests.get(RENABAP_URL)
renabap = gpd.read_file(StringIO(response.text))
renabap_pba = renabap[renabap["provincia"] == "Buenos Aires"]
renabap_pba = renabap_pba.to_crs(USE_CRS)

peligro = gpd.read_file(PELIGRO_PATH)
peligro = peligro.to_crs(USE_CRS)

peligro_bounds = peligro.total_bounds
peligro_bbox = box(*peligro_bounds)

if os.path.exists(PARTIDOS_PATH):
    partidos = gpd.read_file(PARTIDOS_PATH)
else:
    partidos = wfs_to_gdf(
        wfs_url=PARTIDOS_WFS_URL,
        layer_name="idera:Departamento",
        srs="EPSG:5347",
    )

    partidos.to_file(PARTIDOS_PATH, driver="GeoJSON")

partidos = partidos.to_crs(USE_CRS)
la_plata = partidos[partidos["fna"] == "Partido de La Plata"]

# Obtener la geometría principal
main_geom = la_plata.geometry.iloc[0]

# Si es un MultiPolygon, mantener solo el polígono más grande (el partido principal)
# Esto elimina la pequeña isla que aparece en los datos
if main_geom.geom_type == "MultiPolygon":
    # Obtener todos los polígonos y mantener el que tenga mayor área
    largest_polygon = max(main_geom.geoms, key=lambda p: p.area)
    la_plata = la_plata.copy()  # Create a copy to avoid SettingWithCopyWarning
    la_plata.loc[la_plata.index[0], "geometry"] = largest_polygon

la_plata_bbox = la_plata.geometry.iloc[0]

renabap_pba_intersect = renabap_pba[
    renabap_pba.geometry.intersects(la_plata_bbox)
].copy()


if os.path.exists(CUENCAS_PATH):
    cuencas = gpd.read_file(CUENCAS_PATH)
else:
    params = {"where": "1=1", "outFields": "*", "f": "geojson"}

    cuencas_response = requests.get(CUENCAS_API_URL, params=params)
    with open(CUENCAS_PATH, "w", encoding="utf-8") as f:
        f.write(cuencas_response.text)

    cuencas = gpd.read_file(StringIO(cuencas_response.text))

cuencas = cuencas.to_crs(USE_CRS)
cuencas = cuencas.clip(la_plata)

# Map watershed names to axes based on the EJE_MAPPING
cuencas["eje"] = (
    cuencas["Cuenca"]
    .map(
        {
            cuenca: eje
            for eje, cuencas_list in EJE_MAPPING.items()
            for cuenca in cuencas_list
        }
    )
    .fillna("otro")
)

# Calculate total area of RENABAP settlements in hectares (POSGAR projection is in meters)
renabap_total_area_ha = (
    renabap_pba_intersect.geometry.area.sum() / 10000
)  # Convert m² to hectares
la_plata_area_ha = la_plata.geometry.iloc[0].area / 10000
percentage_coverage = (renabap_total_area_ha / la_plata_area_ha) * 100

# Get common bounds for all maps
common_bounds = la_plata.total_bounds

# Intersect settlements with hazard zones
settlement_hazard = gpd.overlay(renabap_pba_intersect, peligro, how="intersection")

settle_hazard_cuencas = gpd.overlay(
    settlement_hazard, cuencas, how="intersection", keep_geom_type=True
)

Hay un total de 166 barrios populares en el Partido de La Plata, que representan 33888 familias. Estos barrios ocupan 1760 hectáreas del Partido de La Plata, o 2.0 por ciento del partido.

3.4 Fuentes de datos

3.4.1 RENABAP

El Registro Nacional de Barrios Populares (RENABAP) es producido por la Subsecretaría de Integración Socio Urbana y proporciona información sobre asentamientos informales en Argentina, incluyendo estimaciones de población y delimitaciones geográficas de estos barrios. Más información sobre el RENABAP está disponible en el Observatorio de Barrios Populares. Los datos fueron obtenidos a través del Mapa de Barrios Populares y están disponibles para descarga como GeoJSON.

3.4.2 GHSL

La Capa Global de Asentamientos Humanos (Global Human Settlement Layer) (Schiavina et al. 2023) es un conjunto de datos de resolución de 100 metros que proporciona estimaciones de población multitemporales (1975-2030) derivadas de datos censales y administrativos, informadas por la distribución y clasificación de áreas construidas. El GHSL ya tiene un uso científico establecido para mapear la exposición poblacional a peligros de inundación (Tellman et al. 2021). Sin embargo, esta fuente presenta limitaciones importantes: estudios sobre modelado de riesgo de inundación con conjuntos de datos globales han demostrado que evaluar la exposición a esta escala de resolución puede llevar a sobreestimaciones de la exposición poblacional en zonas de peligro de inundación en comparación con datos de mayor resolución (Smith et al. 2019).

# Load GHSL data with dask chunking for memory efficiency
ghsl = rioxarray.open_rasterio(
    GHSL_PATH,
    chunks={"x": 1024, "y": 1024},  # Adjust chunk size based on your memory constraints
)

ghsl = ghsl.rio.reproject(dst_crs=USE_CRS)


ghsl_clipped = ghsl.rio.clip(
    [la_plata.geometry.iloc[0]],
    from_disk=True,  # Process from disk to avoid loading entire dataset into memory
)

3.4.3 Peligro de inundación

Los datos de peligro de inundación utilizados en este análisis fueron desarrollados por la Facultad de Ingeniería de la Universidad Nacional de La Plata como parte del Plan de Reducción del Riesgo por Inundaciones en la Región de La Plata (Romanazzi et al. 2019). Estos datos fueron generados mediante la aplicación del modelo hidrológico-hidráulico bidimensional FLO-2D, que permitió simular la dinámica de inundación de todas las cuencas del partido de La Plata para distintos escenarios de eventos pluviométricos extremos. El modelo calcula las principales variables hidráulicas (altura del agua, velocidad y caudal) a lo largo del tiempo, y a partir de estos resultados se generaron los mapas de peligrosidad que combinan el efecto de la profundidad con la velocidad de la corriente, ofreciendo un indicador más completo que los mapas tradicionales de máximas profundidades.

3.4.4 Google-Microsoft-OSM Open Buildings

Los datos de Google-Microsoft-OSM Open Buildings - combined by VIDA (VIDA 2023) representan una forma más precisa de evaluar dónde se ubican los asentamientos humanos. Este conjunto de datos combina Google’s V3 Open Buildings, Microsoft’s GlobalMLFootprints, y OpenStreetMap building footprints, conteniendo más de 2.7 mil millones de huellas de edificios. Estos datos han sido exitosamente aplicados a evaluaciones de riesgo de inundación por empresas globales de riesgo financiero como ICE, demostrando su utilidad para mapear la exposición climática a nivel de huella de edificio individual. Sin embargo, en ausencia de información sobre si los edificios son residenciales o tienen otros usos, y sin datos sobre el número total de unidades en el edificio y habitantes por edificio, solo podemos obtener estimaciones proporcionales aproximadas de dónde se ubican las personas, sin tener una comprensión precisa de quién vive realmente allí y cuántas personas.

def fetch_buildings(geodataframe, temp_file="buildings_filtered.parquet"):
    """Fetch building data for a given GeoDataFrame region"""

    # Get S2 cell and bounds
    center = geodataframe.to_crs(WEB_MERCATOR_CRS).union_all().centroid
    center_wgs84 = (
        gpd.GeoDataFrame(geometry=[center], crs=WEB_MERCATOR_CRS)
        .to_crs(WGS84_CRS)
        .geometry.iloc[0]
    )
    cell = s2sphere.CellId.from_lat_lng(
        s2sphere.LatLng.from_degrees(center_wgs84.y, center_wgs84.x)
    ).parent(10)
    bounds = geodataframe.to_crs(WGS84_CRS).total_bounds

    # Find matching S2 partition
    s3 = boto3.client(
        "s3",
        endpoint_url="https://data.source.coop",
        aws_access_key_id="",
        aws_secret_access_key="",
        config=Config(s3={"addressing_style": "path"}),
    )

    partitions = {
        obj["Key"].split("/")[-1].replace(".parquet", "")
        for obj in s3.list_objects_v2(
            Bucket="vida",
            Prefix="google-microsoft-osm-open-buildings/geoparquet/by_country_s2/country_iso=ARG/",
        ).get("Contents", [])
    }

    parent_id = next(
        str(cell.parent(level).id())
        for level in range(10, 0, -1)
        if str(cell.parent(level).id()) in partitions
    )

    # Setup DuckDB and query
    con = duckdb.connect()
    for cmd in [
        "INSTALL spatial",
        "LOAD spatial",
        "INSTALL httpfs",
        "LOAD httpfs",
        "SET s3_region='us-east-1'",
        "SET s3_endpoint='data.source.coop'",
        "SET s3_use_ssl=true",
        "SET s3_url_style='path'",
    ]:
        con.execute(cmd)

    # Export and read back
    query = f"""
    COPY (SELECT * FROM 's3://vida/google-microsoft-osm-open-buildings/geoparquet/by_country_s2/country_iso=ARG/{parent_id}.parquet'
          WHERE bbox.xmax >= {bounds[0]} AND bbox.xmin <= {bounds[2]} AND
                bbox.ymax >= {bounds[1]} AND bbox.ymin <= {bounds[3]}
    ) TO '{temp_file}' (FORMAT PARQUET);
    """

    con.execute(query)
    df = pd.read_parquet(temp_file)
    df["geometry"] = gpd.GeoSeries.from_wkb(df["geometry"])

    return gpd.GeoDataFrame(df, geometry="geometry", crs=WGS84_CRS)


if os.path.exists(BUILDINGS_PATH):
    buildings = gpd.read_parquet(BUILDINGS_PATH)
else:
    buildings = fetch_buildings(renabap_pba_intersect)


buildings_proj = buildings.to_crs(USE_CRS)

buildings_proj = buildings_proj.clip(la_plata)

fig1, ax1 = create_consistent_map("Asentamientos RENABAP en La Plata", common_bounds)

renabap_pba_intersect_3857 = renabap_pba_intersect.to_crs(WEB_MERCATOR_CRS)

renabap_pba_intersect_3857.plot(
    ax=ax1, facecolor="none", edgecolor="black", linewidth=0.5, legend=False, zorder=10
)

plt.tight_layout()
plt.show()

peligro_clipped = gpd.clip(peligro, la_plata)

peligro_clipped_3857 = peligro_clipped.to_crs(WEB_MERCATOR_CRS)

# Reorder the categories so they map correctly to plasma colormap
peligro_clipped_3857["PELIGROSID_ordered"] = pd.Categorical(
    peligro_clipped_3857["PELIGROSID"],
    categories=["baja", "media", "alta"],
    ordered=True,
)


fig2, ax2 = create_consistent_map("Zonas de Peligro en La Plata", common_bounds)


peligro_clipped_3857.plot(
    ax=ax2,
    column="PELIGROSID_ordered",
    cmap="plasma",
    alpha=0.75,
    legend=True,
    zorder=5,
)

plt.tight_layout()
plt.show()

fig3, ax3 = create_consistent_map("Datos de población GHSL", common_bounds)

ghsl_masked = ma.masked_where(ghsl_clipped.values[0] == 0, ghsl_clipped.values[0])

ghsl_valid = (ghsl_clipped.values[0] != NODATA_VALUE) & (ghsl_clipped.values[0] != 0)
ghsl_valid_data = ghsl_clipped.values[0][ghsl_valid]

plasma_cmap = PLASMA_CMAP

ghsl_clipped_3857 = ghsl_clipped.rio.reproject(WEB_MERCATOR_CRS)

# Mask out zeros AND nodata values
ghsl_masked_3857 = ma.masked_where(
    (ghsl_clipped_3857.values[0] == 0) | (ghsl_clipped_3857.values[0] == NODATA_VALUE),
    ghsl_clipped_3857.values[0],
)

im = ax3.imshow(
    ghsl_masked_3857,
    extent=[
        ghsl_clipped_3857.x.min(),
        ghsl_clipped_3857.x.max(),
        ghsl_clipped_3857.y.min(),
        ghsl_clipped_3857.y.max(),
    ],
    cmap=plasma_cmap,
    alpha=0.75,
)

plt.tight_layout()
plt.show()

fig4, ax4 = create_consistent_map("Huellas de edificios", common_bounds)

buildings_3857 = buildings_proj.to_crs(WEB_MERCATOR_CRS)

buildings_3857.plot(ax=ax4, facecolor="grey", edgecolor="none", alpha=0.7)

plt.tight_layout()
plt.show()

(a) Asentamientos RENABAP en La Plata

(b) Zonas de Peligro en La Plata

(c) Datos de población GHSL

(d) Huellas de edificios

Figure 3.1: Fuentes de datos para análisis de exposición

3.5 Metodología y procesamiento

3.5.1 Interpolación por area

La interpolación areal es un método simple en el que las variables de los datos fuente se ponderan según la superposición entre polígonos fuente y objetivo, luego se reagregan para ajustarse a las geometrías de los polígonos objetivo. En nuestro análisis, esto significa distribuir proporcionalmente la población de cada barrio popular según el área de intersección con diferentes niveles de peligro de inundación. El analasis original de la exposición poblacional a peligros de inundación en la región del Partido de La Plata se realizó utilizando este método.

if renabap_pba_intersect.crs != peligro.crs:
    peligro = peligro.to_crs(renabap_pba_intersect.crs)

hazard_levels = peligro["PELIGROSID"].unique()

renabap_with_porciones = renabap_pba_intersect.copy()
for level in hazard_levels:
    renabap_with_porciones[f"porcion_{level}"] = 0.0

renabap_with_porciones["total_area"] = renabap_with_porciones.geometry.area

for idx, barrio in renabap_with_porciones.iterrows():
    barrio_geom = barrio.geometry
    barrio_total_area = barrio_geom.area

    if barrio_total_area == 0:
        continue

    for level in hazard_levels:
        hazard_subset = peligro[peligro["PELIGROSID"] == level]

        if hazard_subset.empty:
            continue

        intersection_area = 0
        for _, hazard_row in hazard_subset.iterrows():
            try:
                intersection = barrio_geom.intersection(hazard_row.geometry)
                if not intersection.is_empty:
                    intersection_area += intersection.area
            except Exception as e:
                print(
                    f"Error calculating intersection for {barrio.get('nombre_barrio', idx)}: {e}"
                )
                continue

        proportion = (
            intersection_area / barrio_total_area if barrio_total_area > 0 else 0
        )
        renabap_with_porciones.at[idx, f"porcion_{level}"] = proportion

# Create barrio tidy format using consolidated function
barrio_areal_rows = []
for idx, row in renabap_with_porciones.iterrows():
    for level in hazard_levels:
        familias_expuestas = row[f"porcion_{level}"] * row["familias_aproximadas"]
        if familias_expuestas > 0:
            barrio_areal_rows.append(
                {
                    "id_renabap": row["id_renabap"],
                    "PELIGROSID": level,
                    "fam_expuestas": familias_expuestas,
                }
            )

barrio_areal_temp = pd.DataFrame(barrio_areal_rows)
barrio_areal_tidy = create_exposure_tidy_data(
    data=barrio_areal_temp,
    id_column="id_renabap",
    peligrosidad_column="PELIGROSID",
    method_suffix="areal",
    exposure_values=barrio_areal_temp["fam_expuestas"],
    exclude_zero=True,
)

# 2. CUENCA AREAL EXPOSURE - aggregate from barrio level, avoiding double counting
# Get the cuenca for each settlement - but handle settlements that cross cuenca boundaries
settlement_cuenca_mapping = settle_hazard_cuencas[
    ["id_renabap", "Cuenca"]
].drop_duplicates()

# Check if any settlements appear in multiple cuencas
settlement_counts = settlement_cuenca_mapping["id_renabap"].value_counts()
multi_cuenca_settlements = settlement_counts[settlement_counts > 1].index

if len(multi_cuenca_settlements) > 0:
    # For settlements in multiple cuencas, assign to the cuenca with largest intersection
    # For now, just take the first occurrence
    settlement_cuenca_mapping = settlement_cuenca_mapping.drop_duplicates(
        subset=["id_renabap"], keep="first"
    )

cuenca_areal_tidy = barrio_areal_tidy[
    barrio_areal_tidy["peligrosidad"] != NON_HAZARD_VALUE
].merge(settlement_cuenca_mapping, on="id_renabap", how="left")
cuenca_areal_tidy = (
    cuenca_areal_tidy.groupby(["Cuenca", "peligrosidad"])["fam_expuestas_areal"]
    .sum()
    .reset_index()
)

# 3. EJE AREAL EXPOSURE - same fix
cuenca_eje_mapping = settle_hazard_cuencas[["Cuenca", "eje"]].drop_duplicates()

eje_areal_tidy = cuenca_areal_tidy.merge(cuenca_eje_mapping, on="Cuenca")
eje_areal_tidy = (
    eje_areal_tidy.groupby(["eje", "peligrosidad"])["fam_expuestas_areal"]
    .sum()
    .reset_index()
)

3.5.2 Mapeo dasymetrico

El mapeo dasimétrico reorganiza datos cartográficos de una unidad de recolección en áreas más precisas, modificando los límites originales usando datos de apoyo relacionados. Por ejemplo, un atributo de población organizado por tracto censal se vuelve más significativo cuando se eliminan áreas donde es razonable inferir que la gente no vive (cuerpos de agua, terrenos vacíos). En nuestro caso, utilizamos datos GHSL y huellas de edificios como información auxiliar para mejorar la precisión de las estimaciones de distribución poblacional.

### BUILDINGS

# Get ALL buildings per settlement (not just hazard-intersected ones)
buildings_settlement = gpd.overlay(
    buildings_proj, renabap_pba_intersect, how="intersection"
)
total_buildings_per_settlement = (
    buildings_settlement.groupby("id_renabap")
    .size()
    .reset_index(name="total_buildings")
)

# Get buildings intersected with hazard zones
buildings_hazard = gpd.overlay(
    buildings_proj, settle_hazard_cuencas, how="intersection"
)

# 1. Buildings per barrio-hazard (including non-hazard areas)
buildings_barrio_hazard = (
    buildings_hazard.groupby(["id_renabap", "PELIGROSID"])
    .size()
    .reset_index(name="buildings_count")
)

# Calculate ratios using TOTAL buildings per settlement (not just hazard buildings)
barrio_ratios = buildings_barrio_hazard.merge(
    total_buildings_per_settlement, on="id_renabap"
)
barrio_ratios["ratio"] = (
    barrio_ratios["buildings_count"] / barrio_ratios["total_buildings"]
)
barrio_pop = renabap_pba_intersect[
    ["id_renabap", "familias_aproximadas"]
].drop_duplicates()
barrio_exposure = barrio_ratios.merge(barrio_pop, on="id_renabap")
barrio_exposure["fam_expuestas"] = (
    barrio_exposure["ratio"] * barrio_exposure["familias_aproximadas"]
)

# Add non-hazard population for each settlement
settlements_with_hazards = barrio_exposure["id_renabap"].unique()
all_settlements = total_buildings_per_settlement["id_renabap"].unique()

for settlement in all_settlements:
    if settlement in settlements_with_hazards:
        # Calculate non-hazard population
        hazard_pop = barrio_exposure[barrio_exposure["id_renabap"] == settlement][
            "fam_expuestas"
        ].sum()
        total_pop = barrio_pop[barrio_pop["id_renabap"] == settlement][
            "familias_aproximadas"
        ].iloc[0]
        non_hazard_pop = total_pop - hazard_pop

        if non_hazard_pop > 0:
            barrio_exposure = pd.concat(
                [
                    barrio_exposure,
                    pd.DataFrame(
                        [
                            {
                                "id_renabap": settlement,
                                "PELIGROSID": "none",
                                "buildings_count": 0,
                                "total_buildings": total_buildings_per_settlement[
                                    total_buildings_per_settlement["id_renabap"]
                                    == settlement
                                ]["total_buildings"].iloc[0],
                                "ratio": (
                                    total_buildings_per_settlement[
                                        total_buildings_per_settlement["id_renabap"]
                                        == settlement
                                    ]["total_buildings"].iloc[0]
                                    - buildings_hazard[
                                        buildings_hazard["id_renabap"] == settlement
                                    ].shape[0]
                                )
                                / total_buildings_per_settlement[
                                    total_buildings_per_settlement["id_renabap"]
                                    == settlement
                                ]["total_buildings"].iloc[0]
                                if total_buildings_per_settlement[
                                    total_buildings_per_settlement["id_renabap"]
                                    == settlement
                                ]["total_buildings"].iloc[0]
                                > 0
                                else 0,
                                "familias_aproximadas": total_pop,
                                "fam_expuestas": non_hazard_pop,
                            }
                        ]
                    ),
                ],
                ignore_index=True,
            )
    else:
        # Settlement with no hazard intersection - all population is non-hazard
        total_pop = barrio_pop[barrio_pop["id_renabap"] == settlement][
            "familias_aproximadas"
        ].iloc[0]
        barrio_exposure = pd.concat(
            [
                barrio_exposure,
                pd.DataFrame(
                    [
                        {
                            "id_renabap": settlement,
                            "PELIGROSID": "none",
                            "buildings_count": 0,
                            "total_buildings": total_buildings_per_settlement[
                                total_buildings_per_settlement["id_renabap"]
                                == settlement
                            ]["total_buildings"].iloc[0],
                            "ratio": 1.0,
                            "familias_aproximadas": total_pop,
                            "fam_expuestas": total_pop,
                        }
                    ]
                ),
            ],
            ignore_index=True,
        )

# Create buildings tidy dataframe for barrio level using consolidated function
buildings_barrio_tidy = create_exposure_tidy_data(
    data=barrio_exposure,
    id_column="id_renabap",
    peligrosidad_column="PELIGROSID",
    method_suffix="buildings",
    exposure_values=barrio_exposure["fam_expuestas"],
    exclude_zero=False,  # Keep all values including zeros for completeness
)

# 2. Cuenca exposure - using total buildings across all settlements in cuenca
buildings_cuenca_hazard = (
    buildings_hazard.groupby(["Cuenca", "PELIGROSID"])
    .size()
    .reset_index(name="buildings_count")
)
buildings_cuenca_total = (
    buildings_settlement.merge(
        settle_hazard_cuencas[["id_renabap", "Cuenca"]].drop_duplicates(),
        on="id_renabap",
    )
    .groupby("Cuenca")
    .size()
    .reset_index(name="total_buildings_all")
)

cuenca_ratios = buildings_cuenca_hazard.merge(buildings_cuenca_total, on="Cuenca")
cuenca_ratios["ratio"] = (
    cuenca_ratios["buildings_count"] / cuenca_ratios["total_buildings_all"]
)

cuenca_pop = (
    settle_hazard_cuencas.drop_duplicates("id_renabap")
    .groupby("Cuenca")["familias_aproximadas"]
    .sum()
    .reset_index()
)
cuenca_exposure = cuenca_ratios.merge(cuenca_pop, on="Cuenca")
cuenca_exposure["fam_expuestas"] = (
    cuenca_exposure["ratio"] * cuenca_exposure["familias_aproximadas"]
)

# Create buildings tidy dataframe for cuenca level using consolidated function
buildings_cuenca_tidy = create_exposure_tidy_data(
    data=cuenca_exposure,
    id_column="Cuenca",
    peligrosidad_column="PELIGROSID",
    method_suffix="buildings",
    exposure_values=cuenca_exposure["fam_expuestas"],
    exclude_zero=False,
)

# 3. Eje exposure - using total buildings across all settlements in eje
buildings_eje_hazard = (
    buildings_hazard.groupby(["eje", "PELIGROSID"])
    .size()
    .reset_index(name="buildings_count")
)
buildings_eje_total = (
    buildings_settlement.merge(
        settle_hazard_cuencas[["id_renabap", "eje"]].drop_duplicates(), on="id_renabap"
    )
    .groupby("eje")
    .size()
    .reset_index(name="total_buildings_all")
)

eje_ratios = buildings_eje_hazard.merge(buildings_eje_total, on="eje")
eje_ratios["ratio"] = eje_ratios["buildings_count"] / eje_ratios["total_buildings_all"]

eje_pop = (
    settle_hazard_cuencas.drop_duplicates("id_renabap")
    .groupby("eje")["familias_aproximadas"]
    .sum()
    .reset_index()
)
eje_exposure = eje_ratios.merge(eje_pop, on="eje")
eje_exposure["fam_expuestas"] = (
    eje_exposure["ratio"] * eje_exposure["familias_aproximadas"]
)

# Create buildings tidy dataframe for eje level using consolidated function
buildings_eje_tidy = create_exposure_tidy_data(
    data=eje_exposure,
    id_column="eje",
    peligrosidad_column="PELIGROSID",
    method_suffix="buildings",
    exposure_values=eje_exposure["fam_expuestas"],
    exclude_zero=False,
)


### GHSL 

# Convert to the format expected by rasterstats
geometries = [geom for geom in renabap_pba_intersect.geometry]

# Use rasterstats for vectorized zonal statistics
stats = rasterstats.zonal_stats(
    geometries,
    ghsl_clipped.values[0],  # rasterstats expects 2D array
    affine=ghsl_clipped.rio.transform(),
    stats=["sum"],
    nodata=ghsl_clipped.rio.nodata,
)

# Extract the sum values
ghsl_totals = [stat["sum"] if stat["sum"] is not None else 0 for stat in stats]

# Add the GHSL population estimates as a new column
renabap_pba_intersect["ghsl_pop_est"] = ghsl_totals


# Get the reference raster properties from GHSL data
reference_raster = ghsl_clipped
reference_transform = reference_raster.rio.transform()
reference_crs = reference_raster.rio.crs
reference_shape = reference_raster.shape[1:]  # Get 2D shape (height, width)


# Prepare geometries and values for rasterization
geometries_ghsl = [
    (geom, value)
    for geom, value in zip(
        renabap_pba_intersect.geometry, renabap_pba_intersect["ghsl_pop_est"]
    )
]
geometries_familias = [
    (geom, value)
    for geom, value in zip(
        renabap_pba_intersect.geometry, renabap_pba_intersect["familias_aproximadas"]
    )
]

# Create GHSL population raster
ghsl_pop_raster = rasterize(
    geometries_ghsl,
    out_shape=reference_shape,
    transform=reference_transform,
    fill=0,
    dtype=np.float32,
    all_touched=False,
)

# Create familias aproximadas raster
familias_raster = rasterize(
    geometries_familias,
    out_shape=reference_shape,
    transform=reference_transform,
    fill=0,
    dtype=np.float32,
    all_touched=False,
)


# Step 1: Divide original GHSL by the barrio-level GHSL to get fractional population
# Use masking to avoid division on invalid cells
mask = (ghsl_clipped.values[0] != NODATA_VALUE) & (ghsl_pop_raster > 0.1)
ghsl_fractional = np.full_like(ghsl_clipped.values[0], NODATA_VALUE, dtype=np.float64)
ghsl_fractional[mask] = ghsl_clipped.values[0][mask] / ghsl_pop_raster[mask]

# Step 2: Multiply fractional population by familias aproximadas to get downscaled data
mask2 = (ghsl_fractional != NODATA_VALUE) & (familias_raster > 0)
familias_downscaled = np.full_like(
    ghsl_clipped.values[0], NODATA_VALUE, dtype=np.float64
)
familias_downscaled[mask2] = ghsl_fractional[mask2] * familias_raster[mask2]

# Verify the results - exclude NODATA_VALUE from range calculations
ghsl_valid = ghsl_clipped.values[0] != NODATA_VALUE
fractional_valid = ghsl_fractional != NODATA_VALUE
downscaled_valid = familias_downscaled != NODATA_VALUE

# GHSL downscaling for all three levels using the same approach

# 1. BARRIO-HAZARD EXPOSURE using consolidated approach
ghsl_barrio_exposure = []
for idx, row in settlement_hazard.iterrows():
    stats = zonal_stats(
        [row.geometry],
        familias_downscaled,
        affine=reference_transform,
        stats=["sum"],
        nodata=NODATA_VALUE,
    )[0]

    ghsl_barrio_exposure.append(
        {
            "id_renabap": row["id_renabap"],
            "PELIGROSID": row["PELIGROSID"],
            "fam_expuestas": stats["sum"] if stats["sum"] is not None else 0,
        }
    )

ghsl_barrio_temp = pd.DataFrame(ghsl_barrio_exposure)
ghsl_barrio_tidy = create_exposure_tidy_data(
    data=ghsl_barrio_temp,
    id_column="id_renabap",
    peligrosidad_column="PELIGROSID",
    method_suffix="ghsl",
    exposure_values=ghsl_barrio_temp["fam_expuestas"],
    exclude_zero=False,
)

# 2. CUENCA-HAZARD EXPOSURE using consolidated approach
ghsl_cuenca_exposure = []
for cuenca in settle_hazard_cuencas["Cuenca"].unique():
    for peligro in settle_hazard_cuencas["PELIGROSID"].unique():
        # Get all geometries for this cuenca-hazard combination
        geoms = settle_hazard_cuencas[
            (settle_hazard_cuencas["Cuenca"] == cuenca)
            & (settle_hazard_cuencas["PELIGROSID"] == peligro)
        ].geometry.tolist()

        if geoms:
            stats = zonal_stats(
                geoms,
                familias_downscaled,
                affine=reference_transform,
                stats=["sum"],
                nodata=NODATA_VALUE,
            )

            total_pop = sum(
                [stat["sum"] if stat["sum"] is not None else 0 for stat in stats]
            )

            ghsl_cuenca_exposure.append(
                {
                    "Cuenca": cuenca,
                    "PELIGROSID": peligro,
                    "fam_expuestas": total_pop,
                }
            )

ghsl_cuenca_temp = pd.DataFrame(ghsl_cuenca_exposure)
ghsl_cuenca_tidy = create_exposure_tidy_data(
    data=ghsl_cuenca_temp,
    id_column="Cuenca",
    peligrosidad_column="PELIGROSID",
    method_suffix="ghsl",
    exposure_values=ghsl_cuenca_temp["fam_expuestas"],
    exclude_zero=False,
)

# 3. EJE-HAZARD EXPOSURE using consolidated approach
ghsl_eje_exposure = []
for eje in settle_hazard_cuencas["eje"].unique():
    for peligro in settle_hazard_cuencas["PELIGROSID"].unique():
        # Get all geometries for this eje-hazard combination
        geoms = settle_hazard_cuencas[
            (settle_hazard_cuencas["eje"] == eje)
            & (settle_hazard_cuencas["PELIGROSID"] == peligro)
        ].geometry.tolist()

        if geoms:
            stats = zonal_stats(
                geoms,
                familias_downscaled,
                affine=reference_transform,
                stats=["sum"],
                nodata=NODATA_VALUE,
            )

            total_pop = sum(
                [stat["sum"] if stat["sum"] is not None else 0 for stat in stats]
            )

            ghsl_eje_exposure.append(
                {
                    "eje": eje,
                    "PELIGROSID": peligro,
                    "fam_expuestas": total_pop,
                }
            )

ghsl_eje_temp = pd.DataFrame(ghsl_eje_exposure)
ghsl_eje_tidy = create_exposure_tidy_data(
    data=ghsl_eje_temp,
    id_column="eje",
    peligrosidad_column="PELIGROSID",
    method_suffix="ghsl",
    exposure_values=ghsl_eje_temp["fam_expuestas"],
    exclude_zero=False,
)

3.6 Resultados

3.6.1 Comparación de métodos

# 1. BARRIO-LEVEL WIDE DATAFRAME using consolidated function
renabap_info = renabap_pba_intersect[["id_renabap", "nombre_barrio"]].drop_duplicates()

barrio_wide = create_wide_exposure_dataframe(
    areal_data=barrio_areal_tidy,
    ghsl_data=ghsl_barrio_tidy,
    buildings_data=buildings_barrio_tidy,
    id_columns=["id_renabap", "peligrosidad"],
    exclude_hazard_value=NON_HAZARD_VALUE,
)

# Add barrio names and reorder columns
barrio_wide = barrio_wide.merge(renabap_info, on="id_renabap", how="left")
barrio_wide = barrio_wide[
    [
        "id_renabap",
        "nombre_barrio",
        "peligrosidad",
    ]
    + EXPOSURE_COLUMNS
]

# 2. CUENCA-LEVEL WIDE DATAFRAME using consolidated function
settlement_cuenca_mapping = settle_hazard_cuencas[
    ["id_renabap", "Cuenca"]
].drop_duplicates()
settlement_counts = settlement_cuenca_mapping["id_renabap"].value_counts()
multi_cuenca_settlements = settlement_counts[settlement_counts > 1].index
if len(multi_cuenca_settlements) > 0:
    settlement_cuenca_mapping = settlement_cuenca_mapping.drop_duplicates(
        subset=["id_renabap"], keep="first"
    )

# Aggregate buildings data from barrio to cuenca level
cuenca_buildings_wide = (
    buildings_barrio_tidy[buildings_barrio_tidy["peligrosidad"] != NON_HAZARD_VALUE]
    .merge(settlement_cuenca_mapping, on="id_renabap", how="left")
    .groupby(["Cuenca", "peligrosidad"])["fam_expuestas_buildings"]
    .sum()
    .reset_index()
)

# Convert to tidy format for the consolidated function
cuenca_buildings_tidy_for_merge = create_exposure_tidy_data(
    data=cuenca_buildings_wide,
    id_column="Cuenca",
    peligrosidad_column="peligrosidad",
    method_suffix="buildings",
    exposure_values=cuenca_buildings_wide["fam_expuestas_buildings"],
    exclude_zero=False,
)

cuenca_wide = create_wide_exposure_dataframe(
    areal_data=cuenca_areal_tidy,
    ghsl_data=ghsl_cuenca_tidy,
    buildings_data=cuenca_buildings_tidy_for_merge,
    id_columns=["Cuenca", "peligrosidad"],
    exclude_hazard_value=NON_HAZARD_VALUE,
)

cuenca_wide["cuenca"] = cuenca_wide["Cuenca"].str.lower()
cuenca_wide = cuenca_wide[
    [
        "cuenca",
        "peligrosidad",
    ]
    + EXPOSURE_COLUMNS
]

# 3. EJE-LEVEL WIDE DATAFRAME using consolidated function
settlement_eje_mapping = settle_hazard_cuencas[["id_renabap", "eje"]].drop_duplicates()
eje_settlement_counts = settlement_eje_mapping["id_renabap"].value_counts()
multi_eje_settlements = eje_settlement_counts[eje_settlement_counts > 1].index
if len(multi_eje_settlements) > 0:
    settlement_eje_mapping = settlement_eje_mapping.drop_duplicates(
        subset=["id_renabap"], keep="first"
    )

# Aggregate buildings data from barrio to eje level
eje_buildings_wide = (
    buildings_barrio_tidy[buildings_barrio_tidy["peligrosidad"] != NON_HAZARD_VALUE]
    .merge(settlement_eje_mapping, on="id_renabap", how="left")
    .groupby(["eje", "peligrosidad"])["fam_expuestas_buildings"]
    .sum()
    .reset_index()
)

# Convert to tidy format for the consolidated function
eje_buildings_tidy_for_merge = create_exposure_tidy_data(
    data=eje_buildings_wide,
    id_column="eje",
    peligrosidad_column="peligrosidad",
    method_suffix="buildings",
    exposure_values=eje_buildings_wide["fam_expuestas_buildings"],
    exclude_zero=False,
)

eje_wide = create_wide_exposure_dataframe(
    areal_data=eje_areal_tidy,
    ghsl_data=ghsl_eje_tidy,
    buildings_data=eje_buildings_tidy_for_merge,
    id_columns=["eje", "peligrosidad"],
    exclude_hazard_value=NON_HAZARD_VALUE,
)

eje_wide = eje_wide[
    [
        "eje",
        "peligrosidad",
    ]
    + EXPOSURE_COLUMNS
]

barrio_tidy = pd.melt(
    barrio_wide,
    id_vars=["id_renabap", "nombre_barrio", "peligrosidad"],
    value_vars=EXPOSURE_COLUMNS,
    var_name="metodo",
    value_name="fam_expuestas",
)

barrio_tidy["metodo"] = barrio_tidy["metodo"].str.replace(
    COLUMN_MAPPINGS["method_cleanup_prefix"], ""
)

barrio_tidy = barrio_tidy.merge(
    renabap_pba_intersect[["id_renabap", "geometry"]], on="id_renabap", how="left"
)

barrio_tidy["area"] = barrio_tidy.geometry.apply(lambda geom: geom.area)


# Group by id_renabap and peligro, then find which method has the highest fam_expuestas
highest_methods = barrio_tidy.groupby(["id_renabap", "peligrosidad"])[
    "fam_expuestas"
].idxmax()

# Get the method names for the highest estimates
method_counts = barrio_tidy.loc[highest_methods, "metodo"].value_counts()

plt.figure(figsize=(12, 7))
sns.barplot(x=method_counts.index, y=method_counts.values, hue=method_counts.index, palette="viridis", legend=False)
plt.title(
    "Métodos que Más Frecuentemente Devuelven la Mayor Estimación de Familias Expuestas",
    fontsize=14,
    fontweight="bold",
    pad=20,
)
plt.xlabel("Método", fontsize=12, fontweight="bold")
plt.ylabel("Frecuencia", fontsize=12, fontweight="bold")
plt.xticks(rotation=45, ha="right")

plt.tight_layout()
plt.show()


# Group by id_renabap and peligro, then find which method has the lowest fam_expuestas
lowest_methods = barrio_tidy.groupby(["id_renabap", "peligrosidad"])[
    "fam_expuestas"
].idxmin()

# Get the method names for the lowest estimates
lowest_method_counts = barrio_tidy.loc[lowest_methods, "metodo"].value_counts()

# Crear gráfico de barras para métodos con menor estimación
plt.figure(figsize=(12, 7))
sns.barplot(
    x=lowest_method_counts.index, y=lowest_method_counts.values, hue=lowest_method_counts.index, palette="viridis", legend=False
)
plt.title(
    "Métodos que Más Frecuentemente Devuelven la Menor Estimación de Familias Expuestas",
    fontsize=14,
    fontweight="bold",
    pad=20,
)
plt.xlabel("Método", fontsize=12, fontweight="bold")
plt.ylabel("Frecuencia", fontsize=12, fontweight="bold")
plt.xticks(rotation=45, ha="right")

plt.tight_layout()
plt.show()


# First, left join familias_aproximadas from renabap_pba_intersect
final_tidy_with_pop = barrio_tidy.merge(
    renabap_pba_intersect[["id_renabap", "familias_aproximadas"]],
    on="id_renabap",
    how="left",
)

# Calculate the range (highest - lowest) per id_renabap and peligro
range_by_barrio = final_tidy_with_pop.groupby(["id_renabap", "peligrosidad"])[
    "fam_expuestas"
].agg(["max", "min"])
range_by_barrio["range"] = range_by_barrio["max"] - range_by_barrio["min"]

# Merge back to get the total population for each barrio
range_by_barrio = range_by_barrio.reset_index().merge(
    final_tidy_with_pop[["id_renabap", "familias_aproximadas"]].drop_duplicates(),
    on="id_renabap",
)

# Calculate absolute percent difference as fraction of total barrio population
range_by_barrio["abs_pct_diff"] = (
    range_by_barrio["range"] / range_by_barrio["familias_aproximadas"]
) * 100

# Calculate average absolute percent difference per peligro level
avg_pct_diff_by_peligro = range_by_barrio.groupby("peligrosidad")["abs_pct_diff"].mean()

# Calculate absolute percent difference by method
method_errors = final_tidy_with_pop.copy()


# Calculate absolute error as percent of total population for each method
method_errors["abs_error_pct"] = (
    abs(
        method_errors["fam_expuestas"]
        - method_errors.groupby(["id_renabap", "peligrosidad"])[
            "fam_expuestas"
        ].transform("mean")
    )
    / method_errors["familias_aproximadas"]
    * 100
)

# Calculate coefficient of variation for each barrio-peligro combination
barrio_reliability = (
    final_tidy_with_pop.groupby(["id_renabap", "peligrosidad"])
    .agg({"fam_expuestas": ["mean", "std"], "familias_aproximadas": "first"})
    .reset_index()
)

barrio_reliability.columns = [
    "id_renabap",
    "peligrosidad",
    "mean_estimate",
    "std_estimate",
    "familias_aproximadas",
]
barrio_reliability["coefficient_variation"] = (
    barrio_reliability["std_estimate"] / barrio_reliability["mean_estimate"]
)

# Create box plot
plt.figure(figsize=(12, 3))
sns.boxplot(
    data=barrio_reliability,
    y="peligrosidad",
    x="coefficient_variation",
    hue="peligrosidad",
    palette="viridis",
    legend=False,
    width=0.4,
)
plt.title(
    "Variabilidad de Estimaciones entre Métodos", fontsize=16, fontweight="bold", pad=20
)
plt.ylabel("Peligrosidad", fontsize=14, fontweight="bold")
plt.xlabel(
    "Coeficiente de Variación (0 = estimaciones idénticas, 1 = muy variables)",
    fontsize=12,
    fontweight="bold",
)
plt.tight_layout()
plt.show()


# Create scatter plot colored by peligrosidad
plt.figure(figsize=(10, 6))

# Get unique peligrosidad levels
peligrosidad_levels = range_by_barrio["peligrosidad"].unique()

for peligro in peligrosidad_levels:
    # Filter data for this peligrosidad level
    peligro_data = range_by_barrio[range_by_barrio["peligrosidad"] == peligro]

    plt.scatter(
        peligro_data["familias_aproximadas"],
        peligro_data["abs_pct_diff"],
        alpha=0.7,
        label=f"Peligro: {peligro}",
    )

plt.xlabel("Familias Aproximadas (Total Barrio Population)")
plt.ylabel("Absolute Percent Difference (%)")
plt.title("Method Disagreement vs Barrio Size (Colored by Peligrosidad)")

plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

# First, get the area data from final_tidy
area_data = barrio_tidy[["id_renabap", "area"]].drop_duplicates()

# Merge area back into range_by_barrio
range_by_barrio_with_area = range_by_barrio.merge(
    area_data, on="id_renabap", how="left"
)


plt.figure(figsize=(10, 6))

# Get unique peligrosidad levels
peligrosidad_levels = range_by_barrio_with_area["peligrosidad"].unique()

for peligro in peligrosidad_levels:
    # Filter data for this peligrosidad level
    peligro_data = range_by_barrio_with_area[
        range_by_barrio_with_area["peligrosidad"] == peligro
    ]

    plt.scatter(
        peligro_data["area"],
        peligro_data["abs_pct_diff"],
        alpha=0.7,
        label=f"Peligro: {peligro}",
    )

plt.xlabel("Area")
plt.ylabel("Absolute Percent Difference (%)")
plt.title("Method Disagreement vs Area (Colored by Peligrosidad)")

plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()


# Filter for high exposure (alta peligrosidad) using the joined dataframe
alta_data = barrio_tidy[barrio_tidy["peligrosidad"] == "alta"].copy()

# Aggregate by nombre_barrio and sum fam_expuestas for each method
# This handles cases where there are multiple geometries with the same barrio name
alta_aggregated = (
    alta_data.groupby(["nombre_barrio", "metodo"])["fam_expuestas"].sum().reset_index()
)

# Remove cases where the barrio name is "Sin Nombre"
alta_aggregated = alta_aggregated[
    alta_aggregated["nombre_barrio"] != "Sin Nombre"
].copy()

# Calculate total exposure per barrio across all methods
total_exposure = (
    alta_aggregated.groupby("nombre_barrio")["fam_expuestas"]
    .sum()
    .sort_values(ascending=False)
)
top_10_barrios = total_exposure.head(10).index

# Filter aggregated data for top 10 barrios
top_10_data = alta_aggregated[
    alta_aggregated["nombre_barrio"].isin(top_10_barrios)
].copy()

# Create range plot showing min, max, and individual points
plt.figure(figsize=(14, 10))  # Increased height to accommodate longer barrio names


for i, barrio in enumerate(top_10_barrios):
    barrio_data = top_10_data[top_10_data["nombre_barrio"] == barrio]
    if len(barrio_data) > 0:
        values = barrio_data["fam_expuestas"].values
        min_val = values.min()
        max_val = values.max()

        # Plot range line
        plt.plot([min_val, max_val], [i, i], "k-", alpha=0.5, linewidth=2)

        # Plot individual points colored by method
        for _, row in barrio_data.iterrows():
            color = METHOD_COLORS[row["metodo"]]
            plt.plot(row["fam_expuestas"], i, "o", color=color, markersize=6, alpha=0.8)

plt.yticks(range(len(top_10_barrios)), top_10_barrios)
plt.xlabel("Familias Expuestas")
plt.ylabel("Barrio")
plt.title("Range of High Exposure Estimates for Top 10 Barrios", fontsize=14)
plt.grid(True, alpha=0.3)

# Add legend
legend_elements = [
    plt.Line2D(
        [0],
        [0],
        marker="o",
        color="w",
        markerfacecolor=color,
        markersize=8,
        label=method,
    )
    for method, color in METHOD_COLORS.items()
]
plt.legend(handles=legend_elements, title="Método")

plt.tight_layout()
plt.show()

area_data = (
    barrio_tidy[barrio_tidy["metodo"] == "fam_exp_areal"]
    .groupby(["nombre_barrio", "peligrosidad"])["fam_expuestas"]
    .sum()
    .reset_index()
)
area_data = area_data.rename(columns={"fam_expuestas": "fam_exp_areal"})

ghsl_data = (
    barrio_tidy[barrio_tidy["metodo"] == "fam_exp_ghsl"]
    .groupby(["nombre_barrio", "peligrosidad"])["fam_expuestas"]
    .sum()
    .reset_index()
)
ghsl_data = ghsl_data.rename(columns={"fam_expuestas": "fam_exp_ghsl"})

edificios_data = (
    barrio_tidy[barrio_tidy["metodo"] == "fam_exp_edificios"]
    .groupby(["nombre_barrio", "peligrosidad"])["fam_expuestas"]
    .sum()
    .reset_index()
)
edificios_data = edificios_data.rename(
    columns={"fam_expuestas": "fam_exp_edificios"}
)

# Merge all methods together
barrio_summary = area_data.merge(
    ghsl_data, on=["nombre_barrio", "peligrosidad"], how="outer"
)
barrio_summary = barrio_summary.merge(
    edificios_data, on=["nombre_barrio", "peligrosidad"], how="outer"
)

barrio_summary = barrio_summary.fillna(0)


# Sort by nombre_barrio and peligrosidad in descending order
barrio_summary = barrio_summary.sort_values(
    ["nombre_barrio", "peligrosidad"], ascending=True
)

(a) Discrepancia vs población total del barrio

(b) Discrepancia vs área del barrio

(c) Rango de estimaciones para barrios con mayor exposición

(d)

(e)

(f)

Figure 3.2: Análisis comparativo de métodos de estimación por barrio

def plot_method_map(
    method_name,
    final_tidy,
    la_plata,
    PELIGROSIDAD_COLORS,
    common_bounds,
    create_consistent_map,
):
    fig, ax = create_consistent_map(
        f"Barrios por población expuesta estimada - {method_name.capitalize()}",
        common_bounds,
    )

    method_data = final_tidy[
        (final_tidy["metodo"] == method_name)
        & (final_tidy["peligrosidad"].isin(["alta", "media"]))
    ].copy()

    method_gdf = gpd.GeoDataFrame(method_data, geometry="geometry", crs=USE_CRS)
    method_gdf = method_gdf.clip(la_plata)
    method_gdf_3857 = method_gdf.to_crs(WEB_MERCATOR_CRS)

    plotting_order = ["media", "alta"]

    np.random.seed(42)
    for peligrosidad in plotting_order:
        level_data = method_gdf_3857[method_gdf_3857["peligrosidad"] == peligrosidad]
        for _, row in level_data.iterrows():
            centroid = row["geometry"].centroid
            jitter_x = np.random.uniform(-200, 200)
            jitter_y = np.random.uniform(-200, 200)
            x_pos = centroid.x + jitter_x
            y_pos = centroid.y + jitter_y
            color = PELIGROSIDAD_COLORS[row["peligrosidad"]]
            size = max(10, row["fam_expuestas"] * 2 + 15)
            ax.scatter(
                x_pos,
                y_pos,
                s=size,
                color=color,
                alpha=0.9,
                edgecolors="white",
                linewidth=1.0,
            )

    legend_elements = [
        plt.Line2D(
            [0],
            [0],
            marker="o",
            color="w",
            markerfacecolor=color,
            markersize=8,
            label=level.capitalize(),
        )
        for level, color in PELIGROSIDAD_COLORS.items()
    ]
    ax.legend(handles=legend_elements, title="Nivel de Peligrosidad", loc="upper right")
    plt.tight_layout()
    plt.show()


area_data = (
    barrio_tidy[barrio_tidy["metodo"] == "fam_exp_areal"]
    .groupby(["nombre_barrio", "peligrosidad"])["fam_expuestas"]
    .sum()
    .reset_index()
)
area_data = area_data.rename(columns={"fam_expuestas": "fam_exp_areal"})

ghsl_data = (
    barrio_tidy[barrio_tidy["metodo"] == "fam_exp_ghsl"]
    .groupby(["nombre_barrio", "peligrosidad"])["fam_expuestas"]
    .sum()
    .reset_index()
)
ghsl_data = ghsl_data.rename(columns={"fam_expuestas": "fam_exp_ghsl"})

edificios_data = (
    barrio_tidy[barrio_tidy["metodo"] == "fam_exp_edificios"]
    .groupby(["nombre_barrio", "peligrosidad"])["fam_expuestas"]
    .sum()
    .reset_index()
)
edificios_data = edificios_data.rename(
    columns={"fam_expuestas": "fam_exp_edificios"}
)

# Merge all methods together
barrio_summary = area_data.merge(
    ghsl_data, on=["nombre_barrio", "peligrosidad"], how="outer"
)
barrio_summary = barrio_summary.merge(
    edificios_data, on=["nombre_barrio", "peligrosidad"], how="outer"
)

barrio_summary = barrio_summary.fillna(0)


# Sort by nombre_barrio and peligrosidad in descending order
barrio_summary = barrio_summary.sort_values(
    ["nombre_barrio", "peligrosidad"], ascending=True
)

3.6.2 Exposición por barrio

methods = barrio_tidy["metodo"].unique()

for method in methods:
    plot_method_map(
        method,
        barrio_tidy,
        la_plata,
        PELIGROSIDAD_COLORS,
        common_bounds,
        create_consistent_map,
    )

(a) Estimación basada en edificios

(b) Estimación basada en GHSL

(c) Estimación basada en interpolación areal

Figure 3.3: Comparación de métodos de estimación de exposición por barrio

show(round_numeric_columns(barrio_summary))

Loading ITables v2.4.4 from the internet... (need help?)

3.6.3 Exposición por cuenca

cuencas_centroids = cuencas.copy()
cuencas_centroids["geometry"] = cuencas_centroids["geometry"].centroid

# Create a lowercase version of Cuenca for matching
cuencas_centroids["cuenca_lower"] = cuencas_centroids["Cuenca"].str.lower()

cuenca_tidy = pd.melt(
    cuenca_wide,
    id_vars=["cuenca", "peligrosidad"],
    value_vars=EXPOSURE_COLUMNS,
    var_name="metodo",
    value_name="fam_expuestas",
)


cuenca_tidy["metodo"] = cuenca_tidy["metodo"].str.replace(
    COLUMN_MAPPINGS["method_cleanup_prefix"], ""
)

cuenca_tidy_with_geometry = cuenca_tidy.merge(
    cuencas_centroids[["cuenca_lower", "geometry"]],
    left_on="cuenca",
    right_on="cuenca_lower",
    how="left",
)

cuenca_tidy_gdf = gpd.GeoDataFrame(
    cuenca_tidy_with_geometry, geometry="geometry", crs=cuencas.crs
)

methods = cuenca_tidy_gdf["metodo"].unique()

for method in methods:
    plot_method_map(
        method,
        cuenca_tidy_gdf,
        la_plata,
        PELIGROSIDAD_COLORS,
        common_bounds,
        create_consistent_map,
    )

(a) Estimación basada en edificios

(b) Estimación basada en GHSL

(c) Estimación basada en interpolación areal

Figure 3.4: Comparación de métodos de estimación de exposición por cuenca

show(round_numeric_columns(cuenca_wide))

Loading ITables v2.4.4 from the internet... (need help?)

3.6.4 Exposición por eje

eje_tidy = pd.melt(
    eje_wide,
    id_vars=["eje", "peligrosidad"],
    value_vars=EXPOSURE_COLUMNS,
    var_name="metodo",
    value_name="fam_expuestas",
)


eje_tidy["metodo"] = eje_tidy["metodo"].str.replace(
    COLUMN_MAPPINGS["method_cleanup_prefix"], ""
)

# Create eje geodataframe by dissolving cuencas by eje and then taking centroids
ejes = cuencas.dissolve(by="eje").reset_index()
ejes_centroids = ejes.copy()
ejes_centroids["geometry"] = ejes_centroids["geometry"].centroid

eje_tidy_with_geometry = eje_tidy.merge(
    ejes_centroids[["eje", "geometry"]], on="eje", how="left"
)

eje_tidy_gdf = gpd.GeoDataFrame(
    eje_tidy_with_geometry, geometry="geometry", crs=cuencas.crs
)

methods = eje_tidy_gdf["metodo"].unique()

for method in methods:
    plot_method_map(
        method,
        eje_tidy_gdf,
        la_plata,
        PELIGROSIDAD_COLORS,
        common_bounds,
        create_consistent_map,
    )

(a) Estimación basada en edificios

(b) Estimación basada en GHSL

(c) Estimación basada en interpolación areal

Figure 3.5: Comparación de métodos de estimación de exposición por eje de cuenca

show(round_numeric_columns(eje_wide))

Loading ITables v2.4.4 from the internet... (need help?)

3.7 Conclusiones

Romanazzi, Pablo et al. 2019. Plan de Reducción Del Riesgo Por Inundaciones En La Región de La Plata. Edited by Sebastián Guerrini, Pablo Morosi, Eduardo Pablo Spinelli, and Josefina López MacKenzie. 1st ed. La Plata: Universidad Nacional de La Plata. Facultad de Ingeniería; Municipalidad de La Plata. https://sedici.unlp.edu.ar/bitstream/handle/10915/165109/Documento_completo.pdf-PDFA.pdf?sequence=1&isAllowed=y.

Schiavina, M., S. Freire, A. Carioli, and K. MacManus. 2023. “GHS-POP R2023A - GHS Population Grid Multitemporal (1975-2030).” European Commission, Joint Research Centre (JRC). https://doi.org/10.2905/2FF68A52-5B5B-4A22-8F40-C41DA8332CFE.

Smith, A., P. D. Bates, O. Wing, et al. 2019. “New Estimates of Flood Exposure in Developing Countries Using High-Resolution Population Data.” Nature Communications 10: 1814. https://doi.org/10.1038/s41467-019-09282-y.

Tellman, B., J. A. Sullivan, C. Kuhn, et al. 2021. “Satellite Imaging Reveals Increased Proportion of Population Exposed to Floods.” Nature 596: 80–86. https://doi.org/10.1038/s41586-021-03695-w.

VIDA. 2023. “Google-Microsoft-OSM Open Buildings - Combined by VIDA.” https://source.coop/repositories/vida/google-microsoft-osm-open-buildings/access.