---
title: "Análisis de la red vial de Ixil (Yucatán, México) y cambio de la vegetación"
subtitle: "Teoría de grafos con sfnetworks: cómo se relaciona la estructura vial con el cambio medido con NDVI (2001–2025)"
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, porque no esta en py
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: show
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>
```
# Presentación
Este documento lee la red vial de Ixil **como un sistema** y la cruza con el diagnóstico de cambio de la vegetación hecho con la serie de NDVI de 2001 a 2025 (realizado en el anterior codigo). El NDVI resume cuánta vegetación verde hay en cada píxel; aquel diagnóstico entrega, para cada píxel, su NDVI en cada fecha y, si perdió la vegetación de forma permanente, el año en que apareció esa **cicatriz**. Aquí se responde a cuatro preguntas:
1. **¿Qué estructura tiene la red?** Jerarquía, conectividad, centralidad y accesibilidad desde la cabecera.
2. **¿El cambio sigue a los ejes viales?** Se compara qué tipo de cambio ocurre cerca y lejos de las vías.
3. **¿Las zonas mal conectadas concentran algún tipo de cambio?** Se compara el cambio según la accesibilidad por la red.
4. **¿Las cicatrices se acercan o se alejan de las vías con el tiempo?** Se usa el año de aparición de cada cicatriz.
La red se construye y analiza con `sfnetworks` [@vandermeer_sfnetworks] a partir de las vías de OpenStreetMap, descargadas con `osmdata` [@padgham2017].
# Preparación del entorno
```{r}
#| label: librerias
# Al renderizar fuera de RStudio puede no haber un repositorio de CRAN definido
if (identical(unname(getOption("repos")["CRAN"]), "@CRAN@")) {
options(repos = c(CRAN = "https://cloud.r-project.org"))
}
if (!requireNamespace("pacman", quietly = TRUE)) install.packages("pacman")
# sfnetworks depende de sfheaders, que exige Rcpp >= 1.1.0. Si la versión instalada
# es anterior, se actualiza. Si aun así falla la carga, reinicia R y vuelve a ejecutar.
if (!requireNamespace("Rcpp", quietly = TRUE) || packageVersion("Rcpp") < "1.1.0") {
install.packages("Rcpp")
}
pacman::p_load(sf,
sfnetworks,
tidygraph,
tidyverse,
terra,
osmdata,
httr2, # lo necesita osmdata para descargar
spdep,
leaflet,
knitr,
scales
)
terraOptions(progress = 0)
fmt <- function(x, d = 0) format(round(x, d), big.mark = ",", nsmall = d, trim = TRUE)
# 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())
```
## Parámetros del análisis
```{r}
#| label: parametros
# Límite municipal: el mismo shapefile del diagnóstico de cambio
RUTA_SHP_IXIL <- "G:/FLACSO/SEMILLERO_ PAPER EJIDOS MEXICO/DATOS ESPACIALES/ixil_muni.shp"
CAMPO_FILTRO <- NULL
VALOR_FILTRO <- NULL
# Resultado del diagnóstico con NDVI (diagnostico_cambio_ndvi_ixil.qmd)
DIR_NDVI <- "salidas_ixil_ndvi"
RUTA_NDVI <- file.path(DIR_NDVI, "ixil_ndvi_normalizado.tif") # una banda por fecha
RUTA_CIC_ANIO <- file.path(DIR_NDVI, "ixil_cicatrices_anio.tif") # año de cada cicatriz
RUTA_CIC_GPKG <- file.path(DIR_NDVI, "cicatrices_candidatas_ndvi.gpkg") # parches (para el mapa)
# Clases de NDVI: deben ser las mismas que en el diagnóstico con NDVI
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")
# Red vial. NULL = descargar de OpenStreetMap (se guarda una copia en CACHE_OSM
# para no volver a descargar). Para usar una capa propia (p. ej. la Red Nacional
# de Caminos de INEGI) indica su ruta y el campo con el tipo de vía.
# >>> AQUÍ VA EL SHAPEFILE DE VÍAS <<< (y el nombre del campo con el tipo de vía)
RUTA_RED <- NULL
CAMPO_TIPO_RED <- "highway"
# La red se toma con este margen alrededor del municipio. No conviene reducirlo:
# varias partes de Ixil (la franja costera, por ejemplo) solo se conectan entre sí
# por carreteras que pasan por los municipios vecinos, y con un margen pequeño
# quedan como islas sueltas que la limpieza topológica descarta.
BUFFER_RED_M <- 5000
CACHE_OSM <- sprintf("datos_ixil/red_vial_osm_margen_%sm.gpkg", BUFFER_RED_M)
MIN_HA_CATEGORIA <- 10 # categorías de cambio con menos superficie no se leen en
# los índices de concentración (son demasiado ruidosas)
# Origen de las rutas y de la accesibilidad. NULL = el punto de la red con más
# intersecciones a menos de 500 m dentro del municipio, es decir, el centro de la
# traza urbana de la cabecera. Para fijarlo a mano: c(longitud, latitud).
CABECERA_LONLAT <- NULL
# Infraestructura eólica (opcional; si las rutas no existen, la sección se omite)
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
ALFA <- 2 # penalización de las vías poco transitables (0 = solo distancia)
UMBRAL_ACCES_M <- 5000 # umbral de accesibilidad desde la cabecera (m)
BANDAS_VIA_M <- c(0, 100, 250, 500, 1000, Inf) # bandas de distancia a la vía (m)
ANCHO_MORAN_M <- 2000 # distancia máxima entre vecinos en el índice de Moran (m)
MIN_HA_NODO <- 1 # superficie mínima del área servida por un nodo (ha)
# Validación con el índice Kappa: ¿la red vial predice dónde se pierde vegetación?
CATEGORIAS_OBJETIVO <- c("Cicatriz candidata", "Pérdida de vegetación") # lo que se intenta predecir
TAM_BLOQUE_M <- 1000 # lado de los bloques que separan calibración y validación (m)
SEMILLA <- 2026 # fija el sorteo de bloques para que el resultado sea reproducible
CRS_METRICO <- "EPSG:32616" # UTM zona 16 Norte: distancias en metros
# Las salidas llevan el prefijo "red_ndvi_" para no pisar las del diagnóstico
dir_salidas <- "salidas_ixil_ndvi"
dir.create(dir_salidas, showWarnings = FALSE, recursive = TRUE)
dir.create(dirname(CACHE_OSM), showWarnings = FALSE, recursive = TRUE)
```
# Datos de entrada
## Límite municipal y patrón de cambio
```{r}
#| label: datos-limite-cambio
limite <- st_read(RUTA_SHP_IXIL, quiet = TRUE) |> st_make_valid()
if (!is.null(CAMPO_FILTRO)) {
limite <- limite[as.character(limite[[CAMPO_FILTRO]]) == as.character(VALOR_FILTRO), ]
}
limite_m <- st_sf(geometry = st_union(st_geometry(st_transform(limite, CRS_METRICO))))
zona_red <- st_as_sfc(st_bbox(st_buffer(limite_m, BUFFER_RED_M))) # rectángulo con margen
area_mun_km2 <- as.numeric(st_area(limite_m)) / 1e6
if (!file.exists(RUTA_NDVI) || !file.exists(RUTA_CIC_ANIO)) {
stop("No se encontró el resultado del diagnóstico con NDVI en '", DIR_NDVI,
"'. Renderiza primero diagnostico_cambio_ndvi_ixil.qmd.")
}
# Serie de NDVI normalizado (una banda por fecha) y año de aparición de cada cicatriz
ndvi <- rast(RUTA_NDVI)
anios_serie <- as.integer(gsub("\\D", "", names(ndvi)))
anios <- range(anios_serie) # primera y última fecha
r_cambio <- ndvi[[1]] # grilla de referencia del análisis
r_cic_anio <- rast(RUTA_CIC_ANIO)
if (!compareGeom(r_cic_anio, r_cambio, stopOnError = FALSE)) {
stop("La serie de NDVI y el raster de cicatrices no comparten grilla: vuelve a renderizar el diagnóstico.")
}
crs_r <- st_crs(crs(r_cambio))
limite_r <- st_transform(limite_m, crs_r)
# Tabla de píxeles: NDVI en la primera y la última fecha, su diferencia, año de
# cicatriz (NA si no es cicatriz), superficie y coordenadas
v1 <- values(ndvi[[which.min(anios_serie)]])[, 1]
v2 <- values(ndvi[[which.max(anios_serie)]])[, 1]
celdas <- which(!is.na(v1) & !is.na(v2))
px <- data.frame(celda = celdas,
ha = values(cellSize(r_cambio, unit = "ha"))[celdas, 1],
ndvi1 = v1[celdas],
ndvi2 = v2[celdas],
anio_cic = values(r_cic_anio)[celdas, 1])
px$delta <- px$ndvi2 - px$ndvi1
px <- cbind(px, xyFromCell(r_cambio, celdas))
# Categoría de cambio de cada píxel. La cicatriz (pérdida que persiste en todas
# las fechas posteriores) se separa de la simple bajada de clase entre la
# primera y la última fecha, que puede ser pasajera.
clasificar <- function(x) as.integer(cut(x, CORTES_NDVI, labels = CLASES_NDVI,
include.lowest = TRUE, right = FALSE))
c1 <- clasificar(px$ndvi1); c2 <- clasificar(px$ndvi2)
categorias <- c("Persistencia", "Cicatriz candidata", "Pérdida de vegetación",
"Ganancia de vegetación")
paleta_cat <- setNames(c("#e6e6e6", "#d7191c", "#fdae61", "#1a9641"), categorias)
px$cambio <- factor(case_when(
!is.na(px$anio_cic) ~ "Cicatriz candidata",
c2 < c1 ~ "Pérdida de vegetación",
c2 > c1 ~ "Ganancia de vegetación",
TRUE ~ "Persistencia"
), categorias)
px_pts <- st_as_sf(px[, c("x", "y")], coords = c("x", "y"), crs = crs_r) |>
st_transform(CRS_METRICO)
area_total_ha <- sum(px$ha)
# Auxiliar: vector de valores por píxel -> raster sobre la grilla del análisis
a_raster <- function(valores, nombre) {
r <- rast(r_cambio); v <- rep(NA_real_, ncell(r)); v[px$celda] <- valores
values(r) <- v; names(r) <- nombre; r
}
```
```{r}
#| label: mapa-base
# Todos los mapas del documento comparten el mismo fondo: imagen satelital de
# Google (por defecto), Google híbrido, Esri y un mapa base oscuro.
a_ll <- function(x) st_transform(x, 4326)
grupos_base <- c("Google Satélite", "Google Híbrido", "Esri World Imagery", "Mapa oscuro")
mapa_base <- function() {
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 = a_ll(limite_m), fill = FALSE, color = "#000000", weight = 2,
group = "Límite municipal") |>
# Todos los mapas abren encuadrados en el municipio
fitBounds(caja_ll[["xmin"]], caja_ll[["ymin"]], caja_ll[["xmax"]], caja_ll[["ymax"]])
}
caja_ll <- st_bbox(a_ll(limite_m))
# Control para encender y apagar capas y cambiar de fondo
control_capas <- function(m, grupos) {
addLayersControl(m, baseGroups = grupos_base,
overlayGroups = c(grupos, "Límite municipal"),
options = layersControlOptions(collapsed = TRUE))
}
```
```{r}
#| label: tbl-categorias
#| tbl-cap: "Categorías de cambio que se cruzan con la red vial."
px |>
group_by(cambio, .drop = FALSE) |>
summarise(ha = sum(ha), .groups = "drop") |>
mutate(pct = 100 * ha / sum(ha)) |>
mutate(ha = round(ha, 1), pct = round(pct, 2)) |>
setNames(c(paste("Categoría de cambio", anios[1], "→", anios[2]), "ha", "% del municipio")) |>
kable(format.args = list(big.mark = ","))
```
*Persistencia*, *pérdida* y *ganancia* comparan la clase de NDVI de `r anios[1]` con la de `r anios[2]` (`r paste(CLASES_NDVI, collapse = ", ")`). La *cicatriz candidata* viene del diagnóstico: un píxel que tenía vegetación, la perdió y no la recuperó en ninguna fecha posterior.
## Red vial
```{r}
#| label: datos-red
# Descarga desde OpenStreetMap todas las vías (highway) de la zona. El servidor
# público de Overpass a veces rechaza la consulta por saturación (error 429 o
# 504): se reintenta, alternando entre servidores, antes de darse por vencido.
descargar_osm <- function(zona) {
bb <- st_bbox(st_transform(zona, 4326))
servidores <- c("https://overpass-api.de/api/interpreter",
"https://overpass.kumi.systems/api/interpreter")
for (intento in 1:6) {
set_overpass_url(servidores[(intento - 1) %% length(servidores) + 1])
osm <- tryCatch(
opq(bbox = unname(bb), timeout = 120) |>
add_osm_feature(key = "highway") |>
osmdata_sf(),
error = function(e) { message("Intento ", intento, " fallido: ", conditionMessage(e)); NULL })
if (!is.null(osm) && !is.null(osm$osm_lines) && nrow(osm$osm_lines) > 0) {
return(dplyr::select(osm$osm_lines, any_of(c("osm_id", "name")), type = highway))
}
Sys.sleep(20)
}
stop("No se pudo descargar la red de OpenStreetMap. Espera unos minutos y vuelve a ",
"intentar, o indica una capa local en RUTA_RED.")
}
if (!is.null(RUTA_RED)) {
calles <- st_read(RUTA_RED, quiet = TRUE)
calles$type <- as.character(calles[[CAMPO_TIPO_RED]])
fuente_red <- paste("capa local:", basename(RUTA_RED))
} else {
if (!file.exists(CACHE_OSM)) st_write(descargar_osm(zona_red), CACHE_OSM, quiet = TRUE)
calles <- st_read(CACHE_OSM, quiet = TRUE)
fuente_red <- paste("OpenStreetMap, descargada el", format(file.mtime(CACHE_OSM), "%d/%m/%Y"))
}
# CRS proyectado y solo líneas simples. Se conservan enteras las vías que tocan
# la zona (recortarlas rompería conexiones) y se descartan las que están en
# construcción o proyectadas: aún no forman parte de la red.
calles <- calles |>
st_transform(CRS_METRICO) |>
dplyr::filter(!type %in% c("construction", "proposed"), !is.na(type)) |>
st_filter(zona_red)
calles <- suppressWarnings(st_cast(calles, "LINESTRING"))
# Jerarquía vial: agrupa los tipos de OSM en cuatro clases y les asigna una
# transitabilidad (1 = pavimentada y continua; 0.25 = brecha o sendero).
jerarquias <- c("Principal", "Local", "Servicio", "Brecha o sendero", "Otra")
clasificar_via <- function(type) {
factor(case_when(
grepl("^(motorway|trunk|primary|secondary|tertiary)", type) ~ "Principal",
type %in% c("residential", "living_street", "unclassified") ~ "Local",
type == "service" ~ "Servicio",
type %in% c("track", "path", "footway", "cycleway", "bridleway") ~ "Brecha o sendero",
TRUE ~ "Otra"
), jerarquias)
}
transitabilidad_jer <- c("Principal" = 1.00, "Local" = 0.75, "Servicio" = 0.50,
"Brecha o sendero" = 0.25, "Otra" = 0.50)
# Colores claros, legibles sobre la imagen satelital y sobre el fondo negro
paleta_jer <- c("Principal" = "#ff4d4d", "Local" = "#f0f0f0", "Servicio" = "#4db8ff",
"Brecha o sendero" = "#ffc533", "Otra" = "#9e9e9e")
calles$jerarquia <- clasificar_via(calles$type)
# Todas las vías mapeadas dentro del municipio (estén o no conectadas al resto).
# Se usan para longitudes, densidades y distancia a la vía; el grafo, en cambio,
# solo conserva la parte conectada.
vias_mun <- suppressWarnings(st_intersection(calles[, c("type", "jerarquia")], st_geometry(limite_m)))
vias_mun$km <- as.numeric(st_length(vias_mun)) / 1000
km_mun <- sum(vias_mun$km)
# Versiones en lon/lat para los mapas (geometría algo simplificada para aligerar
# el HTML): todas las vías de la zona, y las del municipio y su entorno inmediato
pal_jer_ll <- colorFactor(unname(paleta_jer), levels = names(paleta_jer))
calles_ll <- a_ll(st_simplify(calles[, c("type", "jerarquia")], dTolerance = 3))
vias_ll <- a_ll(st_filter(calles[, c("type", "jerarquia")], st_buffer(limite_m, 1500)))
```
# Construcción del grafo
Las calles (líneas) pasan a ser las **aristas** y sus extremos (puntos) los **nodos**.
```{r}
#| label: grafo
red <- as_sfnetwork(calles, directed = FALSE)
n_bruto <- c(nodos = igraph::vcount(red), aristas = igraph::ecount(red))
km_bruto <- sum(as.numeric(st_length(calles))) / 1000
red
```
# Limpieza topológica
Paso crítico con datos de OSM: se eliminan duplicados y bucles, se parten las aristas donde se cruzan, se quitan los nodos de grado 2 (que no son intersecciones) y se conserva el componente principal, es decir, la parte de la red que está conectada entre sí.
```{r}
#| label: topologia
red <- red |>
activate("edges") |>
filter(!edge_is_multiple(), !edge_is_loop()) |> # duplicados y bucles
convert(to_spatial_subdivision, .clean = TRUE) |> # une cruces interiores
convert(to_spatial_smooth, # quita nodos de grado 2
summarise_attributes = "first", .clean = TRUE) |>
activate("nodes") |>
filter(group_components() == 1) # componente principal
n_limpio <- c(nodos = igraph::vcount(red), aristas = igraph::ecount(red))
```
# Pesos de las aristas
Cada arista recibe dos pesos:
- **longitud** (m): la distancia real;
- **costo de tránsito**: la longitud penalizada por la mala transitabilidad, $costo = longitud \times (1 + \alpha\,(1 - transitabilidad))$. Con $\alpha$ = `r ALFA`, un kilómetro de brecha "cuesta" lo mismo que `r 1 + ALFA * (1 - 0.25)` km de carretera pavimentada.
```{r}
#| label: pesos-aristas
red <- red |>
activate("edges") |>
mutate(
longitud = as.numeric(edge_length()), # metros
jerarquia = clasificar_via(type),
transitabilidad = unname(transitabilidad_jer[as.character(jerarquia)]),
costo_transito = longitud * (1 + ALFA * (1 - transitabilidad))
)
```
# Métricas de centralidad
Se calculan el grado de cada nodo (cuántos tramos llegan a él), su cercanía al resto de la red y la **centralidad de intermediación** de nodos y aristas: cuántas de las rutas más cortas entre todos los pares de nodos pasan por ese elemento [@freeman1977].
```{r}
#| label: metricas
red <- red |>
activate("nodes") |>
mutate(
grado = centrality_degree(),
cercania = centrality_closeness(weights = longitud, normalized = TRUE),
interm = centrality_betweenness(weights = longitud, normalized = TRUE)
) |>
activate("edges") |>
mutate(flujo = centrality_edge_betweenness(weights = longitud))
nodos <- st_as_sf(red, "nodes")
aristas <- st_as_sf(red, "edges")
km_limpio <- sum(aristas$longitud) / 1000
# Aristas en lon/lat para los mapas (geometría algo simplificada)
aristas_ll <- a_ll(st_simplify(aristas[, c("type", "jerarquia", "longitud", "transitabilidad", "flujo")],
dTolerance = 3))
```
```{r}
#| label: fig-red
#| fig-cap: "Red vial analizada (componente principal, por jerarquía) y límite municipal de Ixil, sobre la imagen satelital. En gris punteado, las vías mapeadas que quedaron fuera del grafo por no estar conectadas."
mapa_base() |>
addPolylines(data = calles_ll, color = "#bdbdbd", weight = 1, dashArray = "3", opacity = 0.8,
group = "Vías fuera del grafo") |>
addPolylines(data = aristas_ll, color = ~pal_jer_ll(jerarquia), opacity = 0.95,
weight = ~ifelse(jerarquia == "Principal", 4, 2),
popup = ~paste0(jerarquia, " (", type, ")"), group = "Red analizada") |>
addLegend(pal = pal_jer_ll, values = droplevels(aristas_ll$jerarquia), title = "Jerarquía",
position = "bottomleft") |>
control_capas(c("Red analizada", "Vías fuera del grafo"))
```
## Indicadores de estructura
```{r}
#| label: tbl-jerarquia
#| tbl-cap: "Longitud de las vías mapeadas dentro del municipio por jerarquía vial."
tab_jer <- vias_mun |>
st_drop_geometry() |>
group_by(jerarquia) |>
summarise(km = sum(km), .groups = "drop") |>
mutate(pct = 100 * km / sum(km))
tab_jer |>
mutate(km = round(km, 1), pct = round(pct, 1)) |>
setNames(c("Jerarquía", "km", "% de la red")) |>
kable()
```
```{r}
#| label: tbl-conectividad
#| tbl-cap: "Indicadores de conectividad de la red."
V <- n_limpio[["nodos"]]; E <- n_limpio[["aristas"]]
prop_sin_salida <- mean(nodos$grado == 1)
indices <- c(
beta = E / V, # aristas por nodo
gamma = E / (3 * (V - 2)), # aristas existentes / máximo posible
alfa = (E - V + 1) / (2 * V - 5) # circuitos existentes / máximo posible
)
num <- function(x, d = 0) format(round(x, d), big.mark = ",", nsmall = d, trim = TRUE)
data.frame(
indicador = c("Nodos del grafo (intersecciones y extremos)", "Aristas del grafo (tramos)",
"Longitud del grafo, municipio más margen (km)",
"Vías mapeadas dentro del municipio (km)",
"Densidad vial del municipio (km de vía por km²)",
"Nodos sin salida (grado 1, %)", "Intersecciones de 3 o más tramos (%)",
"Índice beta (aristas por nodo)",
"Índice gamma (0 = mínima conexión, 1 = máxima)",
"Índice alfa (0 = red en árbol, 1 = malla completa)"),
valor = c(num(V), num(E), num(km_limpio, 1), num(km_mun, 1), num(km_mun / area_mun_km2, 2),
num(100 * prop_sin_salida, 1), num(100 * mean(nodos$grado >= 3), 1),
num(indices[["beta"]], 2), num(indices[["gamma"]], 2), num(indices[["alfa"]], 2))
) |> kable(col.names = c("Indicador", "Valor"), align = "lr")
```
Un índice alfa cercano a 0 describe una red **en árbol**, esto significa que como en este caso casi no hay rutas alternativas y cada camino lleva a un solo lugar. Es lo esperable en un municipio rural como Ixil donde las brechas entran al monte y terminan ahí.
```{r}
#| label: fig-betweenness
#| fig-cap: "Centralidad de intermediación de las aristas: cuántas de las rutas más cortas entre pares de nodos pasan por cada tramo. Es un indicador del flujo potencial."
# La intermediación varía en varios órdenes de magnitud: el color y el grosor
# siguen su logaritmo para que también se distingan los tramos poco centrales.
aristas_ll$log_flujo <- log10(aristas_ll$flujo + 1)
pal_flujo <- colorNumeric("YlOrRd", domain = range(aristas_ll$log_flujo))
mapa_base() |>
addPolylines(data = aristas_ll, color = ~pal_flujo(log_flujo), opacity = 0.95,
weight = ~rescale(log_flujo, to = c(1, 6)),
popup = ~paste0("Intermediación: ", format(round(flujo), big.mark = ","),
"<br>", jerarquia),
group = "Intermediación") |>
addLegend(pal = pal_flujo, values = aristas_ll$log_flujo, title = "Intermediación",
labFormat = labelFormat(transform = function(x) round(10^x - 1), big.mark = ","),
position = "bottomleft") |>
control_capas("Intermediación")
```
# Ruta más corta y ruta más transitable
Se comparan dos formas de ir **de la cabecera a la costa**: la ruta más corta en metros y la que minimiza el costo de tránsito (evita brechas aunque el recorrido sea más largo).
```{r}
#| label: origen-destino
nodos_dentro <- which(lengths(st_intersects(nodos, limite_m)) > 0)
# Origen: la cabecera. Por defecto, el nodo del municipio con más intersecciones a
# menos de 500 m: el centro de la traza urbana, donde la red es más densa.
if (is.null(CABECERA_LONLAT)) {
nd <- nodos[nodos_dentro, ]
vecinos_500 <- lengths(st_is_within_distance(nd, nd, dist = 500))
cabecera <- st_geometry(nd)[which.max(vecinos_500)]
} else {
cabecera <- st_transform(st_sfc(st_point(CABECERA_LONLAT), crs = 4326), CRS_METRICO)
}
origen <- st_nearest_feature(cabecera, nodos)
# Destino: el nodo más al norte dentro del municipio (la costa)
destino <- nodos_dentro[which.max(st_coordinates(nodos)[nodos_dentro, 2])]
if (origen == destino) stop("El origen y el destino coinciden: revisa CABECERA_LONLAT.")
# ¿El grafo llega hasta las vías mapeadas más al norte del municipio? Si esas vías
# quedaron sueltas (fuera del componente principal), la ruta no alcanza la costa.
falta_norte_m <- unname(st_bbox(vias_mun)["ymax"]) - st_coordinates(nodos)[destino, 2]
txt_norte <- if (falta_norte_m > 1000) {
paste0("**Atención:** las vías mapeadas más al norte del municipio están ",
format(round(falta_norte_m), big.mark = ","), " m más allá del destino y no están ",
"conectadas al grafo. Con esta capa de vías la ruta no llega a la costa, sino al ",
"último punto conectado con la cabecera.")
} else ""
```
```{r}
#| label: tbl-rutas
#| tbl-cap: "Ruta más corta frente a ruta más transitable, de la cabecera a la costa."
ruta_corta <- st_network_paths(red, from = origen, to = destino, weights = "longitud")
ruta_transi <- st_network_paths(red, from = origen, to = destino, weights = "costo_transito")
e_corta <- aristas[ruta_corta$edge_paths[[1]], ]
e_transi <- aristas[ruta_transi$edge_paths[[1]], ]
d_euclid <- as.numeric(st_distance(nodos[origen, ], nodos[destino, ]))
describir_ruta <- function(e, nombre) {
tibble(ruta = nombre,
longitud_km = sum(e$longitud) / 1000,
transitabilidad = weighted.mean(e$transitabilidad, e$longitud),
pct_brecha = 100 * sum(e$longitud[e$jerarquia == "Brecha o sendero"]) / sum(e$longitud),
circuidad = sum(e$longitud) / d_euclid) # red / línea recta
}
resumen <- bind_rows(describir_ruta(e_corta, "Más corta"),
describir_ruta(e_transi, "Más transitable"))
resumen |>
mutate(across(where(is.numeric), ~ round(.x, 2))) |>
setNames(c("Ruta", "Longitud (km)", "Transitabilidad media", "% del recorrido en brecha",
"Circuidad (red / línea recta)")) |>
kable()
```
`r txt_norte`
```{r}
#| label: fig-rutas
#| fig-cap: "Ruta más corta (blanca, discontinua) y ruta más transitable (celeste) entre la cabecera y la costa, sobre la transitabilidad de cada tramo y la imagen satelital."
pal_transi <- colorNumeric(c("#ff7043", "#ffee58", "#66bb6a"), domain = c(0.25, 1))
mapa_base() |>
addPolylines(data = aristas_ll, color = ~pal_transi(transitabilidad), weight = 1.5, opacity = 0.9,
group = "Transitabilidad de la red") |>
addPolylines(data = aristas_ll[ruta_transi$edge_paths[[1]], ], color = "#00e5ff", weight = 6,
opacity = 0.95, group = "Ruta más transitable") |>
addPolylines(data = aristas_ll[ruta_corta$edge_paths[[1]], ], color = "#ffffff", weight = 3,
dashArray = "8", opacity = 1, group = "Ruta más corta") |>
addMarkers(data = a_ll(nodos[c(origen, destino), ]),
popup = c("Origen: cabecera", "Destino: nodo más al norte del municipio")) |>
addLegend(pal = pal_transi, values = c(0.25, 1), title = "Transitabilidad", position = "bottomleft") |>
control_capas(c("Ruta más corta", "Ruta más transitable", "Transitabilidad de la red"))
```
## Curva de compromiso: ¿cuánto se alarga el viaje por evitar las brechas?
```{r}
#| label: tbl-sensibilidad
#| tbl-cap: "Sensibilidad de la ruta cabecera–costa al parámetro α."
sensibilidad <- lapply(c(0, 0.5, 1, 2, 4, 8), function(a) {
r <- red |>
activate("edges") |>
mutate(c = longitud * (1 + a * (1 - transitabilidad)))
p <- st_network_paths(r, from = origen, to = destino, weights = "c")
e <- aristas[p$edge_paths[[1]], ]
tibble(alfa = a,
longitud_m = sum(e$longitud),
transitabilidad_media = weighted.mean(e$transitabilidad, e$longitud))
}) |>
bind_rows() |>
mutate(desvio_m = longitud_m - first(longitud_m))
sensibilidad |>
mutate(longitud_m = round(longitud_m), desvio_m = round(desvio_m),
transitabilidad_media = round(transitabilidad_media, 3)) |>
setNames(c("α", "Longitud (m)", "Transitabilidad media", "Desvío respecto a la más corta (m)")) |>
kable(format.args = list(big.mark = ","))
```
```{r}
#| label: fig-sensibilidad
#| fig-cap: "Compromiso entre distancia y transitabilidad según α (la etiqueta de cada punto)."
#| fig-height: 4
ggplot(sensibilidad, aes(longitud_m, transitabilidad_media, label = alfa)) +
geom_path(color = "grey60") +
geom_point(size = 3, color = "#00e5ff") +
geom_text(vjust = -1, colour = "white") +
scale_x_continuous(labels = label_number(big.mark = ",")) +
labs(x = "Longitud de la ruta (m)", y = "Transitabilidad media de la ruta") +
tema_oscuro()
```
# Accesibilidad desde la cabecera
```{r}
#| label: accesibilidad
# Costo de llegar desde la cabecera a cada nodo, con los dos pesos
nodos$dist_red <- as.numeric(st_network_cost(red, from = origen, weights = "longitud")[1, ])
nodos$costo_red <- as.numeric(st_network_cost(red, from = origen, weights = "costo_transito")[1, ])
# Enlace de cada píxel del mapa de cambio con la red:
# - la vía mapeada más cercana, conectada o no (distancia y jerarquía);
# - la arista del grafo más cercana (centralidad);
# - el nodo más cercano (por dónde "entra" el píxel a la red).
# Distancia exacta de cada punto a la línea que se le indica (su línea más cercana),
# calculada segmento a segmento con operaciones vectorizadas. Da el mismo resultado
# que st_distance(..., by_element = TRUE), pero en segundos en lugar de minutos.
# 'lineas' debe ser de tipo LINESTRING.
dist_a_linea <- function(xy, lineas, i_linea) {
co <- st_coordinates(lineas) # X, Y y L1 (número de línea)
n <- nrow(co)
seg <- data.frame(x1 = co[-n, "X"], y1 = co[-n, "Y"], x2 = co[-1, "X"], y2 = co[-1, "Y"],
linea = co[-n, "L1"])[co[-n, "L1"] == co[-1, "L1"], ]
seg_por_linea <- split(seq_len(nrow(seg)), seg$linea)
d <- numeric(nrow(xy))
puntos_por_linea <- split(seq_len(nrow(xy)), i_linea)
for (g in names(puntos_por_linea)) {
p <- puntos_por_linea[[g]]; s <- seg[seg_por_linea[[g]], ]
x <- xy[p, 1]; y <- xy[p, 2]
dmin <- rep(Inf, length(p))
for (k in seq_len(nrow(s))) {
vx <- s$x2[k] - s$x1[k]; vy <- s$y2[k] - s$y1[k]; largo2 <- vx^2 + vy^2
# posición (0 a 1) del punto del segmento más cercano a cada píxel
t <- if (largo2 == 0) 0 else pmin(1, pmax(0, ((x - s$x1[k]) * vx + (y - s$y1[k]) * vy) / largo2))
dmin <- pmin(dmin, sqrt((x - (s$x1[k] + t * vx))^2 + (y - (s$y1[k] + t * vy))^2))
}
d[p] <- dmin
}
d
}
xy_px <- st_coordinates(px_pts)
i_via <- st_nearest_feature(px_pts, calles)
px$dist_via <- dist_a_linea(xy_px, calles, i_via)
px$jer_via <- calles$jerarquia[i_via]
principales <- calles[calles$jerarquia == "Principal", ]
px$dist_principal <- if (nrow(principales) > 0) {
dist_a_linea(xy_px, principales, st_nearest_feature(px_pts, principales))
} else NA_real_
px$flujo_via <- percent_rank(aristas$flujo)[st_nearest_feature(px_pts, aristas)] # percentil 0-1
px$nodo <- st_nearest_feature(px_pts, nodos)
xy_nodo <- st_coordinates(nodos)[px$nodo, ]
dist_nodo <- sqrt((xy_px[, 1] - xy_nodo[, 1])^2 + (xy_px[, 2] - xy_nodo[, 2])^2)
# Accesibilidad del píxel = costo por la red hasta su nodo + tramo final fuera de
# la red. En la versión de costo, ese tramo se penaliza como terreno sin camino.
px$acc_dist <- nodos$dist_red[px$nodo] + dist_nodo
px$acc_costo <- nodos$costo_red[px$nodo] + dist_nodo * (1 + ALFA)
```
```{r}
#| label: fig-accesibilidad
#| fig-cap: "Nodos alcanzables desde la cabecera dentro del umbral de distancia, por la red."
alcanzables <- a_ll(dplyr::filter(nodos, dist_red <= UMBRAL_ACCES_M))
pal_dist <- colorNumeric("viridis", domain = c(0, UMBRAL_ACCES_M))
mapa_base() |>
addPolylines(data = vias_ll, color = "#e0e0e0", weight = 1, opacity = 0.7, group = "Red vial") |>
addCircleMarkers(data = alcanzables, radius = 4, stroke = FALSE, fillOpacity = 0.95,
fillColor = ~pal_dist(dist_red),
popup = ~paste0(format(round(dist_red), big.mark = ","), " m por la red"),
group = "Nodos alcanzables") |>
addMarkers(data = a_ll(cabecera), popup = "Cabecera (origen)", group = "Cabecera") |>
addLegend(pal = pal_dist, values = c(0, UMBRAL_ACCES_M), title = "Distancia por la red (m)",
position = "bottomleft") |>
control_capas(c("Nodos alcanzables", "Red vial", "Cabecera"))
```
## Ejercicio 1: accesibilidad con el costo penalizado
> *Repetir el paso 8 con el peso penalizado: ¿cuánto se reduce el área accesible?*
```{r}
#| label: tbl-ejercicio1
#| tbl-cap: "Ejercicio 1. Alcance desde la cabecera con el mismo umbral, midiendo solo distancia o penalizando las vías poco transitables."
ej1 <- tibble(
peso = c("Distancia (longitud)", "Costo de tránsito (penaliza vías poco transitables)"),
nodos_alcanzables = c(sum(nodos$dist_red <= UMBRAL_ACCES_M), sum(nodos$costo_red <= UMBRAL_ACCES_M)),
ha_alcanzables = c(sum(px$ha[px$acc_dist <= UMBRAL_ACCES_M]),
sum(px$ha[px$acc_costo <= UMBRAL_ACCES_M]))
) |>
mutate(pct_municipio = 100 * ha_alcanzables / area_total_ha)
reduccion_ej1 <- 1 - ej1$ha_alcanzables[2] / ej1$ha_alcanzables[1]
ej1 |>
mutate(ha_alcanzables = round(ha_alcanzables), pct_municipio = round(pct_municipio, 1)) |>
setNames(c("Peso", "Nodos alcanzables", "Superficie alcanzable (ha)", "% del municipio")) |>
kable(format.args = list(big.mark = ","))
```
**Respuesta.** Con un umbral de `r format(UMBRAL_ACCES_M, big.mark = ",")` m, penalizar la mala transitabilidad reduce el área accesible desde la cabecera en `r round(100 * reduccion_ej1, 1)` %: de `r format(round(ej1$ha_alcanzables[1]), big.mark = ",")` a `r format(round(ej1$ha_alcanzables[2]), big.mark = ",")` ha. Esa diferencia es el territorio que en el mapa parece cercano pero al que solo se llega por vías de baja jerarquía o caminando fuera de la red.
```{r}
#| label: fig-acc-pixel
#| fig-cap: "Accesibilidad de cada píxel desde la cabecera (costo de tránsito, en km equivalentes de carretera pavimentada), sobre la imagen satelital. Los tonos amarillos son las zonas mal conectadas. La capa se puede apagar para ver el fondo."
pal_acc <- colorNumeric("viridis", domain = range(px$acc_costo / 1000), na.color = "transparent")
mapa_base() |>
addRasterImage(a_raster(px$acc_costo / 1000, "acc"), colors = pal_acc, opacity = 0.75,
group = "Accesibilidad") |>
addPolylines(data = vias_ll, color = "#ffffff", weight = 1, opacity = 0.8, group = "Red vial") |>
addMarkers(data = a_ll(cabecera), popup = "Cabecera (origen)", group = "Cabecera") |>
addLegend(pal = pal_acc, values = range(px$acc_costo / 1000), title = "km equivalentes",
position = "bottomleft") |>
control_capas(c("Accesibilidad", "Red vial", "Cabecera"))
```
# Red vial y patrón de cambio
Esta sección es el cruce con el diagnóstico de NDVI. Para cada zona se calcula un **índice de concentración** de cada categoría de cambio:
$$\text{concentración} = \frac{\text{parte de la categoría que cae en la zona}}{\text{parte del municipio que ocupa la zona}}$$
Un valor de 1 significa que el cambio se reparte igual que el territorio; **mayor que 1**, que ese cambio se concentra en la zona; menor que 1, que la evita. Las categorías con menos de `r MIN_HA_CATEGORIA` ha en todo el municipio aparecen en las tablas pero no en las figuras ni en los hallazgos: con tan poca superficie el índice es ruido.
```{r}
#| label: funciones-cruce
# Superficie, tasa e índice de concentración de cada categoría de cambio por zona
cruce_por_zona <- function(zona) {
d <- data.frame(zona = zona, cambio = px$cambio, ha = px$ha)
totales <- d |> group_by(cambio, .drop = FALSE) |> summarise(ha_cat = sum(ha), .groups = "drop")
d |>
group_by(zona, cambio, .drop = FALSE) |>
summarise(ha = sum(ha), .groups = "drop") |>
group_by(zona) |>
mutate(ha_zona = sum(ha)) |>
ungroup() |>
left_join(totales, by = "cambio") |>
mutate(pct_zona = 100 * ha / ha_zona,
concentracion = (ha / ha_cat) / (ha_zona / area_total_ha))
}
tabla_cruce <- function(cruce, nombre_zona) {
cruce |>
dplyr::filter(cambio != "Persistencia", ha_cat > 0) |>
transmute(zona, ha_zona = round(ha_zona), cambio,
valor = paste0(format(round(ha, 1), big.mark = ",", trim = TRUE),
" (", number(concentracion, 0.01), ")")) |>
pivot_wider(names_from = cambio, values_from = valor) |>
rename(!!nombre_zona := zona, "Superficie de la zona (ha)" = ha_zona) |>
kable(format.args = list(big.mark = ","))
}
figura_cruce <- function(cruce, etiqueta_x) {
cruce |>
dplyr::filter(cambio != "Persistencia", ha_cat >= MIN_HA_CATEGORIA) |>
ggplot(aes(zona, concentracion, colour = cambio, group = cambio)) +
geom_hline(yintercept = 1, linetype = "dashed", colour = "grey70") +
geom_line(linewidth = 0.9) + geom_point(size = 2.5) +
scale_colour_manual(values = paleta_cat, name = NULL) +
labs(x = etiqueta_x, y = "Índice de concentración") +
guides(colour = guide_legend(nrow = 2)) +
tema_oscuro() + theme(legend.position = "bottom")
}
conc <- function(cruce, z, cat) cruce$concentracion[cruce$zona == z & cruce$cambio == cat]
# Lectura continua, que no depende de los cortes de clase: cambio medio de NDVI
# entre la primera y la última fecha en cada zona
delta_por_zona <- function(zona, nombre_zona) {
data.frame(zona = zona, delta = px$delta, ha = px$ha) |>
group_by(zona) |>
summarise(ha = sum(ha), delta_medio = mean(delta), delta_mediano = median(delta),
pct_cae = 100 * mean(delta <= -0.2), .groups = "drop") |>
mutate(ha = round(ha), across(c(delta_medio, delta_mediano), ~ round(.x, 3)),
pct_cae = round(pct_cae, 2)) |>
setNames(c(nombre_zona, "Superficie (ha)", "Cambio medio de NDVI", "Cambio mediano de NDVI",
"% de píxeles con caída de 0.2 o más")) |>
kable(format.args = list(big.mark = ","))
}
```
```{r}
#| label: fig-cambio-red
#| fig-cap: "Categorías de cambio y red vial sobre la imagen satelital. Solo se colorean los píxeles que cambiaron; donde hubo persistencia se ve el fondo. Las vías principales, en blanco grueso."
# La persistencia se deja transparente para que se vea la imagen de fondo
codigo_cat <- as.integer(px$cambio)
codigo_cat[px$cambio == "Persistencia"] <- NA
cat_mapa <- categorias[-1]
pal_cat <- colorFactor(unname(paleta_cat[cat_mapa]), levels = match(cat_mapa, categorias),
na.color = "transparent")
mapa_base() |>
addRasterImage(a_raster(codigo_cat, "categoria"), colors = pal_cat, method = "ngb",
opacity = 0.85, group = "Cambio de vegetación") |>
addPolylines(data = vias_ll, color = "#cfcfcf", weight = 1, opacity = 0.8, group = "Red vial") |>
addPolylines(data = vias_ll[vias_ll$jerarquia == "Principal", ], color = "#ffffff", weight = 4,
opacity = 1, group = "Red vial") |>
addLegend(colors = unname(paleta_cat[cat_mapa]), labels = cat_mapa, opacity = 1,
title = "Cambio de vegetación", position = "bottomleft") |>
control_capas(c("Cambio de vegetación", "Red vial"))
```
## ¿El cambio sigue a los ejes viales?
```{r}
#| label: tbl-bandas
#| tbl-cap: "Cambio según la distancia a la vía más cercana (cualquier jerarquía). En cada celda: hectáreas y, entre paréntesis, índice de concentración."
etiquetas_banda <- paste0(head(BANDAS_VIA_M, -1), "–", tail(BANDAS_VIA_M, -1), " m")
etiquetas_banda[length(etiquetas_banda)] <- paste0("> ", BANDAS_VIA_M[length(BANDAS_VIA_M) - 1], " m")
banda <- function(d) cut(d, BANDAS_VIA_M, right = FALSE, labels = etiquetas_banda)
cruce_via <- cruce_por_zona(banda(px$dist_via))
tabla_cruce(cruce_via, "Distancia a la vía")
```
```{r}
#| label: fig-bandas
#| fig-cap: "Índice de concentración de cada categoría de cambio según la distancia a la vía más cercana. Una curva que baja de izquierda a derecha indica un cambio pegado a las vías."
#| fig-height: 4.5
figura_cruce(cruce_via, "Distancia a la vía más cercana")
```
```{r}
#| label: tbl-delta-bandas
#| tbl-cap: "Cambio de NDVI entre la primera y la última fecha según la distancia a la vía más cercana."
delta_por_zona(banda(px$dist_via), "Distancia a la vía")
```
```{r}
#| label: tbl-bandas-principal
#| tbl-cap: "Cambio según la distancia a la vía PRINCIPAL más cercana (carreteras primarias, secundarias y terciarias). Hectáreas e índice de concentración."
cruce_principal <- cruce_por_zona(banda(px$dist_principal))
tabla_cruce(cruce_principal, "Distancia a la vía principal")
```
```{r}
#| label: tbl-jerarquia-cambio
#| tbl-cap: "Cambio a menos de 250 m de una vía, según la jerarquía de esa vía. Hectáreas e índice de concentración."
cerca <- px$dist_via < 250
zona_jer <- factor(ifelse(cerca, as.character(px$jer_via), "A 250 m o más de cualquier vía"),
c(jerarquias, "A 250 m o más de cualquier vía"))
cruce_jer <- cruce_por_zona(droplevels(zona_jer))
tabla_cruce(cruce_jer, "Vía más cercana")
```
## ¿Las zonas mal conectadas concentran algún tipo de cambio?
El municipio se divide en cinco partes de igual número de píxeles según su accesibilidad desde la cabecera (costo de tránsito): Q1 es el quinto mejor conectado y Q5 el peor.
```{r}
#| label: tbl-quintiles
#| tbl-cap: "Cambio según la accesibilidad desde la cabecera, en quintiles. Hectáreas e índice de concentración."
cortes_q <- quantile(px$acc_costo, 0:5 / 5)
quintil <- cut(px$acc_costo, cortes_q, include.lowest = TRUE,
labels = c("Q1 (mejor conectado)", "Q2", "Q3", "Q4", "Q5 (peor conectado)"))
cruce_acc <- cruce_por_zona(quintil)
tabla_cruce(cruce_acc, "Accesibilidad")
```
```{r}
#| label: fig-quintiles
#| fig-cap: "Índice de concentración de cada categoría de cambio según la accesibilidad desde la cabecera."
#| fig-height: 4.5
figura_cruce(cruce_acc, "Accesibilidad desde la cabecera")
```
```{r}
#| label: tbl-delta-quintiles
#| tbl-cap: "Cambio de NDVI entre la primera y la última fecha según la accesibilidad desde la cabecera."
delta_por_zona(quintil, "Accesibilidad")
```
## Perfil vial de cada categoría de cambio
```{r}
#| label: tbl-perfil
#| tbl-cap: "Perfil vial de cada categoría de cambio (medianas de los píxeles de cada categoría)."
perfil <- px |>
group_by(cambio) |>
summarise(dist_via_m = median(dist_via),
dist_principal_m = median(dist_principal),
acc_km_equiv = median(acc_costo) / 1000,
centralidad = 100 * median(flujo_via),
delta_medio = mean(delta),
ha = sum(ha),
.groups = "drop") |>
relocate(ha, .after = cambio)
perfil |>
mutate(ha = round(ha, 1), across(c(dist_via_m, dist_principal_m), round),
acc_km_equiv = round(acc_km_equiv, 1), centralidad = round(centralidad),
delta_medio = round(delta_medio, 3)) |>
setNames(c("Categoría de cambio", "ha", "Distancia mediana a una vía (m)",
"Distancia mediana a una vía principal (m)",
"Accesibilidad mediana desde la cabecera (km equiv.)",
"Percentil de centralidad de la vía más cercana", "Cambio medio de NDVI")) |>
kable(format.args = list(big.mark = ","))
```
## ¿Las cicatrices se acercan o se alejan de las vías con el tiempo?
Cada cicatriz tiene un año de aparición. Si los desmontes avanzaran a lo largo de los caminos, las cicatrices más recientes quedarían cada vez más lejos de la cabecera pero igual de cerca de una vía; si respondieran a otra causa (inundación, costa), no guardarían relación con la red.
```{r}
#| label: tbl-cicatrices-anio
#| tbl-cap: "Cicatrices candidatas según su año de aparición y su posición respecto a la red vial."
hay_cicatrices <- any(px$cambio == "Cicatriz candidata")
if (hay_cicatrices) {
cic_anio <- px[px$cambio == "Cicatriz candidata", ] |>
group_by(anio = anio_cic) |>
summarise(dist_via_m = median(dist_via),
dist_principal_m = median(dist_principal),
acc_km_equiv = median(acc_costo) / 1000,
pct_a_250 = 100 * mean(dist_via < 250),
ha = sum(ha),
.groups = "drop") |>
relocate(ha, .after = anio)
cic_anio |>
mutate(anio = as.character(anio), ha = round(ha, 1),
across(c(dist_via_m, dist_principal_m), round),
acc_km_equiv = round(acc_km_equiv, 1), pct_a_250 = round(pct_a_250, 1)) |>
setNames(c("Aparece en", "ha", "Distancia mediana a una vía (m)",
"Distancia mediana a una vía principal (m)",
"Accesibilidad mediana desde la cabecera (km equiv.)",
"% a menos de 250 m de una vía")) |>
kable(format.args = list(big.mark = ","))
} else {
cat("El diagnóstico con NDVI no detectó cicatrices candidatas.")
}
```
```{r}
#| label: fig-cicatrices-anio
#| fig-cap: "Distancia a la vía más cercana de las cicatrices candidatas, por año de aparición. El ancho de cada caja es proporcional a la superficie de cicatriz de ese año."
#| fig-height: 4
if (hay_cicatrices) {
ggplot(px[px$cambio == "Cicatriz candidata", ], aes(factor(anio_cic), dist_via)) +
geom_boxplot(varwidth = TRUE, fill = "#d7191c", alpha = 0.6, colour = "grey85",
outlier.size = 0.4, outlier.colour = "grey70") +
geom_hline(yintercept = median(px$dist_via), linetype = "dashed", colour = "grey70") +
scale_y_continuous(labels = label_number(big.mark = ",")) +
labs(x = "Año de aparición", y = "Distancia a la vía más cercana (m)",
caption = "Línea discontinua: mediana de todo el municipio") +
tema_oscuro()
}
```
# Ejercicio 2: ¿las vías principales son también las más centrales?
```{r}
#| label: tbl-ejercicio2
#| tbl-cap: "Ejercicio 2. Intermediación de las aristas según la jerarquía vial."
umbral_top <- quantile(aristas$flujo, 0.9)
ej2 <- aristas |>
st_drop_geometry() |>
group_by(jerarquia) |>
summarise(tramos = n(), km = sum(longitud) / 1000,
flujo_mediano = median(flujo), flujo_medio = mean(flujo),
en_top10 = sum(flujo >= umbral_top), .groups = "drop") |>
mutate(pct_tramos = 100 * tramos / sum(tramos),
pct_top10 = 100 * en_top10 / sum(en_top10))
prueba_kw <- kruskal.test(flujo ~ jerarquia, data = st_drop_geometry(aristas))
rho <- suppressWarnings(cor(aristas$transitabilidad, aristas$flujo, method = "spearman"))
ej2 |>
transmute(jerarquia, tramos, km = round(km, 1), flujo_mediano = round(flujo_mediano),
flujo_medio = round(flujo_medio), pct_tramos = round(pct_tramos, 1),
pct_top10 = round(pct_top10, 1)) |>
setNames(c("Jerarquía", "Tramos", "km", "Intermediación mediana", "Intermediación media",
"% de los tramos", "% del 10 % más central")) |>
kable(format.args = list(big.mark = ","))
```
```{r}
#| label: fig-ejercicio2
#| fig-cap: "Distribución de la intermediación de aristas por jerarquía vial (escala logarítmica)."
#| fig-height: 4
ggplot(st_drop_geometry(aristas), aes(jerarquia, flujo + 1, fill = jerarquia)) +
geom_boxplot(outlier.size = 0.6, show.legend = FALSE, colour = "grey75",
outlier.colour = "grey60") +
scale_fill_manual(values = paleta_jer) +
scale_y_log10(labels = label_number(big.mark = ",")) +
labs(x = NULL, y = "Intermediación de la arista (+1)") +
tema_oscuro()
```
```{r}
#| label: respuesta-ej2
fila_p <- ej2[ej2$jerarquia == "Principal", ]
hay_principal <- nrow(fila_p) == 1
txt_ej2 <- if (!hay_principal) {
"La red no tiene tramos clasificados como principales, así que la comparación no aplica."
} else {
paste0("Las vías principales son el ", round(fila_p$pct_tramos, 1), " % de los tramos y aportan el ",
round(fila_p$pct_top10, 1), " % del 10 % de tramos más centrales. Su intermediación mediana es ",
format(round(fila_p$flujo_mediano), big.mark = ","), ", frente a ",
format(round(median(aristas$flujo[aristas$jerarquia != "Principal"])), big.mark = ","),
" del resto de la red. La diferencia entre jerarquías ",
ifelse(prueba_kw$p.value < 0.05, "es estadísticamente significativa",
"no es estadísticamente significativa"),
" (Kruskal-Wallis, p ", ifelse(prueba_kw$p.value < 0.001, "< 0.001",
paste("=", round(prueba_kw$p.value, 3))),
"), y la correlación de Spearman entre transitabilidad e intermediación es ",
round(rho, 2), ". ",
ifelse(fila_p$pct_top10 > 50,
"Jerarquía y centralidad coinciden en lo esencial: los ejes principales son también los corredores por donde pasaría el flujo.",
"Jerarquía y centralidad coinciden solo en parte: buena parte de los tramos más centrales no son vías principales, sino calles locales o caminos que funcionan como conectores obligados."))
}
```
**Respuesta.** `r txt_ej2`
# Ejercicio 3: Moran con distancia de red frente a distancia euclidiana
> *Se construye una matriz de pesos por distancia de red con `st_network_cost()` y se compara el Moran's I con el euclidiano.*
La unidad de análisis es el **área servida por cada nodo**: los píxeles cuyo nodo más cercano es ese. Se usan dos variables: el **cambio medio de NDVI** de esa área entre la primera y la última fecha, y el porcentaje del área que cambió de clase. El índice de Moran mide si las áreas vecinas se parecen entre sí, y se calcula con `spdep` [@bivand2018]. Se construyen dos matrices de pesos con la misma regla (inverso de la distancia, hasta `r format(ANCHO_MORAN_M, big.mark = ",")` m) y lo único que cambia es cómo se mide la distancia entre nodos: en línea recta o por la red.
```{r}
#| label: ejercicio3
areas_nodo <- px |>
group_by(nodo) |>
summarise(pct_cambio = 100 * sum(ha[cambio != "Persistencia"]) / sum(ha),
delta_medio = mean(delta), # cambio medio de NDVI del área del nodo
ha_nodo = sum(ha),
.groups = "drop") |>
dplyr::filter(ha_nodo >= MIN_HA_NODO)
ids <- areas_nodo$nodo
D_red <- st_network_cost(red, from = ids, to = ids, weights = "longitud")
D_red <- matrix(as.numeric(D_red), nrow = length(ids))
D_euc <- matrix(as.numeric(st_distance(nodos[ids, ])), nrow = length(ids))
# Inverso de la distancia hasta ANCHO_MORAN_M. La distancia se acota por abajo
# en 30 m (un píxel) para que dos nodos casi encimados no dominen los pesos.
matriz_pesos <- function(D) {
W <- ifelse(is.finite(D) & D <= ANCHO_MORAN_M, 1 / pmax(D, 30), 0)
diag(W) <- 0
suppressWarnings(mat2listw(W, style = "W", zero.policy = TRUE))
}
lw_red <- matriz_pesos(D_red)
lw_euc <- matriz_pesos(D_euc)
moran_fila <- function(x, lw, variable, distancia) {
m <- moran.test(x, lw, zero.policy = TRUE)
tibble(variable = variable, distancia = distancia,
I = unname(m$estimate[1]), esperado = unname(m$estimate[2]), p = m$p.value,
vecinos_medios = mean(card(lw$neighbours)),
sin_vecinos = sum(card(lw$neighbours) == 0))
}
ej3 <- bind_rows(
moran_fila(areas_nodo$pct_cambio, lw_euc, "% del área que cambió", "Euclidiana"),
moran_fila(areas_nodo$pct_cambio, lw_red, "% del área que cambió", "Por la red"),
moran_fila(areas_nodo$delta_medio, lw_euc, "Cambio medio de NDVI", "Euclidiana"),
moran_fila(areas_nodo$delta_medio, lw_red, "Cambio medio de NDVI", "Por la red")
)
```
```{r}
#| label: tbl-ejercicio3
#| tbl-cap: "Ejercicio 3. Índice de Moran con pesos por distancia euclidiana y por distancia de red."
ej3 |>
mutate(I = round(I, 3), esperado = round(esperado, 4),
p = ifelse(p < 0.001, "< 0.001", as.character(round(p, 3))),
vecinos_medios = round(vecinos_medios, 1)) |>
setNames(c("Variable", "Distancia", "I de Moran", "I esperado", "p", "Vecinos por nodo (media)",
"Nodos sin vecinos")) |>
kable()
```
```{r}
#| label: respuesta-ej3
I_euc <- ej3$I[1]; I_red <- ej3$I[2]
txt_ej3 <- paste0(
"Con ", length(ids), " áreas de nodo, el I de Moran del porcentaje de cambio es ", round(I_euc, 3),
" con distancia euclidiana y ", round(I_red, 3), " con distancia de red. ",
ifelse(abs(I_red - I_euc) < 0.02,
"Los dos valores son prácticamente iguales: en Ixil, estar cerca en línea recta y estar cerca por la red es casi lo mismo a esta escala, y la red no añade estructura al patrón de cambio.",
ifelse(I_red > I_euc,
"La autocorrelación es mayor por la red: los lugares que se parecen en su cambio son los que están conectados entre sí, no solo los que están cerca. El patrón de cambio está organizado por la red vial.",
"La autocorrelación es mayor en línea recta: el cambio forma manchas continuas en el territorio que no siguen la conectividad vial. El patrón responde más a la vecindad física (tipo de suelo, humedad, cobertura contigua) que a la red.")),
" La matriz por red tiene menos vecinos por nodo (", round(ej3$vecinos_medios[2], 1), " frente a ",
round(ej3$vecinos_medios[1], 1), "), porque la distancia por la red nunca es menor que la recta.")
```
**Respuesta.** `r txt_ej3`
# Caminos e infraestructura eólica
Los caminos de acceso a los aerogeneradores son demasiado angostos para el píxel de 30 m de la serie de NDVI, pero deberían aparecer en la red vial si están mapeados. Aquí se mide cuánta vía mapeada hay dentro del área de influencia eólica, de qué jerarquía, y a qué distancia queda cada aerogenerador de la vía más cercana.
```{r}
#| label: eolico
rutas_eol <- RUTA_EOLICO[file.exists(RUTA_EOLICO)]
hay_eolico <- length(rutas_eol) > 0
if (hay_eolico) {
eolico <- do.call(rbind, lapply(rutas_eol, function(f) {
x <- st_read(f, quiet = TRUE) |> st_make_valid() |> st_transform(CRS_METRICO)
st_sf(archivo = rep(basename(f), nrow(x)), geometry = st_geometry(x))
}))
zona_eol <- st_intersection(st_union(st_buffer(eolico, BUFFER_EOLICO_M)), st_geometry(limite_m))
hay_eolico <- length(zona_eol) > 0 && as.numeric(sum(st_area(zona_eol))) > 0
}
if (hay_eolico) {
area_eol_km2 <- as.numeric(sum(st_area(zona_eol))) / 1e6
dentro <- suppressWarnings(st_intersection(vias_mun[, "jerarquia"], zona_eol))
dentro$km <- as.numeric(st_length(dentro)) / 1000
tab_eol <- tab_jer |>
dplyr::select(jerarquia, km_total = km) |>
left_join(dentro |> st_drop_geometry() |> group_by(jerarquia) |>
summarise(km_dentro = sum(km), .groups = "drop"), by = "jerarquia") |>
mutate(km_dentro = replace_na(km_dentro, 0), km_fuera = km_total - km_dentro)
dens_dentro <- sum(tab_eol$km_dentro) / area_eol_km2
dens_fuera <- sum(tab_eol$km_fuera) / (area_mun_km2 - area_eol_km2)
# Distancia de cada elemento puntual (aerogenerador) a la vía mapeada más cercana
puntos_eol <- eolico[st_geometry_type(eolico) %in% c("POINT", "MULTIPOINT"), ]
puntos_eol <- puntos_eol[lengths(st_intersects(puntos_eol, limite_m)) > 0, ]
dist_aero <- if (nrow(puntos_eol) > 0) {
as.numeric(st_distance(puntos_eol, calles[st_nearest_feature(puntos_eol, calles), ],
by_element = TRUE))
} else numeric(0)
}
```
```{r}
#| label: tbl-eolico
#| tbl-cap: "Vías mapeadas dentro y fuera del área de influencia de la infraestructura eólica, dentro del municipio (km)."
if (hay_eolico) {
tab_eol |>
transmute(jerarquia, km_dentro = round(km_dentro, 1), km_fuera = round(km_fuera, 1)) |>
bind_rows(tibble(jerarquia = "Total", km_dentro = round(sum(tab_eol$km_dentro), 1),
km_fuera = round(sum(tab_eol$km_fuera), 1)),
tibble(jerarquia = "Densidad (km de vía por km²)", km_dentro = round(dens_dentro, 2),
km_fuera = round(dens_fuera, 2))) |>
setNames(c("Jerarquía", "Dentro del área de influencia", "Resto del municipio")) |>
kable()
} else {
cat("No se cargó infraestructura eólica dentro del municipio.")
}
```
```{r}
#| label: tbl-aerogeneradores
#| tbl-cap: "Distancia de los aerogeneradores situados dentro del municipio a la vía mapeada más cercana."
if (hay_eolico && length(dist_aero) > 0) {
tibble(indicador = c("Aerogeneradores dentro del municipio", "Distancia mínima (m)",
"Distancia mediana (m)", "Distancia máxima (m)",
"A menos de 100 m de una vía mapeada"),
valor = c(length(dist_aero), round(min(dist_aero)), round(median(dist_aero)),
round(max(dist_aero)), sum(dist_aero < 100))) |>
kable(col.names = c("Indicador", "Valor"), format.args = list(big.mark = ","))
}
```
```{r}
#| label: fig-eolico
#| fig-cap: "Red vial e infraestructura eólica sobre la imagen satelital. En celeste, el área de influencia dentro del municipio; los puntos son los aerogeneradores y el contorno blanco, los polígonos de parque cargados."
if (hay_eolico) {
tipo_eol <- as.character(st_geometry_type(eolico))
eol_pts <- st_geometry(eolico)[tipo_eol %in% c("POINT", "MULTIPOINT")]
eol_pol <- st_geometry(eolico)[tipo_eol %in% c("POLYGON", "MULTIPOLYGON")]
m <- mapa_base() |>
addPolygons(data = a_ll(zona_eol), stroke = FALSE, fillColor = "#4db8ff", fillOpacity = 0.3,
group = "Área de influencia eólica") |>
addPolylines(data = vias_ll, color = ~pal_jer_ll(jerarquia), opacity = 0.95,
weight = ~ifelse(jerarquia == "Principal", 4, 2),
popup = ~paste0(jerarquia, " (", type, ")"), group = "Red vial") |>
addLegend(pal = pal_jer_ll, values = droplevels(vias_ll$jerarquia), title = "Jerarquía vial",
position = "bottomleft")
if (length(eol_pol) > 0) {
m <- addPolygons(m, data = a_ll(suppressWarnings(st_cast(eol_pol, "MULTIPOLYGON"))), fill = FALSE,
color = "#ffffff", weight = 2, group = "Infraestructura eólica")
}
if (length(eol_pts) > 0) {
m <- addCircleMarkers(m, data = a_ll(suppressWarnings(st_cast(eol_pts, "POINT"))), radius = 5,
color = "#ffffff", fillColor = "#2b83ba", fillOpacity = 1, weight = 1,
group = "Infraestructura eólica")
}
control_capas(m, c("Red vial", "Área de influencia eólica", "Infraestructura eólica"))
}
```
# La red y las cicatrices sobre Google Earth
El mapa reúne la red vial por jerarquía, las cicatrices candidatas del diagnóstico con NDVI (coloreadas por año de aparición), la cabecera usada como origen y la infraestructura eólica, sobre la imagen satelital de Google. Sirve para comprobar a ojo si una cicatriz está junto a un camino y qué es.
```{r}
#| label: fig-mapa-red-cicatrices
#| fig-cap: "Red vial, cicatrices candidatas e infraestructura eólica sobre la imagen satelital de Google."
m <- mapa_base() |>
addPolylines(data = vias_ll, color = ~pal_jer_ll(jerarquia), opacity = 0.9,
weight = ~ifelse(jerarquia == "Principal", 4, 2), popup = ~paste0(jerarquia, " (", type, ")"),
group = "Red vial") |>
addLegend(pal = pal_jer_ll, values = droplevels(vias_ll$jerarquia), title = "Jerarquía vial",
position = "bottomright") |>
addMarkers(data = a_ll(cabecera), popup = "Cabecera (origen de rutas y accesibilidad)",
group = "Cabecera")
grupos_mapa <- c("Red vial", "Cabecera")
if (file.exists(RUTA_CIC_GPKG)) {
cic_ll <- a_ll(st_read(RUTA_CIC_GPKG, quiet = TRUE))
cic_ll$anio_f <- factor(cic_ll$anio, levels = sort(unique(cic_ll$anio)))
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.1,
popup = ~paste0("<b>Cicatriz que aparece en ", anio, "</b><br>", round(ha, 2), " ha"),
group = "Cicatrices candidatas") |>
addLegend(pal = pal_anio, values = cic_ll$anio_f, title = "Cicatriz: aparece en",
position = "bottomleft")
grupos_mapa <- c("Cicatrices candidatas", grupos_mapa)
}
if (hay_eolico && nrow(puntos_eol) > 0) {
m <- addCircleMarkers(m, data = a_ll(suppressWarnings(st_cast(puntos_eol, "POINT"))), radius = 4,
color = "#ffffff", fillColor = "#2b83ba", fillOpacity = 1, weight = 1,
group = "Aerogeneradores")
grupos_mapa <- c(grupos_mapa, "Aerogeneradores")
}
control_capas(m, grupos_mapa)
```
# Validación: índice Kappa {#sec-validacion}
Las secciones anteriores describen dónde ocurre el cambio respecto a la red. Aquí se pone a prueba esa lectura como **modelo**: si la posición de un píxel respecto a la red vial dijera algo sobre la pérdida de vegetación, debería servir para acertar dónde se perdió en una parte del municipio que el modelo no ha visto.
1. **Qué se predice:** si un píxel perdió vegetación entre `r anios[1]` y `r anios[2]` (categorías: `r paste(CATEGORIAS_OBJETIVO, collapse = " y ")`).
2. **Con qué:** tres variables de la red, y nada más: distancia a la vía más cercana, distancia a la vía principal más cercana y accesibilidad desde la cabecera. Se combinan en una regresión logística.
3. **Cómo se valida:** el municipio se divide en bloques de `r format(TAM_BLOQUE_M, big.mark = ",")` m; la mitad, elegida al azar, sirve para calibrar el modelo y la otra mitad para validarlo. Separar por bloques y no por píxeles sueltos evita un acierto inflado: dos píxeles vecinos se parecen tanto que predecir uno con el otro no prueba nada.
4. **Cómo se mide:** con el **índice Kappa**, que indica cuánto del acuerdo entre lo predicho y lo observado supera al esperado por azar, leído con la escala de @landis1977.
Kappa no se reporta solo [@pontius2011]. Se acompaña del desacuerdo de **cantidad** (se predice más o menos pérdida de la que hubo) y de **asignación** (la cantidad es correcta pero está en otro sitio), y de un **modelo nulo** que dice que no se pierde nada: como la pérdida ocupa una parte menor del municipio, ese modelo acierta casi siempre, y su exactitud es la vara con la que hay que medir la del modelo.
```{r}
#| label: validacion-modelo
set.seed(SEMILLA)
# Bloques espaciales: la mitad calibra y la otra mitad valida
xy_m <- st_coordinates(px_pts)
bloque <- paste(floor(xy_m[, 1] / TAM_BLOQUE_M), floor(xy_m[, 2] / TAM_BLOQUE_M))
ids_bloque <- unique(bloque)
en_calib <- bloque %in% sample(ids_bloque, round(length(ids_bloque) / 2))
val <- data.frame(perdio = px$cambio %in% CATEGORIAS_OBJETIVO,
log_via = log1p(px$dist_via), # las distancias, en logaritmo:
log_principal = log1p(px$dist_principal), # importa más pasar de 50 a 500 m que de 5,000 a 5,450
acc_km = px$acc_costo / 1000)
# Calibración: regresión logística solo con los bloques de calibración
modelo <- glm(perdio ~ log_via + log_principal + acc_km, family = binomial, data = val[en_calib, ])
val$prob <- predict(modelo, newdata = val, type = "response")
# Umbral: el que hace que el modelo prediga, en calibración, tanta pérdida como la observada
umbral <- quantile(val$prob[en_calib], 1 - mean(val$perdio[en_calib]))
val$pred <- val$prob >= umbral
evaluar_bin <- function(pred, obs, prob = NULL) {
cm <- table(factor(pred, c(FALSE, TRUE)), factor(obs, c(FALSE, TRUE)))
p <- cm / sum(cm); fila <- rowSums(p); col <- colSums(p)
po <- sum(diag(p)); pe <- sum(fila * col)
cantidad <- sum(abs(fila - col)) / 2
# Área bajo la curva ROC: capacidad de ordenar los píxeles, sin depender del umbral
auc <- if (is.null(prob)) NA else {
n1 <- as.numeric(sum(obs)); n0 <- as.numeric(sum(!obs))
(sum(rank(prob)[obs]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}
data.frame(exactitud = po, kappa = (po - pe) / (1 - pe), cantidad = cantidad,
asignacion = (1 - po) - cantidad,
sensibilidad = cm[2, 2] / sum(cm[, 2]), precision = cm[2, 2] / sum(cm[2, ]), auc = auc)
}
en_valid <- !en_calib
tabla_val <- bind_rows(
"Modelo de red, en calibración" = evaluar_bin(val$pred[en_calib], val$perdio[en_calib], val$prob[en_calib]),
"Modelo de red, en validación" = evaluar_bin(val$pred[en_valid], val$perdio[en_valid], val$prob[en_valid]),
"Modelo nulo (no se pierde nada), en validación" =
evaluar_bin(rep(FALSE, sum(en_valid)), val$perdio[en_valid]),
.id = "modelo")
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: tbl-validacion-bloques
#| tbl-cap: "Reparto del municipio entre calibración y validación."
data.frame(conjunto = ifelse(en_calib, "Calibración", "Validación"), bloque = bloque,
ha = px$ha, perdio = val$perdio) |>
group_by(conjunto) |>
summarise(bloques = n_distinct(bloque), superficie = sum(ha), perdida = sum(ha[perdio]),
.groups = "drop") |>
mutate(pct_perdida = round(100 * perdida / superficie, 1),
superficie = round(superficie), perdida = round(perdida)) |>
setNames(c("Conjunto", "Bloques", "Superficie (ha)", "Con pérdida de vegetación (ha)",
"% con pérdida")) |>
kable(format.args = list(big.mark = ","))
```
```{r}
#| label: tbl-kappa-red
#| tbl-cap: "Validación del modelo de red vial frente a la pérdida de vegetación observada."
tabla_val |>
mutate(across(c(exactitud, cantidad, asignacion, sensibilidad, precision), ~ round(100 * .x, 1)),
across(c(kappa, auc), ~ round(.x, 3)),
lectura = leer_kappa(kappa)) |>
setNames(c("Modelo", "Exactitud global (%)", "Kappa", "Desacuerdo de cantidad (%)",
"Desacuerdo de asignación (%)", "Pérdida real que se detecta (%)",
"Pérdida predicha que es real (%)", "AUC", "Lectura del Kappa")) |>
kable()
```
```{r}
#| label: tbl-confusion-red
#| tbl-cap: "Matriz de confusión en los bloques de validación (ha). Filas: predicho; columnas: observado."
etiq <- c("Sin pérdida", "Con pérdida")
cm_ha <- xtabs(ha ~ pred + obs,
data = data.frame(ha = px$ha[en_valid],
pred = factor(val$pred[en_valid], c(FALSE, TRUE), etiq),
obs = factor(val$perdio[en_valid], c(FALSE, TRUE), etiq)))
names(dimnames(cm_ha)) <- c("Predicho", "Observado")
kable(round(addmargins(cm_ha), 1), format.args = list(big.mark = ","))
```
```{r}
#| label: tbl-coeficientes
#| tbl-cap: "Sentido de cada variable en el modelo. Una razón de probabilidades mayor que 1 indica que, al aumentar la variable, la pérdida de vegetación se vuelve más probable."
coefs <- coef(modelo)[-1]
data.frame(
variable = c("Distancia a la vía más cercana (logaritmo)",
"Distancia a la vía principal más cercana (logaritmo)",
"Accesibilidad desde la cabecera (km equivalentes)"),
razon = round(exp(coefs), 3),
sentido = ifelse(coefs > 0, "más lejos o peor conectado: más pérdida",
"más lejos o peor conectado: menos pérdida")
) |>
kable(col.names = c("Variable", "Razón de probabilidades", "Sentido"), row.names = FALSE)
```
Los coeficientes se leen por su sentido, no por su significancia: con decenas de miles de píxeles vecinos y parecidos entre sí, cualquier valor *p* saldría diminuto sin que eso signifique nada.
```{r}
#| label: fig-mapa-validacion
#| fig-cap: "Probabilidad de pérdida de vegetación según el modelo de red vial y pérdida realmente observada, sobre la imagen satelital. En tono claro, los bloques usados para validar."
pal_prob <- colorNumeric("inferno", domain = c(0, 1), na.color = "transparent")
perdida_obs <- ifelse(val$perdio, 1, NA)
bloques_val <- ifelse(en_valid, 1, NA)
mapa_base() |>
addRasterImage(a_raster(bloques_val, "validacion"),
colors = colorFactor(c("#ffffff", "#ffffff"), levels = 1, na.color = "transparent"),
method = "ngb", opacity = 0.25, group = "Bloques de validación") |>
addRasterImage(a_raster(val$prob, "probabilidad"), colors = pal_prob, opacity = 0.75,
group = "Probabilidad según la red") |>
addRasterImage(a_raster(perdida_obs, "observada"),
colors = colorFactor(c("#00e5ff", "#00e5ff"), levels = 1, na.color = "transparent"),
method = "ngb", opacity = 0.9, group = "Pérdida observada") |>
addPolylines(data = vias_ll, color = "#ffffff", weight = 1, opacity = 0.8, group = "Red vial") |>
addLegend(pal = pal_prob, values = c(0, 1), title = "Probabilidad de pérdida", position = "bottomleft") |>
addLegend(colors = "#00e5ff", labels = "Pérdida observada", opacity = 1, position = "bottomleft") |>
control_capas(c("Pérdida observada", "Probabilidad según la red", "Bloques de validación", "Red vial")) |>
hideGroup("Bloques de validación")
```
```{r}
#| label: lectura-validacion
kv <- tabla_val[2, ]; kc <- tabla_val[1, ]; kn <- tabla_val[3, ]
txt_valid <- paste0(
"Con solo tres variables de la red vial, el modelo predice la pérdida de vegetación en los bloques de validación con un Kappa de ",
fmt(kv$kappa, 3), " (acuerdo ", leer_kappa(kv$kappa), "), frente a ", fmt(kc$kappa, 3),
" en los bloques con que se calibró. ",
ifelse(kc$kappa - kv$kappa > 0.1,
"La caída entre uno y otro indica que parte de lo aprendido es propio de los bloques de calibración y no se traslada al resto del municipio. ",
"Que los dos valores sean parecidos indica que la relación se mantiene en zonas que el modelo no había visto. "),
"Detecta el ", fmt(100 * kv$sensibilidad, 1), " % de la pérdida real, y el ", fmt(100 * kv$precision, 1),
" % de la pérdida que predice es real. Su exactitud global es de ", fmt(100 * kv$exactitud, 1),
" %, frente a ", fmt(100 * kn$exactitud, 1), " % del modelo nulo: ",
ifelse(kv$exactitud > kn$exactitud,
"supera a la hipótesis de que no se pierde nada.",
"no supera en exactitud a la hipótesis de que no se pierde nada, aunque, a diferencia de ella, sí localiza una parte de la pérdida."),
" El error se reparte en ", fmt(100 * kv$cantidad, 1), " % de cantidad y ", fmt(100 * kv$asignacion, 1),
" % de asignación. El área bajo la curva ROC es ", fmt(kv$auc, 3),
" (0.5 equivale a ordenar los píxeles al azar y 1 a ordenarlos sin error).")
```
**Lectura.** `r txt_valid`
Este Kappa valida una **asociación espacial**, no una causa. Un valor alto puede deberse a que la pérdida se concentra en la franja costera, que además es la zona peor conectada: el modelo acertaría por la geografía y no porque los caminos expliquen la pérdida. El resultado depende también del tamaño de los bloques (`TAM_BLOQUE_M`), del sorteo (`SEMILLA`) y del umbral con que la probabilidad se convierte en sí o no; conviene repetirlo con otros valores antes de darlo por firme.
# Principales hallazgos {#sec-hallazgos}
```{r}
#| label: calculos-hallazgos
fmt <- function(x, d = 0) format(round(x, d), big.mark = ",", nsmall = d, trim = TRUE)
pct_jer <- function(j) { p <- tab_jer$pct[tab_jer$jerarquia == j]; if (length(p) == 0) 0 else p }
# Categoría (distinta de persistencia) más concentrada en una zona dada
mas_concentrada <- function(cruce, z) {
d <- cruce |> dplyr::filter(zona == z, cambio != "Persistencia", ha_cat >= MIN_HA_CATEGORIA) |>
arrange(desc(concentracion))
list(cat = as.character(d$cambio[1]), conc = d$concentracion[1], ha = d$ha[1])
}
b1 <- etiquetas_banda[1]; bn <- etiquetas_banda[length(etiquetas_banda)]
cerca_via <- mas_concentrada(cruce_via, b1)
lejos_via <- mas_concentrada(cruce_via, bn)
q5 <- mas_concentrada(cruce_acc, "Q5 (peor conectado)")
q1 <- mas_concentrada(cruce_acc, "Q1 (mejor conectado)")
# Cicatrices candidatas: dónde están respecto a la red y cómo se mueven en el tiempo
cic_total <- sum(px$ha[px$cambio == "Cicatriz candidata"])
txt_cic <- if (hay_cicatrices) {
fila_cic <- perfil[perfil$cambio == "Cicatriz candidata", ]
paste0("Las cicatrices candidatas (", fmt(cic_total, 1), " ha) tienen un índice de concentración de ",
fmt(conc(cruce_via, b1, "Cicatriz candidata"), 2), " en la banda de ", b1, ", de ",
fmt(conc(cruce_via, bn, "Cicatriz candidata"), 2), " en la de ", bn, " y de ",
fmt(conc(cruce_acc, "Q5 (peor conectado)", "Cicatriz candidata"), 2),
" en el quintil peor conectado. Su distancia mediana a una vía es de ",
fmt(fila_cic$dist_via_m), " m, frente a ", fmt(median(px$dist_via)),
" m del conjunto del municipio.")
} else "El diagnóstico con NDVI no detectó cicatrices candidatas."
txt_cic_anio <- if (hay_cicatrices && nrow(cic_anio) >= 3) {
ct <- suppressWarnings(cor.test(px$anio_cic[px$cambio == "Cicatriz candidata"],
px$dist_via[px$cambio == "Cicatriz candidata"], method = "spearman"))
pico <- cic_anio[which.max(cic_anio$ha), ]
paste0("La fecha con más superficie nueva de cicatriz es ", pico$anio, " (", fmt(pico$ha, 1),
" ha), a una distancia mediana de ", fmt(pico$dist_via_m), " m de una vía. La distancia a la vía ",
ifelse(abs(ct$estimate) < 0.1, "no cambia de forma apreciable con el año de aparición",
ifelse(ct$estimate > 0, "aumenta con el año de aparición: las cicatrices recientes están más lejos de las vías",
"disminuye con el año de aparición: las cicatrices recientes están más cerca de las vías")),
" (correlación de Spearman ", fmt(unname(ct$estimate), 2), "). El ",
fmt(100 * mean(px$dist_via[px$cambio == "Cicatriz candidata"] < 250), 1),
" % de la superficie de cicatrices está a menos de 250 m de una vía mapeada.")
} else "No hay suficientes fechas con cicatrices para leer su evolución respecto a la red."
txt_eolico <- if (hay_eolico) {
km_eol <- sum(tab_eol$km_dentro)
paste0(
"Dentro del área de influencia eólica (", fmt(area_eol_km2, 1), " km² del municipio) hay ",
fmt(km_eol, 1), " km de vías mapeadas: ", fmt(dens_dentro, 2), " km/km², frente a ",
fmt(dens_fuera, 2), " km/km² en el resto del municipio.",
if (km_eol > 0) paste0(" El ", fmt(100 * sum(tab_eol$km_dentro[tab_eol$jerarquia %in%
c("Servicio", "Brecha o sendero")]) / km_eol, 1),
" % son caminos de servicio o brechas.") else "",
if (length(dist_aero) > 0) paste0(
" De los ", length(dist_aero), " aerogeneradores situados dentro del municipio, ",
sum(dist_aero < 100), " están a menos de 100 m de una vía mapeada; la distancia mediana es de ",
fmt(median(dist_aero)), " m.",
if (median(dist_aero) > 250) paste0(
" Un aerogenerador construido necesita un camino de acceso a pie de torre, así que esta ",
"distancia indica una de dos cosas: los caminos del parque no están dibujados en ",
"OpenStreetMap") else "") else "")
} else "No se cargó infraestructura eólica."
# Ruta cabecera–costa: ¿hay alternativa a la ruta más corta?
desvio_m <- 1000 * (resumen$longitud_km[2] - resumen$longitud_km[1])
txt_ruta <- if (desvio_m < 1) {
paste0("La ruta más transitable coincide con la más corta: ",
ifelse(resumen$pct_brecha[1] > 0,
"no hay forma de evitar el tramo de brecha sin salir de la red analizada.",
"la ruta más corta ya va por vías transitables."))
} else {
paste0("La ruta más transitable mide ", fmt(resumen$longitud_km[2], 1), " km y recorre un ",
fmt(resumen$pct_brecha[2], 1), " % por brecha: evitar las brechas alarga el viaje ",
fmt(desvio_m), " m.")
}
```
**1. Estructura de la red.** Ixil tiene `r fmt(km_mun, 1)` km de vías mapeadas (`r fmt(km_mun / area_mun_km2, 2)` km por km²). El `r fmt(pct_jer("Principal"), 1)` % son vías principales y el `r fmt(pct_jer("Brecha o sendero"), 1)` % brechas o senderos. El `r fmt(100 * prop_sin_salida, 1)` % de los nodos no tiene salida y el índice alfa es `r fmt(indices[["alfa"]], 2)`: `r ifelse(indices[["alfa"]] < 0.15, "una red en árbol, con muy pocas rutas alternativas", "una red con cierto grado de malla y rutas alternativas")`.
**2. Cabecera–costa.** La ruta más corta mide `r fmt(resumen$longitud_km[1], 1)` km (circuidad `r fmt(resumen$circuidad[1], 2)`) y recorre un `r fmt(resumen$pct_brecha[1], 1)` % por brecha. `r txt_ruta` `r txt_norte`
**3. Accesibilidad (ejercicio 1).** Penalizar la mala transitabilidad reduce el área alcanzable desde la cabecera en `r fmt(100 * reduccion_ej1, 1)` % con el umbral de `r fmt(UMBRAL_ACCES_M)` m.
**4. Cambio y ejes viales.** A menos de `r BANDAS_VIA_M[2]` m de una vía, el cambio más concentrado es *`r cerca_via$cat`* (índice `r fmt(cerca_via$conc, 2)`). En la banda más alejada (`r bn`), el más concentrado es *`r lejos_via$cat`* (índice `r fmt(lejos_via$conc, 2)`).
**5. Zonas mal conectadas.** En el quintil peor conectado, el cambio más concentrado es *`r q5$cat`* (índice `r fmt(q5$conc, 2)`, `r fmt(q5$ha, 1)` ha). En el mejor conectado, *`r q1$cat`* (índice `r fmt(q1$conc, 2)`).
**6. Cicatrices y red vial.** `r txt_cic`
**7. Cicatrices en el tiempo.** `r txt_cic_anio`
**8. Jerarquía y centralidad (ejercicio 2).** `r txt_ej2`
**9. Estructura espacial del cambio (ejercicio 3).** `r txt_ej3`
**10. Infraestructura eólica.** `r txt_eolico`
**11. Validación (índice Kappa).** `r txt_valid`
# Verificaciones {#sec-verificaciones}
```{r}
#| label: tbl-verificaciones
#| tbl-cap: "Controles de consistencia de la lectura de red."
data.frame(
control = c("Fuente de la red vial",
paste0("Tramos y km en la zona (municipio + margen de ", BUFFER_RED_M, " m)"),
"Nodos y aristas tras la limpieza topológica",
"Longitud conservada en el componente principal (%)",
"Píxeles del mapa de cambio enlazados a la red",
"Superficie del mapa de cambio (ha) frente al polígono (ha)",
"Ruta cabecera–costa: metros que faltan hasta la vía mapeada más al norte",
"Píxeles a más de 2 km de cualquier vía (%)"),
resultado = c(fuente_red,
paste0(fmt(nrow(calles)), " tramos; ", fmt(km_bruto, 1), " km"),
paste0(fmt(V), " nodos; ", fmt(E), " aristas"),
fmt(100 * km_limpio / km_bruto, 1),
paste0(fmt(nrow(px)), " de ", fmt(length(celdas))),
paste0(fmt(area_total_ha), " / ", fmt(100 * area_mun_km2)),
fmt(max(falta_norte_m, 0)),
fmt(100 * mean(px$dist_via > 2000), 1))
) |> kable(col.names = c("Control", "Resultado"))
```
Si la longitud conservada baja mucho de 90 %, la capa de vías tiene partes sueltas que el grafo descarta (aparecen punteadas en el mapa de la red). Suele pasar cuando a la capa le faltan los caminos que las unen, por ejemplo las brechas: la accesibilidad, las rutas y la centralidad dejan entonces de representar esas zonas.
# Referencias
Los datos de la red vial son © colaboradores de OpenStreetMap, bajo licencia ODbL.
::: {#refs}
:::
```{r}
#| label: sesion
#| code-fold: true
sessionInfo()
```
