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
