Acceso a QGIS, GDAL, SAGA o GRASS desde R

por

La capacidad de R para interactuar con SIG dedicados le brinda unas capacidades geoespaciales asombrosas, llegando a veces a superar a los SIG en el modelado estadístico espacial y en algunas áreas de geocomputación, como la creación de mapas interactivos.

Como os contábamos en entradas anteriores, podemos emplear el plugin Processing R Provider para  integrar R dentro de QGIS y ejecutar nuestros algoritmos R a través de interfaces gráficas amigables en este SIG. Pero, ¿y si queremos hacer lo contrario? ¿Podemos aprovechar la amplia funcionalidad de QGIS e integrarla en nuestros scripts en R? ¡La respuesta es sí!

En este artículo, vamos a aprender a emplear la librería qgisprocess para establecer ‘puentes’ desde R a tres productos SIG Open Source que todo desarrollador especializado en los sistemas de información geográfica debe conocer: QGIS, GDAL, SAGA y GRASS.

1. Conexión R – QGIS

Desde su versión 3.14 en adelante, QGIS proporciona el API de desarrollo qgis_process, que permite ejecutar algoritmos de QGIS y de terceros desde una línea de comandos. Podéis consultar toda la información referente a este API en la propia documentación de QGIS.

Una API (interfaz de programación de aplicaciones) es un intermediario entre sistemas que permite a diferentes aplicaciones comunicarse entre sí y compartir información y funcionalidades.

La librería de la que os hablamos hoy, qgisprocess, se nutre de este API haciendo posible llamar desde nuestras sesiones de R a la mayoría de algoritmos que usamos desde el propio QGIS: algoritmos nativos de QGIS, GDAL, GRASS, SAGA…

2. Instalación y configuración de qgisprocess

El primer paso, como siempre, es instalar y cargar la librería. Actualmente no hay una versión estable de qgisprocess en CRAN, así que debemos instalar la versión de desarrollo del repositorio GitHub.

## Instalamos la librería qgisprocess ubicada en Github.
# Para ello se requiere el paquete 'remotes', que debe instalarse previamente si no se ha hecho ya.
# install.packages("remotes")
remotes::install_github("paleolimbot/qgisprocess")
library(qgisprocess)


# > Attempting to load the cache ... Success! 
# > QGIS version: 3.28.4-Firenze 
# > Having access to 1090 algorithms from 7 QGIS processing providers.

Cuando carguemos la librería, si la instalación ha ido bien, veremos en la consola nuestra versión de QGIS y el número de proveedores al que tenemos acceso. Si tenemos más de una versión de QGIS instalada o si al cargar el paquete no detecta automáticamente nuestra instalación de QGIS, podemos configurar la ruta a nuestro SIG de la siguiente forma:

# Configuración de ruta de acceso al API de QGIS -> options(qgisprocess.path ="tu_ruta"). Por ejemplo:
options(qgisprocess.path = "/Applications/QGIS-LTR.app/Contents/MacOS/bin/qgis_process")

La librería incluye las siguientes funciones muy útiles:

  • qgis_path() para conocer la ruta configurada de acceso al API de QGIS.
  • qgis_version() para consultar la versión de QGIS configurada.
  • qgis_providers() para obtener información sobre los proveedores y el número de algoritmos que están disponibles para ser utilizados desde nuestra sesión de R.
  • qgis_algorithms() para obtener información detallada de todos los algoritmos a los que tenemos acceso desde nuestra sesión R.
## Ruta al API de QGIS
qgis_path()
# [1] "/Applications/QGIS-LTR.app/Contents/MacOS/bin/qgis_process"

## Versión de QGIS
qgis_version()
# [1] "3.28.4-Firenze"

## Proveedores y el nº algoritmos disponibles por proveedor
qgis_providers()
#   A tibble: 7 × 3
#  provider provider_title    algorithm_count
#  <chr>    <chr>                       <int>
#  1 gdal     GDAL                           56
#  2 grass7   GRASS                         304
#  3 qgis     QGIS                           50
#  4 3d       QGIS (3D)                       1
#  5 native   QGIS (native c++)             242
#  6 r        R                              71
#  7 saga     SAGA                          366

¡Ya estamos listos para empezar a usar los algoritmos que brinda QGIS desde nuestra sesión de R! Para entender cómo funciona qgisprocess, vamos a ver varios ejemplos prácticos en los que usaremos herramientas nativas de QGIS y algoritmos no nativos.

3. Uso de algoritmos nativos de QGIS desde R

En este primer ejemplo, vamos a calcular un área de influencia o buffer sobre un capa vectorial del entrada. Para conocer el nombre de los geoprocesos a los que tenemos acceso, podemos usar el siguiente truco:

## Búsqueda de algoritmos que tengan buffer en su nombre
herramientas = qgis_algorithms()
grep("buffer", herramientas$algorithm, value = TRUE)

> grep("buffer", herramientas$algorithm, value = TRUE) 
# [1] "gdal:buffervectors" "gdal:onesidebuffer" "grass7:r.buffer" 
# [4] "grass7:r.buffer.lowmem" "grass7:v.buffer" "native:buffer" 
# [7] "native:bufferbym" "native:multiringconstantbuffer" "native:singlesidedbuffer" 
# [10] "native:taperedbuffer" "native:wedgebuffers" "qgis:variabledistancebuffer" 
# [13] "saga:fixeddistancebuffer" "saga:rasterbuffer" "saga:rasterproximitybuffer" 
# [16] "saga:thresholdrasterbuffer" "saga:variabledistancebuffer"

En mi caso, he obtenido 18 algoritmos que contienen la palabra «buffer» en su nombre. Los algoritmos propios de QGIS se corresponden con aquellos que se llaman «native:__» . Así que vamos a empezar usando la herramienta nativa de geoprocesamiento de QGIS «native:buffer».

Antes de ejecutar un proceso, es imprescindible obtener información sobre el algoritmo específico qué queremos usar. Para ello usaremos qgis_show_help().

> qgis_show_help("native:buffer")
Buffer (native:buffer)

----------------
Description
----------------
This algorithm computes a buffer area for all the features in an input layer, using a fixed or dynamic distance.

The segments parameter controls the number of line segments to use to approximate a quarter circle when creating rounded offsets.

The end cap style parameter controls how line endings are handled in the buffer.

The join style parameter specifies whether round, miter or beveled joins should be used when offsetting corners in a line.

The miter limit parameter is only applicable for miter join styles, and controls the maximum distance from the offset curve to use when creating a mitered join.

----------------
Arguments
----------------

INPUT: Input layer
    Argument type:	source
    Acceptable values:
        - Path to a vector layer
DISTANCE: Distance
    Default value:	10
    Argument type:	distance
    Acceptable values:
        - A numeric value
SEGMENTS: Segments
    Default value:	5
    The segments parameter controls the number of line segments to use to approximate a quarter circle when creating rounded offsets.
    Argument type:	number
    Acceptable values:
        - A numeric value
END_CAP_STYLE: End cap style
    Default value:	0
    Argument type:	enum
    Available values:
        - 0: Round
        - 1: Flat
        - 2: Square
    Acceptable values:
        - Number of selected option, e.g. '1'
        - Comma separated list of options, e.g. '1,3'
JOIN_STYLE: Join style
    Default value:	0
    Argument type:	enum
    Available values:
        - 0: Round
        - 1: Miter
        - 2: Bevel
    Acceptable values:
        - Number of selected option, e.g. '1'
        - Comma separated list of options, e.g. '1,3'
MITER_LIMIT: Miter limit
    Default value:	2
    Argument type:	number
    Acceptable values:
        - A numeric value
DISSOLVE: Dissolve result
    Default value:	false
    Argument type:	boolean
    Acceptable values:
        - 1 for true/yes
        - 0 for false/no
OUTPUT: Buffered
    Argument type:	sink
    Acceptable values:
        - Path for new vector layer

----------------
Outputs
----------------

OUTPUT: <outputVector>
    Buffered

Este comando nos dará toda la información necesaria para poder ejecutar cualquier proceso: descripción del algoritmo, parámetros de entrada que necesitamos, resultado esperado… etc.

Vamos a ver un primer ejemplo de uso de esta herramienta para generar un área de influencia de 100 metros de las coordenadas de los centros docentes de la ciudad de Cáceres, obtenidos del Portal de Datos Abiertos del Ayuntamiento de Cáceres. 

## Cargamos los datos de los centros docentes de Cáceres
# install.packages("sf")
library(sf)

centros_docentes <- st_read("http://opendata.caceres.es/GetData/GetData?dataset=schema:EducationalOrganization&format=geojson")
centros_docentes_25830 <-  st_transform(centros_docentes, 25830) # Proyectamos los datos a ETRS89 para trabajar con metros en lugar de grados.

# plot(st_geometry(centros_docentes), axes=TRUE, col="red", main="Centros docentes de Cáceres")

La función que nos permite ejecutar cualquier algoritmo es qgis_run_algorithm(). Esta función requiere los siguientes argumentos:

  • algorithm: el nombre de la herramienta a emplear.
  • … : el resto de argumentos dependen del propio algoritmo que vamos a ejecutar. En el caso de «native:buffer», los argumentos que debemos usar, tal y como nos lo ha mostrado qgis_show_help(«native:buffer»), son INPUT, DISTANCE, SEGMENTS, END_CAP_STYLE, JOIN_STYLE, OUTPUT… Es muy importante consultar la información sobre cada algoritmo para emplear los argumentos correctos y no cometer errores.
## Creación de un fichero temporal para guardar el resultados del cálculo del área de influencia
fichero_salida_temp <- file.path(tempdir(), "centros_docentes_bff_100.gpkg")

## Ejecución
qgis_run_algorithm(
  algorithm = "native:buffer",
  INPUT = centros_docentes_25830,
  DISTANCE = 100,
  OUTPUT = fichero_salida_temp
)

## Lectura en nuestra sesión el resultado del geoproceso
centros_docentes_bff_100 <- st_read(fichero_salida_temp)

## Mostramos un colegio con su buffer
colegio_san_jose <- centros_docentes_25830[(centros_docentes_25830$rdfs_label %in% "Colegio San Jose"),]
colegio_san_jose_bff <- centros_docentes_bff_100[(centros_docentes_bff_100$rdfs_label %in% "Colegio San Jose"),]
plot(st_geometry(colegio_san_jose_bff), col="green", main = "Buffer a 100m del Colegio San José")
plot(st_geometry(colegio_san_jose), pch=18, col="red", add=TRUE)

4. Uso de algoritmos no nativos de QGIS desde R: SAGA, GRASS…

Como comentábamos al principio de este artículo, el API qgis_process nos permite acceder a las herramientas nativas de QGIS y a todos los algoritmos de los proveedores a los que este tiene acceso.

Para finalizar, vamos a ver cómo usar un algoritmo no nativo de QGIS, donde el funcionamiento es exactamente el mismo. En este caso, vamos a usar la herramienta «saga:clippointswithpolygons»  de SAGA para obtener el número de bibliotecas que se encuentran a menos de 100 metros de alguno de los centros docentes.

> ## Información sobre el algoritmo a usar
> qgis_show_help("saga:clippointswithpolygons")
Clip points with polygons (saga:clippointswithpolygons)

----------------
Description
----------------

----------------
Arguments
----------------

POINTS: Points
    Argument type:	source
    Acceptable values:
        - Path to a vector layer
POLYGONS: Polygons
    Argument type:	source
    Acceptable values:
        - Path to a vector layer
FIELD: Add Attribute to Clipped Points
    Argument type:	field
    Acceptable values:
        - The name of an existing field
        - ; delimited list of existing field names
METHOD: Clipping Options
    Argument type:	enum
    Available values:
        - 0: [0] one layer for all points
        - 1: [1] separate layer for each polygon
    Acceptable values:
        - Number of selected option, e.g. '1'
        - Comma separated list of options, e.g. '1,3'
CLIPS: Clipped Points
    Argument type:	vectorDestination
    Acceptable values:
        - Path for new vector layer

----------------
Outputs
----------------

CLIPS: <outputVector>
    Clipped Points
## Creación de un fichero temporal para guardar los resultados del recorte
fichero_salida_temp2 <- file.path(tempdir(), "recorte_bff_biblio.geojson")

## Ejecución
qgis_run_algorithm(
  algorithm = "saga:clippointswithpolygons",
  POINTS = bibliotecas_25830,
  POLYGONS = centros_docentes_bff_100,
  METHOD = 0,
  FIELD = "rdfs_label",
  CLIPS = fichero_salida_temp2
)
## Lectura en nuestra sesión el resultado del geoproceso
recorte_bff_biblio <- st_read(fichero_salida_temp2)

Podemos ver que únicamente tres bibliotecas se encuentran a menos de 100 metros de los centros docentes.

> recorte_bff_biblio$rdfs_label
[1] "Zamora Vicente" "Zamora Vicente" "Central UEX"

5. Conclusiones

La librería qgisprocess nos permite el acceso a más de 900 geoalgoritmos y nos facilita considerablemente el trabajo de automatización de nuestros geoprocesos. No obstante, hay que tener en cuenta que actualmente es una versión en desarrollo. Esto implica que aún no existe una versión estable y que podemos encontrarnos con ciertos errores en el camino.

Si te ha gustado esta entrada y quieres aprender a trabajar con datos espaciales en el entrono de R, inscríbete ya en el curso online de R y GIS: usa R como un GIS.

3 comentarios en «Acceso a QGIS, GDAL, SAGA o GRASS desde R»

  2. Estupendo post, muchas gracias por la explicación esto es de gran ayuda a la hora de llamar funciones de QGIS directamente desde R y automatizar y documentar procesos.

Los comentarios están cerrados.