9.2 (R)QGIS

QGIS is one of the most popular open-source GIS (Table 9.1; Graser and Olaya 2015).Its main advantage lies in the fact that it provides a unified interface to several other open-source GIS.This means that you have access to GDAL, GRASS and SAGA through QGIS (Graser and Olaya 2015).To run all these geoalgorithms (frequently more than 1000 depending on your set-up) outside of the QGIS GUI, QGIS provides a Python API.RQGIS establishes a tunnel to this Python API through the reticulate package.Basically, functions set_env() and open_app() are doing this.Note that it is optional to run set_env() and open_app() since all functions depending on their output will run them automatically if needed.Before running RQGIS, make sure you have installed QGIS and all its (third-party) dependencies such as SAGA and GRASS.To install RQGIS a number of dependencies are required, as described in the install_guide vignette, which covers installation on Windows, Linux and Mac.At the time of writing (autumn 2018) RQGIS only supports the Long Term Release (2.18), but support for QGIS 3 is in the pipeline (see RQGIS3).

  1. library(RQGIS)
  2. set_env(dev = FALSE)
  3. #> $`root`
  4. #> [1] "C:/OSGeo4W64"
  5. #> ...

Leaving the path-argument of set_env() unspecified will search the computer for a QGIS installation.Hence, it is faster to specify explicitly the path to your QGIS installation.Subsequently, open_app() sets all paths necessary to run QGIS from within R, and finally creates a so-called QGIS custom application (see http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/intro.html#using-pyqgis-in-custom-applications).

  1. open_app()

We are now ready for some QGIS geoprocessing from within R!The example below shows how to unite polygons, a process that unfortunately produces so-called slivers, tiny polygons resulting from overlaps between the inputs that frequently occur in real-world data.We will see how to remove them.

For the union, we use again the incongruent polygons we have already encountered in Section 4.2.5.Both polygon datasets are available in the spData package, and for both we would like to use a geographic CRS (see also Chapter 6).

  1. data("incongruent", "aggregating_zones", package = "spData")
  2. incongr_wgs = st_transform(incongruent, 4326)
  3. aggzone_wgs = st_transform(aggregating_zones, 4326)

To find an algorithm to do this work, find_algorithms() searches all QGIS geoalgorithms using regular expressions.Assuming that the short description of the function contains the word “union”, we can run:

  1. find_algorithms("union", name_only = TRUE)
  2. #> [1] "qgis:union"        "saga:fuzzyunionor" "saga:union"

Short descriptions for each geoalgorithm can also be provided, by setting name_only = FALSE.If one has no clue at all what the name of a geoalgorithm might be, one can leave the search_term-argument empty which will return a list of all available QGIS geoalgorithms.You can also find the algorithms in the QGIS online documentation.

The next step is to find out how qgis:union can be used.open_help() opens the online help of the geoalgorithm in question.get_usage() returns all function parameters and default values.

  1. alg = "qgis:union"
  2. open_help(alg)
  3. get_usage(alg)
  4. #>ALGORITHM: Union
  5. #>  INPUT <ParameterVector>
  6. #>  INPUT2 <ParameterVector>
  7. #>  OUTPUT <OutputVector>

Finally, we can let QGIS do the work.Note that the workhorse function run_qgis() accepts R named arguments, i.e., you can specify the parameter names as returned by get_usage() in run_qgis() as you would do in any other regular R function.Note also that run_qgis() accepts spatial objects residing in R’s global environment as input (here: aggzone_wgs and incongr_wgs).But of course, you could also specify paths to spatial vector files stored on disk.Setting the load_output to TRUE automatically loads the QGIS output as an sf-object into R.

  1. union = run_qgis(alg, INPUT = incongr_wgs, INPUT2 = aggzone_wgs, 
  2.                  OUTPUT = file.path(tempdir(), "union.shp"),
  3.                  load_output = TRUE)
  4. #> $`OUTPUT`
  5. #> [1] "C:/Users/geocompr/AppData/Local/Temp/RtmpcJlnUx/union.shp"

Note that the QGIS union operation merges the two input layers into one layer by using the intersection and the symmetrical difference of the two input layers (which by the way is also the default when doing a union operation in GRASS and SAGA).This is not the same as st_union(incongr_wgs, aggzone_wgs) (see Exercises)!The QGIS output contains empty geometries and multipart polygons.Empty geometries might lead to problems in subsequent geoprocessing tasks which is why they will be deleted.st_dimension() returns NA if a geometry is empty, and can therefore be used as a filter.

  1. # remove empty geometries
  2. union = union[!is.na(st_dimension(union)), ]

Next we convert multipart polygons into single-part polygons (also known as explode geometries or casting).This is necessary for the deletion of sliver polygons later on.

  1. # multipart polygons to single polygons
  2. single = st_cast(union, "POLYGON")

One way to identify slivers is to find polygons with comparatively very small areas, here, e.g., 25000 m2 (see blue colored polygons in the left panel of Figure 9.1).

  1. # find polygons which are smaller than 25000 m^2
  2. x = 25000
  3. units(x) = "m^2"
  4. single$area = st_area(single)
  5. sub = dplyr::filter(single, area < x)

The next step is to find a function that makes the slivers disappear.Assuming the function or its short description contains the word “sliver”, we can run:

  1. find_algorithms("sliver", name_only = TRUE)
  2. #> [1] "qgis:eliminatesliverpolygons"

This returns only one geoalgorithm whose parameters can be accessed with the help of get_usage() again.

  1. alg = "qgis:eliminatesliverpolygons"
  2. get_usage(alg)
  3. #>ALGORITHM: Eliminate sliver polygons
  4. #>  INPUT <ParameterVector>
  5. #>  KEEPSELECTION <ParameterBoolean>
  6. #>  ATTRIBUTE <parameters from INPUT>
  7. #>  COMPARISON <ParameterSelection>
  8. #>  COMPARISONVALUE <ParameterString>
  9. #>  MODE <ParameterSelection>
  10. #>  OUTPUT <OutputVector>
  11. #>  ...

Conveniently, the user does not need to specify each single parameter.In case a parameter is left unspecified, run_qgis() will automatically use the corresponding default value as an argument if available.To find out about the default values, run get_args_man().

To remove the slivers, we specify that all polygons with an area less or equal to 25000 m2 should be joined to the neighboring polygon with the largest area (see right panel of Figure 9.1).

  1. clean = run_qgis("qgis:eliminatesliverpolygons",
  2.                  INPUT = single,
  3.                  ATTRIBUTE = "area",
  4.                  COMPARISON = "<=",
  5.                  COMPARISONVALUE = 25000,
  6.                  OUTPUT = file.path(tempdir(), "clean.shp"),
  7.                  load_output = TRUE)
  8. #> $`OUTPUT`
  9. #> [1] "C:/Users/geocompr/AppData/Local/Temp/RtmpcJlnUx/clean.shp"

Sliver polygons colored in blue (left panel). Cleaned polygons (right panel).
Figure 9.1: Sliver polygons colored in blue (left panel). Cleaned polygons (right panel).

Other notes

  • Leaving the output parameter(s) unspecified saves the resulting QGIS output to a temporary folder created by QGIS.run_qgis() prints these paths to the console after successfully running the QGIS engine.
  • If the output consists of multiple files and you have set load_output to TRUE, run_qgis() will return a list with each element corresponding to one output file.
    To learn more about RQGIS, please refer to Muenchow, Schratz, and Brenning (2017).