1 Configuración y datos reales

Análisis integrado de PM2.5 en Alemania (año 2022) con datos reales descargados en vivo: estaciones y concentración de PM2.5 del Umweltbundesamt (UBA) —la red oficial alemana que reporta a la EEA—, instalaciones industriales de la EEA (IED/E-PRTR), y regiones NUTS-3 con densidad poblacional de Eurostat. Al knitear en RStudio, el informe reproduce estos resultados reales.

Qué hace este bloque: carga los paquetes espaciales (sf, gstat, spatstat, spdep, eurostat) y fija los parámetros del estudio (país, contaminante, año, sistema de coordenadas métrico).

## Paquetes (CRAN) y parametros. Requiere Internet para descargar los datos reales.
pkgs <- c("sf","gstat","spatstat","spdep","eurostat","giscoR","terra","geodata",
          "httr2","jsonlite","dplyr","ggplot2","viridis","units")
inst <- pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly = TRUE)]
if (length(inst)) install.packages(inst, repos = "https://cloud.r-project.org")
suppressPackageStartupMessages({
  library(sf); library(gstat); library(spatstat); library(spdep); library(eurostat)
  library(terra); library(geodata)
  library(httr2); library(jsonlite); library(dplyr); library(ggplot2); library(viridis); library(units)})
suppressWarnings(sf::sf_use_s2(FALSE))
set.seed(2026); CRS_M <- 3035; YEAR <- 2022
cat("Entorno listo. Alemania (DE) | PM2.5 | año", YEAR, "\n")
## Entorno listo. Alemania (DE) | PM2.5 | año 2022

Qué hace este bloque: descarga en vivo las fuentes reales — PM2.5 por estación (UBA), instalaciones emisoras de PM (EEA), NUTS-3 + densidad (Eurostat) — y además las covariables reales para el kriging con covariables y el modelo de intensidad: elevación (DEM satelital, Copernicus/SRTM) y % de suelo construido (uso del suelo, equivalente a CORINE). Las extrae en las estaciones y en la grilla, y las prepara como imágenes para spatstat.

## ============ DATOS REALES (descarga en vivo) ============
## (A) Estaciones y PM2.5 anual REAL -- Umweltbundesamt (UBA), red oficial alemana (reporta a la EEA)
uba_get <- function(path, ...) {
  r <- request(paste0("https://www.umweltbundesamt.de/api/air_data/v3/", path)) |>
    req_url_query(..., lang = "en") |> req_perform()
  jsonlite::fromJSON(resp_body_string(r), simplifyVector = FALSE)
}
ab <- uba_get("annualbalances/json", component = 9, year = YEAR)$data      # PM2.5 = componente 9
pm <- do.call(rbind, lapply(ab, function(v)
  data.frame(station_id = as.character(v[[1]]), pm25 = as.numeric(v[[2]]))))
stj <- uba_get("stations/json")$data                                # estaciones bajo $data
coords <- do.call(rbind, lapply(names(stj), function(k) { v <- stj[[k]]
  if (length(v) < 9) return(NULL)
  data.frame(station_id = as.character(v[[1]]),
             lon = as.numeric(v[[8]]), lat = as.numeric(v[[9]]),
             entorno = if (length(v) >= 15) as.character(v[[15]]) else NA_character_,
             tipo    = if (length(v) >= 17) as.character(v[[17]]) else NA_character_,
             stringsAsFactors = FALSE) }))
stations <- merge(pm, coords, by = "station_id")
stations <- stations[is.finite(stations$lon) & is.finite(stations$lat) & is.finite(stations$pm25), ]
stations <- st_transform(st_as_sf(stations, coords = c("lon","lat"), crs = 4326), CRS_M)
cat("Estaciones PM2.5 (UBA,", YEAR, "):", nrow(stations), "\n")
## Estaciones PM2.5 (UBA, 2022 ): 253
## (B) NUTS-3 (400) y densidad poblacional REAL -- Eurostat
nuts3 <- tryCatch(
  get_eurostat_geospatial(output_class = "sf", resolution = "01",
                          nuts_level = 3, year = 2021) |>
    dplyr::filter(CNTR_CODE == "DE") |> st_transform(CRS_M),
  error = function(e) {                                     # respaldo: GeoJSON directo de GISCO
    u <- "https://gisco-services.ec.europa.eu/distribution/v2/nuts/geojson/NUTS_RG_01M_2021_3035_LEVL_3.geojson"
    st_read(u, quiet = TRUE) |> dplyr::filter(CNTR_CODE == "DE") |> st_transform(CRS_M)
  })
germany <- st_union(nuts3)
dens <- get_eurostat("demo_r_d3dens", time_format = "num") |>
  dplyr::filter(startsWith(geo, "DE"), TIME_PERIOD == YEAR) |>
  dplyr::transmute(NUTS_ID = geo, dens = values)

## (C) Instalaciones IED/E-PRTR emisoras de PM REAL -- EEA
ied <- "https://air.discomap.eea.europa.eu/arcgis/rest/services/Air/IED_SiteMap/MapServer/0/query"
fj <- request(ied) |> req_url_query(
  where = "countryCode='DE' AND pollutants LIKE '%PM%'",
  outFields = "x_4258,y_4258", returnGeometry = "false", f = "json",
  resultRecordCount = 2000) |> req_perform() |> resp_body_string() |> jsonlite::fromJSON()
facilities <- st_transform(st_as_sf(fj$features$attributes,
                                    coords = c("x_4258","y_4258"), crs = 4258), CRS_M)
cat("NUTS-3:", nrow(nuts3), "| densidad:", nrow(dens),
    "| instalaciones PM:", nrow(facilities), "\n")
## NUTS-3: 401 | densidad: 456 | instalaciones PM: 883
## (D) COVARIABLES REALES para el kriging con covariables (KED) y el modelo de intensidad:
##     elevacion (DEM satelital, Copernicus/SRTM) y % de suelo construido (uso del suelo, ~CORINE).
DE_EXT <- terra::ext(5, 16, 47, 56)
elev  <- geodata::elevation_30s(country = "DEU", path = tempdir())
## covariable de uso del suelo / urbanizacion (con respaldos, por robustez):
built <- tryCatch(
  terra::crop(geodata::landcover(var = "built", path = tempdir()), DE_EXT),
  error = function(e) tryCatch(
    terra::crop(geodata::footprint(2009, path = tempdir()), DE_EXT),
    error = function(e2) terra::crop(geodata::population(2020, res = 2.5, path = tempdir()), DE_EXT)))
extr <- function(sf_pts, r)
  terra::extract(r, terra::vect(st_transform(sf_pts, 4326)), ID = FALSE)[, 1]
stations$elev  <- extr(stations, elev)
stations$built <- extr(stations, built)

## Grilla de prediccion de 10 km, con las covariables extraidas
grid <- st_as_sf(st_make_grid(germany, cellsize = 10000, what = "centers"))
grid <- grid[lengths(st_intersects(grid, germany)) > 0, ]
grid$elev  <- extr(grid, elev)
grid$built <- extr(grid, built)

## Covariables como imagenes (spatstat) para modelar la intensidad (Obj.3)
rast_to_im <- function(r) {
  r <- terra::project(r, "EPSG:3035")
  v <- terra::as.matrix(r, wide = TRUE); v <- v[nrow(v):1, ]
  e <- as.vector(terra::ext(r))
  spatstat.geom::im(v, xrange = e[1:2], yrange = e[3:4])
}
im_elev  <- rast_to_im(elev)
im_built <- rast_to_im(built)
cat("Covariables reales listas (elevacion + % construido) en estaciones, grilla e imagenes\n")
## Covariables reales listas (elevacion + % construido) en estaciones, grilla e imagenes

La siguiente tabla resume los datos reales efectivamente descargados y su cobertura:

Tabla 1. Datos reales descargados y su cobertura.
Fuente N Detalle
Estaciones PM2.5 (UBA) 253 PM2.5 2022: 4.1-15.0 ug/m3 (media 9.6)
Instalaciones emisoras de PM (EEA) 883 IED/E-PRTR que reportan material particulado
Regiones NUTS-3 (Eurostat) 401 unidad areal; densidad demo_r_d3dens

2 Objetivo 1 — Variación continua (geoestadística)

Qué hace este bloque: estima el variograma empírico y ajusta el mejor modelo; luego interpola por kriging ordinario (OK, sin covariables) y por kriging con deriva externa (KED, con las covariables reales elevación y % construido) en una grilla de 10 km, y valida ambos con validación cruzada dejar-uno-fuera. Produce el variograma, los mapas de concentración (OK y KED) y el mapa de error.

## Obj.1 -- Variograma y kriging: ordinario (OK, sin covariables) y con deriva externa (KED, con covariables)
v_emp <- variogram(pm25 ~ 1, stations)
v_fit <- fit.variogram(v_emp, vgm(c("Exp","Sph","Mat")))
rango_km <- round(v_fit$range[nrow(v_fit)] / 1000, 1)
plot(v_emp, v_fit)                                        # Figura: variograma

## (a) Kriging ordinario (SIN covariables)
ok   <- krige(pm25 ~ 1, stations, grid, model = v_fit)
## [using ordinary kriging]
rmse <- round(sqrt(mean(krige.cv(pm25 ~ 1, stations, model = v_fit)$residual^2)), 2)

## (b) Kriging con deriva externa (CON covariables reales: elevacion + % construido)
s2 <- stations[is.finite(stations$elev) & is.finite(stations$built), ]
g2 <- grid[is.finite(grid$elev) & is.finite(grid$built), ]
v_ked <- fit.variogram(variogram(pm25 ~ elev + built, s2), vgm(c("Exp","Sph","Mat")))
ked   <- krige(pm25 ~ elev + built, s2, g2, model = v_ked)
## [using universal kriging]
rmse_ked <- round(sqrt(mean(krige.cv(pm25 ~ elev + built, s2, model = v_ked)$residual^2)), 2)

print(ggplot() + geom_sf(data = ok, aes(color = var1.pred), size = .8) +
  geom_sf(data = st_boundary(nuts3), color = "grey85", linewidth = .1) +
  scale_color_viridis_c(option = "magma", name = "PM2.5") +
  labs(title = "(a) Kriging ordinario (sin covariables)") + theme_minimal())

print(ggplot() + geom_sf(data = ked, aes(color = var1.pred), size = .8) +
  scale_color_viridis_c(option = "magma", name = "PM2.5") +
  labs(title = "(b) Kriging con covariables (KED): elevación + % construido") + theme_minimal())

print(ggplot() + geom_sf(data = ok, aes(color = sqrt(var1.var)), size = .8) +
  scale_color_viridis_c(name = "SE") + labs(title = "Error estándar (OK)") + theme_minimal())

cat("Rango =", rango_km, "km | RMSE OK =", rmse, "| RMSE KED =", rmse_ked, "ug/m3 |",
    "PM2.5 estaciones: min", round(min(stations$pm25),1), "max", round(max(stations$pm25),1), "\n")
## Rango = 44.2 km | RMSE OK = 1.3 | RMSE KED = 0.99 ug/m3 | PM2.5 estaciones: min 4.1 max 15
Tabla 2. Resumen geoestadístico (Obj.1): kriging ordinario vs. con covariables.
Métrica Valor
Modelo de variograma Mat
Rango (km) 44.2
LOO-CV RMSE — OK (ug/m3) 1.3
LOO-CV RMSE — KED con covariables (ug/m3) 0.99
Nº estaciones 253
PM2.5 medio (ug/m3) 9.6

Interpretación. El variograma alcanza la meseta hacia un rango de 44.2 km (correlación espacial del PM2.5). Se comparan dos interpolaciones: el kriging ordinario (sin covariables) y el kriging con deriva externa (KED), que incorpora la elevación y el % de suelo construido como covariables. La validación cruzada da RMSE = 1.3 µg/m³ (OK) frente a 0.99 µg/m³ (KED): las covariables mejoran la predicción. El mapa de error estándar señala dónde la predicción es más incierta (mín 4.1, máx 15 µg/m³ en estaciones).

3 Objetivo 2 — Variación areal (NUTS-3)

Qué hace este bloque: agrega la superficie de PM2.5 a las 400 regiones NUTS-3, calcula el conteo y la densidad de fuentes por región y une la densidad poblacional; luego mide la autocorrelación espacial (índice de Moran) y ajusta una regresión del PM2.5 frente a urbanización y fuentes.

## Obj.2 -- Agregacion a NUTS-3 y relacion con densidad y fuentes (datos reales)
agg <- aggregate(ok["var1.pred"], nuts3, FUN = mean)
areal <- nuts3
areal$pm25     <- agg$var1.pred
areal$n_fac    <- lengths(st_intersects(nuts3, facilities))
areal$area_km2 <- as.numeric(set_units(st_area(nuts3), "km^2"))
areal$fac_dens <- areal$n_fac / areal$area_km2 * 100
areal <- dplyr::left_join(areal, dens, by = "NUTS_ID")
ar <- areal[is.finite(areal$pm25) & is.finite(areal$dens), ]      # completos para modelar
nb <- poly2nb(ar); lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
MI <- moran.test(ar$pm25, lw, zero.policy = TRUE); moran_I <- round(MI$estimate[[1]], 3)
m  <- lm(pm25 ~ log1p(dens) + fac_dens, data = ar); R2 <- round(summary(m)$r.squared, 3)
print(ggplot(areal) + geom_sf(aes(fill = pm25), color = "white", linewidth = .05) +
  scale_fill_viridis_c(option = "magma", name = "PM2.5") +
  labs(title = "PM2.5 medio por NUTS-3") + theme_minimal())

print(ggplot(areal) + geom_sf(aes(fill = fac_dens), color = "white", linewidth = .05) +
  scale_fill_viridis_c(name = "fuentes/100km2") +
  labs(title = "Densidad de fuentes de PM por NUTS-3") + theme_minimal())

cat("Moran I =", moran_I, "| R2 =", R2,
    "| beta(log dens) =", round(coef(m)[2], 3), "| beta(fac_dens) =", round(coef(m)[3], 3), "\n")
## Moran I = 0.876 | R2 = 0.167 | beta(log dens) = 0.155 | beta(fac_dens) = 0.049
Tabla 3. Coeficientes de la regresión areal (PM2.5 ~ log densidad + densidad de fuentes).
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.525 0.145 58.749 0
log1p(dens) 0.155 0.026 6.023 0
fac_dens 0.049 0.011 4.607 0
Tabla 4. Regiones NUTS-3 con mayor PM2.5 estimado.
NUTS-3 PM2.5 Densidad (hab/km2) Fuentes/100km2
DEA12 11.6 2329.0 35.7
DEA17 11.5 2740.3 0.0
DEA32 11.4 2498.9 14.2
DEA31 11.3 1178.2 1.0
DEA13 11.3 2781.8 5.2
DEA14 11.3 1694.2 3.6
DEA16 11.3 1883.5 0.0
DEA51 11.1 2481.6 0.7

Interpretación. El índice de Moran = 0.876 indica la autocorrelación espacial del PM2.5 por región. En la regresión, la urbanización (densidad) y la densidad de fuentes se asocian al PM2.5 (R² = 0.167; β densidad = 0.155, β fuentes = 0.049).

4 Objetivo 3 — Procesos puntuales (fuentes de emisión)

Qué hace este bloque: trata las instalaciones emisoras de PM como un patrón de puntos, estima su intensidad (suavizado kernel), contrasta si el patrón es aleatorio o agregado (función L de Ripley con envolvente + test de Clark-Evans) y ajusta un modelo de intensidad en función de covariables (ppm con elevación y % construido).

## Obj.3 -- Proceso puntual de las fuentes de PM (datos reales EEA)
bb <- st_bbox(facilities)
W  <- owin(c(bb[["xmin"]], bb[["xmax"]]), c(bb[["ymin"]], bb[["ymax"]]))
xy <- st_coordinates(facilities)
pp <- unique(ppp(xy[,1], xy[,2], window = W))
lam <- density(pp, sigma = bw.diggle(pp))
plot(lam, main = "Intensidad kernel de fuentes de PM"); plot(pp, add = TRUE, pch = ".")

Lenv <- envelope(pp, Lest, nsim = 39, verbose = FALSE); plot(Lenv, . - r ~ r, main = "Función L de Ripley")

ce  <- clarkevans.test(pp, correction = "none", alternative = "clustered")
Rce <- round(as.numeric(ce$statistic), 3)
## Intensidad en funcion de COVARIABLES reales (elevacion + % construido) -- Poisson no homogeneo
fit_ppm <- ppm(pp ~ elev + built, covariates = list(elev = im_elev, built = im_built))
print(round(coef(fit_ppm), 4))                        # efecto de cada covariable sobre log-intensidad
## (Intercept)        elev       built 
##    -21.8015     -0.0030      5.5249
cat("Fuentes:", pp$n, "| Clark-Evans R =", Rce, "(R<1 => agregado) | p =", signif(ce$p.value, 3), "\n")
## Fuentes: 175 | Clark-Evans R = 0.384 (R<1 => agregado) | p = 4.49e-55
Tabla 5. Resumen del proceso puntual (Obj.3).
Métrica Valor
Nº de fuentes (patrón de puntos) 175
Índice de Clark-Evans R 0.384
p-valor 4.49e-55
Diagnóstico Agregado (clustered)

Interpretación. El índice de Clark-Evans R = 0.384 (< 1) indica que las 175 fuentes de PM están agregadas (no al azar); la función L por encima de la envolvente lo confirma. El modelo ppm describe la intensidad en función de covariables reales (elevación y % de suelo construido): sus coeficientes cuantifican cómo cambia la densidad de fuentes con el terreno y la urbanización.

5 Objetivos 4 y 5 — Integración y pregunta final

Qué hace este bloque: lleva la intensidad de fuentes (Obj.3) a las estaciones, mide su correlación con la concentración (Obj.1) y, además, ejecuta un kriging con deriva externa (KED) usando la intensidad de fuentes como covariable — la integración formal que pide el Obj.4. Luego descompone la varianza areal del PM2.5 entre urbanización, industria y un residual.

## Obj.4 y 5 -- Integracion (intensidad de fuentes vs concentracion) y descomposicion
sxy <- st_coordinates(stations)
stations$src <- interp.im(lam, sxy[,1], sxy[,2])
corr <- round(cor(stations$pm25, stations$src, use = "complete.obs"), 3)
print(ggplot(stations, aes(src, pm25)) + geom_point(color = "#B22222") +
  geom_smooth(method = "lm", se = TRUE, color = "black") +
  labs(x = "Intensidad de fuentes (Obj.3)", y = "PM2.5 en estación (Obj.1)",
       title = paste0("Concentración vs intensidad de fuentes (r = ", corr, ")")) + theme_minimal())

## KED formal: kriging usando la INTENSIDAD DE FUENTES como covariable (integracion Obj.4)
grid$src <- interp.im(lam, st_coordinates(grid)[,1], st_coordinates(grid)[,2])
s3 <- stations[is.finite(stations$src), ]; g3 <- grid[is.finite(grid$src), ]
v_src   <- fit.variogram(variogram(pm25 ~ src, s3), vgm(c("Exp","Sph")))
ked_src <- krige(pm25 ~ src, s3, g3, model = v_src)
## [using universal kriging]
rmse_src <- round(sqrt(mean(krige.cv(pm25 ~ src, s3, model = v_src)$residual^2)), 2)
print(ggplot() + geom_sf(data = ked_src, aes(color = var1.pred), size = .8) +
  scale_color_viridis_c(option = "magma", name = "PM2.5") +
  labs(title = "PM2.5 por KED con la intensidad de fuentes como covariable (Obj.4)") + theme_minimal())

cat("Obj.4 KED con intensidad de fuentes | RMSE =", rmse_src, "ug/m3\n")
## Obj.4 KED con intensidad de fuentes | RMSE = 1.23 ug/m3
av <- anova(lm(pm25 ~ log1p(dens) + fac_dens, data = ar))
contrib <- round(100 * av$"Sum Sq" / sum(av$"Sum Sq"), 1)
names(contrib) <- c("urbanizacion","industria","residual")
print(contrib)
## urbanizacion    industria     residual 
##         12.1          4.6         83.3
Tabla 6. Descomposición de la varianza areal del PM2.5 (Obj.5).
Factor Contribución a la varianza (%)
Urbanización (densidad poblacional) 12.1
Industria (densidad de fuentes) 4.6
Residual (uso del suelo, meteorología, transporte) 83.3

Interpretación. La correlación concentración–intensidad de fuentes es r = 0.295, y el kriging con la intensidad de fuentes como covariable (KED, Obj.4) produce un mapa guiado por las fuentes con RMSE = 1.23 µg/m³. La descomposición reparte la varianza areal del PM2.5 entre urbanización (12.1%), industria (4.6%) y un residual (83.3%).

Respuesta a la pregunta final. Con datos reales, el nivel regional de PM2.5 en Alemania se explica por una combinación de urbanización e industria (más el residual de uso del suelo, meteorología y transporte). La geoestadística aporta el campo continuo y su incertidumbre; la rama areal, la lectura por región controlando la autocorrelación (Moran = 0.876); los procesos puntuales, el mecanismo y la escala del agrupamiento de las fuentes (Clark-Evans = 0.384); y su integración, el vínculo directo fuentes→concentración (r = 0.295).

6 Conclusión

Este trabajo integró, con datos reales de Alemania (año 2022), las tres ramas de la estadística espacial y respondió a la pregunta del proyecto:

  • Objetivo 1 (geoestadística). A partir de 253 estaciones reales de PM2.5, el variograma mostró correlación espacial hasta un rango de 44.2 km; se interpoló por kriging ordinario (RMSE 1.3) y por kriging con covariables (elevación y % construido; RMSE 0.99). Aporta el campo de exposición y su incertidumbre donde no hay estaciones.
  • Objetivo 2 (areal). Agregando a 401 regiones NUTS-3, el PM2.5 mostró autocorrelación espacial (Moran I = 0.876) y se asoció positivamente a la urbanización y a la densidad de fuentes (R² = 0.167). Aporta la lectura por unidad de gestión, controlando la dependencia espacial.
  • Objetivo 3 (procesos puntuales). Las 175 instalaciones emisoras de PM forman un patrón agregado (Clark-Evans R = 0.384, p = 4.49^{-55}). Aporta el mecanismo y la escala del agrupamiento de las fuentes.
  • Objetivo 4 (integración). La intensidad de fuentes se relaciona con la concentración (r = 0.295) y se incorporó como covariable en un kriging con deriva externa (RMSE 1.23), conectando formalmente el mecanismo (dónde están las fuentes) con el efecto (cuánto PM2.5).
  • Objetivo 5 (síntesis). La varianza regional del PM2.5 se reparte entre urbanización (12.1%), industria (4.6%) y un residual (83.3%). La respuesta a la pregunta final es que el nivel regional se explica por una combinación de urbanización e industria, modulada por el territorio, y que cada rama aporta evidencia que las otras dos no pueden dar por sí solas.

Todo el análisis es reproducible: descarga los datos reales (UBA, EEA, Eurostat) y ejecuta el mismo código de R/ y de los chunks del .Rmd.

7 Referencias

Baddeley, Rubak & Turner (2015); Banerjee, Carlin & Gelfand (2014); Diggle & Ribeiro (2007); Diggle (2013); Gelfand et al. (2010); Rue & Held (2005). Datos: EEA (IED/E-PRTR), Eurostat (NUTS 2021, demo_r_d3dens), Umweltbundesamt (air_data v3, PM2.5). Ver references.bib y docs/DATA_SOURCES.md.