coordinate system - Strange output from RSAGA topographic wetness index?

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)

Answer

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...

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 :)

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.