---
title: "Diagnóstico del cambio de cobertura vegetal en Ixil (Yucatán, México)"
subtitle: "Serie de NDVI Landsat y Sentinel-2 (2001–2025): matriz de transición observada, cicatrices de proyectos eólicos y propiedad ejidal"
author: "Carla Machicado Sanchez"
date: today
lang: es
jupyter: python3
engine: knitr # el código es R: sin esta línea, "jupyter: python3" impediría ejecutarlo
bibliography: referencias_ixil.bib
csl: chicago-author-date.csl
link-citations: true
format:
html:
theme: cyborg
backgroundcolor: "#000000"
toc: true
toc-depth: 3
number-sections: true
code-fold: true
code-tools: true
embed-resources: true
fig-width: 9
fig-height: 6
execute:
warning: false
message: false
cache: false
---
```{=html}
<style>
/* Leyendas y control de capas de los mapas, a juego con el fondo negro */
.leaflet-container .info, .leaflet-container .leaflet-control-layers {
background: rgba(0, 0, 0, 0.8); color: #f0f0f0;
}
.leaflet-container .legend { color: #f0f0f0; }
.leaflet-container .leaflet-control-layers-toggle { filter: invert(1); }
</style>
```
# Preparación del entorno
```{r}
#| label: paquetes
paquetes <- c("terra", "sf", "dplyr", "tidyr", "ggplot2", "knitr", "scales", "leaflet")
faltantes <- paquetes[!paquetes %in% rownames(installed.packages())]
if (length(faltantes) > 0) install.packages(faltantes)
invisible(lapply(paquetes, library, character.only = TRUE))
terraOptions(progress = 0)
fmt <- function(x, d = 0) format(round(x, d), big.mark = ",", nsmall = d, trim = TRUE)
pct <- function(x, d = 1) paste0(fmt(100 * x, d), " %")
# Tema oscuro para las gráficas, a juego con el fondo negro del reporte
tema_oscuro <- function() {
theme_minimal() +
theme(plot.background = element_rect(fill = "#000000", colour = NA),
panel.background = element_rect(fill = "#000000", colour = NA),
legend.background = element_rect(fill = "#000000", colour = NA),
legend.key = element_rect(fill = "#000000", colour = NA),
panel.grid = element_line(colour = "#2b2b2b"),
text = element_text(colour = "#e6e6e6"),
axis.text = element_text(colour = "#cfcfcf"),
strip.text = element_text(colour = "#e6e6e6"))
}
theme_set(tema_oscuro())
```
El procesamiento espacial se hace con `terra` [@hijmans_terra] y `sf` [@pebesma2018], y los mapas interactivos con `leaflet` [@cheng_leaflet].
## Parámetros del análisis
```{r}
#| label: parametros
# Límite municipal de Ixil
RUTA_SHP_IXIL <- "G:/FLACSO/SEMILLERO_ PAPER EJIDOS MEXICO/DATOS ESPACIALES/ixil_muni.shp"
MUNICIPIO <- "Ixil"
# En este paso las imagenes ya tienen la banda calculada de ndvi
DIR_IMAGENES <- "G:/flacso/modulo_5/modelacion_cienciadatos_prospecterr/trabajo final/img"
imagenes <- tribble(
~anio, ~sensor, ~archivo, ~banda_nir, ~banda_roja,
2001, "Landsat 5/7", "ndvi_landsat30m_en_grilla10m_ixil_2001.tif", NA, NA,
2003, "Landsat 5/7", "ndvi_landsat30m_en_grilla10m_ixil_2003.tif", NA, NA,
2004, "Landsat 5/7", "ndvi_landsat30m_en_grilla10m_ixil_2004.tif", NA, NA,
2005, "Landsat 5/7", "ndvi_landsat30m_en_grilla10m_ixil_2005.tif", NA, NA,
2007, "Landsat 5/7", "ndvi_landsat30m_en_grilla10m_ixil_2007.tif", NA, NA,
2008, "Landsat 5/7", "ndvi_landsat30m_en_grilla10m_ixil_2008.tif", NA, NA,
2010, "Landsat 5/7", "ndvi_landsat30m_en_grilla10m_ixil_2010.tif", NA, NA,
2015, "Landsat 8", "ndvi_landsat30m_en_grilla10m_ixil_2015.tif", NA, NA,
2020, "Sentinel-2", "ndvi_sentinel2_10m_ixil_2020.tif", NA, NA,
2025, "Sentinel-2", "ndvi_sentinel2_10m_ixil_2025.tif", NA, NA
)
SIN_DATO <- -32768 # valor "sin dato" de los archivos exportados
# Resolución del análisis (m). Landsat tiene 30 m reales, así que toda la serie se
# analiza a 30 m; Sentinel-2 se agrega a esa grilla. En el mapa interactivo,
# Sentinel-2 se muestra a sus 10 m originales.
RES_ANALISIS <- 30
# Normalización entre fechas (ver la sección correspondiente). El año de
# referencia es la primera fecha disponible salvo que se indique otro.
NORMALIZAR <- TRUE
ANIO_REFERENCIA <- NULL
# Clases de NDVI para la matriz de transición
CORTES_NDVI <- c(-1, 0.2, 0.4, 0.6, 1)
CLASES_NDVI <- c("Sin vegetación", "Vegetación escasa", "Vegetación media", "Vegetación densa")
paleta <- setNames(c("#a50026", "#fdae61", "#a6d96a", "#1a9850"), CLASES_NDVI)
# Definición de cicatriz candidata (sobre el NDVI normalizado)
UMBRAL_VEGETADO <- 0.50 # NDVI mínimo ANTES para considerar que había vegetación
UMBRAL_DESMONTE <- 0.35 # NDVI máximo DESPUÉS, en todas las fechas posteriores
CAIDA_MINIMA <- 0.20 # caída mínima entre antes y después
MIN_HA_PARCHE <- 0.27 # parches menores se descartan (0.27 ha = 3 píxeles de 30 m)
# Validación con el índice Kappa: tres fechas separadas por el mismo intervalo
# (calibración inicial, calibración final y fecha a predecir). NULL = se eligen
# solas: la terna que termina en la fecha más reciente y, entre esas, la de
# intervalo más largo. Para fijarlas: c(2005, 2015, 2025).
ANIOS_VALIDACION <- NULL
# Capas opcionales. Si la ruta no existe, la sección se omite y el resto corre igual.
# - Tenencia: polígonos de tierras ejidales (aquí, tierras de uso común del RAN).
# - Infraestructura eólica: puntos, líneas o polígonos, en una o varias rutas.
RUTA_EJIDOS <- "G:/FLACSO/SEMILLERO_ PAPER EJIDOS MEXICO/DATOS ESPACIALES/Acaparamiento ilegal (empresarios) cut.shp"
CAMPO_EJIDO <- NULL # campo con el nombre de cada polígono (NULL = numerar)
ETIQUETA_EJIDOS <- "Tierras ejidales (capa cargada)"
RUTA_EOLICO <- c("G:/FLACSO/SEMILLERO_ PAPER EJIDOS MEXICO/DATOS ESPACIALES/Aerogeneradores1.shp",
"G:/flacso/SEMILLERO_ PAPER EJIDOS MEXICO/SHAPEFILES-20260609T232032Z-3-001/SHAPEFILES/GEOCOMUNES/_3_1_3_parque_eolico_peninsula_poli_2020 (1)/_3_1_3_parque_eolico_peninsula_poli_2020.shp")
BUFFER_EOLICO_M <- 500 # área de influencia alrededor de la infraestructura (m)
# Línea de transmisión eléctrica (opcional): se mide si las cicatrices se le sobreponen
RUTA_LINEA <- "G:/FLACSO/SEMILLERO_ PAPER EJIDOS MEXICO/DATOS ESPACIALES/Línea de transmisión eléctrica.shp"
BUFFER_LINEA_M <- 100 # corredor a cada lado de la línea (m)
BANDAS_LINEA_M <- c(0, 100, 250, 500, 1000, 3000, Inf) # bandas de distancia a la línea (m)
# Costa. La franja costera es la parte del municipio a ANCHO_COSTA_M o menos de la
# línea de costa. RUTA_COSTA = NULL toma la costa del propio límite municipal (su
# borde norte, que da al mar); para usar otra, indica aquí un shapefile de líneas.
RUTA_COSTA <- NULL
ANCHO_COSTA_M <- 4000
BANDAS_COSTA_M <- c(0, 1000, 2000, 3000, 4000, 5000, Inf) # bandas de distancia a la costa (m)
# Proyección de Markov y escenarios. Cada escenario multiplica por un factor la
# probabilidad de que un píxel pase a cicatriz: 1 = la tendencia observada sigue
# igual; 2 = la pérdida permanente se duplica; 0.5 = se reduce a la mitad.
ESCENARIOS <- c("Tendencial" = 1, "Disruptivo" = 2, "De política" = 0.5, "De política, fuerte" = 0)
PASOS_PROYECCION <- 2 # periodos hacia adelante
SENSIBILIDAD_K <- c(0, 0.25, 0.5, 0.75, 1, 1.5, 2, 3) # factores para el análisis de sensibilidad
CRS_METRICO <- "EPSG:32616" # UTM zona 16 Norte (Yucatán)
dir_salidas <- "salidas_ixil_ndvi"
dir.create(dir_salidas, showWarnings = FALSE, recursive = TRUE)
```
# Límite municipal de Ixil
```{r}
#| label: limite
if (!file.exists(RUTA_SHP_IXIL)) {
stop("No se encontró el shapefile de Ixil en '", RUTA_SHP_IXIL, "'.")
}
limite <- st_read(RUTA_SHP_IXIL, quiet = TRUE) |> st_make_valid()
if (is.na(st_crs(limite))) stop("El shapefile no tiene sistema de referencia (falta el .prj).")
limite_m <- st_sf(municipio = MUNICIPIO,
geometry = st_union(st_geometry(st_transform(limite, CRS_METRICO))))
limite_v <- vect(limite_m)
area_limite_ha <- as.numeric(st_area(limite_m)) / 1e4
bb <- round(st_bbox(st_transform(limite_m, 4326)), 3)
cat("Superficie del polígono:", fmt(area_limite_ha), "ha\n",
"Extensión (lon/lat):", paste(names(bb), bb, collapse = " "))
```
# Datos: imágenes Landsat y Sentinel-2
- **2001 a 2015: Landsat**, con 30 m de resolución [@wulder2019]. Las nubes y sus sombras se eliminaron con la banda de calidad derivada del algoritmo Fmask [@zhu2012].
- **2020 y 2025: Sentinel-2**, con 10 m de resolución [@drusch2012]. Las nubes se eliminaron con Cloud Score+ [@pasquarella2023].
```{r}
#| label: lectura-ndvi
if (!dir.exists(DIR_IMAGENES)) {
stop("La carpeta de imágenes no existe: '", DIR_IMAGENES, "'. Revisa que la ruta esté bien escrita.")
}
imagenes <- imagenes |>
mutate(ruta = ifelse(file.exists(archivo), archivo, file.path(DIR_IMAGENES, archivo)),
existe = file.exists(ruta)) |>
arrange(anio)
faltan <- imagenes[!imagenes$existe, ]
disp <- imagenes[imagenes$existe, ]
if (nrow(disp) < 2) {
stop("Se necesitan al menos dos imágenes. Revisa DIR_IMAGENES y la tabla 'imagenes'.")
}
anios <- disp$anio
anios_chr <- as.character(anios)
# Lee un archivo y devuelve el NDVI (-1 a 1) recortado al municipio, en CRS métrico
leer_ndvi <- function(ruta, nir, rojo) {
r <- rast(ruta)
if (nlyr(r) == 1) {
x <- subst(r, SIN_DATO, NA)
# Los archivos exportados guardan NDVI x 10000 como entero
if (global(abs(x), "max", na.rm = TRUE)[1, 1] > 1.5) x <- x / 10000
} else {
if (is.na(nir) || is.na(rojo)) {
stop("'", basename(ruta), "' tiene ", nlyr(r), " bandas: indica banda_nir y banda_roja.")
}
x <- (r[[nir]] - r[[rojo]]) / (r[[nir]] + r[[rojo]])
}
if (!same.crs(x, CRS_METRICO)) x <- project(x, CRS_METRICO, method = "bilinear")
x <- crop(x, limite_v, snap = "out")
names(x) <- "ndvi"
x
}
nativos <- lapply(seq_len(nrow(disp)), function(i) {
leer_ndvi(disp$ruta[i], disp$banda_nir[i], disp$banda_roja[i])
})
names(nativos) <- anios_chr
disp$pixel_m <- vapply(nativos, function(x) res(x)[1], 0)
# Grilla común de análisis: se agrega cada imagen a RES_ANALISIS con la media
agregar <- function(x) {
f <- max(1, round(RES_ANALISIS / res(x)[1]))
if (f > 1) aggregate(x, fact = f, fun = "mean", na.rm = TRUE) else x
}
plantilla <- agregar(nativos[[1]])
a_grilla <- function(x) {
x <- agregar(x)
if (!compareGeom(x, plantilla, stopOnError = FALSE)) x <- resample(x, plantilla, method = "bilinear")
mask(x, limite_v)
}
ndvi_bruto <- rast(lapply(nativos, a_grilla))
names(ndvi_bruto) <- anios_chr
ha_pixel <- prod(res(plantilla)) / 1e4
en_mun <- !is.na(values(rasterize(limite_v, plantilla))[, 1])
```
```{r}
#| label: aviso-faltantes
#| echo: false
#| results: asis
if (nrow(faltan) > 0) {
cat("::: {.callout-warning}\n## Faltan imágenes\nNo se encontraron los archivos de: **",
paste(faltan$anio, collapse = ", "), "**. El análisis se hizo con las ", nrow(disp),
" fechas disponibles (", paste(anios, collapse = ", "),
"). Cuando descargues las demás a `", DIR_IMAGENES,
"`, vuelve a generar el reporte y entrarán solas.\n:::\n", sep = "")
}
```
```{r}
#| label: tbl-imagenes
#| tbl-cap: "Imágenes usadas y resumen del NDVI dentro del municipio, antes de normalizar."
M_bruto <- values(ndvi_bruto)
resumen_img <- data.frame(
anio = anios,
sensor = disp$sensor,
pixel = disp$pixel_m,
p10 = apply(M_bruto[en_mun, , drop = FALSE], 2, quantile, 0.10, na.rm = TRUE),
mediana = apply(M_bruto[en_mun, , drop = FALSE], 2, median, na.rm = TRUE),
p90 = apply(M_bruto[en_mun, , drop = FALSE], 2, quantile, 0.90, na.rm = TRUE),
sin_dato = 100 * colMeans(is.na(M_bruto[en_mun, , drop = FALSE]))
)
resumen_img |>
mutate(across(c(p10, mediana, p90), ~ round(.x, 3)), sin_dato = round(sin_dato, 2)) |>
setNames(c("Año", "Sensor", "Píxel del archivo (m)", "NDVI percentil 10", "NDVI mediana",
"NDVI percentil 90", "% sin dato")) |>
kable(row.names = FALSE)
```
::: {.callout-note}
## Resolución real
Los archivos de Landsat vienen en una grilla de 10 m para poder superponerlos con Sentinel-2, pero su detalle real es de 30 m. Por eso toda la serie se analiza a `r RES_ANALISIS` m. La resolución de 10 m de Sentinel-2 se conserva en el mapa interactivo.
:::
## Normalización entre fechas
El NDVI de dos fechas no es directamente comparable, aunque sean de la misma temporada, por tres razones: un año puede ser más seco que otro; los sensores de Landsat 5, 7 y 8 no miden exactamente lo mismo [@roy2016]; y Sentinel-2 se usa aquí sin corrección atmosférica, lo que baja su NDVI respecto a Landsat. La @tbl-imagenes muestra cuánto se mueve la mediana de un año a otro por estas causas.
Para que esas diferencias generales no se lean como cambio de cobertura, cada fecha se ajusta linealmente a la fecha de referencia, de modo que compartan mediana y rango intercuartílico dentro del municipio. La consecuencia es doble: las cicatrices, que son cambios **locales**, se detectan mejor; pero un cambio que afectara por igual a **todo** el municipio quedaría oculto.
```{r}
#| label: normalizacion
j_ref <- if (is.null(ANIO_REFERENCIA)) 1L else match(ANIO_REFERENCIA, anios)
if (is.na(j_ref)) stop("ANIO_REFERENCIA no está entre las fechas disponibles.")
iqr <- apply(M_bruto[en_mun, , drop = FALSE], 2, IQR, na.rm = TRUE)
M <- M_bruto
if (NORMALIZAR) {
for (j in seq_along(anios)) {
M[, j] <- (M_bruto[, j] - resumen_img$mediana[j]) * (iqr[j_ref] / iqr[j]) + resumen_img$mediana[j_ref]
}
M[] <- pmin(pmax(M, -1), 1)
}
ndvi <- setValues(ndvi_bruto, M)
```
## Las imágenes sobre Google Earth
```{r}
#| label: capas-opcionales
leer_opcional <- function(rutas, solo_geometria = FALSE) {
rutas <- rutas[file.exists(rutas)]
if (length(rutas) == 0) return(NULL)
capas <- lapply(rutas, function(ruta) {
x <- st_read(ruta, quiet = TRUE) |> st_make_valid() |> st_transform(CRS_METRICO)
if (solo_geometria) st_sf(archivo = rep(basename(ruta), nrow(x)), geometry = st_geometry(x)) else x
})
x <- do.call(rbind, capas)
if (nrow(x) == 0) NULL else x
}
# --- Tenencia ejidal ------------------------------------------------------
ejidos <- leer_opcional(RUTA_EJIDOS)
if (!is.null(ejidos)) {
ejidos <- suppressWarnings(st_intersection(ejidos, st_geometry(limite_m)))
ejidos <- ejidos[st_geometry_type(ejidos) %in% c("POLYGON", "MULTIPOLYGON"), ]
if (nrow(ejidos) == 0) ejidos <- NULL
}
hay_ejidos <- !is.null(ejidos)
if (hay_ejidos) {
ejidos$ejido <- if (!is.null(CAMPO_EJIDO) && CAMPO_EJIDO %in% names(ejidos)) {
as.character(ejidos[[CAMPO_EJIDO]])
} else paste("Polígono", seq_len(nrow(ejidos)))
ejidos$id_ej <- seq_len(nrow(ejidos))
}
# --- Infraestructura eólica -----------------------------------------------
# No se recorta al municipio: un aerogenerador vecino también proyecta su área
# de influencia hacia dentro de Ixil.
eolico <- leer_opcional(RUTA_EOLICO, solo_geometria = TRUE)
hay_eolico <- !is.null(eolico)
if (hay_eolico) {
zona_eolica <- st_sf(dentro = 1L, geometry = st_union(st_buffer(eolico, BUFFER_EOLICO_M)))
# Para el mapa, cada tipo de geometría se dibuja por separado
tipo_geom <- as.character(st_geometry_type(eolico))
de_tipo <- function(tipos, destino) {
x <- st_geometry(eolico)[tipo_geom %in% tipos]
if (length(x) == 0) NULL else st_transform(suppressWarnings(st_cast(x, destino)), 4326)
}
eol_puntos <- de_tipo(c("POINT", "MULTIPOINT"), "POINT")
eol_lineas <- de_tipo(c("LINESTRING", "MULTILINESTRING"), "MULTILINESTRING")
eol_poligonos <- de_tipo(c("POLYGON", "MULTIPOLYGON"), "MULTIPOLYGON")
}
# --- Línea de transmisión eléctrica ---------------------------------------
linea <- leer_opcional(RUTA_LINEA)
if (!is.null(linea)) {
linea <- linea[st_geometry_type(linea) %in% c("LINESTRING", "MULTILINESTRING"), ]
if (nrow(linea) == 0) linea <- NULL
}
hay_linea <- !is.null(linea)
if (hay_linea) {
linea_union <- st_union(st_geometry(linea))
zona_linea <- st_sf(dentro = 1L, geometry = st_buffer(linea_union, BUFFER_LINEA_M))
}
# --- Línea de costa -------------------------------------------------------
# Sin un shapefile propio, la costa se toma del límite municipal: los tramos de
# su borde cuya cara exterior mira al norte (hacia el mar), en la parte norte del
# municipio. Conviene comprobar en los mapas que la línea celeste sigue la playa.
costa_del_limite <- function(poligono, tolerancia = 4000) {
pol <- st_cast(st_geometry(poligono), "POLYGON")
pol <- pol[which.max(st_area(pol))]
co <- st_coordinates(pol); co <- co[co[, "L1"] == 1, c("X", "Y")] # anillo exterior
n <- nrow(co); a <- co[-n, , drop = FALSE]; b <- co[-1, , drop = FALSE]
dx <- b[, 1] - a[, 1]; dy <- b[, 2] - a[, 2]; largo <- sqrt(dx^2 + dy^2)
sentido <- if (sum(a[, 1] * b[, 2] - b[, 1] * a[, 2]) > 0) 1 else -1 # +1: vértices en sentido antihorario
normal_norte <- -sentido * dx / pmax(largo, 1e-9) # cuánto mira al norte la cara exterior
es_costa <- largo > 0 & normal_norte > 0.5 & (a[, 2] + b[, 2]) / 2 > max(co[, 2]) - tolerancia
st_sfc(st_multilinestring(lapply(which(es_costa), function(i) rbind(a[i, ], b[i, ]))),
crs = st_crs(poligono))
}
costa <- if (!is.null(RUTA_COSTA) && file.exists(RUTA_COSTA)) {
x <- st_read(RUTA_COSTA, quiet = TRUE) |> st_make_valid() |> st_transform(CRS_METRICO)
suppressWarnings(st_intersection(st_union(st_geometry(x)), st_buffer(st_geometry(limite_m), 2000)))
} else costa_del_limite(limite_m)
fuente_costa <- if (!is.null(RUTA_COSTA) && file.exists(RUTA_COSTA)) basename(RUTA_COSTA) else "borde norte del límite municipal"
km_costa <- as.numeric(sum(st_length(costa))) / 1000
# Mapa base: imagen satelital de Google, con el límite y las capas opcionales
mapa_base <- function() {
m <- leaflet() |>
addTiles(urlTemplate = "https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}",
attribution = "Imágenes © Google", group = "Google Satélite",
options = tileOptions(maxZoom = 20)) |>
addTiles(urlTemplate = "https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}",
attribution = "Imágenes © Google", group = "Google Híbrido",
options = tileOptions(maxZoom = 20)) |>
addProviderTiles(providers$Esri.WorldImagery, group = "Esri World Imagery") |>
addProviderTiles(providers$CartoDB.DarkMatter, group = "Mapa oscuro") |>
addPolygons(data = st_transform(limite_m, 4326), fill = FALSE, color = "#ffff00",
weight = 2, group = "Límite municipal") |>
# Todos los mapas abren encuadrados en el municipio
fitBounds(bb[["xmin"]], bb[["ymin"]], bb[["xmax"]], bb[["ymax"]])
if (hay_ejidos) {
m <- addPolygons(m, data = st_transform(ejidos, 4326), fill = FALSE, color = "#00e5ff",
weight = 2, dashArray = "6", popup = ~ejido, group = ETIQUETA_EJIDOS)
}
if (hay_eolico) {
if (!is.null(eol_poligonos)) {
m <- addPolygons(m, data = eol_poligonos, fill = FALSE, color = "#ffffff", weight = 2,
group = "Infraestructura eólica")
}
if (!is.null(eol_lineas)) {
m <- addPolylines(m, data = eol_lineas, color = "#ffffff", weight = 2,
group = "Infraestructura eólica")
}
if (!is.null(eol_puntos)) {
m <- addCircleMarkers(m, data = eol_puntos, radius = 4, color = "#ffffff",
fillColor = "#2b83ba", fillOpacity = 1, weight = 1,
group = "Infraestructura eólica")
}
}
if (hay_linea) {
m <- addPolylines(m, data = st_transform(st_cast(st_geometry(linea), "MULTILINESTRING"), 4326),
color = "#ff9800", weight = 3, dashArray = "8", opacity = 1,
group = "Línea de transmisión")
}
m <- addPolylines(m, data = st_transform(costa, 4326), color = "#00e5ff", weight = 3, opacity = 1,
group = "Línea de costa")
m
}
grupos_fijos <- c("Límite municipal", if (hay_ejidos) ETIQUETA_EJIDOS,
if (hay_eolico) "Infraestructura eólica",
if (hay_linea) "Línea de transmisión", "Línea de costa")
grupos_base <- c("Google Satélite", "Google Híbrido", "Esri World Imagery", "Mapa oscuro")
```
```{r}
#| label: fig-mapa-ndvi
#| fig-cap: "NDVI de cada fecha sobre la imagen satelital de Google. Con el control de la derecha se enciende y apaga cada año; al apagar todos queda visible la imagen de fondo."
pal_ndvi <- colorNumeric("RdYlGn", domain = c(0, 1), na.color = "transparent")
m <- mapa_base()
for (a in anios_chr) {
# Sentinel-2 se muestra a su resolución original; Landsat, a la de análisis
capa <- if (grepl("Sentinel", disp$sensor[anios_chr == a])) mask(nativos[[a]], limite_v) else ndvi_bruto[[a]]
m <- addRasterImage(m, clamp(capa, 0, 1), colors = pal_ndvi, opacity = 0.85,
group = paste("NDVI", a), maxBytes = 16 * 1024^2)
}
m |>
addLegend(pal = pal_ndvi, values = c(0, 1), title = "NDVI", position = "bottomleft") |>
addLayersControl(baseGroups = grupos_base,
overlayGroups = c(paste("NDVI", anios_chr), grupos_fijos),
options = layersControlOptions(collapsed = FALSE)) |>
hideGroup(paste("NDVI", head(anios_chr, -1)))
```
# Cruce de las fechas
Para construir una matriz de transición, el NDVI normalizado se agrupa en cuatro clases: `r paste0(CLASES_NDVI, " (", head(CORTES_NDVI, -1), " a ", tail(CORTES_NDVI, -1), ")", collapse = "; ")`. La matriz compara la primera fecha con la última; las fechas intermedias se usan después para seguir la trayectoria y fechar las cicatrices.
```{r}
#| label: cruce
ANIO_T1 <- anios[1]
ANIO_T2 <- anios[length(anios)]
clasificar <- function(x) cut(x, CORTES_NDVI, labels = CLASES_NDVI, include.lowest = TRUE, right = FALSE)
# Píxeles del análisis: dentro del municipio y con dato en la primera y la última fecha
validos <- en_mun & !is.na(M[, 1]) & !is.na(M[, ncol(M)])
px <- data.frame(celda = which(validos), ha = ha_pixel)
N <- M[px$celda, , drop = FALSE] # NDVI normalizado: píxeles x fechas
px$g1 <- clasificar(N[, 1])
px$g2 <- clasificar(N[, ncol(N)])
px <- cbind(px, xyFromCell(plantilla, px$celda))
niveles <- CLASES_NDVI
# Auxiliar: vector de valores por píxel -> raster sobre la grilla de análisis
a_raster <- function(valores, nombre) {
r <- rast(plantilla); v <- rep(NA_real_, ncell(r)); v[px$celda] <- as.numeric(valores)
values(r) <- v; names(r) <- nombre
r
}
```
# Matriz de transición observada
## En hectáreas
Las filas son la clase de origen (`r ANIO_T1`) y las columnas la de destino (`r ANIO_T2`). La diagonal es lo que **persistió**; todo lo demás es cambio.
```{r}
#| label: matriz
M_ha <- xtabs(ha ~ g1 + g2, data = px)
names(dimnames(M_ha)) <- c("desde", "hacia")
# Probabilidades de transición: cada fila dividida por su total
tot_fila <- rowSums(M_ha)
P <- M_ha / tot_fila
P[tot_fila == 0, ] <- NA # clase ausente en T1: su fila no está definida
# Verificación: toda fila con superficie en T1 debe sumar 1
stopifnot(all(abs(rowSums(P)[tot_fila > 0] - 1) < 1e-8))
```
```{r}
#| label: tbl-matriz-ha
#| tbl-cap: "Matriz de transición observada (ha). Filas: origen; columnas: destino."
kable(round(addmargins(M_ha), 1), format.args = list(big.mark = ","))
```
## En probabilidades
```{r}
#| label: tbl-matriz-p
#| tbl-cap: "Matriz de transición observada (% de la superficie de origen). Cada fila suma 100."
kable(round(100 * P, 2))
```
```{r}
#| label: fig-matriz
#| fig-cap: "Matriz de transición: el color indica qué porcentaje de la clase de origen terminó en cada destino; el texto añade las hectáreas."
#| fig-height: 5
df_mat <- as.data.frame(as.table(M_ha), responseName = "ha") |>
left_join(as.data.frame(as.table(P), responseName = "p"), by = c("desde", "hacia")) |>
mutate(etiqueta = ifelse(ha > 0, paste0(number(100 * p, 0.1), "%\n",
number(ha, 1, big.mark = ","), " ha"), ""))
ggplot(df_mat, aes(hacia, desde, fill = p)) +
geom_tile(colour = "#000000") +
geom_text(aes(label = etiqueta, colour = p > 0.5), size = 3, lineheight = 0.9) +
scale_fill_distiller(palette = "YlOrRd", direction = 1, labels = percent,
na.value = "#1a1a1a", name = "% del origen") +
scale_colour_manual(values = c("TRUE" = "white", "FALSE" = "black"), guide = "none") +
scale_y_discrete(limits = rev(niveles), labels = label_wrap(16)) +
scale_x_discrete(position = "top", labels = label_wrap(12)) +
labs(x = paste("Hacia", ANIO_T2), y = paste("Desde", ANIO_T1)) +
tema_oscuro() +
theme(axis.text.x = element_text(size = 9), panel.grid = element_blank())
```
# Lectura de la matriz: qué cambia y hacia dónde
## Balance por clase
Para cada clase se separan cinco cantidades [@pontius2004]:
- **persistencia**: superficie que no cambió (diagonal);
- **pérdida** y **ganancia**: lo que salió de la clase y lo que entró;
- **cambio neto** = ganancia − pérdida: lo que se ve al comparar solo los totales;
- **intercambio** = 2 × mín(ganancia, pérdida): cambio que los totales esconden, porque la clase pierde en un sitio y gana en otro.
```{r}
#| label: balance
balance <- data.frame(
clase = niveles,
ha_t1 = as.numeric(rowSums(M_ha)),
ha_t2 = as.numeric(colSums(M_ha)),
persistencia = as.numeric(diag(M_ha))
) |>
mutate(perdida = ha_t1 - persistencia,
ganancia = ha_t2 - persistencia,
cambio_neto = ganancia - perdida,
intercambio = 2 * pmin(ganancia, perdida))
area_total <- sum(M_ha)
cambio_total <- area_total - sum(diag(M_ha))
cambio_neto_paisaje <- sum(abs(balance$cambio_neto)) / 2
intercambio_paisaje <- cambio_total - cambio_neto_paisaje
# Sentido del cambio: bajar de clase es perder vegetación; subir, ganarla
codigo1 <- as.integer(px$g1); codigo2 <- as.integer(px$g2)
ha_baja <- sum(px$ha[codigo2 < codigo1])
ha_sube <- sum(px$ha[codigo2 > codigo1])
```
```{r}
#| label: tbl-balance
#| tbl-cap: "Balance del cambio por clase de NDVI (ha)."
balance |>
mutate(across(where(is.numeric), ~ round(.x, 1))) |>
setNames(c("Clase", paste("ha", ANIO_T1), paste("ha", ANIO_T2), "Persistencia", "Pérdida",
"Ganancia", "Cambio neto", "Intercambio")) |>
kable(format.args = list(big.mark = ","))
```
```{r}
#| label: fig-balance
#| fig-cap: "Ganancias y pérdidas por clase. El punto negro marca el cambio neto."
#| fig-height: 3.6
balance |>
transmute(clase = factor(clase, rev(niveles)), Ganancia = ganancia, Pérdida = -perdida) |>
pivot_longer(-clase, names_to = "componente", values_to = "ha") |>
ggplot(aes(ha, clase)) +
geom_col(aes(fill = componente)) +
geom_point(data = transmute(balance, clase = factor(clase, rev(niveles)), ha = cambio_neto),
size = 2.5, colour = "white") +
geom_vline(xintercept = 0, colour = "white") +
scale_fill_manual(values = c("Ganancia" = "#4575b4", "Pérdida" = "#d73027"), name = NULL) +
scale_x_continuous(labels = label_number(big.mark = ",")) +
labs(x = "Hectáreas", y = NULL) +
tema_oscuro() + theme(legend.position = "bottom")
```
## Transiciones dominantes
```{r}
#| label: tbl-transiciones
#| tbl-cap: "Transiciones ordenadas por superficie: de qué clase a qué clase pasó el suelo."
trans <- df_mat |>
filter(desde != hacia, ha > 0) |>
arrange(desc(ha)) |>
transmute(desde, hacia, ha,
pct_del_origen = 100 * p,
pct_del_cambio = 100 * ha / cambio_total)
trans |>
head(12) |>
mutate(ha = round(ha, 1), across(starts_with("pct"), ~ round(.x, 2))) |>
setNames(c("Desde", "Hacia", "ha", "% de la clase de origen", "% de todo el cambio")) |>
kable(format.args = list(big.mark = ","))
```
## Trayectoria a lo largo de la serie
```{r}
#| label: trayectoria
# Clase de cada píxel en cada fecha (código entero; NA si no hay dato)
C <- matrix(as.integer(clasificar(N)), nrow = nrow(N))
composicion <- do.call(rbind, lapply(seq_along(anios), function(j) {
data.frame(anio = anios[j], clase = factor(CLASES_NDVI, CLASES_NDVI),
ha = tabulate(C[, j], nbins = length(CLASES_NDVI)) * ha_pixel)
}))
# Cambio entre fechas consecutivas
periodos <- do.call(rbind, lapply(seq_len(length(anios) - 1), function(j) {
a <- C[, j]; b <- C[, j + 1]; ok <- !is.na(a) & !is.na(b)
data.frame(periodo = paste(anios[j], "→", anios[j + 1]),
baja = sum(b[ok] < a[ok]) * ha_pixel,
sube = sum(b[ok] > a[ok]) * ha_pixel,
igual = sum(b[ok] == a[ok]) * ha_pixel)
}))
```
```{r}
#| label: fig-composicion
#| fig-cap: "Superficie de cada clase de NDVI en cada fecha de la serie."
#| fig-height: 4
ggplot(composicion, aes(factor(anio), ha, fill = clase)) +
geom_col() +
scale_fill_manual(values = paleta, name = NULL) +
scale_y_continuous(labels = label_number(big.mark = ",")) +
labs(x = NULL, y = "Hectáreas") +
tema_oscuro() + theme(legend.position = "bottom")
```
```{r}
#| label: tbl-periodos
#| tbl-cap: "Superficie que baja o sube de clase de NDVI entre fechas consecutivas (ha)."
periodos |>
mutate(across(where(is.numeric), ~ round(.x, 1))) |>
setNames(c("Periodo", "Baja de clase (pierde vegetación)", "Sube de clase (gana vegetación)",
"Sin cambio de clase")) |>
kable(format.args = list(big.mark = ","))
```
Una superficie grande que baja de clase en un periodo y vuelve a subir en el siguiente no es desmonte: es la variación entre años que la normalización no alcanza a quitar, o ruido de la propia imagen. Por eso las cicatrices no se definen con una sola pareja de fechas.
# Cicatrices candidatas
Una **cicatriz candidata** es un píxel que cumple tres condiciones a la vez:
1. **tenía vegetación**: la mediana de su NDVI en las fechas anteriores es de al menos `r UMBRAL_VEGETADO`;
2. **la perdió**: su NDVI cae al menos `r CAIDA_MINIMA`;
3. **no la recuperó**: en *todas* las fechas posteriores su NDVI se queda en `r UMBRAL_DESMONTE` o menos.
La tercera condición es la que separa un desmonte permanente (plataforma, camino, banco de material, expansión urbana) de una quema, una sequía o una parcela en descanso, que reverdecen. El **año de aparición** es la primera fecha en que se cumplen las tres. Las cicatrices que aparecen en la última fecha no tienen fechas posteriores que las confirmen, y se señalan como *sin confirmar*.
```{r}
#| label: cicatrices
n_f <- ncol(N)
anio_cic <- rep(NA_integer_, nrow(N))
caida <- rep(NA_real_, nrow(N))
for (k in 2:n_f) {
antes <- if (k == 2) N[, 1] else apply(N[, 1:(k - 1), drop = FALSE], 1, median, na.rm = TRUE)
despues <- suppressWarnings(apply(N[, k:n_f, drop = FALSE], 1, max, na.rm = TRUE))
es <- is.na(anio_cic) & is.finite(despues) & !is.na(antes) &
antes >= UMBRAL_VEGETADO & despues <= UMBRAL_DESMONTE & (antes - despues) >= CAIDA_MINIMA
anio_cic[es] <- anios[k]
caida[es] <- (antes - despues)[es]
}
# Parches: píxeles contiguos de cicatriz. Se descartan los menores de MIN_HA_PARCHE.
hay_cicatrices <- any(!is.na(anio_cic))
if (hay_cicatrices) {
r_anio <- a_raster(anio_cic, "anio")
r_caida <- a_raster(caida, "caida")
parches <- patches(!is.na(r_anio), directions = 8, zeroAsNA = TRUE)
id_parche <- values(parches)[px$celda, 1]
tam <- table(id_parche) * ha_pixel
conservar <- as.integer(names(tam)[tam >= MIN_HA_PARCHE])
fuera <- !is.na(id_parche) & !(id_parche %in% conservar)
anio_cic[fuera] <- NA; caida[fuera] <- NA; id_parche[is.na(anio_cic)] <- NA
hay_cicatrices <- any(!is.na(anio_cic))
}
px$anio_cic <- anio_cic
px$caida <- caida
px$tipo <- ifelse(!is.na(px$anio_cic), "Cicatriz candidata", "Sin cicatriz")
px$vegetado_t1 <- N[, 1] >= UMBRAL_VEGETADO
if (hay_cicatrices) {
r_parche <- a_raster(id_parche, "id_parche")
cicatrices <- st_as_sf(as.polygons(r_parche))
por_parche <- px[!is.na(id_parche), ] |>
mutate(id_parche = id_parche[!is.na(id_parche)]) |>
group_by(id_parche) |>
summarise(ha = sum(ha), anio = min(anio_cic), caida_media = mean(caida), .groups = "drop")
cicatrices <- left_join(cicatrices, por_parche, by = "id_parche")
cicatrices$confirmada <- ifelse(cicatrices$anio < ANIO_T2, "Sí", "Sin confirmar")
if (hay_ejidos) {
cicatrices <- suppressWarnings(st_join(cicatrices, ejidos["ejido"], largest = TRUE))
}
if (hay_eolico) {
cicatrices$dist_eolico_m <- round(as.numeric(st_distance(cicatrices, st_union(eolico))[, 1]))
}
# Distancia de cada parche (su borde más próximo) a la línea eléctrica y a la costa
if (hay_linea) {
cicatrices$dist_linea_m <- round(as.numeric(st_distance(cicatrices, linea_union)[, 1]))
}
cicatrices$dist_costa_m <- round(as.numeric(st_distance(cicatrices, costa)[, 1]))
} else {
cicatrices <- NULL
}
```
```{r}
#| label: tbl-cicatrices-anio
#| tbl-cap: "Cicatrices candidatas según la fecha en que aparecen."
if (hay_cicatrices) {
cic_anio <- px[!is.na(px$anio_cic), ] |>
group_by(anio = anio_cic) |>
summarise(ha = sum(ha), caida_media = mean(caida), .groups = "drop") |>
left_join(st_drop_geometry(cicatrices) |> count(anio, name = "parches"), by = "anio") |>
mutate(parches = replace_na(parches, 0L),
confirmada = ifelse(anio < ANIO_T2, "Sí", "Sin confirmar"))
cic_anio |>
mutate(anio = as.character(anio), ha = round(ha, 1), caida_media = round(caida_media, 2)) |>
setNames(c("Aparece en", "ha", "Caída media de NDVI", "Parches que empiezan ese año", "Confirmada")) |>
kable(format.args = list(big.mark = ","))
} else {
cat("No se detectaron cicatrices candidatas con la definición actual.")
}
```
```{r}
#| label: tbl-parches
#| tbl-cap: "Parches de cicatriz candidata por tamaño."
if (hay_cicatrices) {
cicatrices$tamano <- cut(cicatrices$ha, c(0, 0.5, 2, 10, Inf),
labels = c("< 0.5 ha", "0.5 a 2 ha", "2 a 10 ha", "> 10 ha"))
cicatrices |>
st_drop_geometry() |>
group_by(tamano, .drop = FALSE) |>
summarise(parches = n(), ha = round(sum(ha), 1), .groups = "drop") |>
setNames(c("Tamaño del parche", "Número de parches", "ha")) |>
kable(format.args = list(big.mark = ","))
}
```
## Las cicatrices sobre Google Earth
Este mapa es la herramienta de verificación: cada parche se dibuja solo con su contorno, para que debajo se vea la imagen satelital reciente y se pueda comprobar a ojo si corresponde a una plataforma, un camino, una casa o un falso positivo. Al pulsar un parche aparecen su año de aparición, su superficie y la caída de NDVI.
```{r}
#| label: fig-mapa-cicatrices
#| fig-cap: "Cicatrices candidatas, coloreadas por año de aparición, sobre la imagen satelital de Google."
m <- mapa_base()
if (hay_cicatrices) {
cic_ll <- st_transform(cicatrices, 4326)
cic_ll$anio_f <- factor(cic_ll$anio, levels = anios[-1])
pal_anio <- colorFactor("YlOrRd", levels = levels(cic_ll$anio_f))
m <- m |>
addPolygons(data = cic_ll, color = ~pal_anio(anio_f), weight = 2, opacity = 1,
fillOpacity = 0.08, group = "Cicatrices candidatas",
popup = ~paste0("<b>Aparece en ", anio, "</b> (", confirmada, ")<br>",
round(ha, 2), " ha<br>Caída media de NDVI: ", round(caida_media, 2))) |>
addLegend(pal = pal_anio, values = cic_ll$anio_f, title = "Aparece en", position = "bottomleft")
}
m |>
addLayersControl(baseGroups = grupos_base,
overlayGroups = c(if (hay_cicatrices) "Cicatrices candidatas", grupos_fijos),
options = layersControlOptions(collapsed = FALSE))
```
## Cicatrices y propiedad ejidal
```{r}
#| label: resumen-zonas
valor_raster <- function(capa_sf, campo) {
values(rasterize(vect(capa_sf), plantilla, field = campo))[px$celda, 1]
}
if (hay_ejidos) {
px$ejido <- ejidos$ejido[valor_raster(ejidos, "id_ej")]
px$tenencia <- ifelse(is.na(px$ejido), "Fuera de la capa", ETIQUETA_EJIDOS)
}
if (hay_eolico) {
px$eolico <- ifelse(is.na(valor_raster(zona_eolica, "dentro")),
paste("A más de", BUFFER_EOLICO_M, "m"),
paste("A", BUFFER_EOLICO_M, "m o menos de infraestructura eólica"))
}
# Distancia de cada píxel a la costa y a la línea de transmisión
px_pts <- st_as_sf(px[, c("x", "y")], coords = c("x", "y"), crs = CRS_METRICO)
etiquetar_bandas <- function(cortes) {
e <- paste0(format(head(cortes, -1), big.mark = ",", trim = TRUE), "–",
format(tail(cortes, -1), big.mark = ",", trim = TRUE), " m")
e[length(e)] <- paste0("> ", format(cortes[length(cortes) - 1], big.mark = ",", trim = TRUE), " m")
e
}
en_bandas <- function(d, cortes) cut(d, cortes, right = FALSE, labels = etiquetar_bandas(cortes))
px$dist_costa <- as.numeric(st_distance(px_pts, costa))
px$franja <- ifelse(px$dist_costa <= ANCHO_COSTA_M,
paste0("Franja costera (a ", format(ANCHO_COSTA_M, big.mark = ","), " m o menos del mar)"),
"Resto del municipio")
if (hay_linea) {
px$dist_linea <- as.numeric(st_distance(px_pts, linea_union))
px$corredor <- ifelse(px$dist_linea <= BUFFER_LINEA_M,
paste("A", BUFFER_LINEA_M, "m o menos de la línea"),
paste("A más de", BUFFER_LINEA_M, "m"))
}
# Para cada zona: cuánta superficie vegetada había al inicio y cuánta quedó como
# cicatriz. La "concentración" compara la parte de las cicatrices que cae en la
# zona con la parte de la vegetación inicial que la zona tenía: un valor mayor
# que 1 indica que la zona carga con más cicatriz de la que le correspondería
# si el desmonte se repartiera de forma pareja.
resumen_zona <- function(zona) {
px |>
mutate(zona = zona) |>
group_by(zona) |>
summarise(superficie = sum(ha),
vegetado_t1 = sum(ha[vegetado_t1]),
cicatriz = sum(ha[tipo == "Cicatriz candidata"]),
.groups = "drop") |>
mutate(pct_vegetado_perdido = 100 * cicatriz / vegetado_t1,
pct_de_cicatrices = 100 * cicatriz / sum(cicatriz),
concentracion = (cicatriz / sum(cicatriz)) / (vegetado_t1 / sum(vegetado_t1)))
}
mostrar_zonas <- function(tab, nombre_zona) {
tab |>
mutate(across(c(superficie, vegetado_t1, cicatriz), ~ round(.x, 1)),
across(c(pct_vegetado_perdido, pct_de_cicatrices, concentracion), ~ round(.x, 2))) |>
setNames(c(nombre_zona, "Superficie (ha)", paste("Con vegetación en", ANIO_T1, "(ha)"),
"Cicatriz candidata (ha)", "% de su vegetación convertida en cicatriz",
"% del total de cicatrices", "Índice de concentración")) |>
kable(format.args = list(big.mark = ","))
}
```
```{r}
#| label: tbl-tenencia
#| tbl-cap: "Cicatrices según tenencia de la tierra."
if (hay_ejidos) {
zona_tenencia <- resumen_zona(px$tenencia)
mostrar_zonas(zona_tenencia, "Tenencia")
}
```
```{r}
#| label: tbl-ejidos
#| tbl-cap: "Cicatrices en cada polígono de la capa de tenencia."
if (hay_ejidos) {
zona_ejido <- resumen_zona(ifelse(is.na(px$ejido), "(fuera de la capa)", px$ejido))
mostrar_zonas(arrange(zona_ejido, desc(cicatriz)), "Polígono")
}
```
```{r}
#| label: sin-ejidos
#| echo: false
#| results: asis
if (!hay_ejidos) {
cat("> No se encontró la capa de tenencia en `", RUTA_EJIDOS,
"`. Añádela para obtener el cruce con la propiedad ejidal.", sep = "")
}
```
## Cicatrices e infraestructura eólica
```{r}
#| label: tbl-eolico
#| tbl-cap: "Cicatrices según cercanía a la infraestructura eólica."
if (hay_eolico) {
zona_eolico <- resumen_zona(px$eolico)
mostrar_zonas(zona_eolico, "Zona")
}
```
```{r}
#| label: tbl-eolico-anio
#| tbl-cap: "Cicatrices candidatas (ha) por año de aparición y cercanía a la infraestructura eólica."
if (hay_eolico && hay_cicatrices) {
cruce_ea <- xtabs(ha ~ anio_cic + eolico, data = px[!is.na(px$anio_cic), ])
names(dimnames(cruce_ea)) <- c("Aparece en", "")
kable(round(addmargins(cruce_ea), 1), format.args = list(big.mark = ","))
}
```
```{r}
#| label: tbl-eolico-ejido
#| tbl-cap: "Cicatrices candidatas (ha) según tenencia y cercanía a la infraestructura eólica."
if (hay_eolico && hay_ejidos && hay_cicatrices) {
cruce_ee <- xtabs(ha ~ tenencia + eolico, data = px[px$tipo == "Cicatriz candidata", ])
kable(round(addmargins(cruce_ee), 1), format.args = list(big.mark = ","))
}
```
```{r}
#| label: sin-eolico
#| echo: false
#| results: asis
if (!hay_eolico) {
cat("> No se encontró la capa de infraestructura eólica. Sin ella, las cicatrices ",
"candidatas no pueden atribuirse a un proyecto.", sep = "")
}
```
## Cicatrices y línea de transmisión eléctrica
Una línea de transmisión deja una cicatriz propia: la brecha que se desmonta bajo los cables. Si la línea estuviera construida, a lo largo de su trazo deberían aparecer cicatrices alargadas. Aquí se comprueba si las cicatrices detectadas se sobreponen a la línea, de dos formas: por superficie (cuánta cicatriz cae en un corredor de `r BUFFER_LINEA_M` m a cada lado) y por parches (cuántos tocan la línea o entran en ese corredor).
```{r}
#| label: tbl-linea-datos
#| tbl-cap: "Línea de transmisión cargada."
if (hay_linea) {
km_linea_total <- as.numeric(sum(st_length(linea))) / 1000
km_linea_mun <- as.numeric(sum(st_length(suppressWarnings(
st_intersection(linea_union, st_geometry(limite_m)))))) / 1000
campos <- intersect(c("Nombre", "Tipo", "Empresa", "Fase"), names(linea))
datos_linea <- data.frame(
dato = c(campos, "Longitud total (km)", "Longitud dentro del municipio (km)"),
valor = c(vapply(campos, function(k) paste(unique(na.omit(as.character(linea[[k]]))), collapse = "; "), ""),
fmt(km_linea_total, 2), fmt(km_linea_mun, 2)))
kable(datos_linea, col.names = c("Dato", "Valor"), row.names = FALSE)
}
```
```{r}
#| label: tbl-linea-corredor
#| tbl-cap: "Cicatrices dentro y fuera del corredor de la línea de transmisión."
if (hay_linea) {
zona_corredor <- resumen_zona(px$corredor)
mostrar_zonas(zona_corredor, "Zona")
}
```
```{r}
#| label: tbl-linea-bandas
#| tbl-cap: "Cicatrices según la distancia a la línea de transmisión."
if (hay_linea) {
zona_bandas_linea <- resumen_zona(en_bandas(px$dist_linea, BANDAS_LINEA_M))
mostrar_zonas(zona_bandas_linea, "Distancia a la línea")
}
```
```{r}
#| label: linea-parches
cicatriz_total_ha <- sum(px$ha[px$tipo == "Cicatriz candidata"])
txt_linea <- if (!hay_linea) {
"No se cargó ninguna línea de transmisión."
} else if (!hay_cicatrices) {
"No hay cicatrices candidatas que cruzar con la línea de transmisión."
} else {
n_tocan <- sum(cicatrices$dist_linea_m == 0)
n_corredor <- sum(cicatrices$dist_linea_m <= BUFFER_LINEA_M)
ha_corredor <- sum(px$ha[px$tipo == "Cicatriz candidata" & px$dist_linea <= BUFFER_LINEA_M])
mas_cercano <- min(cicatrices$dist_linea_m)
fase <- if ("Fase" %in% names(linea)) paste(unique(na.omit(as.character(linea$Fase))), collapse = "; ") else ""
paste0(
"De los ", nrow(cicatrices), " parches de cicatriz, ", n_tocan, " se sobreponen a la línea de transmisión y ",
n_corredor, " entran en el corredor de ", BUFFER_LINEA_M, " m (", fmt(ha_corredor, 1),
" ha de cicatriz, el ", fmt(100 * ha_corredor / cicatriz_total_ha, 2), " % del total). ",
"El parche más cercano queda a ", fmt(mas_cercano), " m de la línea. ",
ifelse(n_corredor == 0,
paste0("**Las cicatrices no se sobreponen a la línea de transmisión**: a lo largo de sus ",
fmt(km_linea_mun, 1), " km dentro del municipio no hay ninguna pérdida permanente de vegetación."),
"Hay cicatrices sobre el trazo o junto a él: conviene revisarlas una a una en el mapa para ver si siguen la forma alargada de una brecha."),
ifelse(nzchar(fase), paste0(" La capa registra la línea con la fase «", fase,
"», lo que es coherente con que no haya dejado huella en las imágenes."), ""))
}
```
**Resultado.** `r txt_linea`
## Cicatrices en la costa
La costa se toma de `r fuente_costa` (`r fmt(km_costa, 1)` km, la línea celeste de los mapas). Se llama **franja costera** a la parte del municipio situada a `r format(ANCHO_COSTA_M, big.mark = ",")` m o menos del mar; es la zona de playa, duna, ciénaga y manglar. La tabla por bandas muestra cómo cambia el resultado si la franja se define más angosta o más ancha.
```{r}
#| label: tbl-costa-franja
#| tbl-cap: "Cicatrices en la franja costera y en el resto del municipio."
zona_franja <- resumen_zona(px$franja)
mostrar_zonas(zona_franja, "Zona")
```
```{r}
#| label: tbl-costa-bandas
#| tbl-cap: "Cicatrices según la distancia a la costa."
zona_bandas_costa <- resumen_zona(en_bandas(px$dist_costa, BANDAS_COSTA_M))
mostrar_zonas(zona_bandas_costa, "Distancia a la costa")
```
```{r}
#| label: fig-costa-bandas
#| fig-cap: "Superficie de cicatrices candidatas según la distancia a la costa."
#| fig-height: 3.8
if (hay_cicatrices) {
ggplot(zona_bandas_costa, aes(zona, cicatriz)) +
geom_col(fill = "#d7191c") +
geom_text(aes(label = paste0(fmt(cicatriz, 1), " ha")), vjust = -0.4, colour = "white", size = 3.3) +
scale_y_continuous(labels = label_number(big.mark = ","), expand = expansion(mult = c(0, 0.12))) +
labs(x = "Distancia a la costa", y = "Cicatriz candidata (ha)") +
tema_oscuro()
}
```
```{r}
#| label: tbl-costa-anio
#| tbl-cap: "Cicatrices candidatas (ha) por año de aparición, en la franja costera y en el resto del municipio."
if (hay_cicatrices) {
cruce_ca <- xtabs(ha ~ anio_cic + franja, data = px[!is.na(px$anio_cic), ])
names(dimnames(cruce_ca)) <- c("Aparece en", "")
kable(round(addmargins(cruce_ca), 1), format.args = list(big.mark = ","))
}
```
```{r}
#| label: costa-parches
txt_costa <- if (!hay_cicatrices) {
"No hay cicatrices candidatas que situar respecto a la costa."
} else {
fr <- zona_franja[grepl("^Franja", zona_franja$zona), ]
# Un parche se cuenta en la costa si su punto interior representativo cae en la franja
dist_centro <- as.numeric(st_distance(suppressWarnings(st_point_on_surface(st_geometry(cicatrices))), costa))
n_costa <- sum(dist_centro <= ANCHO_COSTA_M)
ha_parches_costa <- sum(cicatrices$ha[dist_centro <= ANCHO_COSTA_M])
a_1km <- sum(px$ha[px$tipo == "Cicatriz candidata" & px$dist_costa <= 1000])
paste0(
"**", n_costa, " de los ", nrow(cicatrices), " parches de cicatriz (", fmt(100 * n_costa / nrow(cicatrices), 1),
" %) se alojan en la franja costera.** Por superficie, la franja concentra ", fmt(fr$cicatriz, 1),
" ha de las ", fmt(cicatriz_total_ha, 1), " ha de cicatrices (", fmt(fr$pct_de_cicatrices, 1),
" %), aunque solo ocupa el ", fmt(100 * fr$superficie / sum(zona_franja$superficie), 1),
" % del municipio: su índice de concentración es ", fmt(fr$concentracion, 2),
". En esa franja, el ", fmt(fr$pct_vegetado_perdido, 1), " % de la superficie que tenía vegetación en ",
ANIO_T1, " terminó como cicatriz. A menos de 1 km del mar hay ", fmt(a_1km, 1), " ha. ",
"Fuera de la franja quedan ", nrow(cicatrices) - n_costa, " parches (",
fmt(cicatriz_total_ha - fr$cicatriz, 1), " ha).")
}
```
**Resultado.** `r txt_costa`
Estas cicatrices costeras no tienen la forma de una obra (plataformas, caminos o brechas): son manchas amplias en la zona de ciénaga y manglar. Lo más probable es que reflejen un cambio del propio humedal (más agua superficial, pérdida o muerte de manglar, suelo salino descubierto) y no un desmonte. Atribuirlas a una causa exige revisar la imagen de fondo y, de ser posible, verificarlo en campo.
# Validación: índice Kappa {#sec-validacion}
Una matriz de transición describe el pasado, pero también es un modelo: si las probabilidades de cambio se mantuvieran, permitiría anticipar la fecha siguiente. Validarla consiste en ponerla a prueba con una fecha que ya se conoce:
1. **calibrar** la matriz con dos fechas pasadas;
2. **predecir** con ella una tercera fecha, separada por el mismo número de años;
3. **comparar** la predicción con lo que realmente se observó en esa tercera fecha.
La comparación se resume con el **índice Kappa**, que mide cuánto del acuerdo entre dos mapas supera al que se obtendría por azar: 0 es un acuerdo no mejor que el azar y 1 es un acuerdo perfecto. Se interpreta con la escala habitual de @landis1977. Como las clases de NDVI están ordenadas, se añade el **Kappa ponderado**, que castiga menos confundir dos clases vecinas que dos clases alejadas.
Kappa tiene dos límites conocidos, y por eso no se reporta solo [@pontius2011]. Primero, mezcla dos errores distintos, que aquí se separan: el **desacuerdo de cantidad** (el modelo predice una superficie equivocada de cada clase) y el **desacuerdo de asignación** (acierta la superficie pero la pone en otro sitio). Segundo, un modelo que diga "nada cambia" suele obtener un Kappa alto, porque la mayor parte del territorio persiste; por eso la predicción se compara con ese **modelo nulo de persistencia**. Un modelo de cambio solo aporta algo si lo supera.
```{r}
#| label: validacion-kappa
K <- length(CLASES_NDVI)
# Terna de fechas: calibración inicial, calibración final y fecha a predecir
if (is.null(ANIOS_VALIDACION)) {
ternas <- expand.grid(a = anios, b = anios, c = anios) |>
dplyr::filter(a < b, b < c, b - a == c - b) |>
arrange(desc(c), desc(b - a))
terna <- if (nrow(ternas) > 0) as.integer(unlist(ternas[1, ])) else NULL
} else {
if (!all(ANIOS_VALIDACION %in% anios)) stop("ANIOS_VALIDACION contiene fechas que no están en la serie.")
terna <- as.integer(sort(ANIOS_VALIDACION))
if (diff(terna)[1] != diff(terna)[2]) {
warning("Los dos intervalos de ANIOS_VALIDACION no duran lo mismo: la matriz se aplica tal cual.")
}
}
hay_validacion <- !is.null(terna) && length(terna) == 3
if (hay_validacion) {
jv <- match(terna, anios)
ok3 <- complete.cases(C[, jv]) # píxeles con dato en las tres fechas
c_a <- C[ok3, jv[1]]; c_b <- C[ok3, jv[2]]; c_o <- C[ok3, jv[3]]
# 1. Calibración: matriz de transición entre las dos primeras fechas
T_cal <- table(factor(c_a, 1:K), factor(c_b, 1:K))
P_cal <- matrix(as.numeric(T_cal), K, K) / as.numeric(rowSums(T_cal))
sin_origen <- rowSums(T_cal) == 0 # clase ausente al inicio: se supone que persiste
P_cal[sin_origen, ] <- diag(K)[sin_origen, ]
dimnames(P_cal) <- list(desde = CLASES_NDVI, hacia = CLASES_NDVI)
stopifnot(all(abs(rowSums(P_cal) - 1) < 1e-8))
# 2. Predicción de la tercera fecha. La cadena de Markov dice CUÁNTA superficie
# de cada clase pasa a cada otra, pero no QUÉ píxeles. Se prueban dos reglas:
# - clase más probable: cada píxel recibe el destino más probable de su clase;
# - reparto según NDVI: se respeta la superficie que espera la cadena y, dentro
# de cada clase de origen, bajan de clase los píxeles de NDVI más bajo y
# suben los de NDVI más alto.
repartir_segun_ndvi <- function(clase, ndvi_valor, P) {
pred <- integer(length(clase))
for (i in 1:K) {
idx <- which(clase == i)
if (length(idx) == 0) next
cortes <- round(cumsum(P[i, ]) * length(idx)); cortes[K] <- length(idx)
pred[idx[order(ndvi_valor[idx])]] <- rep(1:K, times = diff(c(0, cortes)))
}
pred
}
pred_nulo <- c_b # modelo nulo: nada cambia
pred_map <- max.col(P_cal, ties.method = "first")[c_b]
pred_markov <- repartir_segun_ndvi(c_b, N[ok3, jv[2]], P_cal)
ha_markov <- as.numeric(tabulate(c_b, K) %*% P_cal) * ha_pixel
# 3. Comparación con lo observado
pesos_lin <- 1 - abs(outer(1:K, 1:K, "-")) / (K - 1) # acuerdo parcial entre clases vecinas
evaluar <- function(pred, obs, base) {
cm <- table(factor(pred, 1:K), factor(obs, 1:K))
p <- cm / sum(cm); fila <- rowSums(p); col <- colSums(p)
po <- sum(diag(p)); pe <- sum(fila * col)
po_w <- sum(pesos_lin * p); pe_w <- sum(pesos_lin * outer(fila, col))
cantidad <- sum(abs(fila - col)) / 2
# Figura de mérito: de todo el cambio observado o predicho, qué parte se acertó
cambio_obs <- obs != base; cambio_pred <- pred != base
merito <- sum(cambio_obs & cambio_pred & pred == obs) / max(1, sum(cambio_obs | cambio_pred))
list(matriz = cm,
resumen = data.frame(exactitud = po, kappa = (po - pe) / (1 - pe),
kappa_pond = (po_w - pe_w) / (1 - pe_w),
cantidad = cantidad, asignacion = (1 - po) - cantidad,
merito = merito))
}
val_nulo <- evaluar(pred_nulo, c_o, c_b)
val_map <- evaluar(pred_map, c_o, c_b)
val_markov <- evaluar(pred_markov, c_o, c_b)
tabla_val <- bind_rows("Persistencia (nada cambia)" = val_nulo$resumen,
"Cadena de Markov, clase más probable" = val_map$resumen,
"Cadena de Markov, reparto según NDVI" = val_markov$resumen,
.id = "modelo")
# Escala de Landis y Koch (1977)
leer_kappa <- function(k) {
as.character(cut(k, c(-Inf, 0, 0.2, 0.4, 0.6, 0.8, Inf), right = TRUE,
labels = c("pobre (no mejor que el azar)", "leve", "aceptable", "moderado",
"considerable", "casi perfecto")))
}
}
```
```{r}
#| label: aviso-validacion
#| echo: false
#| results: asis
if (hay_validacion) {
cat("La matriz se calibra con **", terna[1], " → ", terna[2], "** y se usa para predecir **", terna[3],
"**, que se compara con el NDVI observado ese año. Los dos intervalos duran ", diff(terna)[1],
" años, así que la matriz se aplica directamente. Entran en la validación ",
fmt(sum(ok3) * ha_pixel), " ha: los píxeles con dato en las tres fechas.", sep = "")
} else {
cat("> No hay tres fechas separadas por el mismo intervalo, así que no se puede validar la ",
"matriz. Fija `ANIOS_VALIDACION` a mano o añade fechas a la serie.", sep = "")
}
```
```{r}
#| label: tbl-matriz-calibracion
#| tbl-cap: "Matriz de transición de calibración (% de la superficie de origen). Cada fila suma 100."
if (hay_validacion) kable(round(100 * P_cal, 2))
```
```{r}
#| label: tbl-kappa
#| tbl-cap: "Validación de la predicción frente al NDVI observado."
if (hay_validacion) {
tabla_val |>
mutate(across(c(exactitud, cantidad, asignacion), ~ round(100 * .x, 2)),
across(c(kappa, kappa_pond, merito), ~ round(.x, 3)),
lectura = leer_kappa(kappa)) |>
setNames(c("Modelo", "Exactitud global (%)", "Kappa", "Kappa ponderado",
"Desacuerdo de cantidad (%)", "Desacuerdo de asignación (%)", "Figura de mérito",
"Lectura del Kappa")) |>
kable()
}
```
```{r}
#| label: tbl-confusion
#| tbl-cap: "Matriz de confusión de la cadena de Markov con reparto según NDVI (ha). Filas: clase predicha; columnas: clase observada."
if (hay_validacion) {
cm_ha <- unclass(val_markov$matriz) * ha_pixel
dimnames(cm_ha) <- list(Predicha = CLASES_NDVI, Observada = CLASES_NDVI)
kable(round(addmargins(as.table(cm_ha)), 1), format.args = list(big.mark = ","))
}
```
```{r}
#| label: tbl-cantidades
#| tbl-cap: "Superficie de cada clase en la fecha predicha: observada, esperada por la cadena de Markov y según el modelo nulo (ha)."
if (hay_validacion) {
cantidades <- data.frame(
clase = CLASES_NDVI,
observada = tabulate(c_o, K) * ha_pixel,
markov = ha_markov,
nulo = tabulate(c_b, K) * ha_pixel
) |>
mutate(error_markov = markov - observada, error_nulo = nulo - observada)
cantidades |>
mutate(across(where(is.numeric), ~ round(.x, 1))) |>
setNames(c("Clase", paste("Observada", terna[3]), "Esperada por Markov",
paste0("Modelo nulo (= ", terna[2], ")"), "Error de Markov", "Error del modelo nulo")) |>
kable(format.args = list(big.mark = ","))
}
```
```{r}
#| label: fig-mapa-validacion
#| fig-cap: "Dónde falla la predicción de la cadena de Markov con reparto según NDVI, sobre la imagen satelital. Solo se colorean los píxeles en que la clase predicha y la observada no coinciden."
if (hay_validacion) {
# 1 = se observó menos vegetación de la predicha; 2 = se observó más
discrepancia <- rep(NA_real_, nrow(px))
discrepancia[ok3] <- ifelse(pred_markov == c_o, NA, ifelse(c_o < pred_markov, 1, 2))
pal_disc <- colorFactor(c("#ff4d4d", "#4db8ff"), levels = 1:2, na.color = "transparent")
mapa_base() |>
addRasterImage(a_raster(discrepancia, "discrepancia"), colors = pal_disc, method = "ngb",
opacity = 0.85, group = "Discrepancias") |>
addLegend(colors = c("#ff4d4d", "#4db8ff"), opacity = 1, position = "bottomleft",
labels = c("Menos vegetación que la predicha", "Más vegetación que la predicha"),
title = paste("Predicción de", terna[3])) |>
addLayersControl(baseGroups = grupos_base, overlayGroups = c("Discrepancias", grupos_fijos),
options = layersControlOptions(collapsed = FALSE))
}
```
```{r}
#| label: lectura-validacion
txt_valid <- if (hay_validacion) {
km <- val_markov$resumen; kn <- val_nulo$resumen
supera <- km$kappa > kn$kappa + 0.005
peor_cantidad <- cantidades[which.max(abs(cantidades$error_markov)), ]
paste0(
"Calibrada con ", terna[1], "–", terna[2], ", la cadena de Markov con reparto según NDVI predice ", terna[3],
" con una exactitud global de ", pct(km$exactitud), " y un Kappa de ", fmt(km$kappa, 3),
" (acuerdo ", leer_kappa(km$kappa), "; Kappa ponderado ", fmt(km$kappa_pond, 3),
"). Con la regla de la clase más probable el Kappa es ", fmt(val_map$resumen$kappa, 3),
". El modelo nulo de persistencia obtiene un Kappa de ", fmt(kn$kappa, 3), ": el modelo de cambio ",
ifelse(supera, "**supera**", ifelse(abs(km$kappa - kn$kappa) <= 0.005, "**iguala**", "**no supera**")),
" a la hipótesis de que nada cambia. Su error se reparte en ", pct(km$cantidad),
" de cantidad y ", pct(km$asignacion), " de asignación: ",
ifelse(km$asignacion > km$cantidad,
"predomina el error de localización, es decir, el modelo estima razonablemente cuánta superficie hay de cada clase pero no dónde.",
"predomina el error de cantidad, es decir, las tasas de cambio del periodo de calibración no se mantuvieron en el siguiente."),
" La clase peor estimada en superficie es *", peor_cantidad$clase, "* (",
ifelse(peor_cantidad$error_markov > 0, "sobrestimada", "subestimada"), " en ",
fmt(abs(peor_cantidad$error_markov)), " ha). La figura de mérito es ", fmt(km$merito, 3),
": de todo el cambio observado o predicho, esa es la parte que el modelo acertó.")
} else "No hay tres fechas con el mismo intervalo, así que la matriz no se pudo validar."
```
**Lectura.** `r txt_valid`
Al interpretar este resultado conviene recordar tres cosas. La cadena de Markov solo dice cuánta superficie cambia; qué píxeles cambian depende de la regla de reparto. La de la clase más probable manda todos los píxeles de una clase al mismo destino y por eso hace desaparecer clases enteras; la de reparto según NDVI respeta las superficies esperadas, pero supone que dentro de cada clase pierden vegetación los píxeles que ya tenían menos, lo cual es un supuesto y no un hecho. En segundo lugar, las tres fechas provienen de sensores distintos y, aunque están normalizadas, parte del desacuerdo es ruido de las imágenes y no cambio real. Por último, un Kappa alto no valida el modelo si el modelo nulo obtiene uno igual o mayor: indica solo que la mayor parte del territorio no cambió.
# Proyección de Markov y escenarios {#sec-proyeccion}
Si las probabilidades de cambio de la última década se mantuvieran, ¿cuánta vegetación más se perdería? Para responderlo se construye una cadena de Markov con cinco estados: las cuatro clases de NDVI y un quinto estado, **cicatriz**, para los píxeles que perdieron la vegetación de forma permanente. La cicatriz es un estado sin retorno: por definición, un píxel que entra en él no vuelve a salir.
La matriz se calibra con el intervalo más reciente de la serie y se aplica hacia adelante en pasos de la misma duración. Sobre esa base se comparan los escenarios, que se distinguen por un solo parámetro: el factor por el que se multiplica la probabilidad de pasar a cicatriz.
- **Tendencial** (factor 1): la última década se repite.
- **Disruptivo** (factor mayor que 1): la pérdida permanente se acelera, porque el deterioro del humedal costero se intensifica o porque se ejecutan obras que hoy están en proyecto.
- **De política** (factor menor que 1): protección o restauración efectiva que frena la pérdida.
```{r}
#| label: proyeccion
ESTADOS <- c("Cicatriz", CLASES_NDVI)
S <- length(ESTADOS)
# Estado de cada píxel en la fecha j: cicatriz si ya había aparecido; si no, su clase de NDVI
estado_en <- function(j) {
e <- C[, j] + 1L
e[!is.na(px$anio_cic) & px$anio_cic <= anios[j]] <- 1L
e
}
# Calibración: el último intervalo de la terna de validación o, si no la hay, las dos últimas fechas
if (hay_validacion) {
j_a <- match(terna[2], anios); j_b <- match(terna[3], anios)
} else {
j_a <- length(anios) - 1; j_b <- length(anios)
}
paso <- anios[j_b] - anios[j_a]
anio_base <- anios[j_b]
anios_proy <- anio_base + paso * seq_len(PASOS_PROYECCION)
e_a <- estado_en(j_a); e_b <- estado_en(j_b)
matriz_estados <- function(x, y, sel = rep(TRUE, length(x))) {
ok <- sel & !is.na(x) & !is.na(y)
tab <- table(factor(x[ok], 1:S), factor(y[ok], 1:S))
P <- matrix(as.numeric(tab), S, S) / pmax(as.numeric(rowSums(tab)), 1)
vacias <- rowSums(tab) == 0
P[vacias, ] <- diag(S)[vacias, ] # estado ausente: se supone que persiste
dimnames(P) <- list(desde = ESTADOS, hacia = ESTADOS)
P
}
ha_por_estado <- function(x, sel = rep(TRUE, length(x))) tabulate(x[sel & !is.na(x)], S) * ha_pixel
# Escenario: la probabilidad de pasar a cicatriz se multiplica por k, y la
# diferencia se toma de (o se devuelve a) la persistencia del estado de origen
escalar <- function(P, k) {
Q <- P
for (i in 2:S) {
nueva <- min(P[i, 1] * k, P[i, 1] + P[i, i])
Q[i, i] <- P[i, i] - (nueva - P[i, 1]); Q[i, 1] <- nueva
}
Q
}
proyectar <- function(v0, P, pasos) {
v <- v0; filas <- vector("list", pasos)
for (p in seq_len(pasos)) { v <- as.numeric(v %*% P); filas[[p]] <- v }
do.call(rbind, filas)
}
P_proy <- matriz_estados(e_a, e_b)
stopifnot(all(abs(rowSums(P_proy) - 1) < 1e-8), P_proy[1, 1] == 1)
v_base <- ha_por_estado(e_b)
# Proyección de cada escenario: una fila por año proyectado y una columna por estado
proy <- lapply(ESCENARIOS, function(k) {
m <- proyectar(v_base, escalar(P_proy, k), PASOS_PROYECCION)
dimnames(m) <- list(anios_proy, ESTADOS); m
})
vegetado <- CLASES_NDVI[CORTES_NDVI[seq_along(CLASES_NDVI)] >= 0.4] # clases con vegetación media o densa
tabla_esc <- do.call(rbind, lapply(names(proy), function(nm) {
m <- proy[[nm]]
fila <- data.frame(escenario = nm, factor = ESCENARIOS[[nm]], base = v_base[1])
for (a in rownames(m)) fila[[paste0("cic_", a)]] <- m[a, "Cicatriz"]
fila$veg_base <- sum(v_base[match(vegetado, ESTADOS)])
fila$veg_final <- sum(m[nrow(m), vegetado])
fila
}))
```
La calibración usa el intervalo **`r anios[j_a]` → `r anio_base`** (`r paso` años) y la proyección avanza en pasos de `r paso` años hasta `r max(anios_proy)`.
```{r}
#| label: tbl-matriz-proyeccion
#| tbl-cap: "Matriz de transición usada para proyectar (% de la superficie de origen). La cicatriz es un estado sin retorno."
kable(round(100 * P_proy, 2))
```
## Escenarios
```{r}
#| label: tbl-escenarios
#| tbl-cap: "Superficie de cicatriz (pérdida permanente de vegetación) en cada escenario, y vegetación media y densa que quedaría al final (ha)."
tabla_esc |>
mutate(across(where(is.numeric) & !factor, ~ round(.x))) |>
setNames(c("Escenario", "Factor sobre la pérdida", paste("Cicatriz", anio_base),
paste("Cicatriz", anios_proy), paste("Vegetación media y densa", anio_base),
paste("Vegetación media y densa", max(anios_proy)))) |>
kable(format.args = list(big.mark = ","), row.names = FALSE)
```
```{r}
#| label: fig-escenarios
#| fig-cap: "Cicatriz acumulada: lo observado en la serie y lo proyectado en cada escenario."
#| fig-height: 4.5
observado <- data.frame(anio = anios, escenario = "Observado",
ha = vapply(anios, function(a) sum(px$ha[!is.na(px$anio_cic) & px$anio_cic <= a]), 0))
trayectorias <- do.call(rbind, lapply(names(proy), function(nm) {
data.frame(anio = c(anio_base, anios_proy), escenario = nm, ha = c(v_base[1], proy[[nm]][, "Cicatriz"]))
}))
colores_esc <- c("Observado" = "#ffffff",
setNames(c("#ffc533", "#ff4d4d", "#66bb6a", "#4db8ff", "#b39ddb", "#f48fb1")[seq_along(proy)],
names(proy)))
ggplot(mapping = aes(anio, ha, colour = escenario)) +
geom_line(data = trayectorias, linewidth = 1, linetype = "dashed") +
geom_point(data = trayectorias, size = 2.5) +
geom_line(data = observado, linewidth = 1.2) + geom_point(data = observado, size = 2.5) +
scale_colour_manual(values = colores_esc, breaks = names(colores_esc), name = NULL) +
scale_y_continuous(labels = label_number(big.mark = ",")) +
labs(x = NULL, y = "Cicatriz acumulada (ha)") +
tema_oscuro() + theme(legend.position = "bottom")
```
## Sensibilidad
El resultado depende sobre todo del factor aplicado a la pérdida permanente. La tabla muestra cuánta cicatriz habría con distintos valores, sin cambiar nada más.
```{r}
#| label: tbl-sensibilidad
#| tbl-cap: "Sensibilidad de la cicatriz proyectada (ha) al factor sobre la probabilidad de pasar a cicatriz."
sens <- t(vapply(SENSIBILIDAD_K, function(k) proyectar(v_base, escalar(P_proy, k), PASOS_PROYECCION)[, 1],
numeric(PASOS_PROYECCION)))
if (PASOS_PROYECCION == 1) sens <- t(sens)
data.frame(factor = SENSIBILIDAD_K, round(sens)) |>
setNames(c("Factor", paste("Cicatriz", anios_proy))) |>
kable(format.args = list(big.mark = ","))
```
```{r}
#| label: sensibilidad-periodo
# Segunda sensibilidad: calibrar con el intervalo anterior, si existe en la serie
j_ant <- match(anios[j_a] - paso, anios)
hay_periodo_previo <- !is.na(j_ant)
if (hay_periodo_previo) {
P_previa <- matriz_estados(estado_en(j_ant), e_a)
alt_previa <- proyectar(v_base, P_previa, PASOS_PROYECCION)[, 1]
}
tend <- proy[[names(ESCENARIOS)[which.min(abs(ESCENARIOS - 1))]]]
txt_periodo <- if (hay_periodo_previo) {
paste0("Si la matriz se calibrara con el intervalo anterior (", anios[j_ant], " → ", anios[j_a],
"), la tendencia daría ", fmt(alt_previa[PASOS_PROYECCION]), " ha de cicatriz en ", max(anios_proy),
" en lugar de ", fmt(tend[PASOS_PROYECCION, "Cicatriz"]), " ha.")
} else "La serie no tiene un intervalo anterior de la misma duración con el que repetir la calibración."
```
`r txt_periodo`
## Dónde ocurriría
Una sola matriz para todo el municipio reparte la pérdida como si cualquier lugar tuviera el mismo riesgo. Como las cicatrices se concentran en la costa, la cadena se calcula también por separado para la franja costera y para el resto, con la tendencia observada en cada zona.
```{r}
#| label: tbl-proyeccion-zonas
#| tbl-cap: "Proyección tendencial calculada por separado en la franja costera y en el resto del municipio (ha)."
zonas <- sort(unique(px$franja))
P_zona <- lapply(setNames(zonas, zonas), function(z) matriz_estados(e_a, e_b, px$franja == z))
proy_zona <- do.call(rbind, lapply(zonas, function(z) {
v0 <- ha_por_estado(e_b, px$franja == z)
m <- proyectar(v0, P_zona[[z]], PASOS_PROYECCION)
fila <- data.frame(zona = z, superficie = sum(v0), base = v0[1])
for (p in seq_len(PASOS_PROYECCION)) fila[[paste0("cic_", anios_proy[p])]] <- m[p, 1]
fila$veg_base <- sum(v0[match(vegetado, ESTADOS)])
fila$veg_final <- sum(m[PASOS_PROYECCION, match(vegetado, ESTADOS)])
fila
}))
proy_zona |>
mutate(across(where(is.numeric), ~ round(.x))) |>
setNames(c("Zona", "Superficie", paste("Cicatriz", anio_base), paste("Cicatriz", anios_proy),
paste("Vegetación media y densa", anio_base),
paste("Vegetación media y densa", max(anios_proy)))) |>
kable(format.args = list(big.mark = ","), row.names = FALSE)
```
```{r}
#| label: fig-mapa-riesgo
#| fig-cap: "Probabilidad de que cada píxel pase a cicatriz en el próximo periodo, según su estado actual y la tendencia de su zona, sobre la imagen satelital. En blanco, el contorno de las cicatrices que ya existen. Solo se colorean los píxeles con probabilidad apreciable."
# Riesgo de cada píxel que aún no es cicatriz: la probabilidad de su estado en la matriz de su zona
riesgo <- rep(NA_real_, nrow(px))
for (z in zonas) {
sel <- px$franja == z & !is.na(e_b) & e_b > 1
riesgo[sel] <- P_zona[[z]][e_b[sel], 1]
}
riesgo[!is.na(riesgo) & riesgo < 0.01] <- NA # se deja ver el fondo donde el riesgo es casi nulo
pal_riesgo <- colorNumeric("YlOrRd", domain = c(0, max(riesgo, na.rm = TRUE)), na.color = "transparent")
m <- mapa_base() |>
addRasterImage(a_raster(riesgo, "riesgo"), colors = pal_riesgo, method = "ngb", opacity = 0.8,
group = "Riesgo de cicatriz") |>
addLegend(pal = pal_riesgo, values = c(0, max(riesgo, na.rm = TRUE)), position = "bottomleft",
title = paste0("Probabilidad en ", paso, " años"),
labFormat = labelFormat(suffix = " %", transform = function(x) round(100 * x)))
if (hay_cicatrices) {
m <- addPolygons(m, data = st_transform(cicatrices, 4326), fill = FALSE, color = "#ffffff", weight = 1,
group = "Cicatrices actuales")
}
m |> addLayersControl(baseGroups = grupos_base,
overlayGroups = c("Riesgo de cicatriz", if (hay_cicatrices) "Cicatrices actuales", grupos_fijos),
options = layersControlOptions(collapsed = FALSE))
```
## El escenario disruptivo y los proyectos eólicos
```{r}
#| label: proyectos-eolicos
# Si alguna capa eólica trae la superficie de cambio de uso de suelo declarada
# por el proyecto, se usa para dimensionar su huella directa.
etiquetas_proy <- c(Nombre = "Proyecto", Fase = "Fase", Empresa = "Empresa", Municipio = "Municipios",
PotenciaMW = "Potencia (MW)", N_Turbin = "Aerogeneradores",
Camb_Suelo = "Cambio de uso de suelo declarado (ha)",
Area_Ejido = "Superficie sobre ejido (ha)", Fecha = "Año")
proyectos <- NULL
for (f in RUTA_EOLICO[file.exists(RUTA_EOLICO)]) {
x <- st_read(f, quiet = TRUE)
if (!"Camb_Suelo" %in% names(x)) next
x <- st_transform(st_make_valid(x), CRS_METRICO)
x <- x[lengths(st_intersects(x, limite_m)) > 0, ]
if (nrow(x) == 0) next
x$ha_en_municipio <- round(as.numeric(st_area(suppressWarnings(st_intersection(st_geometry(x), st_geometry(limite_m))))) / 1e4, 1)
proyectos <- st_drop_geometry(x)[, c(intersect(names(etiquetas_proy), names(x)), "ha_en_municipio")]
}
hay_proyectos <- !is.null(proyectos)
```
```{r}
#| label: tbl-proyectos
#| tbl-cap: "Proyectos eólicos que tocan el municipio, según la capa cargada."
if (hay_proyectos) {
mostrar <- proyectos
names(mostrar) <- c(etiquetas_proy[names(proyectos)[-ncol(proyectos)]], "Polígono dentro del municipio (ha)")
kable(t(mostrar), col.names = rep("", nrow(mostrar)))
}
```
```{r}
#| label: lectura-proyeccion
nuevas_tend <- tend[1, "Cicatriz"] - v_base[1] # cicatriz nueva en el primer paso tendencial
txt_huella <- if (hay_proyectos) {
cs <- suppressWarnings(sum(as.numeric(proyectos$Camb_Suelo), na.rm = TRUE))
paste0("La capa registra ",
ifelse(nrow(proyectos) == 1, "un proyecto que toca", paste(nrow(proyectos), "proyectos que tocan")),
" el municipio, con ", fmt(cs, 1), " ha de cambio de uso de suelo declaradas en total",
if (hay_linea) paste0(" y una línea de transmisión de ", fmt(km_linea_mun, 1), " km dentro de Ixil") else "",
". Esa huella directa equivale al ", fmt(100 * cs / nuevas_tend, 0),
" % de la cicatriz nueva que la tendencia añade en ", paso, " años (", fmt(nuevas_tend),
" ha), pero caería en el interior del municipio, donde hoy casi no hay cicatrices.")
} else "Las capas eólicas cargadas no traen la superficie de cambio de uso de suelo declarada, así que la huella directa de los proyectos no se puede dimensionar aquí."
fila_esc <- function(nm) tabla_esc[tabla_esc$escenario == nm, ]
ultimo <- paste0("cic_", max(anios_proy))
txt_proy <- paste0(
"Calibrada con ", anios[j_a], "–", anio_base, ", la cadena de Markov lleva la cicatriz de ", fmt(v_base[1]),
" ha en ", anio_base, " a ", fmt(tend[1, "Cicatriz"]), " ha en ", anios_proy[1],
if (PASOS_PROYECCION > 1) paste0(" y ", fmt(tend[PASOS_PROYECCION, "Cicatriz"]), " ha en ", max(anios_proy)) else "",
" si la tendencia se mantiene. En ", max(anios_proy), " los escenarios van de ",
fmt(min(tabla_esc[[ultimo]])), " ha (", tabla_esc$escenario[which.min(tabla_esc[[ultimo]])], ") a ",
fmt(max(tabla_esc[[ultimo]])), " ha (", tabla_esc$escenario[which.max(tabla_esc[[ultimo]])], "). ",
"Calculada por zonas, la franja costera pasaría de ",
fmt(proy_zona$base[grepl("^Franja", proy_zona$zona)]), " a ",
fmt(proy_zona[[ultimo]][grepl("^Franja", proy_zona$zona)]), " ha de cicatriz, y su vegetación media y densa de ",
fmt(proy_zona$veg_base[grepl("^Franja", proy_zona$zona)]), " a ",
fmt(proy_zona$veg_final[grepl("^Franja", proy_zona$zona)]), " ha.")
```
`r txt_huella`
Esto separa dos presiones que conviene no confundir. La pérdida que se proyecta es, casi toda, la del humedal costero, y no depende de los proyectos eólicos, que hasta ahora no han dejado huella. Si esos proyectos se construyeran, su efecto directo sería menor en superficie, pero recaería en otra parte del municipio y sobre otro tipo de tenencia.
**Lectura.** `r txt_proy`
## Cómo leer esta proyección
- **Es una comparación de escenarios, no un pronóstico.** En la validación, la cadena de Markov no superó al modelo que supone que nada cambia. Los números sirven para comparar qué pasaría con más o menos presión, no para anunciar una cifra.
- **Los factores son supuestos.** Duplicar la pérdida o reducirla a la mitad son valores ilustrativos. Un escenario de política bien fundado debería apoyarse en lo que la normativa permite o prohíbe, y uno disruptivo en la superficie que los proyectos declaran.
- **La cicatriz no tiene retorno.** A muy largo plazo, una cadena con un estado así convierte todo en cicatriz. Es una propiedad matemática del modelo y no una predicción; por eso la proyección se limita a `r PASOS_PROYECCION` pasos.
- **La cadena dice cuánto, no dónde.** El mapa de riesgo asigna a cada píxel la probabilidad de su estado y de su zona; no usa la distancia al agua, la elevación ni ninguna otra variable que explique por qué un sitio se pierde y otro no.
- **Las cicatrices de la última fecha no están confirmadas** por una fecha posterior, y forman parte del periodo con que se calibra.
# Principales hallazgos {#sec-hallazgos}
```{r}
#| label: calculos-hallazgos
frase_trans <- function(i) {
if (i > nrow(trans)) return("(no hay más transiciones)")
paste0("*", trans$desde[i], " → ", trans$hacia[i], "* (", fmt(trans$ha[i], 1), " ha, ",
fmt(trans$pct_del_cambio[i], 1), " % de todo el cambio)")
}
con_area <- balance[balance$ha_t1 > 0, ]
persist <- con_area$persistencia / con_area$ha_t1
mas_estable <- con_area$clase[which.max(persist)]
mas_dinamico <- con_area$clase[which.min(persist)]
vegetado_ha <- sum(px$ha[px$vegetado_t1])
cicatriz_ha <- sum(px$ha[px$tipo == "Cicatriz candidata"])
peor_periodo <- periodos[which.max(periodos$baja), ]
txt_cic <- if (hay_cicatrices) {
pico <- cic_anio[which.max(cic_anio$ha), ]
sin_conf <- sum(cic_anio$ha[cic_anio$confirmada != "Sí"])
paste0(fmt(cicatriz_ha, 1), " ha (", pct(cicatriz_ha / vegetado_ha, 2), " de la superficie que tenía vegetación en ",
ANIO_T1, ") perdieron la vegetación y no la recuperaron. Se agrupan en ", nrow(cicatrices),
" parches; el parche mediano mide ", fmt(median(cicatrices$ha), 2), " ha y el mayor ",
fmt(max(cicatrices$ha), 1), " ha. La fecha con más superficie nueva de cicatriz es ", pico$anio,
" (", fmt(pico$ha, 1), " ha). De ese total, ", fmt(sin_conf, 1),
" ha aparecen en la última fecha y están sin confirmar.")
} else "No se detectó ninguna con la definición actual."
txt_ejidos <- if (hay_ejidos && hay_cicatrices) {
ej <- zona_tenencia[zona_tenencia$zona == ETIQUETA_EJIDOS, ]
if (nrow(ej) == 0) {
"Ningún píxel del municipio cae dentro de la capa de tenencia cargada."
} else {
paste0("El ", fmt(ej$pct_de_cicatrices, 1), " % de la superficie de cicatrices (", fmt(ej$cicatriz, 1),
" ha) cae dentro de la capa de tenencia, que concentraba el ", pct(ej$vegetado_t1 / vegetado_ha),
" de la vegetación de ", ANIO_T1, " (índice de concentración: ", fmt(ej$concentracion, 2), ").")
}
} else "Sin capa de tenencia (o sin cicatrices) no es posible cuantificar este cruce."
txt_eolico <- if (hay_eolico && hay_cicatrices) {
eo <- zona_eolico[grepl("o menos", zona_eolico$zona), ]
if (nrow(eo) == 0) {
paste0("Ningún píxel del municipio está a ", BUFFER_EOLICO_M, " m o menos de la infraestructura cargada.")
} else {
paste0("A ", BUFFER_EOLICO_M, " m o menos de la infraestructura eólica se localiza el ",
fmt(eo$pct_de_cicatrices, 1), " % de las cicatrices (", fmt(eo$cicatriz, 1),
" ha); allí se convirtió en cicatriz el ", fmt(eo$pct_vegetado_perdido, 2),
" % de la vegetación que había en ", ANIO_T1, " (índice de concentración: ",
fmt(eo$concentracion, 2), ").")
}
} else "Sin capa de infraestructura eólica las cicatrices no pueden atribuirse a un proyecto."
```
**1. Magnitud del cambio.** Entre `r ANIO_T1` y `r ANIO_T2` cambió de clase de NDVI el `r pct(cambio_total / area_total)` de Ixil (`r fmt(cambio_total)` de `r fmt(area_total)` ha). De ese cambio, `r fmt(cambio_neto_paisaje)` ha son **cambio neto** y `r fmt(intercambio_paisaje)` ha son **intercambio** (la misma clase pierde en un lugar y gana en otro).
**2. Sentido del cambio.** `r fmt(ha_baja)` ha bajaron de clase (menos vegetación) y `r fmt(ha_sube)` ha subieron (más vegetación).
**3. Hacia dónde.** La transición dominante es `r frase_trans(1)`. Le siguen `r frase_trans(2)` y `r frase_trans(3)`.
**4. Persistencia.** La clase más estable es *`r mas_estable`* (`r pct(max(persist))` de su superficie no cambió) y la más dinámica es *`r mas_dinamico`* (`r pct(min(persist))`).
**5. Trayectoria.** El periodo con más superficie que baja de clase es `r peor_periodo$periodo` (`r fmt(peor_periodo$baja)` ha).
**6. Cicatrices candidatas.** `r txt_cic`
**7. Propiedad ejidal.** `r txt_ejidos`
**8. Proyectos eólicos.** `r txt_eolico`
**9. Línea de transmisión eléctrica.** `r txt_linea`
**10. Costa.** `r txt_costa`
**11. Validación (índice Kappa).** `r txt_valid`
**12. Proyección y escenarios.** `r txt_proy` `r txt_huella`
# Verificaciones {#sec-verificaciones}
```{r}
#| label: tbl-verificaciones
#| tbl-cap: "Controles de consistencia del diagnóstico."
data.frame(
control = c("Imágenes usadas / previstas",
"Las filas de la matriz de probabilidades suman 1",
"Superficie analizada frente a la del polígono de Ixil (diferencia %)",
"Píxeles del municipio excluidos por no tener dato en la primera o la última fecha (%)",
"Suma de transiciones + persistencia = superficie analizada",
"Normalización entre fechas",
"Resolución del análisis (m)"),
resultado = c(
paste(nrow(disp), "/", nrow(imagenes)),
ifelse(all(abs(rowSums(P)[tot_fila > 0] - 1) < 1e-8), "Sí", "NO"),
fmt(100 * (area_total - area_limite_ha) / area_limite_ha, 2),
fmt(100 * (1 - sum(validos) / sum(en_mun)), 2),
ifelse(isTRUE(all.equal(sum(trans$ha) + sum(diag(M_ha)), area_total)), "Sí", "NO"),
ifelse(NORMALIZAR, paste("Sí, respecto a", anios[j_ref]), "No"),
as.character(RES_ANALISIS))
) |> kable(col.names = c("Control", "Resultado"))
```
# Exportación
```{r}
#| label: exportar
write.csv(as.data.frame.matrix(round(M_ha, 2)),
file.path(dir_salidas, "ndvi_matriz_transicion_ha.csv"))
write.csv(as.data.frame.matrix(round(P, 6)),
file.path(dir_salidas, "ndvi_matriz_transicion_probabilidad.csv"))
write.csv(balance, file.path(dir_salidas, "ndvi_balance_por_clase.csv"), row.names = FALSE)
write.csv(trans, file.path(dir_salidas, "ndvi_transiciones_ordenadas.csv"), row.names = FALSE)
write.csv(periodos, file.path(dir_salidas, "ndvi_cambio_por_periodo.csv"), row.names = FALSE)
# Serie de NDVI normalizado (una banda por fecha) y año de aparición de cada cicatriz
names(ndvi) <- paste0("ndvi_", anios)
writeRaster(ndvi, file.path(dir_salidas, "ixil_ndvi_normalizado.tif"), overwrite = TRUE)
writeRaster(a_raster(px$anio_cic, "anio_cicatriz"),
file.path(dir_salidas, "ixil_cicatrices_anio.tif"), overwrite = TRUE, datatype = "INT2U")
if (hay_cicatrices) {
ruta_cic <- file.path(dir_salidas, "cicatrices_candidatas_ndvi.gpkg")
if (file.exists(ruta_cic)) file.remove(ruta_cic)
st_write(cicatrices, ruta_cic, quiet = TRUE)
}
```
El archivo `cicatrices_candidatas_ndvi.gpkg` se puede abrir en QGIS.
# Referencias
::: {#refs}
:::
```{r}
#| label: sesion
#| code-fold: true
sessionInfo()
```
