Interfacing QGIS processing from R

# .center[Interfacing QGIS processing from R]

.pull-left[.footnote[<img src="images/FOSS4G 2023 - Logo Main.png" width="30%"/>]]

## .pull-right-wide[.footnote[.right[Floris Vanderhaeghe, Dewey Dunnington, Jan Caha, Jannes Münchow, Jakub Nowosad, Robin Lovelace<br/>2023-06-30]]]

---
class: middle

# .center[My involvement]

.center[
Research Institute for Nature and Forest (INBO), Brussels, Belgium

INBO & Open Science

Helping other scientists to use R in ecology

Contributing to open source applications

Current maintainer of **qgisprocess**

<br/>

[@florisvdh <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)
]

.footnote[.left[<img src="images/flanders.png" width="30%"/>]]

.pull-right[.footnote[.right[<img src="images/inbo.jpg" width="55%"/>]]]

---

# About R

- Free and open-source, cross-platform

- Easy-to-learn interactive programming language

- Data manipulation, calculation and visualization

- Numerical computing & statistics built-in

- Strong ties to scientific research

- Vast package ecosystem: currently 19773 packages on CRAN

- Packages to represent, analyze or visualize geospatial data

- some popular packages: **sf**, **terra**, **mapview**, **tmap**

- many more are availabe: <https://cran.r-project.org/view=Spatial>
  
---
class: inverse, middle, center

# Why interface QGIS processing from R?

---

# 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

<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.32.0-Lima 'Lima' (3.32.0-Lima)
Usage: /usr/bin/qgis_process.bin [--help] [--version] [--json] [--verbose] [--no-python] 
[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/
  
- **qgis** (Jan Caha): functions for each algorithm; uses **qgisprocess** package
  - https://jancaha.github.io/r_package_qgis/

---
class: inverse, center, middle

# A simple example

---
layout: true

# Basic usage

---

---

Load the R package:

```r
library(qgisprocess)
```

```
## Attempting to load the cache ... Success!
## QGIS version: 3.32.0-Lima
## Having access to 1868 algorithms from 13 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/RtmpyYt8uj/file660c36d0a63a/file660c14a85bc2.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: 13 × 3
##    provider         provider_title                         algorithm_count
##    <chr>            <chr>                                            <int>
##  1 cartographytools Cartography tools                                    5
##  2 gdal             GDAL                                                56
##  3 grass7           GRASS                                              306
##  4 latlontools      Lat Lon tools                                        8
##  5 otb              OTB                                                109
##  6 qgis             QGIS                                                50
##  7 3d               QGIS (3D)                                            1
##  8 native           QGIS (native c++)                                  243
##  9 pdal             QGIS (PDAL)                                         17
## 10 qneat3           QNEAT3 - Qgis Network Analysis Toolbox              14
## 11 sagang           SAGA Next Gen                                      509
## 12 visibility       Visibility analysis                                  4
## 13 wbt              WhiteboxTools                                      546
```

---

```r
qgis_plugins()
```

```
## # A tibble: 9 × 2
##   name                    enabled
##   <chr>                   <lgl>  
## 1 QNEAT3                  TRUE   
## 2 ViewshedAnalysis        TRUE   
## 3 cartography_tools       TRUE   
## 4 grassprovider           TRUE   
## 5 latlontools             TRUE   
## 6 otbprovider             TRUE   
## 7 processing              TRUE   
## 8 processing_saga_nextgen TRUE   
## 9 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 8
##  $ MAX             : num 547
##  $ MEAN            : num 348
##  $ MIN             : num 141
##  $ OUTPUT_HTML_FILE: 'qgis_outputHtml' chr "/tmp/RtmpyYt8uj/file660c36d0a63a/file660c1de1cfc4"
##  $ 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("grass7: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/RtmpyYt8uj/file660c36d0a63a/file660c23779d27.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

# Further processing in R

```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
```

---

```r
library(ggplot2)
ggplot(data = zonal_stats, 
       aes(x = distance, y = dem_median, ymin = dem_min, ymax = dem_max)) +
  geom_line() +
  geom_ribbon(fill = "#93b023", alpha = 0.3) +
  facet_wrap(~name, scales = "free_x") +
  labs(y = "elevation [m]", title = plottitle)
```

---
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 !

---

# Future developments

- Include more tutorials
  - options & environment variables
  - accepted R object types
  - copy JSON from QGIS
  
- Expand function documentation

- Adapt to `qgis_process` development

---
class: inverse, center, middle

# Contributors

---
class: inverse, center, middle

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

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