# Altitude Above channel in GRASS 6.4

i'm fairly new to the GIS world and GRASS (come from econometric background), i would like to do some species distribution modeling using R and GRASS.

So i started to read the book of Tomislav Hengl : A practical guide to geostatistical mapping, and specially chapter 8 about species modeling (using the bei dataset in the spatstat R package).

The same analysis is also on the author wiki page : http://spatial-analyst.net/wiki/index.php?title=Species_Distribution_Modelling

Tomislav Hengl, used mainly SAGA-GIS and R, and i would like to do the same analysis using R and GRASS (through spgrass6).

The problem is that i don't know the equivalent of SAGA ta_channel module 3 (altitude above channel) in GRASS, so far i tried r.stream.distance to compute elevation above stream but the result are really different.

This is the part of original code i want to do (problem is wherer there are **):

``````library(gstat)
library(spatstat)
library(splancs)
library(rgeos)
library(maptools)
library(adehabitat)
library(rgdal)
library(lattice)
library(RSAGA)
# General settings for plots:
trellis.par.get("fontsize") -> fontsize
fontsize\$default<-34
fontsize\$points<-22
fontsize\$text<-34
# environnement pour SAGA
monenv <- rsaga.env(path = "/usr/bin/", modules = "/usr/lib/saga")
#-------------------------------
# 1. Data import and preparation
#-------------------------------

# Load the data:
data(bei)
str(bei)
str(bei.extra, max.level=1)
plot(bei.extra\$elev, main="Beilschmiedia (Tropical rain forest trees)")
plot(bei, add=TRUE, pch=16, cex=0.3)
class(bei.extra[[1]])
# Export to a GIS (for visual exploration):
grids <- as(bei.extra[[1]], "SpatialGridDataFrame")
names(grids)[1] <- "elev"
grids\$grad <- as(bei.extra[[2]], "SpatialGridDataFrame")\$v
writeGDAL(grids[1], "dem.mpr", "ILWIS")
write.asciigrid(grids["elev"], "dem.asc")
write.asciigrid(grids["grad"], "grad.asc")
bei.pnt <- data.frame(x=bei\$x, y=bei\$y, no=rep(1, length(bei\$x)))
coordinates(bei.pnt) <- ~x+y
writeOGR(bei.pnt, "bei.shp", "bei", "ESRI Shapefile")
# reformat the point pattern so it fits the grids:
bei.ppp <- ppp(coordinates(bei.pnt)[,1], coordinates(bei.pnt)[,2], marks=bei.pnt\$no, window=as(grids[1], "owin"))
# Generate the wetness index and vertical distance from the channel network in SAGA and import back to R:
rsaga.esri.to.sgrd(in.grids="dem.asc", out.sgrds="dem.sgrd", in.path=getwd())
rsaga.geoprocessor(lib="ta_hydrology", module=15, param=list(DEM="dem.sgrd", C="catharea.sgrd", GN="catchslope.sgrd", CS="modcatharea.sgrd", SB="twi.sgrd", T=10))
rsaga.geoprocessor(lib="ta_channels", module=0, param=list(ELEVATION="dem.sgrd", CHNLNTWRK="chnlntwrk.sgrd", CHNLROUTE="channel_route.sgrd", SHAPES="channels.shp", INIT_GRID="dem.sgrd", DIV_CELLS=5, MINLEN=10))
**rsaga.geoprocessor(lib="ta_channels", module=3, param=list(ELEVATION="dem.sgrd", CHANNELS="chnlntwrk.sgrd", ALTITUDE="achan.sgrd", THRESHOLD=0.1, NOUNDERGROUND=TRUE))
rsaga.sgrd.to.esri(in.sgrds=c("twi.sgrd","achan.sgrd"), out.grids=c("twi.asc","achan.asc"), prec=1, out.path=getwd())**
grids\$achan <- readGDAL("achan.asc")\$band1
grids\$twi <- readGDAL("twi.asc")\$band1
``````

Here is my code (the problem is where i put **):

``````library(gstat)
library(spatstat)
library(splancs)
library(rgeos)
library(maptools)
library(adehabitat)
library(rgdal)
library(spgrass6)
library(ggExtra)
# environnement pour GRASS
loc <- initGRASS("/usr/lib/grass64", home=tempdir(), override = TRUE)
#-------------------------------
# 1. Data import and preparation
#-------------------------------
# Load the data:

data(bei)
bei.data <- as.data.frame(bei)
bei.elev.data <- as.data.frame(bei.extra\$elev)
ggplot(bei.elev.data, aes(x, y, z = value)) + geom_tile(aes(fill = value)) + geom_point(aes(x, y, z = NULL), data = bei.data, size = 2, pch = "+") + scale_fill_gradientn(colour = terrain.colors(4)) + opts(title = "Beilschmiedia (Tropical rain forest trees)")
class(bei.extra[[1]])
# Export to a GIS (for visual exploration):
grids <- as(bei.extra[[1]], "SpatialGridDataFrame")
names(grids)[1] <- "elev"
grids\$grad <- as(bei.extra[[2]], "SpatialGridDataFrame")\$v
writeRAST6(grids["elev"], "dem", zcol = "elev",  ignore.stderr = TRUE)
writeRAST6(grids["grad"], "grad", zcol = "grad", ignore.stderr = TRUE)
bei.pnt <- data.frame(x=bei\$x, y=bei\$y, no=rep(1, length(bei\$x)))
coordinates(bei.pnt) <- ~x+y
# reformat the point pattern so it fits the grids:
bei.ppp <- ppp(coordinates(bei.pnt)[,1], coordinates(bei.pnt)[,2], marks=bei.pnt\$no, window=as(grids[1], "owin"))
# Generate the wetness index and vertical distance from the channel network in GRASS and import back to R:
# Work on  GRASS:
execGRASS("g.region", parameters=list(rast="dem"))
execGRASS("g.region", parameters=list(rast="grad"))
execGRASS("r.info", parameters = list(map = "dem"))
execGRASS("r.info", parameters = list(map = "grad"))
gmeta6()
# wetness index
execGRASS("r.topidx", flags = "overwrite", parameters = list(input = "dem", output = "twi"))

# r.stream.extract
execGRASS("g.list", parameters =  list(type = "rast"))
execGRASS("r.stream.extract", flags =  "overwrite", parameters = list(elevation = "dem", stream_rast = "stream", direction = "stream_dir", threshold = as.integer(10)))
#execGRASS("r.watershed", flags = c("m", "overwrite"), parameters = list(elevation = "dem", stream = "stream", threshold = as.integer(10)))
# altitude above channel
**execGRASS("r.stream.distance", flags = "overwrite", parameters = list(stream = "stream", dir = "stream_dir",  dem = "dem", method = "downstream", elevation = "achan"))**
# spline to fill
execGRASS("r.fillnulls", flags = "overwrite" ,parameters = list(input = "achan", output = "achan"))
# exportation
grids\$achan <- readRAST6("achan")@data\$achan
grids\$twi <- readRAST6("twi")@data\$twi
``````

So my question is what is really Saga altitude above channel (if possible in plain english) and how to compute this using GRASS ? I'm made several researchs (for weeks now), and may be it's right place here. Thank you in advance and sorry for my poor english

-
Hmm, from my half-minute of research, it looks like Altitude above channel as a function name is pretty specific to Saga GIS. I do not know what it is called in other software packages. – jvangeld Mar 11 '11 at 19:46
For a (vague) definition and a sketch of Saga's algorithm see books.google.com/… – whuber Mar 11 '11 at 20:09
@whuber: Thank's for the link, i hope that it's still possible to do this using GRASS. I'll try to understand the algorithm – dickoa Mar 11 '11 at 20:15

## 2 Answers

The algorithm that SAGA uses is explained in this paper(p17): http://downloads.sourceforge.net/saga-gis/hbpl19_02.pdf

-
Very good article to understand without math the algorithm used...i'll try to check the C++ code later to understand the math and the implementation...nevertheless i'm new to grass, so it will take time for me to port code this into grass – dickoa Mar 13 '11 at 17:57
I have no experience with saga, but somebody wrote about calling saga from grass. May be more interesting than recoding the algorithm. sourceforge.net/apps/trac/saga-gis/wiki/… Apart from that, we (saga, grass, sextante) should really cooperate more on these algorithms and their documentation. – johanvdw Mar 13 '11 at 20:18
i agree with you about compatibility with those algorithms, and i would also like to thank you for the sagacvs repository for ubuntu, there was a seg_fault bug in 2.05 and now the 2.06 version runs fine. So i'll use grass as my main gis tool, and saga when it will be more appropriate. – dickoa Mar 14 '11 at 19:31

I contacted the author of r.stream.distance (and other modules, Jarek Jasiewicz) and he answered (reported here with permission):

"Very complex approach. But I wonder what means really different. The main reason could be that grass uses MFD apporach and generate different channel network and different overland flow distance structure.

[Approach] a) is SAGA, results are based on SFD flow direction but is GRASS based on MFD. In fact SAGA and GRASS nad TAUDEM uses the same algorithm, for elevation above channels and as I tested the results are very similar of course with overhead on different flow routing algorithms.

Hope it is the reason."

Furthermore he told me:

"BTW GRASS wetness index and stream distance won't be compatible due to different flow routing algorithm. GRASS wetness index is based on SFD flow routing algorithm. But as far as I know SAGA or TAS offer wetness index based on MFD"

Note: MFD based wetness index is planned to be added to GRASS.

-
Thank you very much for this information..so i'll stick with that result waiting for new feature in GRASS 7. BTW i know that you are really involved into GRASS development, and just want to thank you, it's really a wonderful piece of software. – dickoa Mar 13 '11 at 2:26
Glad to be helpful :) – markusN Mar 13 '11 at 11:58