Capítulo 6 Analizar Información Geográfica desde QGIS y R

Disponemos de la información geográfica de nuestro municipio almacenada como archivos (obtenida en el capítulo 4) y en BD SpatiaLite, GeoPackage y PostGIS (almacenada en el capítulo 5). En este capítulo, vamos a aplicar varias técnicas de análisis a esta información desde QGIS o desde R.

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

Objetivos del capítulo

  • Presentar diversas técnicas para realizar análisis vectorial desde QGIS y R, en concreto:
    • Realizar consultas espaciales básicas mediante SQL sobre PostGIS desde QGIS.
    • Realizar consultas espaciales básicas mediante SQL sobre PostGIS desde R.
  • Presentar distintas técnicas de análisis de datos ráster.
  • Generar datos vectoriales a partir de datos ráster, en particular para el MDT.
  • Generar el MDT a partir de datos vectoriales.
  • Utilizar los datos ráster generados para presentar los datos.

A continuación, se irán presentando los pasos a realizar y cómo llevarlos a cabo con las distintas opciones consideradas. Los ejercicios propuestos están asociados a cada apartado.

6.1 Consultas espaciales sobre BD desde QGIS

QGIS ofrece varias posibilidades para analizar información geográfica. En este caso, el análisis se realiza en la BD (PostGIS), mediante la consulta que ejecuta: QGIS actúa como medio de acceso a la BD y de representación gráfica del resultado. Las operaciones a realizar sobre BD están basadas en Marquez (2015).

Para formular una consulta en SQL desde QGIS sobre una BD espacial, en primer lugar, debemos acceder al administrador de BD de QGIS (pulsando sobre Base de datos > Administrador de bases de datos…) y seleccionar una BD para conectarnos a ella. Previamente hemos de haber definido la conexión, como se ha hecho en el apartado 5.4.2.1.3.

Una vez con acceso a la BD, debemos abrir la pestaña para formular consultas: pulsando sobre la opción Base de datos > Ventana SQL. El resultado es que se abre una pestaña de formulación de consultas asociada a la BD seleccionada en el momento de abrir la ventana. Se pueden tener varias pestañas de consulta abiertas al mismo tiempo, una por cada BD distinta.

6.1.1 Consulta sobre una capa

La consulta que vamos a definir se muestra a continuación.

SELECT ST_Union(ST_Buffer(geometry, 100)) as zona_rio
FROM "lanjaron-jsamos-red-hidrografica";

Obtenemos un buffer (ST_Buffer) a partir de las líneas de una capa. Además, en la misma consulta, los componentes del resultado se fusionan (ST_Union). La consulta se ha definido en PostGIS para la capa correspondiente a la red hidrográfica del municipio.

Consulta y resultado sobre PostGIS desde QGIS.

Figura 6.1: Consulta y resultado sobre PostGIS desde QGIS.

En la figura 6.1, se puede ver la formulación de la consulta en la ventana de consultas de PostGIS. Al pulsar sobre el botón Ejecutar, PostGIS ejecuta la consulta y devuelve el resultado en forma de tabla. Si seleccionamos la opción Cargar como capa nueva, podemos indicar la columna que identifica los registros (si hay alguna) y la columna del resultado que contiene los datos geométricos (en la figura, la columna zona_rio); si pulsamos sobre el botón Cargar, carga los resultados como una capa más de QGIS, sobre la que podemos operar normalmente.

Resultado de la consulta de PostGIS como capa en QGIS.

Figura 6.2: Resultado de la consulta de PostGIS como capa en QGIS.

La capa resultado de la consulta cargada en QGIS se muestra en la figura 6.2. Partiendo de una capa vectorial compuesta por líneas, hemos obtenido una capa compuesta por polígonos mediante las operaciones de transformación que hemos utilizado.

6.1.2 Consulta sobre dos capas

Vamos a realizar un consulta sobre datos geográficos en la que intervienen dos capas. A continuación, se muestra la definición de la consulta considerada para PostGIS.

SELECT p.*
FROM "lanjaron-jsamos-portal" p
WHERE EXISTS (
   SELECT id
   FROM "lanjaron-jsamos-red-hidrografica"
   WHERE ST_DWithin(p.geometry, geometry, 100) 
);

Se trata de obtener los portales (de la capa de CartoCiudad) que se encuentren a menos de una distancia determinada de un elemento de la red hidrográfica. Definiremos la distancia de manera que se obtenga un resultado no vacío.

Consulta sobre dos capas y resultado obtenido.

Figura 6.3: Consulta sobre dos capas y resultado obtenido.

En la figura 6.3, se puede ver la formulación de la consulta. En este caso, podemos indicar la columna que identifica los registros y la columna del resultado que contiene los datos geométricos.

Resultados de las consultas como capas en QGIS.

Figura 6.4: Resultados de las consultas como capas en QGIS.

La capa resultado de la consulta (junto con la capa resultado de la consulta anterior) cargada en QGIS se muestra en la figura 6.4.

Ejercicio 6.1 Realiza dos consultas sobre una BD con datos de tu municipio desde QGIS (pueden ser equivalentes a las que se han mostrado en este apartado, la BD puede ser cualquiera):

  1. utilizando una sola capa,
  2. utilizando dos capas.
Para documentar la realización del ejercicio, para cada consulta, incluye una captura de pantalla de la ventana Administrador de BBDD donde se vea la BD y la consulta, y otra de QGIS donde se presente el resultado de la consulta como una capa del proyecto.

6.2 Consultas espaciales sobre BD desde R

En este apartado vamos a realizar desde R, sobre una BD de nuestro municipio, las mismas consultas desarrolladas en el apartado 6.1 desde QGIS.

6.2.1 Consulta sobre una capa

Para realizar una consulta sobre una BD podemos utilizar la función st_read del paquete sf. En el parámetro dsn indicamos la cadena de conexión, en este caso, la indicada para acceder a la BD PostGIS.

library(sf)

dsn = paste(
  "PG:dbname='lanjaron-jsamos'",
  "host='localhost'",
  "port='5432'",
  "user='postgres'",
  "password='postgres'",
  sep = " "
)

query = paste(
  'SELECT ST_Union(ST_Buffer(geometry, 100)) as zona_rio',
  'FROM "lanjaron-jsamos-red-hidrografica"',
  sep = " "
)

resultado = st_read(dsn, query = query)
plot(st_geometry(resultado))

La consulta también se define como una cadena que se pasa como parámentro a la función st_read. El resultado es una capa que podemos tratar normalmente

Consulta sobre una capa y resultado obtenido en RStudio.

Figura 6.5: Consulta sobre una capa y resultado obtenido en RStudio.

En la figura 6.5, se muestra la consulta y el resultado obtenido en RStudio.

Podemos obtener el mismo resultado utilizando la BD exclusivamente para obtener la capa original y transformándola adecuadamente mediante funciones del paquete sf, como se muestra a contiuación.

redHidro <- st_read(dsn, layer="lanjaron-jsamos-red-hidrografica")

resultado <- st_union(st_buffer(redHidro, dist = 100))
plot(st_geometry(resultado))

Con QGIS podríamos haber realizado una transformación similar a partir de la capa original, mediante las operaciones disponibles en los menús para datos vectoriales.

6.2.2 Consulta sobre dos capas

query = paste(
  'SELECT p.*',
  '  FROM "lanjaron-jsamos-portal" p',
  'WHERE EXISTS (',
  '  SELECT id',
  '  FROM "lanjaron-jsamos-red-hidrografica"',
  '  WHERE ST_DWithin(p.geometry, geometry, 100)',
  ')',
  sep = " "
)

resultado2 = st_read(dsn, query = query)
plot(st_geometry(resultado2), col = 'red',
     add = TRUE)

En el fragmento de código anterior, se muestran las operaciones para realizar la consulta en la que intervienen dos capas y la presentación del resultado junto a al resultado obtenido en el último apartado. Se usa la variable con la cadena de conexión definida anteriormente.

Consulta sobre dos capas y resultado obtenido en RStudio.

Figura 6.6: Consulta sobre dos capas y resultado obtenido en RStudio.

En la figura 6.6, se muestra la consulta sobre dos capas y el resultado obtenido en RStudio.

Ejercicio 6.2 Realiza desde R, mediante SQL, las mismas dos consultas desarrolladas en el ejercicio 6.1 sobre una BD con datos de tu municipio (la BD puede ser cualquiera).

Para documentar la realización del ejercicio, para cada consulta, incluye una captura de pantalla de RStudio donde se vea el código R y el resultado de la consulta.

6.3 Operaciones sobre bandas de satélite

Partiendo de las bandas de satélite para nuestro municipio, podemos realizar diversas composiciones de color. En este caso, vamos a obtener el llamado color natural. También vamos a obtener a partir de ellas un índice de vegetación de la zona.

6.3.1 Color natural

Para obtener la composición color natural necesitamos las bandas 2, 3 y 4 de Landsat-8 o bien de Sentinel-2 (en ambos casos, esas bandas representan la misma información).

6.3.1.1 Obtención en QGIS

En primer lugar, tenemos que crear un ráster virtual con esas bandas (Ráster > Miscelánea > Construir ráster virtual…). En la ventana de selección múltiple de las bandas de entrada, hemos de tener presente el orden de las bandas que seleccionamos ya que, posteriormente, se nombrarán como Banda 1, Banda 2 y Banda 3, en lugar de sus nombres originales. Tenemos que saber qué banda es cada una.

Definir un ráster virtual multibanda.

Figura 6.7: Definir un ráster virtual multibanda.

En la ventana de definición del ráster virtual (figura 6.7), es muy importante dejar seleccionada la opción Place each input file into a separate band. Esta opción es la diferencia entre un ráster virtual multibanda (como es el caso) y otro compuesto por otros ráster que se solapan.

Definir un ráster virtual multibanda.

Figura 6.8: Definir un ráster virtual multibanda.

Una vez creada la nueva capa, tenemos que indicar qué bandas se corresponden con las bandas roja, verde y azul. Esto se define en las propiedades de la nueva capa, en el apartado Simbología (figura 6.8). Partimos de las bandas 2, 3 y 4, al incorporarlas en el ráster virtual son, respectivamente, las bandas Banda 1, Banda 2 y Banda 3. Si miramos en cualquiera de las tablas de descripción de bandas de Landsat-8 o de Sentinel-2, la banda 4 (Banda 3 en el ráster virtual) es la banda roja, la banda 3 (Banda 2 en el ráster virtual) es la verde y la banda 2 (Banda 1 en el ráster virtual) es la azul. Esto es lo que se define en la figura 6.8.

De igual forma, creando ráster virtuales con las bandas adecuadas, podemos generar las combinaciones que deseemos.

6.3.1.2 Cálculo en R

library(raster)

b2 <- raster("lanjaron-jsamos-T30SVF_20210816T105621_B02_10m.tif")
b2 <- stretch(x = b2, minv = 0, maxv = 255)

b3 <- raster("lanjaron-jsamos-T30SVF_20210816T105621_B03_10m.tif")
b3 <- stretch(x = b3, minv = 0, maxv = 255)

b4 <- raster("lanjaron-jsamos-T30SVF_20210816T105621_B04_10m.tif")
b4 <- stretch(x = b4, minv = 0, maxv = 255)

trueColor <- stack(b4, b3, b2)

Para obtener el color verdadero en R, leemos las bandas que necesitamos mediante la función raster del paquete raster. En cada celda de las bandas tenemos valores de radiaciones en determinadas longitudes de onda. Para obtener los colores los valores han de estar en el intervalo [0, 255]. Escalamos los valores mediante la función stretch del mismo paquete.

Por último, obtenemos un nuevo ráster multibanda a partir de las bandas en el orden adecuado: si miramos en cualquiera de las tablas de descripción de bandas de Landsat-8 o de Sentinel-2, la banda 4 es la banda roja, la banda 3 es la verde y la banda 2 es la azul.

6.3.2 Índice de vegetación

El Índice de Vegetación de Diferencia Normalizada (NDVI) se utiliza para estimar la cantidad y desarrollo de la vegetación de una zona. Se considera la medición de la intensidad de la radiación de algunas bandas del espectro electromagnético que la vegetación emite o refleja53.

Se define a partir de la banda roja (\(Red\)), que es la banda 4 de Landsat-8 y de Sentinel-2, y la banda correspondiente al infrarrojo cercano (\(NIR\)), banda 5 de Landsat-8 y banda 8 de Sentinel-2, mediante la siguiente expresión:

\[\begin{equation} NDVI = \frac{NIR - Red}{NIR + Red} \end{equation}\]

6.3.2.1 Obtención en QGIS

Para calcular el índice de vegetación, podemos usar los datos de las bandas ráster conjuntamente en la calculadora ráster, accesible desde Ráster > Calculadora ráster…

Cálculo del índice de vegetación en la calculadora ráster.

Figura 6.9: Cálculo del índice de vegetación en la calculadora ráster.

En la figura 6.9, se muestran las bandas de todos los ráster cargados en el proyecto, seleccionamos las bandas adecuadas según la expresión del índice (pulsando Doble-clic sobre el nombre de la banda en el apartado Bandas ráster, o sobre el operador en el apartado Operadores). La expresión definida se muestra en el apartado Expresión de la calculadora ráster. Podemos indicar el nombre de la capa de salida y el formato, o bien que se genere una capa temporal, en el apartado Capa de resultado.

Selección de la rampa de color.

Figura 6.10: Selección de la rampa de color.

Para mostrar los resultados de una forma más vistosa, podemos cambiar la rampa de color del ráster obtenido. Esa operación se realiza en el apartado Simbología de las propiedades del ráster (figura 6.10). En este caso, se ha seleccionado la rampa de color RdYlGn, en la que corresponden los colores verdes a los valores más altos del índice.

6.3.2.2 Cálculo en R

red <- raster("lanjaron-jsamos-T30SVF_20210816T105621_B04_10m.tif")
nir <- raster("lanjaron-jsamos-T30SVF_20210816T105621_B08_10m.tif")

ndvi <- (nir - red) / (nir + red)

Para calcular el índice de vegetación en R, leemos las bandas que intervienen mediante la función raster y operamos con ellas directamente definiendo la expresión necesaria ya que estamos operando con matrices. El resultado es una banda que podemos presentar como cualquier otra.

Ejercicio 6.3 Obtén el color natural y el índice de vegetación para tu municipio (mediante QGIS o R, solo uno de ellos). Almacena los resultados en formato GeoTIFF, en archivos cuyo nombre tenga 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, incluye una captura de pantalla del entorno donde hayas realizado las operaciones donde se muestre cada uno de los resultados obtenidos.

6.4 Análisis de superficies

En este apartado trabajaremos con el MDT del municipio y realizaremos las operaciones desde QGIS.

Las operaciones que vamos a realizar sobre el MDT se pueden llevar a cabo de varias formas. Podemos acceder a la caja de herramientas pulsado sobre Procesos > Caja de herramientas.

Caja de herramientas de Procesos.

Figura 6.11: Caja de herramientas de Procesos.

Dentro de la caja de herramientas, tenemos dos posibilidades, acceder al apartado Análisis del terreno ráster o bien al apartado GDAL > Análisis ráster (figura 6.11). Adicionalmente, podemos acceder a algunas de las operaciones directamente desde el menú de QGIS, pulsando sobre Ráster > Análisis.

6.4.1 Mapa de sombras

Para calcular el mapa de sombras sobre el MDT desde el menú de QGIS, pulsamos sobre Ráster > Análisis > Mapa de sombras (Hillshade)…

Definición del mapa de sombras.

Figura 6.12: Definición del mapa de sombras.

En la figura 6.12, se muestra la ventana de definición de los parámetros de la operación. Se puede indicar la posición de la luz en cuanto a altitud y azimut, en este caso se han dejado los valores por defecto. Para guardar el resultado, podemos indicar el nombre de un archivo en el campo Mapa de sombras (Hillshade). En el apartado Llamada a la consola de GDAL/OGR se puede ver el comando que se va formando en función de las opciones que seleccionemos.

Al ejecutar la operación, obtenemos el mapa de sombras. Este mapa se suele usar para presentar datos (lo usaremos con ese fin). Por ejemplo, últimamente, en el apartado “`El Tiempo”’ de RTVE suelen usar el mapa de sombras de España como fondo para presentar las previsiones meteorológicas.

6.4.2 Orientación

Para obtener la orientación del terreno desde el menú de QGIS, pulsamos sobre Ráster > Análisis > Orientación… En algunas versiones de QGIS, en lugar de Orientación… aparece Aspecto…

Es frecuente expresar la orientación mediante el azimut, por lo que podemos dejar los valores por defecto de los parámetros. Por la configuración seleccionada para la operación y los datos del terreno, podemos comprobar que los valores están comprendidos en el intervalo [0, 360] para la nueva capa generada (podemos verlo en el apartado Información de las propiedades de la capa o en la ventana Capas).

6.4.3 Pendiente

Para obtener la pendiente del terreno desde el menú de QGIS, pulsamos sobre Ráster > Análisis > Pendiente…

Por defecto expresa la pendiente en grados, aunque también se puede indicar que la exprese mediante un porcentaje. Los parámetros de configuración de esta operación son adecuados. Si hemos indicado que exprese el resultado en grados, los valores deberían estar en el intervalo [0, 90].

Ejercicio 6.4 Utilizando el MDT rectangular de tu municipio, obtén las capas siguientes:

  • mapa de sombras,
  • orientación,
  • pendiente.

Almacena los resultados en formato GeoTIFF, en archivos cuyo nombre tenga 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, incluye una captura de pantalla del entorno donde hayas realizado las operaciones donde se muestre cada uno de los resultados obtenidos.

6.5 Transformación de datos ráster para análisis

A partir del MDT vamos a generar capas vectoriales de curvas (curvas de nivel) y puntos aleatorios a los que asociaremos su correspondiente altitud.

6.5.1 Curvas de nivel

Para obtener las curvas de nivel, pulsamos sobre Ráster > Extracción > Curvas de nivel…

En la ventana de definición de parámetros de la operación, es muy importante el campo Intervalo entre curvas de nivel, que determina la separación y, por tanto, el número de estas. Podemos dejar el valor por defecto o cambiarlo según nuestro propio criterio (según lo accidentado que sea el terreno).

Si miramos la tabla de atributos de las curvas generadas, para cada una tenemos un atributo ELEV con la altitud.

6.5.2 Puntos aleatorios

Para obtener un conjunto de puntos aleatorios a partir del MDT, pulsamos sobre Vectorial > Herramientas de investigación > Puntos aleatorios en la extensión…

Obtención de puntos aleatorios.

Figura 6.13: Obtención de puntos aleatorios.

En la ventana de definición de parámetros de la operación (figura 6.13), en el campo Extensión de entrada seleccionamos la opción Calcular a partir de la capa… y seleccionamos el MDT; podemos indicar el número de puntos (en este caso se ha indicado el valor 1000); también se puede definir una restricción en cuanto a la distancia mínima entre puntos en metros, que en este caso no se ha indicado el valor 50.

Como resultado, se han generado los puntos aleatorios en una nueva capa vectorial, dentro de los límites definidos por el MDT, sin embargo, los puntos no tienen asociada la altitud definida en su posición respecto al MDT. Para ello, necesitamos utilizar el complemento Point sampling tool.

Obtención del valor asociado a los puntos aleatorios.

Figura 6.14: Obtención del valor asociado a los puntos aleatorios.

Una vez instalado, accedemos a él mediante su icono en la barra de herramientas. En la ventana de definición de parámetros de la operación (figura 6.14), definimos la capa que contiene los puntos aleatorios y la lista de capas de las que tomar valores. También indicamos el archivo donde almacenar la nueva capa. En la nueva capa generada, los puntos tienen asociado el valor tomado del MDT.

Ejercicio 6.5 Utilizando el MDT rectangular de tu municipio, obtén las capas siguientes:

  • curvas de nivel,
  • puntos aleatorios con altitud.

Almacena los resultados vectoriales en formato GeoPackage (todos en el mismo archivo), en un archivo cuyo nombre tenga 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, incluye una captura de pantalla del entorno donde hayas realizado las operaciones donde se muestre cada uno de los resultados obtenidos.

Ejercicio 6.6 A partir de las capas vectoriales disponibles para tu municipio busca cómo se hace la operación con QGIS o con R y genera el MDT o MDE.

Para documentar la realización del ejercicio, describe los pasos realizados e incluye las capturas de pantalla del entorno donde hayas realizado las operaciones que consideres necesarias para documentar la operación. Referencia adecuadamente las fuentes que hayas utilizado.

6.6 Presentación de datos

En este apartado, vamos a utilizar el MDT y el mapa de sombras para realizar una representación del terreno con una rampa de color adecuada.

En el apartado Simbología de las propiedades del MDT, en el campo Tipo de renderizador seleccionamos Pseudocolor monobanda, en el campo Rampa de color seleccionamos la opción Crear nueva rampa de color… Se abre la ventana Tipo de rampa de color y seleccionamos el valor Catálogo: cpt:city. Para esta rampa de color, seleccionamos el nombre Topography y la paleta cd-a.

De nuevo en la ventana del apartado Simbología, en el apartado Representación de capas, en el campo Modo de mezcla, elegimos la opción Multiplicar.

Presentación con el MDT y el mapa de sombras.

Figura 6.15: Presentación con el MDT y el mapa de sombras.

El resultado obtenido se puede apreciar en la figura 6.15. El orden de las capas en el proyecto es importante.

Ejercicio 6.7 Utilizando el MDT rectangular de tu municipio y la capa sombras obtenida en el ejercico 6.4 obtén una representación similar a la de la figura 6.15.

Para documentar la realización del ejercicio, incluye una captura de pantalla del entorno donde hayas realizado las operaciones donde se muestre el resultado obtenido.

Bibliografía

Marquez, Angel. 2015. PostGIS Essentials. Packt Publishing.

  1. Se pueden encontrar más detalles en https://eos.com/es/make-an-analysis/ndvi/.↩︎