9.3 (R)SAGA

The System for Automated Geoscientific Analyses (SAGA; Table 9.1) provides the possibility to execute SAGA modules via the command line interface (saga_cmd.exe under Windows and just saga_cmd under Linux) (see the SAGA wiki on modules).In addition, there is a Python interface (SAGA Python API).RSAGA uses the former to run SAGA from within R.

Though SAGA is a hybrid GIS, its main focus has been on raster processing, and here particularly on digital elevation models (soil properties, terrain attributes, climate parameters).Hence, SAGA is especially good at the fast processing of large (high-resolution) raster datasets (Conrad et al. 2015).Therefore, we will introduce RSAGA with a raster use case from Muenchow, Brenning, and Richter (2012).Specifically, we would like to compute the SAGA wetness index from a digital elevation model.First of all, we need to make sure that RSAGA will find SAGA on the computer when called.For this, all RSAGA functions using SAGA in the background make use of rsaga.env().Usually, rsaga.env() will detect SAGA automatically by searching several likely directories (see its help for more information).

  1. library(RSAGA)
  2. rsaga.env()

However, it is possible to have ‘hidden’ SAGA in a location rsaga.env() does not search automatically.linkSAGA searches your computer for a valid SAGA installation.If it finds one, it adds the newest version to the PATH environment variable thereby making sure that rsaga.env() runs successfully.It is only necessary to run the next code chunk if rsaga.env() was unsuccessful (see previous code chunk).

  1. library(link2GI)
  2. saga = linkSAGA()
  3. rsaga.env()

Secondly, we need to write the digital elevation model to a SAGA-format.Note that calling data(landslides) attaches two objects to the global environment - dem, a digital elevation model in the form of a list, and landslides, a data.frame containing observations representing the presence or absence of a landslide:

  1. data(landslides)
  2. write.sgrd(data = dem, file = file.path(tempdir(), "dem"), header = dem$header)

The organization of SAGA is modular.Libraries contain so-called modules, i.e., geoalgorithms.To find out which libraries are available, run:

  1. rsaga.get.libraries()

We choose the library ta_hydrology (ta is the abbreviation for terrain analysis).Subsequently, we can access the available modules of a specific library (here: ta_hydrology) as follows:

  1. rsaga.get.modules(libs = "ta_hydrology")

rsaga.get.usage() prints the function parameters of a specific geoalgorithm, e.g., the SAGA Wetness Index, to the console.

  1. rsaga.get.usage(lib = "ta_hydrology", module = "SAGA Wetness Index")

Finally, you can run SAGA from within R using RSAGA’s geoprocessing workhorse function rsaga.geoprocessor().The function expects a parameter-argument list in which you have specified all necessary parameters.

  1. params = list(DEM = file.path(tempdir(), "dem.sgrd"),
  2.               TWI = file.path(tempdir(), "twi.sdat"))
  3. rsaga.geoprocessor(lib = "ta_hydrology", module = "SAGA Wetness Index", 
  4.                    param = params)

To facilitate the access to the SAGA interface, RSAGA frequently provides user-friendly wrapper-functions with meaningful default values (see RSAGA documentation for examples, e.g., ?rsaga.wetness.index).So the function call for calculating the ‘SAGA Wetness Index’ becomes as simple as:

  1. rsaga.wetness.index(in.dem = file.path(tempdir(), "dem"), 
  2.                     out.wetness.index = file.path(tempdir(), "twi"))

Of course, we would like to inspect our result visually (Figure 9.2).To load and plot the SAGA output file, we use the raster package.

  1. library(raster)
  2. twi = raster::raster(file.path(tempdir(), "twi.sdat"))
  3. # shown is a version using tmap
  4. plot(twi, col = RColorBrewer::brewer.pal(n = 9, name = "Blues"))

SAGA wetness index of Mount Mongón, Peru.
Figure 9.2: SAGA wetness index of Mount Mongón, Peru.

You can find an extended version of this example in vignette("RSAGA-landslides") which includes the use of statistical geocomputing to derive terrain attributes as predictors for a non-linear Generalized Additive Model (GAM) to predict spatially landslide susceptibility (Muenchow, Brenning, and Richter 2012).The term statistical geocomputation emphasizes the strength of combining R’s data science power with the geoprocessing power of a GIS which is at the very heart of building a bridge from R to GIS.