Interfacing QGIS spatial processing algorithms from R

# .center[Interfacing QGIS spatial processing algorithms from R]

.pull-left[.footnote[
<img src="images/useR Salzburg 2024 logo.png" width="20%" hspace="20vw"/>
<img src="images/flanders.png" width="10%" hspace="20vw"/>
<img src="images/inbo.jpg" width="15%"/>
]]

## .pull-right[.footnote[.right[Floris Vanderhaeghe<br/>Dewey Dunnington<br/>Nyall Dawson<br/>Jan Caha<br/>Jannes Muenchow<br/>Jakub Nowosad<br/>Robin Lovelace<br/><br/>[html slides](https://florisvdh.github.io/user-2024-qgisprocess/)<br/>[<svg aria-label="GitHub" role="img" viewBox="0 0 496 512" style="height:1em;width:0.97em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><title>GitHub</title><path d="M165.9 397.4c0 2-2.3 3.6-5.2 3.6-3.3.3-5.6-1.3-5.6-3.6 0-2 2.3-3.6 5.2-3.6 3-.3 5.6 1.3 5.6 3.6zm-31.1-4.5c-.7 2 1.3 4.3 4.3 4.9 2.6 1 5.6 0 6.2-2s-1.3-4.3-4.3-5.2c-2.6-.7-5.5.3-6.2 2.3zm44.2-1.7c-2.9.7-4.9 2.6-4.6 4.9.3 2 2.9 3.3 5.9 2.6 2.9-.7 4.9-2.6 4.6-4.6-.3-1.9-3-3.2-5.9-2.9zM244.8 8C106.1 8 0 113.3 0 252c0 110.9 69.8 205.8 169.5 239.2 12.8 2.3 17.3-5.6 17.3-12.1 0-6.2-.3-40.4-.3-61.4 0 0-70 15-84.7-29.8 0 0-11.4-29.1-27.8-36.6 0 0-22.9-15.7 1.6-15.4 0 0 24.9 2 38.6 25.8 21.9 38.6 58.6 27.5 72.9 20.9 2.3-16 8.8-27.1 16-33.7-55.9-6.2-112.3-14.3-112.3-110.5 0-27.5 7.6-41.3 23.6-58.9-2.6-6.5-11.1-33.3 2.6-67.9 20.9-6.5 69 27 69 27 20-5.6 41.5-8.5 62.8-8.5s42.8 2.9 62.8 8.5c0 0 48.1-33.6 69-27 13.7 34.7 5.2 61.4 2.6 67.9 16 17.7 25.8 31.5 25.8 58.9 0 96.5-58.9 104.2-114.8 110.5 9.2 7.9 17 22.9 17 46.4 0 33.7-.3 75.4-.3 83.6 0 6.5 4.6 14.4 17.3 12.1C428.2 457.8 496 362.9 496 252 496 113.3 383.5 8 244.8 8zM97.2 352.9c-1.3 1-1 3.3.7 5.2 1.6 1.6 3.9 2.3 5.2 1 1.3-1 1-3.3-.7-5.2-1.6-1.6-3.9-2.3-5.2-1zm-10.8-8.1c-.7 1.3.3 2.9 2.3 3.9 1.6 1 3.6.7 4.3-.7.7-1.3-.3-2.9-2.3-3.9-2-.6-3.6-.3-4.3.7zm32.4 35.6c-1.6 1.3-1 4.3 1.3 6.2 2.3 2.3 5.2 2.6 6.5 1 1.3-1.3.7-4.3-1.3-6.2-2.2-2.3-5.2-2.6-6.5-1zm-11.4-14.7c-1.6 1-1.6 3.6 0 5.9 1.6 2.3 4.3 3.3 5.6 2.3 1.6-1.3 1.6-3.9 0-6.2-1.4-2.3-4-3.3-5.6-2z"/></svg>](https://github.com/florisvdh) [<svg aria-label="Mastodon" role="img" viewBox="0 0 448 512" style="height:1em;width:0.88em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><title>Mastodon</title><path d="M433 179.11c0-97.2-63.71-125.7-63.71-125.7-62.52-28.7-228.56-28.4-290.48 0 0 0-63.72 28.5-63.72 125.7 0 115.7-6.6 259.4 105.63 289.1 40.51 10.7 75.32 13 103.33 11.4 50.81-2.8 79.32-18.1 79.32-18.1l-1.7-36.9s-36.31 11.4-77.12 10.1c-40.41-1.4-83-4.4-89.63-54a102.54 102.54 0 0 1-.9-13.9c85.63 20.9 158.65 9.1 178.75 6.7 56.12-6.7 105-41.3 111.23-72.9 9.8-49.8 9-121.5 9-121.5zm-75.12 125.2h-46.63v-114.2c0-49.7-64-51.6-64 6.9v62.5h-46.33V197c0-58.5-64-56.6-64-6.9v114.2H90.19c0-122.1-5.2-147.9 18.41-175 25.9-28.9 79.82-30.8 103.83 6.1l11.6 19.5 11.6-19.5c24.11-37.1 78.12-34.8 103.83-6.1 23.71 27.3 18.4 53 18.4 175z"/></svg>](https://fosstodon.org/@florisvdh)  [<svg aria-label="ORCID" role="img" viewBox="0 0 512 512" style="height:1em;width:1em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><title>ORCID</title><path d="M294.75 188.19h-45.92V342h47.47c67.62 0 83.12-51.34 83.12-76.91 0-41.64-26.54-76.9-84.67-76.9zM256 8C119 8 8 119 8 256s111 248 248 248 248-111 248-248S393 8 256 8zm-80.79 360.76h-29.84v-207.5h29.84zm-14.92-231.14a19.57 19.57 0 1 1 19.57-19.57 19.64 19.64 0 0 1-19.57 19.57zM300 369h-81V161.26h80.6c76.73 0 110.44 54.83 110.44 103.85C410 318.39 368.38 369 300 369z"/></svg>](https://orcid.org/0000-0002-6378-6229)]]]

---
class: inverse, middle, center

# About QGIS

---
layout: false

# QGIS

Popular open-source geographic information system (GIS) program to:

- view geospatial data interactively

- create and edit spatial data layers

- perform advanced geoprocessing with algorithms of QGIS, GRASS GIS, GDAL ...

- make sophisticated maps, atlases and reports

---

---
class: inverse, middle, center

# Why interface QGIS processing from R?

---
layout: false

# Main advantages

- Expand geospatial processing abilities in R

- Still have all other processing capabilities of R

- Reproducibility

- No QGIS project needed

- Unified interface to QGIS, GRASS GIS, SAGA, GDAL and other processing providers (1000+ geospatial algorithms)

<br/>

.center[
<img src="images/qgis.svg" width="10%" hspace="15vw"/>
<img src="images/grass.svg" width="10%" hspace="15vw"/>
<img src="images/saga.png" width="10%" hspace="15vw"/>
<img src="images/gdal.svg" width="10%"/>
]

---
layout: true

# Interfacing QGIS processing from R

---

- Formerly: packages **RQGIS** and  **RQGIS3** (Muenchow & Schratz) for QGIS 2 and 3

- use the QGIS Python API

- set QGIS environment variables
  
  - **RQGIS3** was hard to maintain wrt provider changes in QGIS; also, crashes were observed in RStudio IDE
  
  - no longer developed
  
---

- Since QGIS 3.16: packages **qgisprocess** & **qgis**

- use the recent standalone `qgis_process` shell command from QGIS = a unified entry to all providers and algorithms!

- no more QGIS environment variables needed
  
  - actively developed

---

```
$ qgis_process
QGIS Processing Executor - 3.38.0-Grenoble 'Grenoble' (3.38.0-Grenoble)
Usage: /usr/bin/qgis_process.bin [--help] [--version] [--json] [--verbose]
          [--no-python] [--skip-loading-plugins]
          [command] [algorithm id, path to model file, or path to Python script] [parameters]

...
```

`$ qgis_process run <algorithm> [parameters]`

`$ qgis_process run <algorithm> -`

`$ qgis_process plugins`

---
layout: false

# Current R-packages

- **qgisprocess** (Dewey Dunnington et al.): direct interface to `qgis_process`
  - https://r-spatial.github.io/qgisprocess/

``` r
install.packages("qgisprocess")
```
  
- **qgis** (Jan Caha): functions for each algorithm; uses **qgisprocess** package
  - https://jancaha.github.io/r_package_qgis/

``` r
remotes::install_github("jancaha/qgis")
```

---
class: inverse, center, middle

# A simple example

---
layout: true

# Basic usage

---

---

Load the R package:

``` r
library(qgisprocess)
```

```
#> Attempting to load the package cache ... Success!
#> QGIS version: 3.38.0-Grenoble
#> Having access to 2059 algorithms from 18 QGIS processing providers.
#> Run `qgis_configure(use_cached_data = TRUE)` to reload cache and get more details.
```

---

``` r
seine_path <- "data/seine.gpkg"
```

``` r
(seine <- sf::read_sf(seine_path))
```

```
#> Simple feature collection with 3 features and 1 field
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 518344.7 ymin: 6660431 xmax: 879955.3 ymax: 6938864
#> Projected CRS: RGF93 v1 / Lambert-93
#> # A tibble: 3 × 2
#>   name                                                                      geom
#>   <chr>                                                    <MULTILINESTRING [m]>
#> 1 Marne ((879955.3 6755725, 878440.9 6755688, 876653.8 6756227, 874212.2 675791…
#> 2 Seine ((828893.6 6713873, 828216.3 6715450, 827937.9 6716999, 828199.2 671851…
#> 3 Yonne ((773482.1 6660431, 771342.9 6665712, 771043 6667566, 770931.7 6669151,…
```

---

Run algorithm:

``` r
result <- qgis_run_algorithm(
  algorithm = "native:pointsalonglines",
  INPUT = seine_path,
  DISTANCE = 1e4
)
#> Argument `START_OFFSET` is unspecified (using QGIS default value).
#> Argument `END_OFFSET` is unspecified (using QGIS default value).
#> Using `OUTPUT = qgis_tmp_vector()`
```

or:

``` r
result <- qgis::qgis_pointsalonglines(
  INPUT = seine_path,
  DISTANCE = 1e4
)
```
 
---

Extract output (defaults to the `"OUTPUT"` element of `result`).

``` r
qgis_extract_output(result)
```

```
#> [1] "/tmp/RtmpUtuVpc/file5288578598f1/file52883dccf462.gpkg"
#> attr(,"class")
#> [1] "qgis_outputVector"
```

---

---
layout: false
class: inverse, center, middle

# Convenient functions

---
layout: true

# Finding algorithms

---

``` r
qgis_search_algorithms(algorithm = "point.*line")
```

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> provider </th>
   <th style="text-align:left;"> algorithm </th>
   <th style="text-align:left;"> algorithm_title </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> gdal </td>
   <td style="text-align:left;"> gdal:pointsalonglines </td>
   <td style="text-align:left;"> Points along lines </td>
  </tr>
  <tr>
   <td style="text-align:left;"> native </td>
   <td style="text-align:left;"> native:interpolatepoint </td>
   <td style="text-align:left;"> Interpolate point on line </td>
  </tr>
  <tr>
   <td style="text-align:left;"> native </td>
   <td style="text-align:left;"> native:pointsalonglines </td>
   <td style="text-align:left;"> Points along geometry </td>
  </tr>
  <tr>
   <td style="text-align:left;"> native </td>
   <td style="text-align:left;"> native:randompointsonlines </td>
   <td style="text-align:left;"> Random points on lines </td>
  </tr>
  <tr>
   <td style="text-align:left;"> qgis </td>
   <td style="text-align:left;"> qgis:generatepointspixelcentroidsalongline </td>
   <td style="text-align:left;"> Generate points (pixel centroids) along line </td>
  </tr>
  <tr>
   <td style="text-align:left;"> qgis </td>
   <td style="text-align:left;"> qgis:randompointsalongline </td>
   <td style="text-align:left;"> Random points along line </td>
  </tr>
  <tr>
   <td style="text-align:left;"> sagang </td>
   <td style="text-align:left;"> sagang:convertpointstolines </td>
   <td style="text-align:left;"> Convert points to line(s) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> sagang </td>
   <td style="text-align:left;"> sagang:pointtolinedistances </td>
   <td style="text-align:left;"> Point to line distances </td>
  </tr>
  <tr>
   <td style="text-align:left;"> sagang </td>
   <td style="text-align:left;"> sagang:snappointstolines </td>
   <td style="text-align:left;"> Snap points to lines </td>
  </tr>
</tbody>
</table>

---

``` r
qgis_providers()
```

```
#> # A tibble: 18 × 3
#>    provider         provider_title                         algorithm_count
#>    <chr>            <chr>                                            <int>
#>  1 cartographytools Cartography tools                                    5
#>  2 gdal             GDAL                                                56
#>  3 grass            GRASS                                              307
#>  4 latlontools      Lat Lon tools                                       10
#>  5 model            Models                                               3
#>  6 NetworkGT        NetworkGT                                           33
#>  7 Networks         Networks                                            48
#>  8 pcraster         PCRaster                                           102
#>  9 qgis             QGIS                                                42
#> 10 3d               QGIS (3D)                                            1
#> 11 native           QGIS (native c++)                                  263
#> 12 pdal             QGIS (PDAL)                                         17
#> 13 qneat3           QNEAT3 - Qgis Network Analysis Toolbox              14
#> 14 sagang           SAGA Next Gen                                      587
#> 15 script           Scripts                                              1
#> 16 Valhalla         Valhalla                                            20
#> 17 visibility       Visibility analysis                                  4
#> 18 wbt              WhiteboxTools                                      546
```

---

``` r
qgis_plugins()
```

```
#> # A tibble: 12 × 2
#>    name                    enabled
#>    <chr>                   <lgl>  
#>  1 QNEAT3                  TRUE   
#>  2 ViewshedAnalysis        TRUE   
#>  3 cartography_tools       TRUE   
#>  4 grassprovider           TRUE   
#>  5 latlontools             TRUE   
#>  6 network_gt              TRUE   
#>  7 networks                TRUE   
#>  8 pcraster_tools          TRUE   
#>  9 processing              TRUE   
#> 10 processing_saga_nextgen TRUE   
#> 11 valhalla                TRUE   
#> 12 wbt_for_qgis            TRUE
```

---
layout:false

# Algorithm documentation

``` r
qgis_show_help("native:pointsalonglines")
```

``` r
## Points along geometry (native:pointsalonglines)
## 
## ----------------
## Description
## ----------------
## Creates regularly spaced points along line features.
## This algorithm creates a points layer, with points distributed along the lines of an input vector layer. The distance between points (measured along the line) is defined as a parameter.
## 
## Start and end offset distances can be defined, so the first and last point will not fall exactly on the line's first and last nodes. These start and end offsets are defined as distances, measured along the line from the first and last nodes of the lines.
## 
## ----------------
## Arguments
## ----------------
## 
## INPUT: Input layer
##     Argument type:    source
##     Acceptable values:
##         - Path to a vector layer
## DISTANCE: Distance
##     Default value:    1
## .....
```

---
layout: true

# Supports various R objects as input arguments

---

Spatial QGIS argument types

<br/>

---

Non-spatial QGIS argument types:

<br/>

---

``` r
library(terra)
elev <- rast(system.file("ex/elev.tif", package = "terra"))
class(elev)
```

```
#> [1] "SpatRaster"
#> attr(,"package")
#> [1] "terra"
```

``` r
qgis_run_algorithm("native:rasterlayerstatistics", INPUT = elev, BAND = 1)
```

```
#> <Result of `qgis_run_algorithm("native:rasterlayerstatistics", ...)`>
#> List of 9
#>  $ COUNT           : num 4608
#>  $ MAX             : num 547
#>  $ MEAN            : num 348
#>  $ MIN             : num 141
#>  $ OUTPUT_HTML_FILE: 'qgis_outputHtml' chr "/tmp/RtmpUtuVpc/file5288578598f1/file52887ef8bb06"
#>  $ RANGE           : num 406
#>  $ STD_DEV         : num 80.2
#>  $ SUM             : num 1605135
#>  $ SUM_OF_SQUARES  : num 29646349
```

---
layout: true

# Result handling

---

## Extracting elements from the result

By default, the `OUTPUT` element is selected by `qgis_extract_output()`.

- typically contains a file path

---

## Extracting elements from the result

Which output elements are generated by an algorithm?

``` r
qgis_get_output_specs("grass:r.flow")
```

```
#> # A tibble: 3 × 3
#>   name             description       qgis_output_type
#>   <chr>            <chr>             <chr>           
#> 1 flowaccumulation Flow accumulation outputRaster    
#> 2 flowlength       Flow path length  outputRaster    
#> 3 flowline         Flow line         outputVector
```

So, sometimes you need to specify the output name:

``` r
qgis_extract_output(result, "flowline")
```

---

## Coercing output to R objects

Result object or output element can be coerced to an R object:

``` r
sf::st_as_sf(result) # takes OUTPUT by default
```

``` r
result |> qgis_extract_output("flowline") |> sf::st_as_sf()
```

``` r
qgis_as_terra(result)
```

``` r
qgis_as_raster(result)
```

``` r
stars::st_as_stars(result)
```

---
layout:false
class: inverse, center, middle

# Extra!<br/><br/>Extra!

---

# Algorithm piping

``` r
seine_points_buffer <- seine |> 
  qgis_run_algorithm_p("native:pointsalonglines", DISTANCE = 1e4) |> 
  qgis_run_algorithm_p("native:buffer", DISTANCE = 1000,
                       OUTPUT = "data/buffer.gpkg")
seine_points_buffer
```

```
#> <Result of `qgis_run_algorithm("native:buffer", ...)`>
#> List of 1
#>  $ OUTPUT: 'qgis_outputVector' chr "data/buffer.gpkg"
```

Alternative with **qgis** package:

``` r
seine |> 
  qgis::qgis_pointsalonglines(DISTANCE = 1e4) |> 
  qgis_extract_output() |> 
  qgis::qgis_buffer(DISTANCE = 1000)
```

```
#> <Result of `qgis_run_algorithm("native:buffer", ...)`>
#> List of 1
#>  $ OUTPUT: 'qgis_outputVector' chr "/tmp/RtmpUtuVpc/file5288578598f1/file528896a5b20.gpkg"
```

---
layout:true

# Taking parameters from the QGIS GUI

---

---

``` r
jsonlite::prettify(json_from_qgis)
```

```
#> {
#>     "area_units": "m2",
#>     "distance_units": "meters",
#>     "ellipsoid": "EPSG:7030",
#>     "inputs": {
#>         "COLUMN_PREFIX": "dem_",
#>         "INPUT": "/home/floris/git_repositories/foss4g-2023-qgisprocess/data/buffer.gpkg|layername=buffer",
#>         "INPUT_RASTER": "/home/floris/git_repositories/foss4g-2023-qgisprocess/data/dem_clip.tif",
#>         "OUTPUT": "ogr:dbname='/home/floris/git_repositories/foss4g-2023-qgisprocess/data/zonal_stats.gpkg' table=\"seine_zonal\" (geom)",
#>         "RASTER_BAND": 1,
#>         "STATISTICS": [
#>             2,
#>             3,
#>             4,
#>             5,
#>             6
#>         ]
#>     }
#> }
#> 
```

---

``` r
zonal_stats_result <- qgis_run_algorithm(
  "native:zonalstatisticsfb", 
* .raw_json_input = json_from_qgis
)
```

``` r
zonal_stats_result
```

```
#> <Result of `qgis_run_algorithm("native:zonalstatisticsfb", ...)`>
#> List of 1
#>  $ OUTPUT: 'qgis_outputVector' chr "/home/floris/git_repositories/foss4g-2023-qgisprocess/data/zonal_stats.gpkg|layername=seine_zonal"
```

---
layout:false

# Result

``` r
zonal_stats <- sf::st_as_sf(zonal_stats_result) |> 
  sf::st_drop_geometry() |> 
  mutate(distance = set_units(distance, "m") |> set_units("km"))

zonal_stats
```

```
#> # A tibble: 123 × 8
#>    name  distance  angle dem_mean dem_median dem_stdev dem_min dem_max
#>    <chr>     [km]  <dbl>    <dbl>      <dbl>     <dbl>   <dbl>   <dbl>
#>  1 Marne        0 269.       372.       366.      22.0    344.    412.
#>  2 Marne       10 322.       340.       338       21.9    309.    416.
#>  3 Marne       20 327.       331.       334.      23.9    287.    382.
#>  4 Marne       30 332.       296.       280.      29.5    264.    359.
#>  5 Marne       40   5.76     286.       295.      24.3    238.    317.
#>  6 Marne       50 358.       269.       254.      45.4    221.    378.
#>  7 Marne       60  22.4      235.       220.      32.7    201.    312.
#>  8 Marne       70 350.       202.       188.      24.8    183.    306.
#>  9 Marne       80 335.       196.       184.      26.9    167.    276 
#> 10 Marne       90 334.       182.       184.      19.6    150.    209.
#> # ℹ 113 more rows
```

---
class: inverse, center, middle

# Online documentation

---
class: center

# <https://r-spatial.github.io/qgisprocess>

---
class: center

# <https://jancaha.github.io/r_package_qgis>

---
class: center

# Cheat sheet !

---
class: center

# Section in book 'Geospatial Computation with R'

.footnote[### <https://r.geocompx.org>]

---
class: inverse, center, middle

# Recent developments

---

# Recent developments

- Speed-up since QGIS 3.36 & **qgisprocess** 0.4.0

- Added vector support for R package **terra**

- Vignettes:

- which R objects are accepted as algorithm arguments?
  - how to configure behaviour with options or environment variables?
  - how to pass QGIS expressions?
  
--

- Improved handling of deprecated algorithms

---
class: inverse, center, middle

# Questions?<br/><br/>Ideas?

### .pull-right-wide[.footnote[.right[Slides: <https://florisvdh.github.io/user-2024-qgisprocess/>]]]