Capítulo 4 Transformar e Integrar Información Geográfica con QGIS y R

En este capítulo vamos a trabajar con datos obtenidos en el capítulo 3. Por tanto, es imprescindible tener disponibles los datos obtenidos en los ejercicios 3.2 a 3.5 para poder continuar.

Objetivos del capítulo

  • Repasar los conceptos básicos sobre transformación e integración de información geográfica explicados en clase de teoría materializados en herramientas específicas, de uso generalizado.
  • Usar GDAL/OGR desde QGIS y desde R mediante ejemplos.
  • Transformar nuestro conjunto de datos:
    • tanto para datos ráster como para datos vectoriales, nos focalizaremos en la zona de trabajo;
    • para datos vectoriales, los transformaremos para que todos tengan el mismo SRC.

4.1 Datos de partida y transformaciones a realizar

En este apartado se listan los datos de partida y se detalla (en forma de ejercicio) las transformaciones a realizar sobre ellos.

En el resto de este capítulo se muestra cómo realizar las operaciones necesarias para llevar a cabo las transformaciones indicadas en este único ejercicio del capítulo. Se presentan las operaciones en QGIS y R, aunque solo hay que realizarlas una vez: con QGIS o R, o combinándolos. En operaciones donde el volumen de datos es grande, tanto para ráster como para vectorial, funciona más rápido QGIS (por ejemplo, para la ortofoto o para las capas de CORINE Land Cover).

4.1.1 Datos de partida

Para el municipio seleccionado, los datos de partida que necesitamos, obtenidos en el capítulo 3, son los siguientes:

  • Capas vectoriales del término municipal (en un GeoPackage):
    • término municipal,
    • mínimo rectángulo que contiene al término municipal,
  1. Capas vectoriales temáticas:
    • capa de vías de comunicación (líneas),
    • capa de red hidrográfica (líneas),
    • una capa de puntos,
    • una capa de polígonos.
  2. Ortofotografía reciente (opcional42).
  3. MDT.
  4. Capas de CartoCiudad que incluyan al municipio (portal y manzana).
  5. Dos mapas de uso/ocupación del suelo de CORINE Land Cover de fechas distintas.
  6. Bandas de satélite para Sentinel-2 (B1 a B12)43 o Landsat 8 (B1 a B11), con uno de los dos es suficiente.

4.1.2 Transformaciones sobre los datos

Ejercicio 4.1 Las transformaciones que realizaremos son:

  • Recortar todas las capas mediante la capa del mínimo rectángulo que contiene al término municipal.
  • También recortar todas las capas vectoriales y la capa MDT mediante el polígono de la capa del término municipal.
  • Para las capas ráster, mantendremos el CRS que tenían cuando las hemos obtenido y, preferentemente, las almacenaremos en formato TIF.
  • Para las capas vectoriales:
    • utilizaremos para todas el mismo CRS (el CRS para la proyección UTM de la zona del municipio, para el Datum ETRS89), indica cuál es para tu municipio,
    • comprobaremos que la codificación de los campos texto sea correcta (buscar palabras con tilde) y, si es necesario, la cambiaremos adecuadamente,
    • almacenaremos las capas en formato Shapefile.
  • El nombre del archivo de cada capa resultado tendrá como prefijo el nombre del municipio y tu usuario de correo de la UGR (p.e., en mi caso lanjaron-jsamos-).

Para documentar la realización del ejercicio, en un documento, incluye un apartado para cada elemento numerado de la lista de datos de partida (sección 4.1.1):

  • el título del apartado será el nombre que aparece en la lista de la sección 4.1.1,
  • para cada archivo transformado, contendrá un subapartado con los datos siguientes:
    • una indicación sobre qué capa es (puede ser el título del subapartado),
    • el nombre del archivo de partida, y
    • una captura de pantalla para cada una de las capas obtenidas donde se muestre el contenido de la capa y el entorno donde se ha realizado la transformación (la pantalla completa de QGIS o RStudio). Para las bandas de satélite es suficiente la captura de pantalla solo para una de las bandas.

Los datos obtenidos en este capítulo deberán conservarse y estar disponibles para ser utilizados en las actividades que desarrollaremos en los capítulos siguientes.

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

4.2.1 Leer una capa

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

library(sf)

cuadricula25 <-
  st_read("MTN25_ED50_Peninsula_Baleares.shp",
          options = "ENCODING=WINDOWS-1252")

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

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

4.2.1.2 Leer una capa ráster

Se pueden cargar todas las bandas de un ráster conjuntamente mediante la función brick del paquete raster. Las presentamos combinándolas con la función plotRGB.

library(raster)

ortofoto <- brick("ortofotoLineA.tif")
plotRGB(ortofoto, r = 1, g = 2, b = 3)

También podemos leer las bandas una a una mediante la función raster y, con la función stack, construir un ráter multibanda partiendo de las bandas leídas. Aunque la estructura de datos stack es distinta a brick, el resultado obtenido en ambos casos es el mismo.

ortofotoB1 <- raster("ortofotoLineA.tif")
ortofotoB2 <- raster("ortofotoLineA.tif", band = 2)
ortofotoB3 <- raster("ortofotoLineA.tif", band = 3)

ortofotoStack <- stack(ortofotoB1, ortofotoB2, ortofotoB3)
plotRGB(ortofotoStack, r = 1, g = 2, b = 3)

4.2.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 gdalbuildvrt del paquete gdalUtils, a partir de los archivos de una carpeta, también permite indicarlos explícitamente mediante una lista de nombres de la forma c("nombre1","nombreN").

library(gdalUtils)

gdalbuildvrt(gdalfile = "LineA/*.jp2", output.vrt = "raster.vrt")

rb <- brick("raster.vrt")
plotRGB(rb, r = 1, g = 2, b = 3)

Como resultado, crea el archivo que se indica en el parámetro output.vrt, 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.

4.2.2 Reproyectar una capa

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

clc1990 <-
  st_read("CLC1990_ES.gpkg", layer="CLC90_ES")
e <- st_transform(lanjaron_jsamos, crs = st_crs(clc1990))

4.2.2.2 Reproyectar una capa ráster

Un ráster se puede reproyectar mediante la función projectRaster del paquete raster. Hay que tener presente que en la operación se pierde información ya que la correspondencia no suele ser directa entre las celdas del ráster inicial y del resultado.

4.2.3 Guardar una capa

4.2.3.1 Guardar una capa vectorial

Para guardar una capa vectorial, indicamos la capa y el nombre del archivo (el formato se toma de este nombre). También podemos indicar si queremos que se borre el archivo si ya existe, como se hace en el ejemplo siguiente.

st_write(obj = lanjaron_jsamos_clc1990,
         dsn = "lanjaron-jsamos-bbox-clc1990.shp",
         delete_dsn = TRUE)

4.2.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 datatype se indica el tipo de los datos de las bandas del ráster44, aunque la opción por defecto de este parámetro suele ser adecuada.

writeRaster(ortofotoStack, "lanjaron-jsamos-ortofoto.tif", overwrite = TRUE)

4.2.4 Recortar una capa

4.2.4.1 Recortar una capa vectorial

He podido recortar todas las capas vectoriales sin ningún problema mediante las funciones que se indican a continuación. 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. Lo he solucionado convirtiendo la geometría de la capa antes de recortarla, mediante la función siguiente.

clc2018 <- st_cast(clc2018, "MULTIPOLYGON")
4.2.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.

clc <- st_crop(clc1990, lanjaron_jsamos)
4.2.4.1.2 Mediante cualquier polígono

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

clc2 <- st_intersection(clc, lanjaron_jsamos)

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

4.2.4.2 Recortar una capa ráster

4.2.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 <- crop(rb, lanjaron_jsamos)
4.2.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 <- 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.

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

4.2.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 <- "satelite/LC08_L1TP_200034_20210716_20210721_01_T1/"
lf <- list.files(path=path, pattern="*.TIF")
for (f in lf) {
  r1 <- brick(paste(path, f, sep=""))
  # ...
}

4.2.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 <- "satelite/L2A_T30SVF_A032123_20210816T110306/IMG_DATA/R10m/"
res <- c("R10m", "R20m", "R60m")
for (r in res) {
  path_t <- stringr::str_replace(path, "R10m", r)
  lf <- list.files(path=path_t, pattern="*.jp2")
  for (f in lf) {
    r1 <- brick(paste(path_t, f, sep=""))
    # ...
    # para guardar en formato TIF
    f2 <- stringr::str_replace(f, "jp2", "tif")
    # ...
  }
}

4.3 Transformaciones básicas con QGIS

4.3.1 Codificación de los caracteres en una capa

Podemos comprobar si la codificación de los caracteres de una capa es correcta observando palabras que lleven tilde, por ejemplo abriendo la tabla de atributos de la capa y viendo las columnas con datos de nombres.

Hay muchas codificaciones posibles pero, generalmente, en caso de observar que no se muestran correctamente palabras con tilde, el problema se debe a que la capa tiene codificación UTF-8 y se ha considerado que tiene la del sistema (System) o viceversa. Si ocurre esto, podemos cambiar la codificación en el apartado Fuente de la ventana Layer Properties de la capa (se abre desde el menú contextual de la capa o bien pulsando sobre la opción Capa > Propiedades de la capa…).

4.3.2 Reproyectar una capa

Para cambiar el SRC de una capa (reproyectarla), lo que hacemos es guardar una nueva capa a partir de la original. Como vimos en el apartado 1.5.5, seleccionamos la capa y pulsamos sobre la opción Capa > Guardar como…, también se puede hacer la misma operación desde el menú contextual de la capa, pulsamos sobre Exportar > Guardar objetos como…

Seleccionar SRC al guardar una capa en un archivo.

Figura 4.1: Seleccionar SRC al guardar una capa en un archivo.

En la ventana que se abre (figura 4.1), definimos el formato de la nueva capa, su nombre y ubicación, y, en el campo SRC, pulsando sobre el icono Seleccionar SRC que hay a su derecha, seleccionamos el nuevo SRC45.

Selección del SRC.

Figura 4.2: Selección del SRC.

La selección del SRC se lleva a cabo en la ventana de la figura 4.2, si conocemos el código EPSG de la nueva capa, lo podemos introducir en el campo Filtrar y seleccionamos el SRC correspondiente.

Esta operación se puede aplicar tanto a capas vectoriales como a capas ráster.

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

Si trabajamos con una zona relativamente amplia, seguramente hayamos obtenido varias capas ráster que, conjuntamente, cubren la zona de trabajo, probablemente cubran incluso una zona más amplia que la que necesitamos. Hay varias posibilidades para operar conjuntamente con todas las capas, por ejemplo, podríamos fusionarlas en un solo ráster (se genera una nueva capa con todos los datos). Otra posibilidad más sencilla es generar un ráster virtual a partir de los componentes (ocupa poco espacio nuevo porque se basa en las otras capas). Esta última posibilidad es adecuada sobre todo si después vamos a recortar el ráster resultante, por tanto, es la que vamos a considerar.

Construir un ráster virtual.

Figura 4.3: Construir un ráster virtual.

Para definir un ráster virtual a partir de varias capas, cargamos todas las capas en el visor para poder apreciar el ráster que forman conjuntamente. A continuación, pulsamos sobre Ráster > Miscelaneous > Construir ráster virtual… (figura 4.3).

Opciones para construir un ráster virtual.

Figura 4.4: Opciones para construir un ráster virtual.

En la ventana que se abre (figura 4.4), en el campo Input layers, pulsando sobre el botón con puntos suspensivos situado a la derecha, seleccionamos las capas que compondrán el ráster: si solo hemos cargado las capas del ráster (recomendable), podemos seleccionar todas las capas pulsando sobre el botón Seleccionar todo. Dado que todas las capas tendrán la misma resolución, no afecta el criterio que elijamos en el campo Resolution (puede obtener la resolución media, la más fina o la más gruesa). Es muy importante tener deseleccionada la opción Place each input file into a separate band porque, en caso de estar seleccionada, el resultado sería un ráster multibanda y no es lo que necesitamos. En la parte inferior de la ventana, en Llamada a la consola de GDAL/OGR, se muestra el código que se genera en función de los parámetros que vamos definiendo. Para generar esta nueva capa, QGIS se basa en la funcionalidad que ofrecen las librerías GDAL/OGR. Si vamos a usar la capa virtual en otros proyectos, podemos salvarla indicando el nombre del archivo en el campo Virtual; si solo queremos operar con la nueva capa, podemos dejar que la guarde en un archivo temporal. Una vez definidos todos los parámetros, pulsamos sobre el botón Ejecutar para llevar a cabo el proceso de definición de la nueva capa. En la pestaña Registro de la ventana nos indica que lo va haciendo y también avisa cuando haya terminado.

Una vez generada la nueva capa (le asigna el nombre Virtual), para simplificar, podemos eliminar las capas componentes del visor de capas y trabajar solo con la nueva capa virtual (recomendable).

Si al realizar operaciones sobre el raster virtual, podemos probar a salvarlo como una capa, iniciar un nuevo proyecto con la nueva capa y realizar las operaciones sobre ella.

4.3.4 Recortar una capa ráster

Para recortar una capa ráster, vamos a contemplar dos posibilidades: recortarla mediante un rectángulo orientado de Norte a Sur o hacerlo mediante cualquier polígono. En ambos casos, vamos a usar una capa vectorial como máscara de la operación de recorte.

Cuando vayamos a aplicar operaciones en las que intervengan varias capas, es importante que todas las capas tengan el mismo SRC46. Por otro lado, evitaremos reproyectar capas ráster: si la capa vectorial no tiene el mismo SRC que la capa ráster, generaremos una nueva capa vectorial reproyectando la capa original.

4.3.4.1 Mediante un rectángulo

Visualizamos la capa ráster y la capa vectorial (rectángulo orientado de Norte a Sur) conjuntamente, y pulsamos sobre Ráster > Extraction > Cortar ráster por capa de máscara…

Se abre una ventana para definir los parámetros de la operación. Seleccionamos en el campo Capa de entrada la capa ráster y en el campo Capa de máscara la capa vectorial (aquí podemos comprobar que las capas tengan el mismo SRC). Comprobamos que la única opción seleccionada sea Match the extent of the clipped raster to the extent of the mask layer (Ajustar la extensión del ráster cortado a la extensión de la capa de máscara). En el campo Cortado (máscara), pulsando sobre el botón con puntos suspensivos situado a su derecha, elegimos la opción Guardar a archivo… e indicamos la ubicación, el nombre y el tipo: el tipo debe ser TIF files. También marcamos la opción Abrir el archivo de salida después de ejecutar el algoritmo. El comando que ha generado con la configuración que hemos indicado, se puede ver en el apartado Llamada a la consola GDAL/OGR. Para ejecutar la operación, pulsamos sobre el botón Ejecutar. El resultado nos lo muestra en la pestaña Registro de la misma ventana de definición de parámetros. Podemos cerrar la ventana y desactivar las otras capas para ver el resultado en la capa que llama con el nombre que hayamos indicado para el archivo de salida.

4.3.4.2 Mediante cualquier polígono

Para recortar un ráster por cualquier polígono, al igual que para un rectángulo, visualizamos la capa ráster y la capa vectorial conjuntamente, y pulsamos sobre Ráster > Extracción > Cortar ráster por capa de máscara…

Para mantener los bordes de la capa vectorial, en la definición de los parámetros de la operación de recorte, hemos de seleccionar la opción Crear una banda alfa de salida47. El resto de parámetros se definen igual que para el recorte mediante un rectángulo.

4.3.5 Recortar una capa vectorial

Para recortar una capa vectorial, vamos a presentar la operación básica de recorte y, en casos donde la capa original tenga un estilo predefinido, adicionalmente vamos a copiar el estilo de la capa.

Al igual que para las capas ráster, cuando vayamos a aplicar operaciones en las que intervengan varias capas, es importante que todas las capas tengan el mismo SRC. En el caso de trabajar con capas vectoriales, se podría reproyectar cualquiera de ellas. Si vamos a recortar una capa, generalmente será de mayor tamaño que la capa que vamos a recortar, entonces reproyectaremos la capa de máscara para la operación de recorte.

4.3.5.1 Recortar cualquier capa vectorial

Visualizamos las dos capas vectoriales conjuntamente, y pulsamos sobre Vector > Geoprocessing Tools > Cortar…

Se abre una ventana para definir los parámetros de la operación. Seleccionamos en el campo Capa de entrada la capa vectorial a recortar y en el campo Capa de superposición la capa vectorial a usar como plantilla (aquí podemos comprobar que las capas tengan el mismo SRC). En el campo Cortado, pulsando sobre el botón con puntos suspensivos situado a su derecha, elegimos la opción Guardar a archivo… e indicamos la ubicación, el nombre y el tipo. También marcamos la opción Abrir el archivo de salida después de ejecutar el algoritmo. El comando que ha generado con la configuración que hemos indicado, se puede ver en el apartado Llamada a la consola GDAL/OGR. Para ejecutar la operación, pulsamos sobre el botón Ejecutar. El resultado nos lo muestra en la pestaña Registro de la misma ventana de definición de parámetros. Podemos cerrar la ventana y desactivar las otras capas para ver el resultado.

4.3.5.2 Copiar y pegar el estilo de una capa vectorial

Si comparamos el estilo de la capa resultado con el estilo de la capa original, observamos que la capa resultado lo ha perdido. Se puede copiar el estilo de una capa vectorial y aplicárselo a otra, en este caso, a la capa resultado.

Para copiar el estilo de una capa, en el menú contextual de la capa, pulsamos sobre Estilos > Copiar estilo > Totas las categorías de estilo. De igual forma, en el menú contextual de la capa destino, pulsamos sobre Estilos > Pegar estilo > Totas las categorías de estilo.

Para guardar ese estilo asociado a la capa, lo hacemos desde el menú contextual de la capa, pulsando sobre la opción Exportar > Guardar como archivo de definición de capa… e indicamos la carpeta y el nombre del archivo.

Guardar el estilo de una capa.

Figura 4.5: Guardar el estilo de una capa.

Si usáramos un GeoPackage (o una BD PostGIS) en lugar de un Shapefile para guardar la capa, podríamos guardar el estilo dentro del propio Geopackage. Accedemos a la ventana de propiedades de la capa almacenada en un Geopackage y, en el apartado Simbología, pulsamos sobre la flecha del botón Estilo y seleccionamos la opción Guardar estilo… (figura 4.5). En la ventana que se abre, seleccionamos GeoPackage como destino del estilo, le damos un nombre al estilo e indicamos que se use como estilo predeterminado para la capa.


  1. Ocupa mucho espacio y tiempo de proceso. Se puede intentar, si da problemas se indica en el documento y no se tendrá en cuenta a la hora de corregir. Para prácticas posteriores no habrá problema si no se tiene, se sustituirá por otra capa.↩︎

  2. La banda B10 es posible que no esté disponible en la descarga, si es este el caso, lo indicaremos en el documento.↩︎

  3. Los posibles tipos de datos se pueden consultar en https://www.rdocumentation.org/packages/raster/versions/3.0-7/topics/dataType.↩︎

  4. Seleccionamos el SRC correspondiente a la zona de nuestro municipio (no el de la figura).↩︎

  5. En principio deberían funcionar igual con distinto SRC, la experiencia es que se evitan errores y las operaciones se ejecutan más rápido si en cada una solo se realiza una transformación.↩︎

  6. Si no seleccionamos esta opción, recorta el ráster por el polígono pero, adicionalmente, lo incluye en un rectángulo orientado de Norte a Sur y al resto de celdas de este le asigna el valor 0.↩︎