6

This question already has an answer here:

I am trying to generate a topographic wetness index raster with the RSAGA package in R. The output I am able to generate does not make sense to me, as it has a large number of NA cells and I wonder if I am doing something wrong. I have put together a simple dummy example to illustrate what I mean.

Session info: I am running R 3.2.3 on Ubuntu 15.10. I have SAGA v2.2.0 installed, as well as RSAGA v0.94-5.

library(RSAGA)
library(raster)
library(maptools)

#get DEM for Switzerland
alt <- getData('alt', country='CHE')

# write to SAGA format
tempInfile <- 'saga_temp_twi_in'
writeRaster(alt, filename=tempInfile, format='SAGA')

outfile <- 'saga_temp_twi_out.sgrd'
call <- rsaga.wetness.index(in.dem='saga_temp_twi_in.sgrd', out.wetness.index=outfile, env=rsaga.env(modules='/usr/lib/x86_64-linux-gnu/saga/'))

#read resulting file back in
res <- raster('saga_temp_twi_out.sdat')

par(mfrow=c(1,2))
plot(alt)
plot(res)

enter image description here

marked as duplicate by Andre Silva, whyzar, xunilk, BERA, Kersten Jan 3 '18 at 14:45

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • The duplicate is about QGIS, but the source error and the alternatives to address it are the same. – Andre Silva Jan 3 '18 at 1:26
8
+50

I think I see the problem.

I could reproduce this, both using your script, and SAGA itself. The output from SAGA matches the output from your script. So it's not RSAGA.

But look at the slope output from SAGA... vast areas are near-vertical, with slopes around 1.57 radians (nearly 90 degrees) with some flat areas. You can confirm this looking at the histogram of slope in SAGA...

histogram of slop

You're only seeing water accumulate in the flat areas... where the slope is around 0. The slope is acting like a mask.

I think the problem is the file isn't georeferenced, so SAGA thinks that Switzerland is a few meters across (the size of your raster, which of course is in degrees). It's assuming vertical distances in same units as horizontal.. it's not aware the units are degrees. But Switzerland still a few thousand meters high! The end result is a massive vertical exaggeration.

As a quick sanity check, I ran Fill Sinks, then used Grid calculator to divide elevations by 110000 (approx meters to degrees conversion at latitude 45N) Running SAGA Wetness Index on this, gives a much more sensible result :)

final result of topographic wetness

Not sure how to do this in RSAGA, though. Your best bet is to use a raster with a CRS assigned to it, or assign a CRS in SAGA itself. My division trick isn't going to give results as accurate as you'd get with a properly projected raster.

  • 2
    You're correct! I find that I can't write a projected raster to SAGA format with raster's writeRaster function, but if I use RSAGA's rsaga.import.gdal function to import the projected raster, it worked great! – Pascal Jan 9 '16 at 21:08

Not the answer you're looking for? Browse other questions tagged or ask your own question.