💡 OBJETIVOS DE LA ACTIVIDAD:
► Conocer diferentes procedimientos para delimitar un área de trabajo (Region of interest, ROI).
► Familiarizarse con el tratamiento de objetos espaciales con formato vectorial bajo R.
💡 MATERIALES PARA LA ACTIVIDAD:
Los materiales para el desarrollo de esta actividad son los mismos de la actividad pasada.
► Imagen Sentinel 2. Corrresponde a una imagen del satélite Sentinel 2 de la ciudad de Vancouver, en Canadá.
Las escenas disponibles de Landsat ocupan unos 170 km x 185 km, es decir, unoas 31450 \(km^2\)), mientras que los gránulos de Sentinel 2 abarcan alrededor de 100 km2. Trabajar con las imágenes originales de cualquiera de estas plataformas supone manejar ficheros que pueden llegar a superar 1GB. Este gran tamaño ralentiza los análisis y dificulta el manejo de la información. Por esta razón, los productos originales suelen modificarse, bien recortando las imágenes originales, bien agregando varias en una única.
En Teledetección se utiliza el concepto de “área/región de interés a una porción de una imagen o escena que se quiere filtrar para operar con ella posteriormente. Normalmente responde al acrónimo ROI, que será utilizado de ahora en adelante.
En el siguiente apartado se revisarán algunos procedimientos para la identificación y extracción de una ROI. Esta ROI se puede usar para recortar, reproyectar o ajustar objetos espaciales posteriormente. Estos procedimientos implican el uso de una imagen ráster sobre la cual se desea identificar esa región de interés, y un objeto vectorial que servirá para delimitar esa región de interés.
► Un objeto SpatRaster
obtenido a partir de una imagen
de satélite.
library(terra)
## Warning: package 'terra' was built under R version 4.3.3
## terra 1.8.21
imagen <- rast("D:/G174_2025/LABORATORIO_3_Preprocesamiento_imagenes_satelite/datos/delimitacion_roi/sentinel.tif")
plotRGB(imagen, r=5, g=4, b=2, stretch="lin")
terra::ext()
.Tanto los objetos Raster
como los
vectoriales
contienen información sobre su extensión
espacial, que adquiere denominaciones diferentes
extensión (terra::ext
) en objetos ráster.
cuadro delimitador (sf::bbox
) en objetos
vectoriales.
El método más simple y directo para delimitar una posible ROI es
utilizar la función terra::ext()
[https://rspatial.github.io/terra/reference/ext.html].
Este procedimiento crea un objeto con el mismo Sistema de Coordenadas de
referencia, por lo que no es necesario comprobar su proyección. Con ella
se define un espacio que corresponde a un rectángulo definido por cuatro
valores, que corresponden a las cuatro esquinas del rectángulo de
cierre:
las coordenadas mínima y máxima de la dimensión x
(xmin
, xmax
)
las coordenadas minima y maxima de la dimensión y (ymin, ymax),
A continuación, se construirá el área de interés con las coordenadas (en este caso WGS 84 / UTM zone 10N -EPSG:32610-) que delimitan el área de trabajo, según el orden adecuado (xmin, xmax, ymin, ymax)
extension <- ext(490000, 500000, 5451992, 5456626)
plotRGB(imagen, r=5, g=4, b=2, stretch="lin")
plot(extension, add = TRUE)
Una alternativa al procedimiento anterior es crear una lista de
coordenadas que serán convertidas en un fichero vectorial. Este último
puede ser un objeto sf
o un objeto
SpatVector
.
Este procedimiento se inicia creando dos vectores con las coordenadas del polígono que constituye el área de interés, uno con las longitudes y otro con las latitudes (siempre en este orden). Por simplicidad, se diseña también un rectángulo, aunque se podría utilizar formas poligonales de mayor complejidad.
utm_x <- c(486077.676950998, 496242.577132486, 496242.577132486, 486077.676950998,486077.676950998)
utm_y <- c(5451992.17785844, 5451992.17785844, 5456626.95099819, 5456626.95099819,5451992.17785844)
Ambos dos vectores deben convertirse en una lista
lista <- list(cbind(utm_x, utm_y))
sf
.El procedimiento es sencillo y requiere cargar en memoria el paquete sf.
library(sf)
## Warning: package 'sf' was built under R version 4.3.3
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
El objeto lista
se convierte en un objeto
sf
.
pol_sf <- st_polygon(lista)
El último paso es añadir la proyección al objeto.
roi_sf <- st_sf(st_sfc(pol_sf),
crs= "EPSG:32610")
A continuación, se superpone la imagen con el nuevo polígono.
plotRGB(imagen, r=5, g=4, b=2, stretch="lin")
plot(st_geometry(roi_sf), axes = TRUE, add = TRUE)
Si no se desean mantener estos objetos, deben eliminarse del Global Environment
rm(list=ls())
SpatVector
.Dado que ya está creada la matriz con coordenadas y la lista, sólo
tenemos que aplicar la función terra::vect()
[https://rspatial.github.io/terra/reference/vect.html]
añadiendo como argumentos el tipo de objeto vectorial (`type =
“polygons”) y el sistema de coordenadas de referencia (en este caso,
EPSG:32610).
roi_spatvector <- vect(lista,
type="polygons",
crs="EPSG:32610")
plot(roi_spatvector)
A continuación, se superpone la imagen con el nuevo polígono.
plotRGB(imagen, r=5, g=4, b=2, stretch="lin")
plot(roi_spatvector, add = TRUE)
También es posible apoyarse en Google Earth para generar un objeto vectorial que sirve como ROI. Esto tiene la ventaja de diseñar el área de interés sobre una imagen real de dicha zona, ajustándolo a las características de la zona de interés. El procedimiento consiste en delimitar el área de trabajo con Google Earth, grabar el resultado como fichero .kml, y transformarlo en un objeto vectorial.
A continuación debe importarse el fichero *.kml con la función
sf::st_read
. Esto crea un objeto sf
.
poligono_google_earth <- st_read("D:/G174_2025/LABORATORIO_3_Preprocesamiento_imagenes_satelite/datos/delimitacion_roi/delimitacion_ROI_google_earth/Poligono_Google_Earth.kml")
## Reading layer `Poligono_Google_Earth.kml' from data source
## `D:\G174_2025\LABORATORIO_3_Preprocesamiento_imagenes_satelite\datos\delimitacion_roi\delimitacion_ROI_google_earth\Poligono_Google_Earth.kml'
## using driver `KML'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension: XYZ
## Bounding box: xmin: -123.1494 ymin: 49.24982 xmax: -123.0736 ymax: 49.27519
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84
poligono_google_earth
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension: XYZ
## Bounding box: xmin: -123.1494 ymin: 49.24982 xmax: -123.0736 ymax: 49.27519
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84
## Name Description geometry
## 1 vancouver POLYGON Z ((-123.1243 49.27...
Esta opción importa también una dimensión, denominada z
cuyo valor es 0. Por lo tanto, debe eliminarse para mantener las
dimensiones x e y. La función sf::st_zm
elimina una o
varias dimensiones de un objeto sf
.
roi_sf_google_earth <- st_zm(poligono_google_earth[1],
drop=T,
what='ZM')
plot(roi_sf_google_earth)
roi_sf_google_earth
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -123.1494 ymin: 49.24982 xmax: -123.0736 ymax: 49.27519
## Geodetic CRS: WGS 84
## Name geometry
## 1 vancouver POLYGON ((-123.1243 49.275,...
El último paso es añadir el Sistema de Cordenadas de Referencia
proyectadas al objeto sf
, ya que la proyección original es
geográfica.
roi_sf_google_earth <- st_transform(roi_sf_google_earth, "EPSG:32610")
A continuación, se superpone la imagen con el nuevo polígono.
plotRGB(imagen, r=5, g=4, b=2, stretch="lin")
plot(st_geometry(roi_sf_google_earth), axes = TRUE, add = TRUE)
A continuación se exporta como shapefile con la función
sf::st_write
. El argumento append
es útil para
sobreescribir nuevas versiones del mismo fichero.
st_write(roi_sf_google_earth, "D:/G174_2025/LABORATORIO_3_Preprocesamiento_imagenes_satelite/datos/delimitacion_roi/delimitacion_ROI_google_earth/Poligono_Google_Earth.shp", append = FALSE)
Para definir las dimensiones del área de interés existen otras alternativas, desde el propio entorno R.
terra::draw()
y
terra::sel()
.Una de ellas es la función terra::draw()
[https://rspatial.github.io/terra/reference/draw.html],
que crea un objeto SpatVector
o SpatExtent
(es
decir, rectangular) de manera interactiva. Este procedimiento requiere
varios pasos.
::: {.alert .alert-success}
💡 PROCEDIMIENTO INTERACTIVO:
► Representación gráfica de la imagen completa, usando las funciones
plot
(caso de objetos unibanda) o
raster::plotRGB
(caso de objetos multibanda).
plotRGB(imagen, r=5, g=4, b=2, stretch="lin")
► A continuación activa la función y se traza sobre la imagen
original la región de interés. Para ello, se ubica el cursor sobre la
imagen y con un click se identifica el vértice superior izquierdo de la
zona de interés. Posteriormente se define el límite completo de la zona
de interés seleccionando, con un click, el vértice inferior derecho
(xmax, ymin) en el caso de un polígono.Una vez dibujado el objeto, se
vuelve a R pulsando ESC
. También puede preestablecer el
número máximo de clicks con el argumento n = 1000
(número
máximo de clicks -no se aplica cuando x==“extent” -en este caso
n es 2) y xpd = TRUE
(si TRUE, se puede dibujar
fuera de la zona de dibujo.
borrame <- draw(imagen ="ext", # Tipo de objeto a dibujar: "extent", "polygon", "line", or "points"
col="red", # Color usado
lwd=2, # Anchura de la líneas
id= TRUE) # Si TRUE, un ID numérico se dibuja sobre el mapa
La función terra::sel()
puede crear tanto un polígono
como un rectángulo.
mapview
y mapedit
Otra opción interactiva permite la elaboración de esa región de
interés sin salir de R, mediante los paquetes mapview
y
mapedit
. Éstos ofrecen un entorno gráfico para definir de
manera interactiva el área de interés sobre un mapa, añadiendo cierto
número de herramientas de dibujo de gran utilidad. De nuevo, este
procedimiento requiere
Primero se deben instalar y cargar los citados paquetes:
# install.packages("mapview")
# install.packages("mapedit")
library(mapview)
## Warning: package 'mapview' was built under R version 4.3.2
library(mapedit)
Es muy importante definir qué capas cartográficas servirán de
referencia para el trazado de los límites administrativos con el
argumento basemaps
, desde un gran número de opciones
disponibles en sitio web leaflet providers preview
(https://leaflet-extras.github.io/leaflet-providers/preview/).
En este caso, se seleccionarán tres bases cartográficas diferentes, dos
de ESRI y una del proyecto OSM. Además se configura el color y
transparencia de los elementos trazados (argumento
vector.palette
) y se incluirá un argumento para ubicar en
la ventana el cuadro de herramientas con parámetro (argumento
layers.control.pos
).
mapviewOptions(basemaps = c('Esri.WorldImagery',"Esri.WorldShadedRelief","OpenStreetMap.DE"),
vector.palette = rgb(50,220,230,alpha = .90,maxColorValue = 255),
layers.control.pos = "topright")
Finalmente se define el nombre del objeto que recibirá el producto de la sesión interáctiva realizada en la pestaña Viewer de RStudio.
area_vista <- editMap()
Para terminar la sesión y recuparar los elementos trazados se debe
clicar en el botón Done
en la esquina inferior derecha.
Igual que en los casos anteriores, estos objetos vectoriales pueden grabarse en el disco duro.
Existen varias posibilidades para importar bases de datos cartográficas desde paquetes de R los límites administrativos oficiales de diferentes entidades.
Como es habitual, se comienza instalando el paquete (si no se dispone de él) y activándolo.
# install.packages("geodata")
# Alternativa
# install.packages("remotes")
# remotes::install_github("rspatial/geodata")
library(geodata)
## Warning: package 'geodata' was built under R version 4.3.3
A continuación, deben definirse los parámetros de configuración antes proceder a la descarga. Existen varias posibilidades:
► Identificar el área de interés mediante el Código ISO: la
Organización Internacional de Estandarización (ISO) representa países,
territorios dependientes y áreas de interés especial bajo el código
ISO_3166-1_alpha-3. Este se compone de tres caracteres alfabéticos
que representan a cada uno de los estados. En el caso de España, el
código ISO es ESP, mientras el de Nueva Zelanda es NZL. El argumento
fileNames
es el que incorpora este código.
► Seleccionar el nivel de detalle requerido. La información se
organiza por niveles o jerarquía política y de acuerdo con el país
reciben distintos nombres. Para ello se incluye el argumento
level
,
Nivel | España | Francia |
---|---|---|
0 | Pais | País |
1 | CA | Région |
2 | Provincia | Département |
3 | Comarca | Arrondissement |
4 | Municipio | Commune |
Ejemplo
españa_geodata <- gadm(country="ESP", # Código ISO de 3 letras o nombre completo del país. En caso de varias entidades simultáneas deben proporcionarse todos los nombres y convertirse en un vector con rbind.
level= 3, # Subdivisión administrativa según niveles.
version="latest", # La última versión, son posibles otras.
path=tempdir()) # Ruta dónde almacenar los datos descargados.
plot(españa_geodata)
class(españa_geodata)
## [1] "SpatVector"
## attr(,"package")
## [1] "terra"
El resultado es un objeto “SpatVector”.
rgeoboundaries
El paquete
rgeoboundaries es un cliente para la API de geoBoundaries, que proporciona los límites administrativos políticos de los países. A través de este tutorial se descargarán los límites administrativos de algunos paises mediante este paquete, que serán visualizados con el paquete
ggplot2. El paquete
rgeoboundariesse puede descargar desde GitHub usando el paquete
remotes`
que permite una fácil instalación de paquetes R desde repositorios
remotos como GitHub.
Primero se instalará y cargará dicho paqeute desde GitHub
# install.packages("remotes")
# library(remotes)
# remotes::install_github("wmgeolab/rgeoboundaries")
A continuación se activa el paquete.
library(rgeoboundaries)
## Warning: package 'rgeoboundaries' was built under R version 4.3.3
Para descargar los límites de los países usamos la función
geoboundaries::geoboundaries()
s.
En este caso descargaremos los límites de España.
españa_rgeo <- gb_adm0("Spain")
class(españa_rgeo)
## [1] "sf" "data.frame"
En este caso, el objeto vectorial es un objeto de la clase
sf
. Podría ser transformado en objeto vectorial de la clase
SpatVector
(lo que permitirá su integración con los objetos
SpatRáster) de esta manera:
españa_rgeo_spatvector <- vect(españa_rgeo)
El paquete ggplot2 se puede usar para trazar los
límites administrativos descargados usando la función
ggplot2::geom_sf()
# install.packages("ggplot2")
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
ggplot(data = españa_rgeo) +
geom_sf()
También podemos descargar los límites de varios países juntos al incluir los nombres de los países como un vector y dibujarlos a continuación.
iberia_rgeo <- gb_adm0(c("Spain", "Portugal"))
ggplot(data = iberia_rgeo) +
geom_sf()
También se pueden descargar los límites administrativos de los paises
según niveles. Para ello, se debe incluir el nivel administrativo como
argumento en la función geoboundaries()
. El nivel más bajo
es el 1 (“adm1”), mientras que el nivel administrativo 5 (“adm5”) es el
más bajo.
► Descargar límites administrativos de nivel 1
iberia_nivel_ad1 <- gb_adm1(c("spain", "portugal"))
ggplot(data = iberia_nivel_ad1) +
geom_sf()
► Descargar límites administrativos de nivel 2
iberia_nivel_ad2 <- gb_adm2(c("spain", "portugal"))
ggplot(data = iberia_nivel_ad2) +
geom_sf()
Si imprimimos el límite administrativo descargado, veremos que es una
colección de simple feature
con 1 característica (un único
objeto, la entidad geográfica) y 5 campos. Los campos son los atributos
geométricos incluidos en cada feature
:
shapeGroup
shapeType.
shapeName.
shapeISO.
shapeID.
canonical
geometry
Este es el caso:
españa_rgeo
## Simple feature collection with 1 feature and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -18.16131 ymin: 27.63784 xmax: 4.327785 ymax: 43.79238
## Geodetic CRS: WGS 84
## shapeGroup shapeType shapeName shapeISO shapeID
## 1 ESP ADM0 Spain ESP 43220588B53825046256317
## shapeCanonical geometry
## 1 Unknown MULTIPOLYGON (((-18.00495 2...
El objeto iberia_boundaries
es una colección
simple features
con 2 features
(los dos
países, España y Portugal) ), y 5 campos.
iberia_rgeo
## Simple feature collection with 2 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -31.26882 ymin: 27.63784 xmax: 4.327785 ymax: 43.79238
## Geodetic CRS: WGS 84
## shapeGroup shapeType shapeName shapeISO shapeID
## 1 ESP ADM0 Spain ESP 43220588B53825046256317
## 2 PRT ADM0 Portugal PRT 79933184B91516589295848
## shapeCanonical geometry
## 1 Unknown MULTIPOLYGON (((-18.00495 2...
## 2 Unknown MULTIPOLYGON (((-16.05158 3...
El campo shapeName
contiene los nombres de las
divisiones administrativas. Aprovechando ésto, se puede crear un mapa
con los nombres de las divisiones administrativas utilizando la función
ggplot2::geom_sf_label()
y configurando
label = shapeName
.
ggplot(data = iberia_rgeo) +
geom_sf() +
geom_sf_label(aes(label = shapeName))
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
El etiquetado de las divisiones administrativas también crea
automáticamente las etiquetas de los ejes como x e y. Para cambiar esas
etiquetas se utilizan las funciones ggplot2::xlab()
e
ggplot2::ylab()
, acompañada de la función
ggplot2::ggtitle()
para agregar un título a la trama.
ggplot(data = iberia_rgeo) +
geom_sf() +
geom_sf_label(aes(label = shapeName)) +
xlab("Longitude") +
ylab("Latitude") +
ggtitle("Límites de España y Portugal")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
Una de las mayores ventajas de trabajar con el paquete
sf
es que sus objetos son geometrías asociadas a tablas,
pudiendo usarse el mismo lenguaje de selección que se utiliza para el
análisis de data frames, por ejemplo con el paquete dplyr
.
Se mostrarán a continuación varias posibilidades.
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.2
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:terra':
##
## intersect, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
Comprobar la estructura del objeto sf
.
View(iberia_nivel_ad1)
► Para seleccionar del mapa de CCAA todas menos Canarias.
seleccion_1 <- iberia_nivel_ad1 %>% # Objeto `sf` inicial
select(shapeName)%>% # Selección de una única columna, en este caso shapeName
filter(shapeName != "Canarias") %>% # Selección de la fila que no contiene la etiqueta "Canarias"
plot()
► Para seleccionar del mapa de CCAA sólo Cantabria.
seleccion_1 <- iberia_nivel_ad1 %>% # Objeto `sf` inicial
select(shapeName)%>% # Selección de una única columna, en este caso shapeName
filter(shapeName == "Cantabria") %>% # Selección de la fila que contiene la etiqueta "Cantabria"
plot()
rnaturalearth
Una última opción es el paquete rnaturalearth. Este es un paquete que también proporciona datos cartográficos en formato vectorial que se puede visualizar utilizando otros paquetes R. Contiene los siguientes elementos:
► Un conjunto de mapas vectoriales (funciones):
Mapas de países (admin-0): ne_countries()
.
Mapas con los límites dentro de los países (admin-1) :
ne_states()
.
Mapas de líneas costeras: ne_coastline()
.
► Una función ne_download()
para la descarga de los
mapas vectoriales.
El paquete se activa
library(rnaturalearth)
## Warning: package 'rnaturalearth' was built under R version 4.3.2
Si se deseara descargar España como objeto sf
.
españa <- ne_countries(scale = "large", # Escala del mapa one of ‘110‘, ‘50‘, ‘10‘ or ‘small‘, ‘medium‘, ‘large‘.
type = "countries", # Alternativas ’map_units’, ’sovereignty’, ’tiny_countries’.
country = "spain", # Sustituir por continent para descargar todos los países de un continente.
returnclass = "sv") # Tipo de objeto espacial, "sf" (por defecto) o ‘SpatVector‘
class(españa)
## [1] "SpatVector"
## attr(,"package")
## [1] "terra"
Si se deseara descargar todo África
africa <- ne_countries(continent = "africa")
Una vez delimitada la región de interés, es necesario grabarla para
su uso en ulteriores trabajos. Tanto los objetos sf
como
los objetos SpatVector
pueden ser convertidos en ficheros
vectoriales. Por su uso frecuente, se recomienda el formato shapefile
(.shp) o el formato .geojson.
writeVector(roi_spatvector,
"D:/G174_2025/LABORATORIO_3_Preprocesamiento_imagenes_satelite/datos/delimitacion_roi/ventana_espacial.shp",
overwrite=TRUE)
Por ejemplo, el objeto roi_sf
con formato geojson,
st_write(roi_sf,
"D:/G174_2025/LABORATORIO_3_Preprocesamiento_imagenes_satelite/datos/delimitacion_roi/ventana_espacial.geojson",
overwrite=TRUE)