Using QGIS from R through the Network API plugin
Kevin Stadler
Introduction
The qgisremote package allows programmatic interaction with a running (local or remote) instance of QGIS Desktop. Unlike with the Python API, all communication and calls to QGIS are passed through the HTTP interface provided by the Network API plugin which you’ll need to install (and turn on) first.
By connecting to a running instance of QGIS, you can freely combine batch-processing of spatial data in R with the user-friendly interaction, manual filtering and rendering available through QGIS’ graphical interface. The possibility to control QGIS running on other machines might also be useful for educational purposes.
Setup
If you’re interacting with your local QGIS only, no setup beyond loading this R package should be necessary. By default, all requests are passed on to port 8090 on localhost.
library(qgisremote)
# inspect default connection settings
qgisremote.options('host')## [1] "http://localhost:8090/"If you want to set a different host or password for all future calls in the current R session, you can set new values as follows:
# don't run this unless you've actually changed the port in the Network API plugin config dialog
qgisremote.options(host = 'http://localhost:12345')All the functions demonstrated below allow for arbitrary named parameters to be passed (by using R’s ... construct). If you want to override the API target host or password for one request only, you can simply pass them to the function call in question.
Adding map data
First, make sure we are starting with a blank canvas/project:
## NULLAdd a tile layer as a background.
tilelayer <- iface.addTileLayer('http://a.tile.openstreetmap.org/{z}/{x}/{y}.png')
print(tilelayer)## QGIS map layer with the following properties:
## crs: CRS arguments:
## +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0
## +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs
## valid: [1] TRUE
## name: [1] "a.tile.openstreetmap.org tiles"
## extent: [1] -20037508 -20037508 20037508 20037508
## publicSource: [1] "type=xyz&url=http%3A%2F%2Fa.tile.openstreetmap.org%2F%7Bz%7D%2F%7Bx%7D%2F%7By%7D.png"
## type: [1] 1
## id: [1] "a_tile_openstreetmap_org_tiles20170831124138140"
## isEditable: [1] FALSENow add some vector data using iface.addVectorLayer(). The method is overloaded so data can be provided in a number of different ways.
The first option is to simply pass a QGIS provider string. In its simplest possible form this can be a link to a publicly accessible data set somewhere on the web, but it can also be a specification to some locally restricted resource such as a PostGIS server, including login credentials (see the corresponding documentation in the PyQGIS cookbook). The method always returns summary information about the layer that was added, which is a list of class qgislayer.
remotevectorlayer <- iface.addVectorLayer('https://d2ad6b4ur7yvpq.cloudfront.net/naturalearth-3.3.0/ne_50m_populated_places.geojson', baseName = 'providerstring-layer')
print(remotevectorlayer)## QGIS map layer with the following properties:
## crs: CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## valid: [1] TRUE
## name: [1] "providerstring-layer OGRGeoJSON Point"
## extent: [1] -175.22056 -90.00000 179.21665 78.21668
## publicSource: [1] "https://d2ad6b4ur7yvpq.cloudfront.net/naturalearth-3.3.0/ne_50m_populated_places.geojson"
## type: [1] 0
## id: [1] "providerstring_layer20170831124138639"
## isEditable: [1] FALSEThe second option is by passing one of R’s geospatial data objects to the method, in particular instances of class sf and geo_json. The objects are automatically converted to geojson and transferred to QGIS. (Support for transferring files in other formats directly from disk might be added later, but this will require implementing multi-part messages to be able to process formats such as shapefiles which consist of more than one file.)
# you will normally want to read in local files, but st_read also works on remote urls
data <- sf::st_read('https://d2ad6b4ur7yvpq.cloudfront.net/naturalearth-3.3.0/ne_50m_populated_places.geojson')## Reading layer `OGRGeoJSON' from data source `https://d2ad6b4ur7yvpq.cloudfront.net/naturalearth-3.3.0/ne_50m_populated_places.geojson' using driver `GeoJSON'
## Simple feature collection with 1249 features and 92 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -175.2206 ymin: -90 xmax: 179.2166 ymax: 78.21668
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defssfdatavectorlayer <- iface.addVectorLayer(data, baseName = 'sf-layer')## Converting 1249 features to GeoJSON...## Transferring data...print(sfdatavectorlayer)## QGIS map layer with the following properties:
## crs: CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## valid: [1] TRUE
## name: [1] "sf-layer OGRGeoJSON Point"
## extent: [1] -175.22056 -90.00000 179.21665 78.21668
## publicSource: [1] "/tmp/tmpz1mvp6"
## type: [1] 0
## id: [1] "sf_layer20170831124204583"
## isEditable: [1] FALSEReading map and layer information
Get the number of layers currently added.
## [1] 3Note that this is not necessarily the same as the number of layers currently visible on the map canvas:
## [1] 3x <- mapLayers()
# returns a list of qgislayer objects
class(x)## [1] "list"class(x[[1]])## [1] "qgislayer"# print some information about all layers
print(x)## $providerstring_layer20170831124138639
## QGIS map layer with the following properties:
## crs: CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## valid: [1] TRUE
## name: [1] "providerstring-layer OGRGeoJSON Point"
## extent: [1] -175.22056 -90.00000 179.21665 78.21668
## publicSource: [1] "https://d2ad6b4ur7yvpq.cloudfront.net/naturalearth-3.3.0/ne_50m_populated_places.geojson"
## type: [1] 0
## id: [1] "providerstring_layer20170831124138639"
## isEditable: [1] FALSE
##
## $a_tile_openstreetmap_org_tiles20170831124138140
## QGIS map layer with the following properties:
## crs: CRS arguments:
## +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0
## +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs
## valid: [1] TRUE
## name: [1] "a.tile.openstreetmap.org tiles"
## extent: [1] -20037508 -20037508 20037508 20037508
## publicSource: [1] "type=xyz&url=http%3A%2F%2Fa.tile.openstreetmap.org%2F%7Bz%7D%2F%7Bx%7D%2F%7By%7D.png"
## type: [1] 1
## id: [1] "a_tile_openstreetmap_org_tiles20170831124138140"
## isEditable: [1] FALSE
##
## $sf_layer20170831124204583
## QGIS map layer with the following properties:
## crs: CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## valid: [1] TRUE
## name: [1] "sf-layer OGRGeoJSON Point"
## extent: [1] -175.22056 -90.00000 179.21665 78.21668
## publicSource: [1] "/tmp/tmpz1mvp6"
## type: [1] 0
## id: [1] "sf_layer20170831124204583"
## isEditable: [1] FALSERetrieve attribute data from a layer
Access to layer-specific data and operations normally requires QGIS’ unique layer id. However, the R implementation of the methods support simply passing previously returned qgislayer objects as well:
features <- mapLayer.getFeatures(sfdatavectorlayer)
# identical to mapLayer.getFeatures(sfdatavectorlayer$id)
class(features)## [1] "data.frame"nrow(features)## [1] 1249colnames(features)## [1] "POP1995" "UN_FID" "ADM1NAME" "ISO_A2" "MIN_AREAMI"
## [6] "GTOPO30" "NAMEALT" "SOV0NAME" "MAX_POP310" "DIFFNOTE"
## [11] "POP1985" "POP1980" "MAX_BBXMAX" "LONGITUDE" "LS_MATCH"
## [16] "POP1990" "MAX_POP50" "MAX_BBXMIN" "MAX_AREAKM" "UN_LONG"
## [21] "MAX_POP300" "GEONAMEID" "CHANGED" "RANK_MAX" "ADMIN1_COD"
## [26] "MAX_PERMI" "MIN_BBYMAX" "SCALERANK" "MAX_POP20" "MEGANAME"
## [31] "POP2010" "POP2015" "MAX_NATSCA" "UN_ADM0" "ELEVATION"
## [36] "POP_MIN" "MIN_PERKM" "WORLDCITY" "RANK_MIN" "POP_MAX"
## [41] "CAPIN" "POP2000" "GN_ASCII" "POP2005" "CAPALT"
## [46] "MEAN_BBXC" "MAX_PERKM" "ADM0CAP" "TIMEZONE" "FEATURECLA"
## [51] "POP1950" "NAMEPAR" "POP1955" "LABELRANK" "POP_OTHER"
## [56] "UN_LAT" "GEONAMESNO" "FEATURE_CO" "FEATURE_CL" "DIFFASCII"
## [61] "MIN_BBYMIN" "NAMEDIFF" "NATSCALE" "MIN_PERMI" "MEGACITY"
## [66] "NAME" "MAX_POP10" "SOV_A3" "GN_POP" "POP1975"
## [71] "POP1970" "POP2020" "POP2025" "MIN_BBXMAX" "ADM0_A3"
## [76] "MAX_BBYMIN" "MIN_AREAKM" "CITYALT" "MAX_BBYMAX" "CHECKME"
## [81] "ADM0NAME" "POP1960" "POP1965" "LS_NAME" "MEAN_BBYC"
## [86] "NAMEASCII" "POP2050" "COMPARE" "MIN_BBXMIN" "NOTE"
## [91] "LATITUDE" "MAX_AREAMI"Retrieve geodata from a layer
mapLayer.getFeatures() mimics the different feature filter options available in QGIS.
# select area to get features from
rectangle <- c(10, -10, 30, 0)
layerdata <- mapLayer.getFeatures(sfdatavectorlayer, rect = rectangle, geometry = TRUE)
class(layerdata)## [1] "sf" "tbl_df" "tbl" "data.frame"# a lot of attribute columns, so suppress printing of the actual data
print(layerdata, n = 0)## Simple feature collection with 16 features and 92 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 11.88004 ymin: -9.54 xmax: 29.36001 ymax: -1.6333
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs# check out some of the columns of the filtered data set
print(layerdata[, c('NAME', 'POP1950', 'POP2000', 'POP2050')])## Simple feature collection with 16 features and 4 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 11.88004 ymin: -9.54 xmax: 29.36001 ymax: -1.6333
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 16 x 5
## NAME POP1950 POP2000 POP2050 geometry
## <chr> <int> <int> <int> <S3: sfc_POINT>
## 1 Kananga 24 557 1698 <S3: sfc_POINT>
## 2 Brazzaville 83 986 2150 <S3: sfc_POINT>
## 3 Malanje 0 0 0 <S3: sfc_POINT>
## 4 Kasongo 0 0 0 <S3: sfc_POINT>
## 5 Kikwit 0 0 0 <S3: sfc_POINT>
## 6 Mbuji-Mayi 70 932 2851 <S3: sfc_POINT>
## 7 Kinshasa 202 5485 16762 <S3: sfc_POINT>
## 8 Kalemie 0 0 0 <S3: sfc_POINT>
## 9 Bujumbura 0 0 0 <S3: sfc_POINT>
## 10 Goma 0 0 0 <S3: sfc_POINT>
## 11 Matadi 0 0 0 <S3: sfc_POINT>
## 12 Kayes 0 0 0 <S3: sfc_POINT>
## 13 Luanda 138 2591 8236 <S3: sfc_POINT>
## 14 Franceville 0 0 0 <S3: sfc_POINT>
## 15 Pointe-Noire 0 0 0 <S3: sfc_POINT>
## 16 Bandundu 0 0 0 <S3: sfc_POINT>Map canvas operations
Retrieve information about the map canvas.
## $layers
## [1] "sf_layer20170831124204583"
## [2] "providerstring_layer20170831124138639"
## [3] "a_tile_openstreetmap_org_tiles20170831124138140"
##
## $mapUnitsPerPixel
## [1] 75954.45
##
## $scale
## [1] 287071952
##
## $outputSize
## [1] 647 554
##
## $extent
## [1] -21039384 -21039384 21039384 21039384
##
## $visibleExtent
## [1] -24571266 -21039384 24571266 21039384
##
## $fullExtent
## [1] -20037508 -129099728 20037508 20037508
##
## $rotation
## [1] 0
##
## $mapUnits
## [1] 0
##
## $destinationCrs
## $destinationCrs$postgisSrid
## [1] 3857
##
## $destinationCrs$proj4
## [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
##
## $destinationCrs$srsid
## [1] 3857
##
## $destinationCrs$description
## [1] "WGS 84 / Pseudo Mercator"Manipulating the map canvas
The various zoom and centering functions always return the new/updated zoom level or map canvas center, respectively.
# the Amundsen/Scott data point falls outside the typical Web Mercator range,
# so zooming to the full extent will put the map canvas way off center
mapCanvas.zoomToFullExtent()## [1] 1068324384## [1] 0 -54531110# go back to the geographical center on the equator and pick a good zoom level
mapCanvas.setCenter(0, 0)## [1] 0 0## [1] -91440778 -78297049 91440778 78297049mapCanvas.zoomScale(200000000)## [1] 2e+08## [1] -17118542 -14657917 17118542 14657917Reading out the canvas
x <- mapCanvas.saveAsImage()
# the resolution of the map is determined by the QGIS map canvas size
print(x)## class : RasterBrick
## dimensions : 554, 647, 358438, 4 (nrow, ncol, ncell, nlayers)
## resolution : 52916.67, 52916.67 (x, y)
## extent : -17118542, 17118542, -14657917, 14657917 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs
## data source : /tmp/RtmpE7QnyZ/file2f0c7e97305f
## names : file2f0c7e97305f.1, file2f0c7e97305f.2, file2f0c7e97305f.3, file2f0c7e97305f.4
## min values : 0, 0, 0, 0
## max values : 255, 255, 255, 255raster::plotRGB(x)
Manipulating layer rendering
mapLayer.visible(remotevectorlayer)## [1] TRUEmapLayer.visible(remotevectorlayer) <- FALSE
# the map canvas' layer count will now be different
mapCanvas.layerCount()## [1] 2# whereas the number of layers known to QGIS proper is still the same
mapLayers.count()## [1] 3Accessing and manipulating the renderer
renderer <- mapLayer.renderer(sfdatavectorlayer)
# returned as an XML chunk
class(renderer)## [1] "xml_document" "xml_node"# convert to plain character string (only used much later, see below)
rendererstring <- as.character(renderer)The qgisxml() function allows more handy manipulation by scaffolding access to settings stored in QGIS XML chunks:
qgisxml(renderer)## name attrname value
## 1 renderer-v2 forceraster 0
## 2 renderer-v2 symbollevels 0
## 3 renderer-v2 type singleSymbol
## 4 renderer-v2 enableorderby 0
## 5 symbol alpha 1
## 6 symbol clip_to_extent 1
## 7 symbol type marker
## 8 symbol name 0
## 9 layer pass 0
## 10 layer class SimpleMarker
## 11 layer locked 0
## 12 layer angle 0
## 13 layer color 107,9,134,255
## 14 layer horizontal_anchor_point 1
## 15 layer joinstyle bevel
## 16 layer name circle
## 17 layer offset 0,0
## 18 layer offset_map_unit_scale 0,0,0,0,0,0
## 19 layer offset_unit MM
## 20 layer outline_color 0,0,0,255
## 21 layer outline_style solid
## 22 layer outline_width 0
## 23 layer outline_width_map_unit_scale 0,0,0,0,0,0
## 24 layer outline_width_unit MM
## 25 layer scale_method diameter
## 26 layer size 2
## 27 layer size_map_unit_scale 0,0,0,0,0,0
## 28 layer size_unit MM
## 29 layer vertical_anchor_point 1
## 30 sizescale scalemethod diameterThe name/attrname pairs can be used to query and set specific settings:
qgisxml(renderer, 'layer', 'name')## [1] "circle"# the 'simple marker' renderer supports many different marker types (whose names
# can be found in the 'Style' tab of QGIS)
qgisxml(renderer, 'layer', 'name') <- 'star'If the attribute name is unique, the setting can also be accessed by providing this second name alone:
oldcolor <- qgisxml(renderer, 'color')
print(oldcolor)## [1] "107,9,134,255"Types and valid values for the different settings vary wildly. Color specifications for example accept string specifications in r,g,b,a and #rrggbb format, as well as any English color name understood by Qt’s QColor class (see here for a non-exhaustive list).
qgisxml(renderer, 'color') <- 'red'So far we have only manipulated the renderer XML chunk locally. We still need to call mapLayer.setRenderer() (or its equivalent arrow notation assignment to mapLayer.renderer()) to push the changes to QGIS:
mapLayer.renderer(sfdatavectorlayer) <- renderer
# equivalent
## mapLayer.setRenderer(sfdatavectorlayer, renderer)
raster::plotRGB(mapCanvas.saveAsImage())
Note that it is not required to work with the qgisxml() function and xml_node objects to set/update the renderer. XML chunks can be manipulated by any other means and also applied directly from character strings obtained elsewhere.
# replace old colour with new
newrendererstring <- gsub(oldcolor, '#0000FF', rendererstring)
mapLayer.renderer(sfdatavectorlayer) <- newrendererstring
raster::plotRGB(mapCanvas.saveAsImage())