3 Transformaciones básicas con R

A continuación se presentan las operaciones de transformación básicas que podemos realizar sobre las capas con R. Para trabajar con las capas vectoriales utilizaremos el paquete sf, para las capas ráster usaremos los paquetes terra y satres. Este último paquete se ha desarrollado a raíz de esta asignatura.

Adicionalmente, solo en los casos en los que se indique, usaremos las funciones del paquete de la asignatura sigugr².

3.1 Leer una capa

3.1.1 Leer una capa vectorial

Cuando leamos una capa vectorial, para comprobar que la codificación de los caracteres es la adecuada, tenemos que comprobar que las palabras con tilde se vean correctamente. Por defecto, la función st_read() lee los datos en UTF-8. Si no se leen correctamente deberemos cambiar la codificación, generalmente a WINDOWS-1252, como se muestra en el ejemplo siguiente.

file <- "cnig/cuadriculas/MTN25_ETRS89_Peninsula_Baleares_Canarias.shp"
cuadricula25 <- sf::st_read(file, options = "ENCODING=WINDOWS-1252")

Al leer un formato de archivo que pueda contener varias capas, podemos consultar las capas disponibles mediante la función siguiente:

file <- "lanjaron-jsamos.gpkg"
sf::st_layers(file)

Adicionalmente, al leer un formato de archivo que pueda contener varias capas, podemos indicar la capa que queremos obtener.

file <- "lanjaron-jsamos.gpkg"
layer <- "lanjaron-jsamos"
lanjaron_jsamos <-
  sf::st_read(file, 
              layer = layer)

Podemos representar la geometría de la capa leída mediante las funciones siguientes:

terra::plot(sf::st_geometry(lanjaron_jsamos))

3.1.2 Leer una capa ráster

Se pueden leer todas las bandas de un ráster conjuntamente mediante la función rast() del paquete terra. Las presentamos combinándolas con la función plotRGB().

file <- "linea/ortofoto/102643_ortofoto2016_50cm.jp2"
ortofoto <- terra::rast(file)

terra::plotRGB(ortofoto)

3.1.3 Componer un ráster virtual a partir de varios ráster

Si hemos obtenido varias capas ráster que, conjuntamente, cubren nuestra zona de trabajo, podemos definir un ráster virtual mediante la función vrt(), a partir de un vector de archivos que podemos obtener directamente leyendo el contenido de una carpeta, como se muestra a continuación.

folder <- "cnig/mdt05/"
files <- list.files(path = folder, pattern = "*.tif", full.names = TRUE)
mdt05 <- terra::vrt(files, filename= "tmp/mdt05.vrt", overwrite = TRUE)

terra::plot(mdt05)

Como resultado, crea el archivo que se indica en el parámetro filename, correspondiente al ráster virtual, que podemos tratar como cualquier otro ráster. Es un archivo de pequeño tamaño porque se limita a referenciar al resto de capas que lo componen.

Esta función es adecuada tanto para capas de una sola banda como para capas multibanda.

3.1.4 Obtener un ráster multibanda a partir de varias bandas

Si tenemos las bandas almacenadas en archivos distintos (como ocurre con el caso de las bandas obtenidas al descargar los datos de los satélites), para obtener un ráster multibanda, debemos leerlas independientemente y generar a partir de ellas la nueva estructura.

folder <- "esa/IMG_DATA_01/R10m/"
file1 <- "T30SVF_20230905T105621_B02_10m.jp2"
b2 <- terra::rast(paste0(folder, file1))
file2 <- "T30SVF_20230905T105621_B03_10m.jp2"
b3 <- terra::rast(paste0(folder, file2))
file3 <- "T30SVF_20230905T105621_B04_10m.jp2"
b4 <- terra::rast(paste0(folder, file3))

bands <- c(b4, b3, b2)
terra::plotRGB(bands, stretch = 'lin')

El ráster multibanda se genera al definir un vector de variables ráster. En este caso, se han definido en el orden adecuado para que la función plotRGB() presente una fotografía en color.

3.1.5 Obtener un ráster con las bandas de satélite

Para obtener un ráster a partir de las bandas de satélite, podemos usar el paquete satres que, utilizando las funciones descritas en los apartados anteriores, a partir de la carpeta que almacena los archivos ráster (independientemente de su anidamiento), genera una estructura con todas las bandas obtenidas de satélite, clasificadas por su resolución espacial, generando estructuras de ráster virtual si es necesario.

esa <- satres::satres("esa")

A partir de la estructura generada, seleccionando la resolución espacial, podemos obtener objetos de la clase SpatRaster (definida en el paquete terra), y trabajar normalmente con ellos.

esa |>
  satres::get_spatial_resolution()

r10m <- esa |>
  satres::as_SpatRaster("r10m")

bands <- c(r10m[['B04']], r10m[['B03']], r10m[['B02']])
terra::plotRGB(bands, stretch = 'lin')

3.2 Reproyectar una capa

3.2.1 Reproyectar una capa vectorial

Proyectaremos una capa vectorial mediante la función st_transform(). En el parámetro crs se puede indicar directamente el código EPSG de la capa o bien el CRS obtenido a partir de otra capa mediante las funciones que lo devuelven.

lanjaron_jsamos_c25 <-
  sf::st_transform(lanjaron_jsamos, crs = sf::st_crs(cuadricula25))

La función utilizada para obtener el CRS es adecuada para una capa vectorial y también para una ráster. En ambos casos también podemos usar la función terra::crs().

3.2.2 Reproyectar una capa ráster

Un ráster se puede reproyectar mediante la función project() del paquete terra. Hay que tener presente que la operación suele ser bastante costosa en tiempo de ejecución y es posible que se pierda información ya que la correspondencia entre las celdas del ráster inicial y del resultado puede no ser directa.

ortofoto_lan_sat <-
  terra::project(ortofoto, terra::crs(lanjaron_jsamos_sat))

Por eso, según la operación que vayamos a hacer, debemos valorar podemos trabajar con el CRS del ráster y reproyectar las capas vectoriales.

3.3 Guardar una capa

3.3.1 Guardar una capa vectorial

Para guardar una capa vectorial, indicamos la variable con la capa, el nombre del archivo (el formato se toma de este nombre) y el nombre que le queremos dar a la capa. También podemos indicar si queremos que se borre el archivo si ya existe, mediante el parámetro delete_dsn.

sf::st_write(obj = lanjaron_jsamos,
             dsn = "tmp/lanjaron-jsamos.gpkg",
             layer = "lanjaron-jsamos",
             delete_dsn = TRUE)
sf::st_write(obj = lanjaron_jsamos_sat,
             dsn = "tmp/lanjaron-jsamos.gpkg",
             layer = "lanjaron_jsamos_sat")

3.3.2 Guardar una capa ráster

Independientemente del número de bandas de ráster y de la estructura de datos de este, lo podemos guardar en un archivo mediante la función writeRaster(). Mediante el parámetro filetype se indica el tipo de los datos de las bandas del ráster, aunque para el formato tif la opción por defecto de este parámetro es adecuada.

terra::writeRaster(mdt05,
                   "tmp/mdt05.tif",
                   filetype = "GTiff",
                   overwrite = TRUE)

terra::writeRaster(mdt05,
                   "tmp/mdt05.tif",
                   overwrite = TRUE)

En particular, por su tamaño, para el caso de la ortofotografía, se puede reducir el espacio que ocupa mediante los parámetros siguientes:

terra::writeRaster(ortofoto,
                   "tmp/ortofoto.tif",
                   filetype = "GTiff",
                   overwrite = TRUE,
                   gdal = c("COMPRESS=ZSTD"),
                   datatype = 'INT1U')

En lo que se refiere al tipo de datos, al recortarla usando la función terra::mask, lo cambia a FLT4S (float con signo y 4 bytes), cuando lo que necesita es INT1U (int sin signo y 1 byte).

3.3.3 Guardar las bandas de satélite

Si hemos creado un objeto satres a partir de las bandas de satélite, podemos guardarlas en archivos según su resolución espacial.

esa |>
  satres::save_by_resolution(out_dir = "tmp/esa")

3.4 Recortar una capa

3.4.1 Recortar una capa vectorial

Las capas vectoriales contenidas en un GeoPackage pueden tener asociado un estilo dentro del propio archivo. Por ejemplo, este es el caso de CORINE Land Cover: en los estilos se define la leyenda y los colores asociados a cada código.

Aunque recortemos y almacenemos las capas de interés en un nuevo GeoPackage, los estilos no se verán afectados ni tampoco se incluirán automáticamente en el nuevo archivo.

Por este motivo, hay que copiarlos de un GeoPackage a otro, como se indica en el apartado 3.4.1.4.

3.4.1.1 Mediante un rectángulo

La función st_crop() recorta una capa vectorial mediante el mínimo rectángulo envolvente orientado de Norte a Sur de la capa vectorial que se indica.

layer <- sf::st_crop(layer, lanjaron_jsamos)

3.4.1.2 Mediante cualquier polígono

La función st_intersection() obtiene la intersección de las capas vectoriales que se indican.

layer <- sf::st_intersection(layer, lanjaron_jsamos)

3.4.1.3 Caso especial de CORINE Land Cover

He podido recortar todas las capas vectoriales sin ningún problema mediante las funciones anteriores. En el caso de CORINE Land Cover del año 2018, al recortar de cualquier forma, producía un error relacionado con el tipo de datos (ParseException: Unknown WKB type 12). Está solucionado mediante la función sigugr::clip_multipoligon().

clc <- sigugr::clip_multipoligon(clc, lanjaron_jsamos)

3.4.1.4 Copiar estilos entre archivos GeoPackage

Una vez tengamos los dos archivos GeoPackage disponibles, podemos copiar los estilos de uno a otro mediante la función siguiente:

from <- "cnig/clc/CLC2018_ES.gpkg"
to  <- "tmp/clc_lanjaron_jsamos.gpkg"
sigugr::copy_styles(from = from, to = to)

3.4.2 Recortar una capa ráster

3.4.2.1 Mediante un rectángulo

La función crop() recorta una capa ráster mediante el mínimo rectángulo envolvente orientado de Norte a Sur de la capa vectorial que se indica.

rb1 <- terra::crop(rb, lanjaron_jsamos)

3.4.2.2 Mediante cualquier polígono

La función mask() recorta una capa ráster mediante el polígono de la capa vectorial que se indica.

rb2 <- terra::mask(rb1, lanjaron_jsamos)

En el ejemplo se usa el resultado de recortar una capa mediante la función crop() ya que contiene al polígono que usamos en la función mask().

3.4.3 Recortar las bandas de satélite

Si hemos creado un objeto satres a partir de las bandas de satélite, podemos recortar todas las bandas conjuntamente mediante la función clip_bands().

esa_lan <- esa |>
  satres::clip_bands(lanjaron_jsamos)

3.5 Tratar los archivos de una o varias carpetas

Para realizar la transformación en R podemos desarrollar código que permita tratar todos los archivos de una o varias carpetas conjuntamente. A continuación se presentan algunas de ellas.

Adicionalmente, para simplificar el proceso de transformación, podemos definir una función con los pasos necesarios para transformar un tipo de archivo y llamarla repetidas veces pasándole los parámetros adecuados.

3.5.1 Obtener y tratar los archivos de una carpeta

Mediante la función list.files() obtenemos los archivos contenidos en la una carpeta, del tipo indicado mediante un patrón. A continuación, podemos iterar sobre ellos para tratarlos convenientemente.

path <- "cnig/mdt05/"
lf <- list.files(path = path, pattern = "*.tif")
for (f in lf) {
  paste0(path, f)
  # ...
}

3.5.2 Obtener y tratar los archivos de varias carpetas

Si necesitamos tratar los archivos de varias carpetas, podemos definir una lista de las carpetas (o de la forma cómo obtenerlas) e ir recorriéndola para tratarlas.

path <- "usgs/"
lf <-
  list.files(
    path = path,
    pattern = "*.tif",
    recursive = TRUE,
    full.names = TRUE
  )
for (f in lf) {
  f
  # ...
}

3.5.3 Descomprimir todos los archivos de una carpeta

A veces tenemos que descargar varios archivos en formato zip o tar sobre un mismo tema para tratar su contenido conjuntamente. Podemos descomprimirlos indicando la carpeta en la que se encuentra mediante la función sat_untarzip().

path <- "originales/cnig/mdt"
r <- satres::sat_untarzip(path, out_dir = "tmp/mdt")

Se puede usar en primer lugar con el parámetro only_show_files = TRUE para mostrar los archivos a descomprimir y las carpetas de salida.

---

2. Disponible solo en GitHub. Se instala mediante devtools::install_github("josesamos/sigugr")↩︎