1 INTRODUCCIÓN

Es muy habitual que las escenas con imágenes de satélite sean mucho mayores que el área bajo estudio. Esto supone un gasto inútil de tiempo cuando se realizan diferentes procedimientos sobre esas imágenes originales. La solución es recortar las imágenes originales utilizando para ello un fichero shapefile.

Existen diferentes opciones para definir ese área de estudio.

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

El método más simple y directo es crear el área de estudio mediante una lista de coordenadas que serán convertidas en un shapefile mediante el paquete sf. El procedimiento es sencillo.

Se carga en memoria el paquete sf.

library(sf)

Creación de dos vectores con las coordenadas del polígono que corresponde al área de interés, uno con las longitudes y otro con las latitudes (siempre en este orden.

lon <- c(-3.950244286968632,-3.399265032492788, -3.399397758698454,-3.949447474688914, -3.950244286968632) 
lat <- c(40.25054369881899, 40.24978981099207, 40.64984182476227, 40.64910947255029, 40.25054369881899)

Estos dos vectores deben convertirse en una lista

lista <- list(cbind(lon, lat))

Finalmente, este objeto lista se convierte en un objeto sf.

pol <- st_polygon(lista)

El último paso es añadir la proyección al objeto.

poligono <- st_sf(st_sfc(pol), crs= 4326)
plot(st_geometry(poligono), axes = TRUE)

Si no se desean mantener estos objetos, deben eliminarse del Global Environment

rm(list=ls())

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

Cuando se desea obtener un fichero shapefile sin límites administrativos definidos (un simple cuadrado, un rectágulo) una posible opción es delimitar ese polígono con Google Earth y grabar el resultado como fichero *.kml. Algunos sitios web gratuitos pueden convertir ese archivo kml en un archivo shape, pero también es posible replicar el mismo procedimiento en R.

A continuación debe importarse el fichero *.kml con la función sf::st_read.

poligono <- st_read("D:/G174_OCW/Lab_3_preprocesamiento_imagenes_satelite/Crear_shapefile_limites_ROI/Poligono_Google_Earth.kml")
poligono

Esta opción importa también la dimensión z (aunque son su valor es 0), por lo que se elimina para mantener las dimensiones x e y. La función sf::st_zm elimina una o varias dimensiones de un objeto sf.

poligono <- st_zm(poligono[1], 
                  drop=T, 
                  what='ZM') 
plot(poligono)

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(poligono, "D:/G174_OCW/Trabajo_personal/Crear_shapefile_limites_ROI/Poligono_Google_Earth.shp", append = TRUE)

4 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. Una de ellas es la función raster::drawExtent, que permite el dibujo de una zona con un diseño rectangular de manera interactiva. Este procedimiento requiere varios pasos.

Primero se debe activar el packete raster.

library(raster)

Como en capítulos precedentes, se creará un objeto ráster, ahora formado por 625 píxeles, a las que se incorporarán las coordenadas geográficas.

raster <- raster(ncol= 25,                 # Número de columnas del ráster
               nrow= 25,                   # Número de filas
               xmx= -2,xmn= -5,            # Coordenadas del eje X (considerémoslas la longitud)
               ymx= 42,ymn= 39)            # Coordenadas del eje y (idem latitud)
values(raster) <- as.integer(runif(ncell(raster),
                                 min = 1,
                                 max= 10))
raster

A continuación se representa gráficamente el nuevo objeto ráster.

plot(raster, axis = TRUE)

El siguiente paso consiste en el dibujo del área de trabajo sobre el objeto raster.

roi_bis <- drawExtent(show = TRUE,                          # Nos muestra la región de interés
                      col = "yellow")                       # Color amarillo del borde del aárea de interés
roi_bis

Por último, se recortará la imagen original con las dimensiones de la nueva region de interés obtenida de la función anterior.

vista3 <- crop(raster,                                            # Objeto raster inicial.
               roi_bis)                                           # Objeto extent para recorte.

De la misma manera, se puede construir un polígono irregular con la función raster::drawPoly del paquete raster y su posterior conversión a un objeto sfcon la función sf::st_as_sf y su grabación como fichero shapefile con sf_st_write.

Como en el caso anterior, la digitalización de estos bordes requiere la ejecución desde la consola de R

poligono <- drawPoly(sp = TRUE, 
                     col='yellow', 
                     lwd=2)
poligono <- st_as_sf(poligono)
st_write(poligono, "D:/G174_Teledeteccion_2022/OCW/Lab_3_preprocesamiento_imagenes_satelite/datos/recorte/poligono.shp")

Otra opción interactiva permite la elaboración de esa región de interés mediante el uso de Internet. Para ello se pueden utilizar los paquetes mapview y mapedit, que ofrecen un entorno gráfico para definir de manera interactiva el área de interés sobre un mapa en la web, añadiendo cierto número de herramientas de dibujo de gran utilidad.

Primero se deben instalar y cargar los citados paquetes:

# install.packages("mapview")
# install.packages("mapedit")
library(mapview)
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.

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

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

5.1 GADM

Esta base de datos (https://gadm.org/download_country.html) dispone de los límites administrativos de diferentes entidades del globo con una alta resolución espacial, pudiendo se acceder a ella mediante el paquete GADMTools. Este permite descargar directamente un objeto sf que puede ser utilizado para recortar las imágenes de satélite. A fecha de hoy ha pasado al archivo CRAN https://cran.r-project.org/web/packages/GADMTools/index.html

Como es habitual, se comienza instalando el paquete (si no se dispone de él) y activándolo.

#install.packages(GADMTools)
library(GADMTools)

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

Ejemplos

españa_nivel_1 <- gadm_sf_loadCountries("ESP", 
                                        level = 1, 
                                        basefile = "D:/Descargas/")

Una vez descargado, se puede representar gráficamente mediante la función

gadm_plot(españa_nivel_1)

Dada la complejidad de los niveles superiores, es conveniente revisar los nombres de las unidades políticas mediante el argumento listNames() acompañado del nivel correspondiente

listNames(españa_nivel_1, level = 1)

Por último, se puede filtrar y seleccionar la unidad que nos interesa mediante la función gadm_subset().

cantabria <- gadm_subset(españa_nivel_1,level = 1, regions = c("Cantabria","Asturias"))

Si seleccionamos uno de los slots, se puede crear un mapa que luego se exporta como *.shp

mapa.cantabria <- cantabria$sf[ , c(3, 6)]
names(mapa.cantabria) <- c("nombre", "geometry")

5.2 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 rgeoboundaries se 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")

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

library(rgeoboundaries)
españa_boundary <- gb_adm0("Spain")

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)
ggplot(data = españa_boundary) +
  geom_sf()

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

iberia_boundaries <- gb_adm0(c("Spain", "Portugal")) 
                             
ggplot(data = iberia_boundaries) +
  geom_sf()

5.2.3 Diferentes niveles de límites administrativos

También se pueden descargar los inferiores de los límites administrativos de los paises. 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"), type = "sscgs")

ggplot(data = iberia_nivel_ad1) +
  geom_sf()

Descargar límites administrativos de nivel 2

iberia_nivel_ad2 <- gb_adm2(c("spain", "portugal"), type = "sscgs")

ggplot(data = iberia_nivel_ad2) +
  geom_sf()

5.2.4 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:

  • shapeName.

  • shapeISO.

  • shapeID.

  • shapeGroup

  • shapeType.

españa_boundary

El objeto iberia_boundaries es una colección simple features con con 2 features (los dos países, España y Portugal) ), y 5 campos.

iberia_boundaries

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_boundaries) +
  geom_sf() +
  geom_sf_label(aes(label = shapeName))

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_boundaries) +
  geom_sf() +
  geom_sf_label(aes(label = shapeName)) +
  xlab("Longitude") +
  ylab("Latitude") +
  ggtitle("Límites de España y Portugal")