💡 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á.

script_laboratorio_3_delimitacion_ROI

INTRODUCCIÓN

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")

DELIMITACIÓN DE LA REGIÓN DE INTERÉS MEDIANTE LA FUNCIÓN terra::ext().

Tanto los objetos Raster como los vectoriales contienen información sobre su extensión espacial, que adquiere denominaciones diferentes

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:

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)

DELIMITACIÓN DE LA REGIÓN DE INTERÉS A PARTIR UN VECTOR DE COORDENADAS

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.

Creación del vector de coordenadas

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))

Como objeto 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())

Como objeto 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)

DELIMITACIÓN DE LA REGIÓN DE INTERÉS A PARTIR DE UN FICHERO *. kml

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.

Polígono sobre Google Earth
Polígono sobre Google Earth

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)

DELIMITACIÓN DE LA REGIÓN DE INTERÉS MEDIANTE UN PROCEDIMIENTO INTERACTIVO

Para definir las dimensiones del área de interés existen otras alternativas, desde el propio entorno R.

Uso de las funciones 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.

Uso de los paquetes 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.

Ventana map_edit
Ventana map_edit

Igual que en los casos anteriores, estos objetos vectoriales pueden grabarse en el disco duro.

DELIMITACIÓN DE LA REGIÓN DE INTERÉS MEDIANTE EL ACCESO A ENTIDADES ADMINISTRATIVAS DESDE R

Existen varias posibilidades para importar bases de datos cartográficas desde paquetes de R los límites administrativos oficiales de diferentes entidades.

Con el paquete GEODATA

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”.

Con el paquete 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 paqueteggplot2. El paquetergeoboundariesse puede descargar desde GitHub usando el paqueteremotes` 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

Límites administrativos de un solo país

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()

Límites administrativos de varios países

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()

Diferentes niveles de límites administrativos

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()

Información sobre los límites administrativos

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

Selección de una unidad espacial concreta

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()

Con el paquete 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")

GRABAR ROI COMO ARCHIVO VECTORIAL

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)