Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question
    
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
1  
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 2

up vote 1 down vote accepted

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

share|improve this answer
    
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.

share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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