r-spatial / qgisprocess

R package to use QGIS processing algorithms
https://r-spatial.github.io/qgisprocess/
GNU General Public License v3.0
202 stars 20 forks source link

Run grass algorithm that is 100+ times faster than base R solution #26

Closed Robinlovelace closed 4 years ago

Robinlovelace commented 4 years ago

Apologies for raising two issues in one 'issue' but think they may be related.

(OK, the first issue is fixed with the following updated algorithm - I think this is a great use case):

output = qgis_run_algorithm(
  algorithm = "grass7:v.split",
  input = rnet_10,
  length = 500
)

That works but I get this rather verbose error message, despite the fact that the code runs:

Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 334, in _startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "/home/robin/.local/share/QGIS/QGIS3/profiles/default/python/plugins/shapetools/__init__.py", line 8, in classFactory
    return ShapeTools(iface)
  File "/home/robin/.local/share/QGIS/QGIS3/profiles/default/python/plugins/shapetools/shapeTools.py", line 23, in __init__
    self.canvas = iface.mapCanvas()
AttributeError: 'NoneType' object has no attribute 'mapCanvas'

Full reprex below...

Reproducible example that I hope is of use:

# Aim: test break-up of linestrings
# https://github.com/saferactive/saferactive/issues/54
# https://github.com/paleolimbot/qgisprocess

remotes::install_github("itsleeds/pct")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'pct' from a github remote, the SHA1 (4f909cc2) has not changed since last install.
#>   Use `force = TRUE` to force installation
remotes::install_github("paleolimbot/qgisprocess")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'qgisprocess' from a github remote, the SHA1 (7393a470) has not changed since last install.
#>   Use `force = TRUE` to force installation

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(qgisprocess)
#> Using 'qgis_process' at 'qgis_process'.
#> Run `qgis_configure()` for details.

rnet = pct::get_pct_rnet("isle-of-wight")
nrow(rnet)
#> [1] 1311
rnet$length = as.numeric(sf::st_length(rnet))
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
summary(rnet$length)
#>      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#>     1.318    64.825   160.337   424.649   414.868 12143.957

rnet = sf::st_transform(rnet, 27700)
rnet_10 = rnet %>% top_n(n = 10, wt = length)

qgis_configure()
#> getOption('qgisprocess.path') was not found.
#> Sys.getenv('R_QGISPROCESS_PATH') was not found.
#> Trying 'qgis_process' on PATH
#> Success!
algorithms = qgis_algorithms()
algorithms %>% filter(grepl(pattern = "split", x = algorithm, ignore.case = TRUE))
#> # A tibble: 10 x 5
#>    provider provider_title  algorithm       algorithm_id    algorithm_title     
#>    <chr>    <chr>           <chr>           <chr>           <chr>               
#>  1 grass7   GRASS           grass7:v.split  v.split         v.split             
#>  2 native   QGIS (native c… native:antimer… antimeridiansp… Geodesic line split…
#>  3 native   QGIS (native c… native:splitfe… splitfeaturesb… Split features by c…
#>  4 native   QGIS (native c… native:splitli… splitlinesbyle… Split lines by maxi…
#>  5 native   QGIS (native c… native:splitve… splitvectorlay… Split vector layer  
#>  6 native   QGIS (native c… native:splitwi… splitwithlines  Split with lines    
#>  7 saga     SAGA            saga:splitline… splitlinesatpo… Split lines at poin…
#>  8 saga     SAGA            saga:splitline… splitlineswith… Split lines with li…
#>  9 saga     SAGA            saga:splitrgbb… splitrgbbands   Split RGB bands     
#> 10 saga     SAGA            saga:splitshap… splitshapeslay… Split shapes layer …
qgis_show_help("native:splitlinesbylength")
#> Split lines by maximum length (native:splitlinesbylength)
#> 
#> ----------------
#> Description
#> ----------------
#> Splits lines into parts which are no longer than a specified length.
#> This algorithm takes a line (or curve) layer and splits each feature into multiple parts, where each part is of a specified maximum length.
#> 
#> Z and M values at the start and end of the new line substrings are linearly interpolated from existing values.
#> 
#> ----------------
#> Arguments
#> ----------------
#> 
#> INPUT: Input layer
#>  Argument type:  source
#>  Acceptable values:
#>      - Path to a vector layer
#> LENGTH: Maximum line length
#>  Argument type:  distance
#>  Acceptable values:
#>      - A numeric value
#> OUTPUT: Split
#>  Argument type:  sink
#>  Acceptable values:
#>      - Path for new vector layer
#> 
#> ----------------
#> Outputs
#> ----------------
#> 
#> OUTPUT: <outputVector>
#>  Split
output = qgis_run_algorithm(
  algorithm = "native:splitlinesbylength",
  INPUT = rnet_10,
  LENGTH = 500
  )
#> Using `OUTPUT = qgis_tmp_vector()`
#> Running qgis_process run 'native:splitlinesbylength' \
#>   '--INPUT=/tmp/RtmpRjhvRS/file1b1d5e9bbb68/file1b1d213bdac2.gpkg' \
#>   '--LENGTH=500' \
#>   '--OUTPUT=/tmp/RtmpRjhvRS/file1b1d5e9bbb68/file1b1d398b83f7.gpkg'
#> Traceback (most recent call last):
#>   File "/usr/lib/python3/dist-packages/qgis/utils.py", line 334, in _startPlugin
#>     plugins[packageName] = package.classFactory(iface)
#>   File "/home/robin/.local/share/QGIS/QGIS3/profiles/default/python/plugins/shapetools/__init__.py", line 8, in classFactory
#>     return ShapeTools(iface)
#>   File "/home/robin/.local/share/QGIS/QGIS3/profiles/default/python/plugins/shapetools/shapeTools.py", line 23, in __init__
#>     self.canvas = iface.mapCanvas()
#> AttributeError: 'NoneType' object has no attribute 'mapCanvas'
#> error starting plugin: shapetools
#> Traceback (most recent call last):
#>   File "/usr/lib/python3/dist-packages/qgis/utils.py", line 334, in _startPlugin
#>     plugins[packageName] = package.classFactory(iface)
#>   File "/home/robin/.local/share/QGIS/QGIS3/profiles/default/python/plugins/QuickOSM/__init__.py", line 12, in classFactory
#>     return QuickOSMPlugin(iface)
#>   File "/home/robin/.local/share/QGIS/QGIS3/profiles/default/python/plugins/QuickOSM/quick_osm.py", line 57, in __init__
#>     self.toolbar = self.iface.addToolBar('QuickOSM')
#> AttributeError: 'NoneType' object has no attribute 'addToolBar'
#> 
#> error starting plugin: QuickOSM
#> Inputs
#> ----------------
#> 
#> INPUT:   /tmp/RtmpRjhvRS/file1b1d5e9bbb68/file1b1d213bdac2.gpkg
#> LENGTH:  500
#> OUTPUT:  /tmp/RtmpRjhvRS/file1b1d5e9bbb68/file1b1d398b83f7.gpkg
#> 
#> 
#> 0...10
#> ERROR 1: failed to execute insert : UNIQUE constraint failed: file1b1d398b83f7.fid
#> ...
#> ERROR 1: failed to execute insert : UNIQUE constraint failed: file1b1d398b83f7.fid
#> 
#> ----------------
#> Results
#> ----------------
#> 
#> OUTPUT:  /tmp/RtmpRjhvRS/file1b1d5e9bbb68/file1b1d398b83f7.gpkg
output
#> <Result of `qgis_run_algorithm("native:splitlinesbylength", ...)`>
#> List of 1
#>  $ OUTPUT: 'qgis_outputVector' chr "/tmp/RtmpRjhvRS/file1b1d5e9bbb68/file1b1d398b83f7.gpkg"
output$OUTPUT
#> [1] "/tmp/RtmpRjhvRS/file1b1d5e9bbb68/file1b1d398b83f7.gpkg"
#> attr(,"class")
#> [1] "qgis_outputVector"
rnet_10_split = sf::st_read(output[[1]][1])
#> Reading layer `file1b1d398b83f7' from data source `/tmp/RtmpRjhvRS/file1b1d5e9bbb68/file1b1d398b83f7.gpkg' using driver `GPKG'
#> Simple feature collection with 10 features and 8 fields
#> geometry type:  LINESTRING
#> dimension:      XY
#> bbox:           xmin: 433285.2 ymin: 82229.69 xmax: 461638.6 ymax: 97672.46
#> projected CRS:  OSGB 1936 / British National Grid

sf::st_length(rnet_10)
#> Units: [m]
#>  [1]  9751.946  5227.855 11632.518 12139.753  3948.035  5810.891  5464.340
#>  [8]  8982.441  3916.680  3972.333
sf::st_length(rnet_10_split)
#> Units: [m]
#>  [1] 500 500 500 500 500 500 500 500 500 500
plot(rnet_10$geometry, lwd = 9, col = "grey")
plot(rnet_10_split$geometry, add = TRUE)

image

Created on 2020-10-27 by the reprex package (v0.3.0)

Robinlovelace commented 4 years ago

Final thing to say on this: It did that for 500k lines in ~5 min 🤯 (the R implementation was going to take around 10 hr). So your amazing work on this project is already yielding real world benefits. Many thanks @paleolimbot (and @jannes-m for first demonstrating access to GRASS/SAGA/other algorithms via QGIS from R)!

paleolimbot commented 4 years ago

Sorry I've been MIA on this! I too got unrelated errors for a bit because Python or some plugins weren't loading, which I think was fixed in QGIS. It looks like that plugin is making the assumption that there is a mapCanvas, which in headless mode, there won't be. I don't know if there's a way to address that here...I think that QGIS output will always be noisy and should stay that way (or be totally silenced using silent = TRUE). I suppose there could be another level of verbosity that suppresses the stderr output but leaves the stdout bit?

Robinlovelace commented 4 years ago

I suppose there could be another level of verbosity that suppresses the stderr output but leaves the stdout bit?

Yes maybe. It's a minor issue and was wondering if it's more a local install issue or even an issue for the QGIS issue tracker. Clearly not a priority and a minor issue, just reporting out of interest more than anything, can leave open or close as a "won't fix" but useful for others to know about this 'feature' that may affect different QGIS installations differently, seemingly depending on the plugins installed. Related question: can qgis_run_algorithm() be used to run plugins that are subsequently installed by the user? I imagine so, could be one to document if anyone gets it working - the sDNA plugin would be my use case but AFAIK the plugin only works on QGIS 2 still (heads-up collaborator and developer of sDNA spatial network analysis software @fiftysevendegreesofrad FYI this may be a way for more people to easily access sDNA from a command line environment and integrate with R).

paleolimbot commented 4 years ago

I think that this particular issue is a plugin-specific (shapetools plugin) issue rather than a QGIS issue if you were to file one somewhere.

I'm pretty sure that starting with QGIS 3.16, the list of plugins loaded when you open the QGIS GUI will be the same as when you run qgis_process. It probably worked before too, but I think there's a PR that made it more reliable in headless mode recently.

I'll close the issue but as always open a new one if it shows up again and there's something that we can change here to make it better!

Robinlovelace commented 4 years ago

Reasonable solution. On the question of running functionality in plugins with qgisprocess, any ideas? Would be good to have a concrete use case which I currently lack.

paleolimbot commented 4 years ago

I think that all providers except "native" and "3d" are already Python plugins whose algorithms can be run using qgis_process. You could find another plugin that adds algorithms and see if they show up in qgis_algorithms()?

Robinlovelace commented 4 years ago

You could find another plugin that adds algorithms and see if they show up in qgis_algorithms()?

Good plan, will aim to document what I find.