🗺️ Curso GIS en R: Localización Óptima de Vertedero RSU
Manual de análisis espacial paso a paso con Diagnósticos Visuales
Adolfo Diz | Guía de aprendizaje interactiva
2026
📖 1. Introducción y Fundamentos Teóricos
Este manual está diseñado para guiarte en el proceso de toma de decisiones geoespaciales utilizando R. El objetivo es identificar localizaciones idóneas para un vertedero de Residuos Sólidos Urbanos (RSU) en el entorno metropolitano de Zaragoza, basándonos en estrictos criterios normativos, ambientales y geológicos.
A diferencia del software de menús visuales tradicionales, aquí utilizaremos un enfoque basado en código reproducible. La regla de oro de un analista GIS es el diagnóstico visual continuo: nunca guardaremos un dato en el disco duro sin haber representado e inspeccionado previamente su geometría y sus coordenadas en pantalla.
A lo largo de este documento, dejaremos de lado el uso de menús visuales tradicionales para construir un flujo de trabajo reproducible a través de código. Aprenderás a razonar como un programador de sistemas de información geográfica, entendiendo qué ocurre detrás de cada línea de código.
🛠️ FASE 0: Inicialización y Preparación del Entorno
Antes de realizar cualquier operación espacial, es obligatorio estructurar el entorno de RStudio y cargar las herramientas de análisis.
Paso 0.1: Definición de la ruta física de trabajo
¿Por qué hacemos esto?
R necesita saber de qué carpeta de tu disco duro debe leer los archivos
espaciales y dónde debe escribir los resultados que generemos.
Nota para el alumno: Modifica la ruta de abajo para que coincida con la ubicación de tu carpeta local.
Paso 0.2: Instalación y carga de librerías
¿Por qué hacemos esto?
R es un entorno modular. Para habilitar las operaciones de álgebra de
mapas, proyecciones e interactividad, debemos instalar y cargar
librerías externas. La principal es sf (Simple Features),
que dota a R de la capacidad de leer y manipular vectores geométricos de
la misma manera que lo haría una base de datos espacial o un software de
escritorio.
Cargamos la librería sf (Simple Features) para el motor
vectorial, dplyr para gestionar las tablas de atributos,
ggplot2 para mapas estáticos y leaflet para el
visor interactivo web.
# Comprobamos los paquetes instalados y descargamos los que falten
paquetes <- c("sf", "dplyr", "ggplot2", "maptiles", "tidyterra", "leaflet")
instalados <- paquetes %in% installed.packages()[,"Package"]
if(any(!instalados)) {
install.packages(paquetes[!instalados])
}
# Cargamos cada librería en la memoria activa
library(sf) # Motor de análisis vectorial espacial
library(dplyr) # Herramientas de ordenación y filtrado de tablas
library(ggplot2) # Sistema estándar para cartografía estática
library(maptiles) # Motor de descarga de teselas cartográficas
library(tidyterra) # Conector para dibujar coberturas ráster
library(leaflet) # Librería para visores web interactivosPaso 0.3: Configuración geométrica y constantes cartográficas
¿Por qué hacemos esto?
Por defecto, la librería sf calcula distancias y
relaciones utilizando coordenadas esféricas (sobre un globo
tridimensional). Al trabajar con distancias fijas en metros (como
buffers de 500 metros), las matemáticas planas son más precisas y
rápidas. Por ello, desactivamos el motor esférico S2.
Desactivamos el motor esférico S2 para forzar a R a
trabajar con geometría plana bidimensional.
Además, definiremos la constante EPSG:25830 que corresponde a ETRS89 / UTM Huso 30N, el sistema de referencia oficial en la península ibérica que proyecta el globo sobre un plano en metros.
# Forzamos a 'sf' a trabajar con geometría plana bidimensional
sf_use_s2(FALSE)
# Definimos el código de proyección EPSG para nuestro proyecto (UTM Huso 30N)
UTM_30N <- 25830
# Definimos las rutas relativas de lectura y escritura
ruta_datos <- "datos/"
gpkg_analisis <- "vertedero/Analisis/Analisis.gpkg"
# Creamos las carpetas físicas si no existen en nuestro ordenador
dir.create("vertedero/Analisis", recursive = TRUE, showWarnings = FALSE)
dir.create("plots", showWarnings = FALSE)Paso 0.4: Diseño visual para mapas estáticos
¿Por qué hacemos esto?
Para evitar repetir configuraciones estéticas en cada gráfico de
ggplot2, crearemos un objeto de configuración visual con
fuentes ajustadas, fondos limpios y sin ejes de coordenadas visibles (ya
que en los mapas estáticos estos pueden saturar la vista).
# Definición del estilo visual (Theme) de nuestros mapas estáticos en ggplot2
tema_mapa <- theme_minimal() +
theme(
plot.title = element_text(size = 13, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 9, hjust = 0.5, color = "gray40"),
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_line(color = "gray95"),
legend.position = "bottom"
)📥 FASE 1: Carga de Información e Inspección de Formatos
Antes de operar, debemos entender de dónde vienen nuestros archivos y qué contienen.
💡 Concepto clave: El formato histórico .e00 (PAL vs
LAB)
Los archivos con extensión .e00 son formatos de
exportación antiguos de los sistemas de la familia ESRI ArcInfo (el
predecesor de ArcGIS). Este formato empaquetaba diferentes tipos de
información espacial en “subcapas” estructuradas. Al importarlos, R
necesita saber qué subcapa específica queremos leer: *
PAL (Polygon Attribute Table): Almacena la
topología de polígonos. Si vas a cargar geología, parcelas o lagos,
debes indicarle a R que busque la capa PAL. *
LAB (Label Attribute Table): Contiene
etiquetas de puntos. Si vas a cargar pozos, puntos de agua, vértices o
yacimientos puntuales, debes indicarle a R que busque la capa
LAB.
Paso 1.1: Importación de los Puntos de Agua (IPA)
Paso 1.1.1: Inspección de las subcapas del archivo de puntos de agua (IPA)
¿Por qué hacemos esto?
Comprobaremos de forma interactiva qué subcapas componen el archivo
ipa.e00 antes de proceder a cargarlo.
# Consultamos la lista de subcapas almacenadas en el archivo .e00 de puntos de agua
st_layers(paste0(ruta_datos, "ipa.e00"))## Driver: AVCE00
## Available layers:
## layer_name geometry_type features fields crs_name
## 1 LAB Point 45488 8 unnamedPaso 1.1.2: Cargamos la subcapa LAB del archivo de puntos de agua (IPA)
¿Por qué de esta forma?
Puesto que los puntos de agua son elementos puntuales, seleccionamos la
subcapa LAB. Como el archivo original
carece de definición de coordenadas integrada, le asignamos el EPSG de
forma manual con st_set_crs().
Paso 1.2: Importación de la Geología
Paso 1.2.1: Inspección de las subcapas del archivo de geología
¿Por qué hacemos esto?
Haremos la misma comprobación para el archivo de geología.
# Consultamos la estructura del archivo .e00 de geología
st_layers(paste0(ruta_datos, "geologia.e00"))## Driver: AVCE00
## Available layers:
## layer_name geometry_type features fields crs_name
## 1 ARC Line String 397 5 unnamed
## 2 CNT Point 153 1 unnamed
## 3 LAB Point 153 8 unnamed
## 4 PAL Polygon 152 7 unnamedPaso 1.2.2: Cargamos la subcapa PAL del archivo de geología
¿Por qué de esta forma?
Dado que las rocas se distribuyen sobre superficies, seleccionamos la
subcapa de polígonos PAL. Como el archivo
original carece de definición de coordenadas integrada, le asignamos el
EPSG de forma manual con st_set_crs().
# Cargamos la geología seleccionando la subcapa PAL y asignándole el CRS inicial
geologia_raw <- st_read(paste0(ruta_datos, "geologia.e00"), layer = "PAL", quiet = TRUE) %>%
st_set_crs(UTM_30N)(Nota: En la exploración comprobarás que la columna clave para
clasificar el tipo de roca en la capa geológica se denomina
LITOLOGIA y que, efectivamente, sus coordenadas ya se
encuentran proyectadas en metros UTM).
Paso 1.3: Importación de los Espacios Protegidos (LICs)
¿Por qué de esta forma?
Esta información viene almacenada en un archivo de formato Shapefile
clásico (.shp). A diferencia de los archivos
.e00, este formato ya incorpora un archivo de proyección
asociado (.prj). Por tanto, R ya conoce sus coordenadas
originales y podemos utilizar la función st_transform() de
forma segura para proyectar y recalcular el archivo a nuestro huso de
estudio (UTM Huso 30N).
Paso 1.4: Importación de Límites Provinciales
# Cargamos los recintos provinciales oficiales del Instituto Geográfico Nacional
provincias <- st_read(paste0(ruta_datos, "recintos_provinciales_inspire_peninbal_etrs89.shp"), quiet = TRUE)
# Muestra las columnas en vertical con su tipo de datos
glimpse(provincias)## Rows: 51
## Columns: 10
## $ INSPIREID <chr> "ES.IGN.BDDAE.34205400000", "ES.IGN.BDDAE.34195200000", "ES…
## $ COUNTRY <chr> "ES", "ES", "ES", "ES", "ES", "ES", "ES", "ES", "ES", "ES",…
## $ NATLEV <chr> "https://inspire.ec.europa.eu/codelist/AdministrativeHierar…
## $ NATLEVNAME <chr> "Provincia", "Provincia", "Provincia", "Provincia", "Provin…
## $ NATCODE <chr> "34205400000", "34195200000", "34185100000", "34025000000",…
## $ NAMEUNIT <chr> "Territorio no asociado a ninguna provincia", "Melilla", "C…
## $ CODNUT1 <chr> "ES6", "ES6", "ES6", "ES2", "ES4", "ES2", "ES4", "ES5", "ES…
## $ CODNUT2 <chr> "ES64", "ES64", "ES63", "ES24", "ES41", "ES21", "ES41", "ES…
## $ CODNUT3 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-4.297876 3..., MULTIPOLYGON (…## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -7.034043 ymin: 35.17045 xmax: 0.3855578 ymax: 43.45725
## Geodetic CRS: ETRS89
## INSPIREID COUNTRY
## 1 ES.IGN.BDDAE.34205400000 ES
## 2 ES.IGN.BDDAE.34195200000 ES
## 3 ES.IGN.BDDAE.34185100000 ES
## 4 ES.IGN.BDDAE.34025000000 ES
## 5 ES.IGN.BDDAE.34074900000 ES
## 6 ES.IGN.BDDAE.34164800000 ES
## NATLEV
## 1 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
## 2 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
## 3 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
## 4 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
## 5 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
## 6 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
## NATLEVNAME NATCODE NAMEUNIT CODNUT1
## 1 Provincia 34205400000 Territorio no asociado a ninguna provincia ES6
## 2 Provincia 34195200000 Melilla ES6
## 3 Provincia 34185100000 Ceuta ES6
## 4 Provincia 34025000000 Zaragoza ES2
## 5 Provincia 34074900000 Zamora ES4
## 6 Provincia 34164800000 Bizkaia ES2
## CODNUT2 CODNUT3 geometry
## 1 ES64 <NA> MULTIPOLYGON (((-4.297876 3...
## 2 ES64 <NA> MULTIPOLYGON (((-2.952638 3...
## 3 ES63 <NA> MULTIPOLYGON (((-5.382061 3...
## 4 ES24 <NA> MULTIPOLYGON (((0.3224038 4...
## 5 ES41 <NA> MULTIPOLYGON (((-6.10931 41...
## 6 ES21 <NA> MULTIPOLYGON (((-3.089153 4...## [1] "INSPIREID" "COUNTRY" "NATLEV" "NATLEVNAME" "NATCODE"
## [6] "NAMEUNIT" "CODNUT1" "CODNUT2" "CODNUT3" "geometry"🗺️ FASE 2: Delimitación de Áreas mediante Consultas de Atributos; Proyecciones y Disolución
En esta fase utilizaremos las tablas asociadas a las capas geográficas para aislar la provincia de Zaragoza y los municipios que conforman nuestra zona metropolitana de análisis, verificando visualmente cada transformación.
Paso 2.0: Filtrado de la Provincia de Zaragoza
¿Por qué hacemos esto?
La capa provincias que acabamos de importar contiene toda
la división provincial española. Utilizaremos el campo de atributos
NAMEUNIT (que contiene los nombres de las
provincias) para buscar mediante una consulta de texto la palabra
“Zaragoza”. Una vez aislada su geometría, la proyectaremos a nuestro
sistema UTM Huso 30N.
Paso 2.1: Filtrado y Diagnóstico Visual de Cambio de Proyección (CRS)
¿Por qué hacemos esto?
Para entender la diferencia entre un sistema de coordenadas geográfico
(grados decimales) y uno proyectado en metros (UTM), filtraremos la
provincia de Zaragoza de la capa nacional (provincias) y
compararemos visualmente su estado antes y después de transformarla.
💡 Concepto clave: Persistencia del CRS transformado
¿La transformación de coordenadas con st_transform() es
permanente?
Sí, pero solo en el objeto de R resultante. Al
guardarlo en una nueva variable (ej. zaragoza_utm), la
tabla y su columna geométrica recalculan físicamente todos sus valores a
metros. Si posteriormente escribes este objeto con
st_write(), el archivo final quedará proyectado en metros
de forma permanente.
💻 Diagnóstico Visual 2.1: El impacto del CRS en los Ejes de Coordenadas
Representamos la provincia de Zaragoza en pantalla partida. Presta atención a los valores numéricos de los márgenes (ejes X e Y):
# 1. Filtramos Zaragoza en su estado original (Coordenadas geográficas - Grados)
zaragoza_geo <- provincias %>%
filter(grepl("Zaragoza", NAMEUNIT, ignore.case = TRUE))
# 2. Transformamos Zaragoza al sistema proyectado (Metros UTM)
zaragoza_utm <- st_transform(zaragoza_geo, UTM_30N)
# 3. Representamos ambos mapas frente a frente para inspeccionar sus ejes
par(mfrow = c(1, 2)) # Dividimos la pantalla gráfica en 1 fila y 2 columnas
plot(st_geometry(zaragoza_geo), axes = TRUE, main = "Geográficas (Grados Decimales)")
# Nota los valores de los ejes: entre -2° y -0.5° Longitud | 41° y 42.5° Latitud
plot(st_geometry(zaragoza_utm), axes = TRUE, main = "Proyectadas UTM 30N (Metros)")
Paso 2.2: Filtrado de los municipios metropolitanos
¿Por qué hacemos esto?
Cargaremos la capa de municipios de España
(recintos_municipales...) y utilizaremos el campo
NAMEUNIT para filtrar los 7 municipios que
conforman el área de estudio metropolitana, transformándolos
directamente a UTM.
# Cargamos el archivo de límites municipales de España
ttmm <- st_read(paste0(ruta_datos, "recintos_municipales_inspire_peninbal_etrs89.shp"), quiet = TRUE)
# Definimos el listado de municipios objetivo
municipios_objetivo <- c("Zaragoza", "Utebo", "Pastriz", "La Puebla de Alfindén",
"Villanueva de Gállego", "Villamayor de Gállego", "Cuarte de Huerva", "Cadrete")
# Filtramos los 7 municipios seleccionados y los proyectamos
ZONA <- ttmm %>%
filter(NAMEUNIT %in% municipios_objetivo) %>%
st_transform(UTM_30N)Paso 2.3: Creación de la capa de Ámbito Unificado (Dissolve / Union)
¿Por qué hacemos esto?
La variable ZONA almacena 7 polígonos independientes (uno
por municipio). Para nuestro análisis ambiental, necesitamos fusionar
todos estos términos municipales en una sola geometría exterior
unificada que actúe como nuestra máscara de análisis.
💡 Concepto clave: ¿Es un st_union() lo mismo que un
Dissolve?
- En ArcGIS: El comando
Dissolvetoma varios polígonos adyacentes y elimina los límites comunes entre ellos. Si no se selecciona un campo de atributos de agrupación, el resultado es un único polígono continuo. - En R (
sf): La funciónst_union()realiza exactamente este proceso físico sobre las geometrías. Al combinarlo const_sf(), reestructuramos el resultado de geometría pura a una tabla espacial estructurada, creando el objetoAMBITO.
💻 Diagnóstico Visual 2.3: Verificación del proceso de Disolución
Antes de guardar la información, comprobaremos si la función
st_union() ha eliminado correctamente los límites políticos
internos de los 7 términos municipales:
# Creamos la capa unificada aplicando la fusión espacial de límites
AMBITO <- ZONA %>%
st_union() %>%
st_sf()
# Representamos ambas capas frente a frente para comprobar el resultado
par(mfrow = c(1, 2))
plot(st_geometry(ZONA), col = "lightblue", border = "blue", main = "ZONA: 7 Municipios Independientes")
# Se aprecian las líneas divisorias internas de los términos municipales
plot(st_geometry(AMBITO), col = "lightgreen", border = "darkgreen", lwd = 2, main = "ÁMBITO: Polígono Unificado (Dissolve)")
Una vez que hemos verificado visualmente que ambas capas se han generado de forma correcta, procedemos a guardarlas en nuestra base de datos local:
# Guardamos los límites individuales de los municipios en el GeoPackage
st_write(ZONA, gpkg_analisis, layer = "ZONA", delete_layer = TRUE, quiet = TRUE)
# Guardamos el polígono de máscara unificado (Ámbito) en el GeoPackage
st_write(AMBITO, gpkg_analisis, layer = "AMBITO", delete_layer = TRUE, quiet = TRUE)
cat("✅ Capas 'ZONA' y 'AMBITO' verificadas y escritas correctamente en disco.")## ✅ Capas 'ZONA' y 'AMBITO' verificadas y escritas correctamente en disco.# Filtrado de Zaragoza para el mapa estático
zaragoza <- provincias %>%
filter(grepl("Zaragoza", NAMEUNIT, ignore.case = TRUE)) %>%
st_transform(UTM_30N)
# Generación del mapa de delimitación
ggplot() +
geom_sf(data = zaragoza, fill = "gray93", color = "gray75") +
geom_sf(data = ZONA, fill = "#A6CEE3", color = "#1F78B4", size = 0.5) +
geom_sf(data = AMBITO, fill = NA, color = "red", size = 1.2, linetype = "dashed") +
labs(
title = "Delimitación del Área de Estudio",
subtitle = "Municipios metropolitanos seleccionados y contorno unificado (ÁMBITO)"
) +
tema_mapa
✂️ FASE 3: Recorte Espacial (Clip) e Integridad Geométrica
Nuestras capas ambientales originales abarcan áreas regionales muy
extensas. Debemos recortarlas para quedarnos solo con lo que se sitúa
físicamente dentro de nuestro AMBITO.
Paso 3.1: Saneamiento de campos estructurados como listas en geología
# Detectamos qué columnas de la tabla de atributos de geología son de tipo 'list'
cols_lista <- names(geologia_raw)[sapply(geologia_raw, is.list)]
# Si existen, las eliminamos para evitar errores de escritura en el GeoPackage
if(length(cols_lista) > 0) {
geologia_raw <- geologia_raw %>% select(-all_of(cols_lista))
}Paso 3.2: Declaración de nuestra función de recorte (Clip)
¿Por qué hacemos esto?
Para no repetir código, crearemos la función recortar().
Esta función elimina las dimensiones Z o M sobrantes de los datos con
st_zm(), repara geometrías corruptas con
st_make_valid() y ejecuta el recorte espacial estricto con
st_intersection() (el equivalente al comando
Clip en ArcGIS).
# Programamos nuestra función para un procesamiento espacial seguro
recortar <- function(capas, ambito, nombre_salida) {
capas <- st_zm(capas, drop = TRUE, what = "ZM")
if (!all(st_is_valid(capas))) {
capas <- st_make_valid(capas)
}
clip <- st_intersection(capas, ambito)
st_write(clip, gpkg_analisis, layer = nombre_salida, delete_layer = TRUE, quiet = TRUE)
return(clip)
}Paso 3.3: Ejecución y Diagnóstico Visual del Recorte de Geología
¿Por qué hacemos esto?
Aplicaremos la función de recorte a la capa de geología. Antes de
continuar, comprobaremos si la geología regional se ha limitado de forma
estricta a la forma de nuestro ámbito.
💻 Diagnóstico Visual 3.3: Inspección de Geología antes y después del Recorte
Representamos la geología original frente al resultado de la función
recortar(). Añadiremos el límite del AMBITO en
rojo para comprobar la precisión del ajuste:
# Ejecutamos el recorte de la capa de geología regional
geologia_clip <- recortar(geologia_raw, AMBITO, "geologia_clip")
# Representamos ambos estados en pantalla partida para verificar el ajuste espacial
par(mfrow = c(1, 2))
plot(st_geometry(geologia_raw), axes = TRUE, border = "gray80", main = "Geología Regional (Original)")
plot(st_geometry(AMBITO), border = "red", lwd = 2, add = TRUE) # Superponemos el límite de estudio
plot(st_geometry(geologia_clip), axes = TRUE, main = "Geología Acotada (Clipped)")
plot(st_geometry(AMBITO), border = "red", lwd = 1.5, add = TRUE) # Comprobamos el ajuste milimétrico
Una vez verificado que la geología se ajusta al molde de nuestro ámbito de estudio, recortamos el resto de las capas:
# Recortamos los espacios protegidos (LICs) e infraestructuras de agua (IPA) al Ámbito
lics_clip <- recortar(lics, AMBITO, "lics_clip")
ipa_clip <- recortar(ipa, AMBITO, "ipa_clip")
cat("✅ Procesamiento de recortes completado y verificado.")## ✅ Procesamiento de recortes completado y verificado.📝 FASE 4: Reclasificación de Atributos y Disolución de Formaciones
En esta fase depuraremos la tabla de atributos de la geología y la simplificaremos basándonos en su grado de impermeabilidad.
Paso 4.1: Corrección litológica y asignación de permeabilidad
¿Por qué hacemos esto?
La tabla original presenta nombres de rocas inconsistentes. Utilizaremos
mutate() y case_when() para unificar los
términos. Posteriormente, asignaremos de forma binaria el parámetro de
permeabilidad: NO (Sustrato impermeable -
Apto para albergar el vertedero) o SI
(Sustrato permeable - Excluido del análisis).
El archivo original presentaba inconsistencias textuales (por ejemplo, polígonos intermedios denominados “gravas y calizas” que generaban confusión de clasificación).
Utilizaremos la función case_when() de
dplyr (el equivalente directo y moderno a la calculadora de
campos o Field Calculator de ArcGIS con scripts en Python) para
reclasificar la litología bajo criterios limpios y asignarles de forma
binaria el parámetro de permeabilidad (SI o
NO).
# Corregimos los nombres de las rocas y asignamos la impermeabilidad de forma binaria
geologia_edit <- geologia_clip %>%
mutate(LITOLOGIA = case_when(
grepl("gravas", LITOLOGIA, ignore.case = TRUE) & grepl("calizas", LITOLOGIA, ignore.case = TRUE) ~ "Gravas",
grepl("lutitas", LITOLOGIA, ignore.case = TRUE) & grepl("areniscas", LITOLOGIA, ignore.case = TRUE) ~ "Lutitas",
grepl("margas", LITOLOGIA, ignore.case = TRUE) & grepl("areniscas", LITOLOGIA, ignore.case = TRUE) ~ "Margas",
grepl("margas", LITOLOGIA, ignore.case = TRUE) & grepl("lutitas", LITOLOGIA, ignore.case = TRUE) ~ "Margas",
TRUE ~ LITOLOGIA
)) %>%
mutate(PERME = case_when(
LITOLOGIA %in% c("Gravas", "Calizas", "Areniscas") ~ "SI",
LITOLOGIA %in% c("Yesos", "Margas", "Lutitas") ~ "NO",
TRUE ~ "INDEFINIDO"
))Paso 4.2: Disolución espacial por campos de atributos (Dissolve)
¿Por qué hacemos esto?
Para simplificar el análisis cartográfico, colapsaremos todas las
fronteras geográficas de los polígonos que compartan la misma litología
e impermeabilidad. Agruparemos los registros por PERME y
LITOLOGIA utilizando la tubería
group_by() %>% summarize().
💻 Diagnóstico Visual 4.2: Comparación de Estructura de Atributos y Geometría
Para entender qué ha hecho esta disolución por atributos, compararemos la capa segmentada original frente a la capa simplificada por permeabilidad:
# Ejecutamos la disolución espacial por atributos
geologia_perme <- geologia_edit %>%
group_by(PERME, LITOLOGIA) %>%
summarize() %>%
ungroup()
# Visualizamos ambas capas comparando la distribución de sus atributos
par(mfrow = c(1, 2))
# Mapa 1: Distribución litológica original
plot(geologia_edit["LITOLOGIA"], border = NA, main = "Litología Detallada (Segmentada)")
# Mapa 2: Distribución por permeabilidad (Aptitud Geológica)
plot(geologia_perme["PERME"], border = "gray40", main = "Permeabilidad Unificada (SI / NO)")
Una vez comprobado que la capa de geología muestra de manera clara los polígonos permeables (rojos) e impermeables (verdes), procedemos a guardarla en disco:
# Guardamos la geología simplificada en nuestro GeoPackage
st_write(geologia_perme, gpkg_analisis, layer = "geologia_perme", delete_layer = TRUE, quiet = TRUE)
cat("✅ Capa 'geologia_perme' validada y escrita correctamente.")## ✅ Capa 'geologia_perme' validada y escrita correctamente.⭕ FASE 5: Áreas de Influencia y Amortiguación (Buffers)
Estableceremos las distancias de separación obligatorias alrededor de los elementos geográficos protegidos o vulnerables.
Paso 5.1: Carga y adecuación de elementos auxiliares
# Cargamos los ríos, núcleos y yacimientos, proyectándolos y recortándolos al Ámbito en un solo paso
rios <- st_read(paste0(ruta_datos, "rios50.shp"), quiet = TRUE) %>%
st_transform(UTM_30N) %>%
st_intersection(AMBITO)
nucleos <- st_read(paste0(ruta_datos, "nucleos.shp"), quiet = TRUE) %>%
st_transform(UTM_30N) %>%
st_intersection(AMBITO)
arqueo <- st_read(paste0(ruta_datos, "arqueo.shp"), quiet = TRUE) %>%
st_transform(UTM_30N) %>%
st_intersection(AMBITO)Paso 5.2: Cálculo de Buffers multidistancia
Generaremos las zonas de amortiguación radiales de seguridad en metros planos sobre cada una de nuestras capas críticas.
# Calculamos las áreas de influencia (Buffers)
rios_buf <- st_buffer(rios, 500) # 500 metros alrededor de ríos
ipa_buf <- st_buffer(ipa_clip, 100) # 100 metros alrededor de puntos de agua
nucleos_buf <- st_buffer(nucleos, 1000) # 1000 metros alrededor de núcleos urbanos
lics_buf <- st_buffer(lics_clip, 1500) # 1500 metros alrededor de espacios protegidos (LICs)
arqueo_buf <- st_buffer(arqueo, 500) # 500 metros alrededor de yacimientos arqueológicos💻 Diagnóstico Visual 5.3: El concepto geométrico del Buffer
¿Por qué hacemos esto?
Para que comprendas visualmente cómo actúa un buffer, representaremos de
forma detallada un elemento lineal (los ríos con su buffer de 500m) y un
elemento puntual (los núcleos urbanos con su buffer de 1000m):
# Visualizamos el efecto geométrico de amortiguación (halo) de los buffers
par(mfrow = c(1, 2))
# 1. Detalle del buffer lineal (Ríos)
plot(st_geometry(rios_buf), col = "#A6CEE3", border = "#1F78B4", main = "Amortiguación Fluvial (500m)")
plot(st_geometry(rios), col = "blue", lwd = 1.5, add = TRUE) # Superponemos el río original
# 2. Detalle del buffer puntual (Poblaciones)
plot(st_geometry(nucleos_buf), col = "#FDBF6F", border = "#FF7F00", main = "Amortiguación Urbana (1000m)")
plot(st_geometry(nucleos), col = "red", pch = 20, cex = 1.5, add = TRUE) # Superponemos los puntos originales
Una vez que hemos verificado que los buffers se adaptan correctamente a la forma de sus geometrías base, los guardamos en el GeoPackage:
st_write(rios_buf, gpkg_analisis, "rios_buf", delete_layer = TRUE, quiet = TRUE)
st_write(ipa_buf, gpkg_analisis, "ipa_buf", delete_layer = TRUE, quiet = TRUE)
st_write(nucleos_buf, gpkg_analisis, "nucleos_buf", delete_layer = TRUE, quiet = TRUE)
st_write(lics_buf, gpkg_analisis, "lics_buf", delete_layer = TRUE, quiet = TRUE)
st_write(arqueo_buf, gpkg_analisis, "arqueo_buf", delete_layer = TRUE, quiet = TRUE)
cat("✅ Todos los polígonos de amortiguación han sido validados y exportados.")## ✅ Todos los polígonos de amortiguación han sido validados y exportados.🧩 FASE 6: Modelo de Superposición Espacial (Álgebra de Mapas)
En esta fase combinaremos los polígonos de exclusión y los restaremos de nuestro territorio para obtener las parcelas candidatas idóneas.
Paso 6.0: Visualización de las Zonas de Exclusión Acumuladas
Antes de consolidar las restricciones en un único objeto, podemos analizar espacialmente el solapamiento de todas las áreas de seguridad calculadas en la fase anterior.
ggplot() +
geom_sf(data = AMBITO, fill = "white", color = "black", size = 1) +
geom_sf(data = lics_buf, fill = "#33A02C", alpha = 0.35, color = NA) +
geom_sf(data = nucleos_buf, fill = "#E31A1C", alpha = 0.35, color = NA) +
geom_sf(data = rios_buf, fill = "#1F78B4", alpha = 0.35, color = NA) +
geom_sf(data = arqueo_buf, fill = "#FF7F00", alpha = 0.35, color = NA) +
geom_sf(data = ipa_buf, fill = "#A6CEE3", alpha = 0.5, color = NA) +
labs(
title = "Zonas de Exclusión Acumuladas",
subtitle = "Polígonos de afección y amortiguación legalmente restringidos"
) +
tema_mapa
Paso 6.1: Consolidación de restricciones y diagnóstico de la mancha de exclusión
¿Por qué hacemos esto?
Uniremos todos los buffers en una única mancha de restricción mediante
st_union(). Posteriormente, recortaremos esa mancha de
exclusión al contorno de nuestro ámbito con
st_intersection().
💻 Diagnóstico Visual 6.1: Inspección de la superficie prohibida consolidada
Representaremos de forma conjunta la silueta del AMBITO
de estudio y superpondremos en color naranja la mancha de restricciones
consolidada. Esto nos permitirá comprobar de un vistazo qué zonas del
mapa quedan vetadas por normativa ambiental:
# 1. Unificamos todos los buffers en una única cobertura geométrica (Union)
exclusion_union <- st_union(c(st_geometry(rios_buf), st_geometry(ipa_buf), st_geometry(nucleos_buf),
st_geometry(lics_buf), st_geometry(arqueo_buf)))
# Estructuramos la geometría resultante y definimos su nivel de Aptitud = 2 (Zona Excluida)
RESULT1 <- st_sf(Aptitud = 2, geometry = st_sfc(exclusion_union, crs = st_crs(AMBITO)))
AMBITO <- AMBITO %>% mutate(Aptitud = 1) # Ámbito inicial con Aptitud = 1 (Zona Viable)
# 2. Recortamos las restricciones consolidadas a los límites del Ámbito (Clip)
RESULT2 <- st_intersection(RESULT1, st_geometry(AMBITO)) %>% st_collection_extract("POLYGON")
# Representamos la distribución territorial de las restricciones en nuestra zona de estudio
plot(st_geometry(AMBITO), border = "black", lwd = 2, main = "Superficie de Restricción Consolidada (RESULT2)")
plot(st_geometry(RESULT2), col = "orange", border = "darkorange", add = TRUE)
Paso 6.2: Resta de áreas de exclusión y cruce geológico definitivo
¿Por qué hacemos esto?
Sustraeremos la mancha de restricciones de nuestro territorio completo
mediante st_difference(), obteniendo la capa
RESULT3 (Zonas viables por distancia).
Finalmente, cruzaremos esas áreas viables con las zonas donde la
geología es impermeable (GEOIMPER) mediante
st_intersection(), aislando las candidatas óptimas
definitivas.
💻 Diagnóstico Visual 6.2: El flujo deductivo del Análisis de Superposición
Para verificar si la lógica espacial del algoritmo es coherente, representaremos el proceso de deducción territorial en tres pasos: 1. RESULT3: Zonas que cumplen con todas las distancias mínimas (zonas libres de restricciones). 2. GEOIMPER: Zonas donde el terreno es de material arcilloso o impermeable. 3. ZONAS_APTAS: Las parcelas finales óptimas obtenidas tras cruzar ambas condiciones (intersección espacial).
# 1. Restamos las restricciones del ámbito (Erase / Difference)
RESULT3 <- st_difference(AMBITO, st_geometry(RESULT2)) %>% st_collection_extract("POLYGON")
# 2. Aislamos la geología impermeable (Apta)
GEOIMPER <- geologia_perme %>% filter(PERME == "NO")
# 3. Calculamos la intersección final (Cruce óptimo)
ZONAS_APTAS <- st_intersection(RESULT3, GEOIMPER) %>%
st_collection_extract("POLYGON") %>%
select(LITOLOGIA, PERME, Aptitud)
# Representamos de forma secuencial la deducción para comprobar la coherencia espacial
par(mfrow = c(1, 3))
plot(st_geometry(RESULT3), col = "yellow", main = "1. Viable por Distancias (RESULT3)")
plot(st_geometry(GEOIMPER), col = "lightgreen", main = "2. Geología Impermeable (GEOIMPER)")
plot(st_geometry(ZONAS_APTAS), col = "red", main = "3. Candidatos Óptimos (Intersección)")
Una vez que hemos verificado que los polígonos finales (rojos) se sitúan estrictamente dentro de las zonas libres de restricciones (amarillas) y sobre terreno impermeable (verdes), los escribimos en nuestra base de datos:
st_write(RESULT1, gpkg_analisis, "RESULT1", delete_layer = TRUE, quiet = TRUE)
st_write(RESULT2, gpkg_analisis, "RESULT2", delete_layer = TRUE, quiet = TRUE)
st_write(RESULT3, gpkg_analisis, "RESULT3", delete_layer = TRUE, quiet = TRUE)
st_write(GEOIMPER, gpkg_analisis, "GEOIMPER", delete_layer = TRUE, quiet = TRUE)
st_write(ZONAS_APTAS, gpkg_analisis, "ZONAS_APTAS", delete_layer = TRUE, quiet = TRUE)
cat(sprintf("✅ Análisis de superposición finalizado y validado. Se han guardado %d polígonos candidatos.", nrow(ZONAS_APTAS)))## ✅ Análisis de superposición finalizado y validado. Se han guardado 2 polígonos candidatos.🔍 Chunk de Exploración y Deducción 6
🧠 Deducción de resultados: Observa cómo se ha ido reduciendo la superficie disponible.
RESULT3contiene amplias zonas aptas si solo atendiéramos a las distancias a carreteras, núcleos urbanos o ríos. Sin embargo, al cruzar esa información con las arcillas y yesos (GEOIMPER), obtenemos las verdaderas áreas seguras representadas enZONAS_APTAS.
# Comparación del recuento de elementos candidatos
nrow(RESULT3) # Zonas aptas exclusivamente por distancias de seguridad
nrow(ZONAS_APTAS) # Zonas definitivas óptimas (distancia idónea + geología impermeable)Generemos la cartografía estática final del análisis de superposición ponderada:
ggplot() +
geom_sf(data = AMBITO, fill = "gray90", color = "gray60") +
geom_sf(data = RESULT3, fill = "#FEE08B", alpha = 0.7, color = NA) +
geom_sf(data = ZONAS_APTAS, fill = "#D73027", color = "black", size = 0.5) +
labs(
title = "Resultados del Análisis de Localización Óptima",
subtitle = "Amarillo: Zonas aptas por distancias de seguridad | Rojo: Polígonos de Aptitud Óptima (Sustrato Impermeable)"
) +
tema_mapa
🛰️ FASE 7: Cartografía Estática sobre Imagen Satelital
Pintaremos nuestro resultado final sobre una ortofoto aérea real obtenida mediante sensores de satélite.
Paso 7.1: Reproyección a sistema de coordenadas geográficas (WGS84)
¿Por qué hacemos esto?
Las coberturas y servidores de mapas web consumen imágenes estructuradas
bajo el sistema global WGS84 (Latitud/Longitud -
EPSG:4326). Por tanto, reproyectamos todas nuestras capas de
análisis antes de solicitar la descarga de imágenes.
Paso 7.2: Descarga de teselas de satélite
Paso 7.3: Renderización del mapa estático final mediante
ggplot2
¿Por qué de esta forma?
Superpondremos la imagen de satélite como base, dibujaremos el perímetro
del ámbito de estudio en blanco, las zonas viables genéricas en color
naranja y los candidatos óptimos finales en rojo semitransparente.
# Generamos la composición cartográfica estática final
p_final <- ggplot() +
geom_spatraster_rgb(data = satelite) +
geom_sf(data = AMBITO_4326, fill = NA, color = "white", linewidth = 0.8) +
geom_sf(data = RESULT3_4326, fill = "#FF8C00", alpha = 0.35, color = NA) +
geom_sf(data = ZONAS_APTAS_4326, fill = "#D73027", alpha = 0.7, color = "white", size = 0.4) +
labs(
title = "Localización Óptima de Vertedero RSU",
subtitle = "Filtros de exclusión legal y aptitud litológica sobre ortofoto satelital",
caption = "Fuentes de información cartográfica: Esri World Imagery & Datos Abiertos"
) +
tema_mapa +
theme(panel.background = element_rect(fill = "black"))
# Visualizamos la composición gráfica en pantalla
print(p_final)
🗺️ FASE 8: Compilación de Visores Web Interactivos (Leaflet)
Paso 8.1: Extracción de límites de encuadre (Bounding Box)
Paso 8.2: Construcción y visualización del mapa web interactivo
¿Por qué de esta forma?
Utilizaremos la librería leaflet para integrar el mapa
interactivo final. Añadiremos mapas base intercambiables, definiremos
los polígonos de análisis y configuraremos ventanas emergentes dinámicas
(popups) que mostrarán las propiedades del terreno al hacer clic en
pantalla.
# Renderizamos el visor web interactivo final
mapa_interactivo <- leaflet(height = 650, options = leafletOptions(preferCanvas = TRUE)) %>%
# Añadimos mapas base intercambiables
addTiles(group = "Callejero (OpenStreetMap)") %>%
addProviderTiles(providers$CartoDB.Positron, group = "Mapa Claro (CartoDB)") %>%
addProviderTiles(providers$Esri.WorldImagery, group = "Satélite (ESRI)") %>%
# Añadimos la capa perimetral del Ámbito de estudio
addPolygons(
data = AMBITO_4326,
group = "Ámbito de Estudio",
color = "black", weight = 2.5, fillOpacity = 0.02, dashArray = "5, 5"
) %>%
# Añadimos las superficies viables por distancia de seguridad
addPolygons(
data = RESULT3_4326,
group = "Aptas por Distancia",
color = "#FF8C00", weight = 1, fillOpacity = 0.3, fillColor = "yellow",
popup = "<b>Zona Viable por Distancias</b><br>Cumple con todas las separaciones de seguridad."
) %>%
# Añadimos los polígonos correspondientes a los candidatos definitivos óptimos
addPolygons(
data = ZONAS_APTAS_4326,
group = "Zonas Óptimas (Impermeables)",
color = "white", weight = 1.5, fillOpacity = 0.75, fillColor = "#D73027",
popup = ~paste0("<b>💎 CANDIDATO ÓPTIMO DETECTADO</b><hr>",
"<b>Tipo de Roca:</b> ", LITOLOGIA, "<br>",
"<b>Impermeable:</b> ", PERME, "<br>",
"<b>Aptitud Geotécnica:</b> Excelente para RSU")
) %>%
# Añadimos el panel de control de capas en pantalla
addLayersControl(
baseGroups = c("Mapa Claro (CartoDB)", "Callejero (OpenStreetMap)", "Satélite (ESRI)"),
overlayGroups = c("Ámbito de Estudio", "Aptas por Distancia", "Zonas Óptimas (Impermeables)"),
options = layersControlOptions(collapsed = FALSE)
) %>%
# Forzamos al visor interactivo a encuadrarse sobre nuestro territorio
fitBounds(
lng1 = bb[["xmin"]],
lat1 = bb[["ymin"]],
lng2 = bb[["xmax"]],
lat2 = bb[["ymax"]]
) %>%
# Agregamos la leyenda explicativa
addLegend(
position = "bottomleft",
colors = c("yellow", "#D73027"), #c("#FF8C00", "#D73027"),
labels = c("Apta por distancia de seguridad", "Óptima (Sustrato impermeable)"),
opacity = 0.8,
title = "Candidatas RSU"
)
# Mostramos el mapa web resultante
mapa_interactivo📚 Resumen de Equivalencias de Comandos: ArcGIS vs R
sf
Para facilitar tu transición al desarrollo espacial con código, consulta este cuadro que muestra la correspondencia entre las herramientas típicas de ArcGIS y los comandos de R:
| Herramienta clásica en ArcGIS | Función correspondiente en R (sf /
dplyr) |
Explicación conceptual del proceso |
|---|---|---|
| Add Field / Calculate Field | mutate() + case_when() |
Crea nuevas columnas o recalcula datos de las filas basadas en reglas lógicas. |
| Select / Select By Attributes | filter() |
Filtra filas de la tabla de datos espaciales según condiciones temáticas. |
| Clip (Recortar) | st_intersection() |
Utiliza un polígono frontera como molde para recortar cualquier capa de entrada. |
| Buffer (Áreas de influencia) | st_buffer() |
Genera polígonos de amortiguación radial en metros de radio alrededor de puntos, líneas o polígonos. |
| Dissolve (Disolver) | group_by() %>% summarize() |
Fusiona geometrías que comparten un mismo atributo temático, eliminando fronteras internas. |
| Merge / Append | bind_rows() o st_union() |
Agrupa diferentes registros geográficos independientes en una única colección espacial. |
| Erase (Borrar) | st_difference() |
Resta o recorta una capa de la superficie de otra (operación de diferencia simétrica). |
| Intersect (Intersección) | st_intersection() |
Conserva únicamente las zonas espaciales compartidas donde dos o más capas se solapan. |
| Project / Define Projection | st_transform() /
st_set_crs() |
Gestiona y reescribe los sistemas de coordenadas y proyecciones geográficas. |
🎓 Conclusiones y Guía de Entregables del Proyecto
🧠 ¿Qué hemos aprendido a lo largo de esta práctica?
- Diagnóstico visual continuo: Ahora comprendes que en el trabajo GIS con código es imprescindible representar de manera continua los datos espaciales antes de grabarlos en disco. Esto nos permite comprobar de inmediato la corrección de los sistemas de coordenadas, de los límites territoriales y de la lógica de los algoritmos.
- Metodología y Lógica Reproducible: El gran valor
del análisis espacial mediante programación en R frente al clásico
entorno de ventanas es que este script es una infraestructura viva. Si
el ayuntamiento cambia los límites de protección ambiental o se
actualizan las coordenadas de los yacimientos, basta con volver a
ejecutar este documento (
Knit) y toda la cartografía intermedia e interactiva se actualizará automáticamente en pocos segundos, sin necesidad de rehacer decenas de pasos manuales. - Robustez y Saneamiento Espacial: Hemos aprendido a
programar funciones que limpian las imperfecciones de los datos
vectoriales brutos (topología rota, polígonos inconsistentes, columnas
de tipo lista complejas de los archivos históricos
.e00), preparándolos de forma segura para formatos modernos estandarizados por la OGC como el GeoPackage. - Visualización Multicanal: El uso de código te
capacita para exportar tanto cartografía técnica estática de alta
calidad sobre imágenes satelitales (
ggplot2) listas para incluir en informes en papel, como visores cartográficos web interactivos listos para producción (leaflet) y listos para compartir con terceros de manera dinámica.
📂 Entregables del Alumno
Para la evaluación final de esta práctica de localización óptima, el
alumno deberá entregar un único archivo comprimido .zip
titulado con su nombre y apellidos que contenga la siguiente estructura
de carpetas y archivos:
Curso_GIS_RSU.Rmd: Este archivo de código fuente de Markdown documentado con tus comentarios y apuntes de clase.Curso_GIS_RSU.html: El documento web interactivo generado tras pulsar el botón Knit en RStudio.vertedero/Analisis/Analisis.gpkg: El archivo de base de datos geográfica (GeoPackage) que contiene todas las capas vectoriales generadas a lo largo de las distintas fases (los buffers individuales, la geología corregida y la capa final de áreas óptimas).plots/Resultado_Final_Satelital.png: El mapa estático final sobre ortofoto de satélite exportado en alta resolución dentro de la carpetaplots/.
¡Enhorabuena! Has completado con éxito la migración de un flujo metodológico de SIG de escritorio a la potencia analítica de la programación de R moderno.