Testing GSIF and RQGIS packages
Can we access QGIS spatial algorithms from R?
That’s exactly what the RQGIS package is promising. RQGIS is available from http://jannes-m.github.io/RQGIS/index.html
According to Jannes Muenchow and Patrick Schratz, RQGIS’s authors, its main advantages are:
It provides access to QGIS functionalities. Thereby, it calls Python from the command line (QGIS API) but R users can stay in their programming environment of choice without having to touch Python.
It offers a broad suite of geoalgorithms making it possible to solve virtually any GIS problem.
R users can just use one package (RQGIS) instead of using RSAGA and spgrass to access SAGA and GRASS functions.
This document explores several functionalities of RQGIS. As we are interested in digital soil mapping, we will also test several libraries used for such purposes, specifically the GISF (Global Soil Information Facilities) package.
# GISF software installation instructions
# By: T. Hengl
# http://gsif.isric.org
# These examples based on the MRO version of R (Microsoft R Open 3.2.4)
## Check that all packages have been installed:
#list.of.packages = c("GSIF", "plotKML", "nnet", "plyr", "ROCR", "randomForest", "plyr", "parallel", "psych", "mda", "h2o", "dismo", "grDevices", "snowfall", "hexbin", "lattice", "ranger", "xgboost", "doParallel", "caret", "RQGIS")
#new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
#if(length(new.packages)) install.packages(new.packages)## (optional) install from R-Forge:
#install.packages("GSIF", repos=c("http://R-Forge.R-project.org"), type="source", dependencies=TRUE)
#install.packages("devtools")## (optional) install from github:
#library(devtools)
#devtools::install_github("jannes-m/RQGIS")
## http://stackoverflow.com/questions/8229859/sourcing-an-r-script-from-github-for-global-session-use-from-within-a-wrapper
#source_https <- function(url, ...) {
# load package
# require(RCurl)
# download:
# cat(getURL(url, followlocation = TRUE, cainfo = system.file("CurlSSL", "cacert.pem", package = "RCurl")), file = basename(url))
# source(basename(url))
#}
#source_https("https://raw.githubusercontent.com/cran/GSIF/master/R/OCSKGM.R")
#OCSKGM#install.packages("quantregForest")Let’s test the GSIF package:
library(GSIF)## GSIF version 0.5-3 (2016-07-18)## URL: http://gsif.r-forge.r-project.org/library(sp)
library(boot)
library(aqp)## This is aqp 1.10library(plyr)
library(rpart)
library(splines)
library(gstat)
library(quantregForest)## Loading required package: randomForest## randomForest 4.6-12## Type rfNews() to see new features/changes/bug fixes.## Loading required package: RColorBrewerlibrary(plotKML)## plotKML version 0.5-6 (2016-05-02)## URL: http://plotkml.r-forge.r-project.org/demo(meuse, echo=FALSE)
omm <- fit.gstatModel(meuse, om~dist+ffreq, meuse.grid, method="quantregForest")## Fitting a Quantile Regression Forest model...## Fitting a 2D variogram...## Saving an object of class 'gstatModel'...om.rk <- predict(omm, meuse.grid)## Subsetting observations to fit the prediction domain in 2D...## Prediction error for 'randomForest' model estimated using the 'quantreg' package.## Generating predictions using the trend model (RK method)...## [using ordinary kriging]
##
100% done## Running 5-fold cross validation using 'krige.cv'...## Creating an object of class "SpatialPredictions"om.rk## Variable : om
## Minium value : 1
## Maximum value : 17
## Size : 153
## Total area : 4964800
## Total area (units) : square-m
## Resolution (x) : 40
## Resolution (y) : 40
## Resolution (units) : m
## Vgm model : Exp
## Nugget (residual) : 2.24
## Sill (residual) : 3.77
## Range (residual) : 1750
## RMSE (validation) : 1.755
## Var explained : 73.7%
## Effective bytes : 1221
## Compression method : gzipplotKML(om.rk)## KML file opened for writing...## Reprojecting to +proj=longlat +datum=WGS84 ...## Warning in paths(show.paths = TRUE): Could not locate SAGA GIS! Install
## program and add it to the Windows registry. See http://www.saga-gis.org/en/
## for more info.## Writing to KML...
## Reprojecting to +proj=longlat +datum=WGS84 ...## Writing to KML...## Closing om.rk.kmlNow, let’s test one of the SAGA GIS modules:
## SAGA GIS (https://sourceforge.net/projects/saga-gis/):
if(.Platform$OS.type == "windows"){
saga_cmd = shortPathName(normalizePath("C:/SAGA-GIS/saga_cmd.exe"))
} else {
saga_cmd = "/usr/local/bin/saga_cmd"
}
system(paste(saga_cmd))## Warning: comando ejecutado 'C:\SAGA-GIS\saga_cmd.exe' tiene estatus 1library(rgdal)## rgdal: version: 1.2-5, (SVN revision 648)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
## Path to GDAL shared files: C:/Rlibs/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: C:/Rlibs/rgdal/proj
## Linking to sp version: 1.2-4library(raster)##
## Attaching package: 'raster'## The following objects are masked from 'package:aqp':
##
## metadata, metadata<-data("eberg_grid")
gridded(eberg_grid) <- ~x+y
proj4string(eberg_grid) <- CRS("+init=epsg:31467")
writeGDAL(eberg_grid["DEMSRT6"], "DEMSRT6.sdat", "SAGA")
system(paste(saga_cmd, 'ta_lighting 0 -ELEVATION "DEMSRT6.sgrd" -SHADE "hillshade.sgrd" -EXAGGERATION 2'))
plot(raster("hillshade.sdat"), col=SAGA_pal[[3]])
Now, the turn is for GDAL:
## GDAL (https://trac.osgeo.org/gdal/wiki/DownloadingGdalBinaries):
if(.Platform$OS.type == "windows"){
gdal.dir <- shortPathName("C:/Program files/GDAL")
gdal_translate <- paste0(gdal.dir, "/gdal_translate.exe")
gdalwarp <- paste0(gdal.dir, "/gdalwarp.exe")
} else {
gdal_translate = "gdal_translate"
gdalwarp = "gdalwarp"
}
system(paste(gdalwarp))## Warning: comando ejecutado 'C:\PROGRA~1\GDAL/gdalwarp.exe' tiene estatus 1system(paste(gdalwarp, ' DEMSRT6.sdat DEMSRT6_ll.tif -t_srs \"+proj=longlat +datum=WGS84\"'))
plot(raster("DEMSRT6_ll.tif"))
Finally, let´s test RQGIS:
## RQGIS (https://www.r-bloggers.com/rqgis-0-1-0-release/)
#install.packages("RQGIS")
library(RQGIS)
library(raster)
#env <- set_env()
env <- set_env("C:/OSGeo4W64")
find_algorithms("grass", name_only=TRUE, qgis_env=env)## [1] "C:\\Users\\IVNLIZ~1\\AppData\\Local\\Temp\\RtmpkPXXXu"
## [2] "grass7:nviz"
## [3] "grass7:r.out.gridatb"
## [4] "grass7:r.resample"
## [5] "grass7:v.in.dxf"
## [6] "grass7:v.in.geonames"
## [7] "grass7:v.in.lidar"
## [8] "grass7:v.in.mapgen"
## [9] "grass7:v.mkgrid"
## [10] "grass7:v.out.dxf"
## [11] "grass7:v.out.pov"
## [12] "grass7:v.pack"
## [13] "grass7:v.perturb"
## [14] "grass:r.out.gridatb"
## [15] "grass:r.resample"
## [16] "grass:v.in.dxf"
## [17] "grass:v.mkgrid"
## [18] "grass:v.out.dxf"
## [19] "grass:v.out.pov"
## [20] "grass:v.perturb"find_algorithms("slope", name_only=TRUE, qgis_env=env)## [1] "C:\\Users\\IVNLIZ~1\\AppData\\Local\\Temp\\RtmpkPXXXu"
## [2] "taudem:dinfinityupslopedependence"
## [3] "taudem:d8extremeupslopevalue"
## [4] "taudem:slopeareacombination"
## [5] "taudem:slopeaveragedown"
## [6] "taudem:slopeoverarearatio"
## [7] "grass7:r.slope"
## [8] "grass7:r.slope.aspect"
## [9] "saga:downslopedistancegradient"
## [10] "saga:dtmfilterslopebased"
## [11] "saga:relativeheightsandslopepositions"
## [12] "saga:slopelength"
## [13] "saga:slopeaspectcurvature"
## [14] "saga:upslopearea"
## [15] "saga:vegetationindexslopebased"
## [16] "gdalogr:slope"
## [17] "grass:r.flow"
## [18] "grass:r.slope"
## [19] "grass:r.slope.aspect"args <- get_args_man(alg="gdalogr:slope", options=TRUE, qgis_env=env)
args## $INPUT
## [1] "None"
##
## $BAND
## [1] "1"
##
## $COMPUTE_EDGES
## [1] "False"
##
## $ZEVENBERGEN
## [1] "False"
##
## $AS_PERCENT
## [1] "False"
##
## $SCALE
## [1] "1.0"
##
## $OUTPUT
## [1] "None"# load data into R
data("dem", package = "RQGIS")
# define input
args$INPUT <- dem
# specify output path
args$OUTPUT <- file.path('c:/rlibs/tmp',"slope.tif")
#args$ASPECT <- file.path('c:/rlibs/tmp',"aspect.sdat")targs <- get_args_man("qgis:creategrid", options = TRUE, qgis_env = env)
targs$TYPE <- 1
targs$EXTENT <- "0, 1, 0, 1"
targs$HSPACING <- "0.5"
targs$VSPACING <- "0.5"
targs$CRS <- "EPSG:4326"
targs$OUTPUT <- file.path('c:/rlibs/tmp',"grid.shp")
run_qgis("qgis:creategrid", params = targs, load_output = targs$OUTPUT, qgis_env = env)## class : SpatialPolygonsDataFrame
## features : 4
## extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 4
## names : left, top, right, bottom
## min values : 0, 0.5, 0.5, 0
## max values : 0.5, 1, 1, 0.5slope <- run_qgis(alg = "gdalogr:slope",
params = args, check_params = TRUE, show_msg = TRUE,
load_output = args$OUTPUT,
qgis_env = env)# visualize the result
hs <- hillShade(terrain(dem), terrain(dem, "aspect"), 40, 270)
pal <- RColorBrewer::brewer.pal(6, "Blues")
spplot(slope, col.regions = pal, alpha.regions = 0.8,
scales = list(tck = c(1, 0)),
colorkey = list(space = "bottom",
width = 1, height = 0.5,
axis.line = list(col = "black")),
cuts = length(pal) - 1) +
latticeExtra::as.layer(spplot(hs, col.regions = gray(0:100 / 100)),
under = TRUE)
find_algorithms("terrain", name_only=TRUE, qgis_env=env)## [1] "C:\\Users\\IVNLIZ~1\\AppData\\Local\\Temp\\RtmpkPXXXu"
## [2] "grass7:r.param.scale"
## [3] "saga:randomterraingeneration"
## [4] "saga:terrainruggednessindextri"
## [5] "gdalogr:triterrainruggednessindex"
## [6] "grass:r.param.scale"args <- get_args_man(alg="gdalogr:triterrainruggednessindex", options=TRUE, qgis_env=env)
args## $INPUT
## [1] "None"
##
## $BAND
## [1] "1"
##
## $COMPUTE_EDGES
## [1] "False"
##
## $OUTPUT
## [1] "None"# define input
args$INPUT <- dem
# specify output path
args$OUTPUT <- file.path('c:/rlibs/tmp',"rugged.tif")
#args$ASPECT <- file.path('c:/rlibs/tmp',"aspect.sdat")rugged <- run_qgis(alg = "gdalogr:triterrainruggednessindex",
params = args, check_params = TRUE, show_msg = TRUE,
load_output = args$OUTPUT,
qgis_env = env)# visualize the result
hs <- hillShade(terrain(dem), terrain(dem, "aspect"), 40, 270)
pal <- RColorBrewer::brewer.pal(10, "Reds")## Warning in RColorBrewer::brewer.pal(10, "Reds"): n too large, allowed maximum for palette Reds is 9
## Returning the palette you asked for with that many colorsspplot(rugged, col.regions = pal, alpha.regions = 0.8,
scales = list(tck = c(1, 0)),
colorkey = list(space = "bottom",
width = 1, height = 0.5,
axis.line = list(col = "black")),
cuts = length(pal) - 1) +
latticeExtra::as.layer(spplot(hs, col.regions = gray(0:100 / 100)),
under = TRUE)
One final note: We weren’t able to use SAGA functions within RQGIS. Several attempts to obtain terrain derivatives weren’t unsuccesful. RQGIS error messages simply stated that the package was not able to produce the expected outputs (i.e. slope and aspect). To be continued…