1  Exposición

Mapeo dasimétrico con datos GHSL

1.1 Introducion

En este flujo de trabajo aplicamos mapeo dasimétrico para redistribuir proporcionalmente datos censales argentinos de 2022 desde el nivel de radio censal hasta una resolución de 100 metros para la localidad de Esperanza en la provincia de Santa Fe, Argentina, utilizando datos auxiliares de la Capa Global de Asentamientos Humanos (GHSL) de 2023. Este método representa una forma razonablemente efectiva de reducir la escala de datos censales a mayor resolución para su uso en modelado de riesgo climático con fuentes de datos abiertos.

import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import rasterstats
from rasterio.features import rasterize
from io import BytesIO
from owslib.wfs import WebFeatureService

import rioxarray
import contextily as ctx
from shapely.geometry import box
import xarray as xr


USE_CRS = "EPSG:5347"  # posgar para esperanza
WEB_MERCATOR_CRS = "EPSG:3857"

DEFAULT_FIGSIZE = (12, 10)
MAP_PADDING = 500
PLASMA_CMAP = plt.cm.plasma


def setup_base_map(
    figsize=None, bounds=None, boundary_gdf=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 and boundary_gdf is not None:
        bounds = boundary_gdf.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."""
    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_boundary_outline(ax, boundary_gdf, crs="EPSG:3857"):
    """Add the outline of a boundary geodataframe to a map."""
    boundary_3857 = boundary_gdf.to_crs(crs)
    boundary_3857.plot(
        ax=ax,
        facecolor="none",
        edgecolor="black",
        linewidth=0.5,
        linestyle="--",
        legend=False,
        zorder=5,
    )


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

    add_basemap(ax)

    add_north_arrow(ax)

    add_boundary_outline(ax, boundary_gdf)

    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


base_url = "https://wms.ign.gob.ar/geoserver/ign/ows"

munis = wfs_to_gdf(wfs_url=base_url, layer_name="ign:municipio", srs="EPSG:4326")

esperanza = munis[munis["nam"] == "Esperanza"]
esperanza = esperanza.to_crs(USE_CRS)

1.2 Qué es el mapeo dasimétrico?

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.

1.3 Fuentes de datos

1.3.1 Censo Argentino

Los datos censales provienen del Censo Nacional de Población, Hogares y Viviendas 2022 de Argentina, realizado por el Instituto Nacional de Estadística y Censos (INDEC) (Instituto Nacional de Estadística y Censos (INDEC) 2024). Los datos fueron descargados desde la plataforma Redatam y procesados en R en un repositorio separado, ya que resultó más eficiente trabajar con esa herramienta para el procesamiento inicial. Posteriormente, estos datos fueron unidos con el conjunto de datos de radios censales de Argentina y guardados localmente. Se planea hacer disponible una copia de los datos en formato geoespacial como archivo Parquet almacenado en la nube, pero este trabajo aún está en progreso.

datos_censales = gpd.read_parquet(
    "/home/nissim/Documents/dev/fulbright/ciut-redatam/datos_censales_2022_geo.parquet"
)

datos_censales = datos_censales.to_crs(USE_CRS)

geometria_esperanza = esperanza.geometry.iloc[0]
centroides_dentro = datos_censales.geometry.centroid.within(geometria_esperanza)
completamente_dentro = datos_censales.within(geometria_esperanza)
datos_censales_esperanza = datos_censales[
    completamente_dentro | centroides_dentro
].copy()

fig, ax = create_consistent_map(
    "Población Total por Radio Censal - Censo 2022", esperanza
)

datos_censales_esperanza_3857 = datos_censales_esperanza.to_crs(WEB_MERCATOR_CRS)

datos_censales_esperanza_3857.plot(
    column="POB_TOT_P", ax=ax, cmap=PLASMA_CMAP, legend=False, alpha=0.8, zorder=2
)

plt.title("Población Total por Radio Censal - Censo 2022", fontsize=16, fontweight="bold", pad=20)
plt.show()

1.3.2 Datos 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. Para este análisis se utilizan los datos GHSL de 2023, que son los más recientes disponibles y los más cercanos temporalmente al Censo Argentino 2022. 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).

ghsl = rioxarray.open_rasterio(
    "/home/nissim/Documents/dev/fulbright/ciut-riesgo/notebooks/data/GHS_POP_E2025_GLOBE_R2023A_54009_100_V1_0_R13_C13/GHS_POP_E2025_GLOBE_R2023A_54009_100_V1_0_R13_C13.tif",
    chunks={"x": 1024, "y": 1024},
)

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

ghsl_clipped = ghsl.rio.clip(
    [esperanza.geometry.iloc[0]],
    from_disk=True,
)

ghsl_masked = ghsl_clipped.where(ghsl_clipped > 0)

fig, ax = create_consistent_map(
    "Estimaciones de Población GHSL 2023 - Esperanza, Santa Fe", esperanza
)

ghsl_masked_3857 = ghsl_masked.rio.reproject(WEB_MERCATOR_CRS)

ghsl_masked_3857 = ghsl_masked_3857.where(ghsl_masked_3857 > 0)

ghsl_masked_3857.plot(ax=ax, cmap=PLASMA_CMAP, alpha=0.75, add_colorbar=False, add_labels=False, zorder=2)

plt.title("Estimaciones de Población GHSL 2023 - Esperanza, Santa Fe", fontsize=16, fontweight="bold", pad=20)
plt.show()

1.4 Resultados

geometrias = [geom for geom in datos_censales_esperanza.geometry]

estadisticas = rasterstats.zonal_stats(
    geometrias,
    ghsl_clipped.values[0],
    affine=ghsl_clipped.rio.transform(),
    stats=["sum"],
    nodata=-200,
)

totales_ghsl = [stat["sum"] if stat["sum"] is not None else 0 for stat in estadisticas]

datos_censales_esperanza["estimacion_pob_ghsl"] = totales_ghsl


raster_referencia = ghsl_clipped
transformacion_referencia = raster_referencia.rio.transform()
crs_referencia = raster_referencia.rio.crs
forma_referencia = raster_referencia.shape[1:]
raster_ghsl = ghsl_clipped.values[0]


geometrias_ghsl = [
    (geom, valor)
    for geom, valor in zip(
        datos_censales_esperanza.geometry,
        datos_censales_esperanza["estimacion_pob_ghsl"],
    )
]
geometrias_pob = [
    (geom, valor)
    for geom, valor in zip(
        datos_censales_esperanza.geometry, datos_censales_esperanza["POB_TOT_P"]
    )
]

raster_pob_ghsl = rasterize(
    geometrias_ghsl,
    out_shape=forma_referencia,
    transform=transformacion_referencia,
    fill=0,
    dtype=np.float32,
    all_touched=True,
)

raster_pob_censo = rasterize(
    geometrias_pob,
    out_shape=forma_referencia,
    transform=transformacion_referencia,
    fill=0,
    dtype=np.float32,
    all_touched=True,
)


mascara = (raster_ghsl > 0) & (raster_ghsl != -200) & (raster_pob_ghsl > 0.1)
ghsl_fraccional = np.full_like(raster_ghsl, 0, dtype=np.float64)
ghsl_fraccional[mascara] = raster_ghsl[mascara] / raster_pob_ghsl[mascara]

mascara2 = (ghsl_fraccional > 0) & (raster_pob_censo > 0)
pob_redistribuida = np.full_like(raster_ghsl, 0, dtype=np.float64)
pob_redistribuida[mascara2] = ghsl_fraccional[mascara2] * raster_pob_censo[mascara2]

pob_redistribuida_da = xr.DataArray(
    pob_redistribuida,
    coords={"y": ghsl_clipped.y, "x": ghsl_clipped.x},
    dims=["y", "x"],
    attrs=ghsl_clipped.attrs.copy(),
)

pob_redistribuida_da = pob_redistribuida_da.rio.write_crs(USE_CRS)

pob_redistribuida_enmascarada = pob_redistribuida_da.where(pob_redistribuida_da > 0)

total_redistribuido = pob_redistribuida[pob_redistribuida > 0].sum()
total_censo = datos_censales_esperanza["POB_TOT_P"].sum()

fig, ax = create_consistent_map(
    "Población Redistribuida a 100m - Esperanza, Santa Fe", esperanza
)

pob_redistribuida_enmascarada_3857 = pob_redistribuida_enmascarada.rio.reproject(
    WEB_MERCATOR_CRS
)

pob_redistribuida_enmascarada_3857 = pob_redistribuida_enmascarada_3857.where(
    pob_redistribuida_enmascarada_3857 > 0
)

pob_redistribuida_enmascarada_3857.plot(
    ax=ax, cmap=PLASMA_CMAP, alpha=0.75, add_colorbar=False, add_labels=False, zorder=2
)

plt.title("Población Censal Redistribuida a 100m - Esperanza, Santa Fe", fontsize=16, fontweight="bold", pad=20)
plt.show()

Este flujo de trabajo demuestra la capacidad de redistribuir exitosamente datos de población censal a resolución de 100 metros utilizando datos GHSL para mapeo dasimétrico. El proceso logra una conservación de población del 97.6%, redistribuyendo 45,616 personas de un total censal de 46,757 habitantes. Esta metodología permite una representación espacial más precisa de la distribución poblacional, manteniendo la consistencia con los datos censales oficiales.

Instituto Nacional de Estadística y Censos (INDEC). 2024. “Bases de Datos.” https://www.indec.gob.ar/indec/web/Institucional-Indec-BasesDeDatos-6.

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.