Geomorphometry in R + SAGA + ILWIS + GRASS

Published: Jun 30, 2009 by Tom Hengl

This is a sample code to process Baranja hill DEM using a combination of R, SAGA, ILWIS and Google Earth (under MS Windows machines). All four packages are available as open source or freeware (GE). You first need to obtain and install R (including the maptools, gstat, rgdal and RSAGA packages), ILWIS 3.5, SAGA GIS and GRASS GIS. After you finished installing R, SAGA, ILWIS, GRASS, copy the code down-below and start running it line by line. If you experience problems, send as a post via the R-sig-geo or ILWIS mailing list. Note: in order to use the functionality of ILWIS, SAGA, GRASS, you need to install them first. These are not internal packages in R!

Fig: Screenshots of R + SAGA/GRASS/ILWIS integration

######## R script ###############
# DEM processing and extraction of channel networks;
# ! first download and install SAGA GIS [], ILWIS GIS [] and GRASS GIS [] to your machine!
ILWIS <- "C:\Progra~1\N52\Ilwis35\IlwisClient.exe -C"
MGI_Z6 <- "+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000 +y_0=0 +ellps=bessel +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824 +units=m"
# The coordinate system is "MGI Zone 6"; see []
# obtain the data:
download.file("", destfile=paste(getwd(), "", sep="/"))
fname <- zip.file.extract(file="DEM25m.asc", zipname="")
file.copy(fname, "./DEM25m.asc", overwrite=TRUE)
dem25m <- readGDAL("DEM25m.asc")
proj4string(dem25m) <- CRS(MGI_Z6)
# view the data in ILWIS:
writeGDAL(dem25m[1], "dem25m.mpr", "ILWIS")
# set the right coordinate system!
download.file("", destfile=paste(getwd(), "gk_6.csy", sep="/"))
shell(cmd=paste(ILWIS, "setcsy dem25m.grf gk_6.csy -force"), wait=F)
shell(cmd=paste(ILWIS, "open dem25m.mpr -noask"), wait=F)
# make a relief view in ILWIS:
shell(cmd=paste(ILWIS, "run C:\Progra~1\Ilwis3\Scripts\Hydro-DEM\dem_visualization dem25m.mpr dem25m_c.mpr"), wait=F)
# load to SAGA and derive drainage network:"dem25m.asc", out.sgrds="dem25m.sgrd", in.path=getwd())
# First, fill the spurious sinks:
rsaga.get.usage("ta_preprocessor", 2)
rsaga.geoprocessor(lib="ta_preprocessor", module=2, param=list(DEM="dem25m.sgrd", RESULT="dem25m_f.sgrd", MINSLOPE=0.05))
# Second, extract the channel network:
rsaga.get.usage("ta_channels", 0)
rsaga.geoprocessor(lib="ta_channels", module=0, param=list(ELEVATION="dem25m.sgrd", CHNLNTWRK="channel_ntwrk.sgrd", CHNLROUTE="channel_route.sgrd", SHAPES="channels2.shp", INIT_GRID="dem25m.sgrd", DIV_CELLS=3, MINLEN=40))
channels <- readOGR("channels.shp", "channels")
spplot(channels["LENGTH"], col.regions=bpy.colors())
# derive Topographic Wetness Index:
rsaga.geoprocessor(lib="ta_hydrology", module=15, param=list(DEM="dem25m.sgrd", C="catharea.sgrd", GN="catchslope.sgrd", CS="modcatharea.sgrd", SB="TWI.sgrd", T=10))
# Export of grids to Google Earth using SAGA GIS:
rsaga.geoprocessor(lib="pj_proj4", 2, param=list(SOURCE_PROJ=paste('"', proj4string(dem25m), '"', sep=""), TARGET_PROJ=""+proj=longlat +datum=WGS84"", SOURCE="TWI.sgrd", TARGET="TWI_ll.sgrd", TARGET_TYPE=0, INTERPOLATION=0))
# export to PNG:
rsaga.geoprocessor(lib="io_grid_image", 0, param=list(GRID="TWI_ll.sgrd", FILE="TWI.png"))
# read back to R:"TWI_ll.sgrd", out.grids="TWI_ll.asc", prec=1, out.path=getwd())
TWI.ll <- readGDAL("TWI_ll.asc")
proj4string(TWI.ll) <- CRS("+proj=longlat +datum=WGS84")
TWI.kml <- GE_SpatialGrid(TWI.ll)
kmlOverlay(TWI.kml, kmlfile="TWI.kml", imagefile="TWI.png", name="Topographic Wetness Index")
# Optional: export to Google Earth using ILWIS:
dem25m.ll <- spTransform(dem25m[1], CRS("+proj=longlat +datum=WGS84"))
corrf <- (1+cos((dem25m.ll@bbox[2,2]+dem25m.ll@bbox[2,1])/2*pi/180))/2
geogrd.cell <- corrf*(dem25m.ll@bbox[1,2]-dem25m.ll@bbox[1,1])/dem25m@grid@cells.dim[1]
geoarc <- spsample(dem25m.ll, type="regular", cellsize=c(geogrd.cell,geogrd.cell))
gridded(geoarc) <- TRUE
# resample the map (Bilinear) to the new geographic grid:
shell(cmd=paste(ILWIS, " crgrf geoarc.grf ",geoarc@grid@cells.dim[[2]]," ",geoarc@grid@cells.dim[[1]]," -crdsys=LatlonWGS84 -lowleft=(",geoarc@grid@cellcentre.offset[[1]],",",geoarc@grid@cellcentre.offset[[2]],") -pixsize=",geoarc@grid@cellsize[[1]],sep=""), wait=F)
shell(cmd=paste(ILWIS, "dem25m_ll_c.mpr = MapResample(dem25m_c.mpr, geoarc, BiLinear)"), wait=F)
shell(cmd=paste(ILWIS, "open dem25m_ll_c.mpr -noask"))
shell(cmd=paste(ILWIS, "export tiff(dem25m_ll_c.mpr, dem25m_c.tif)"), wait=F)
# generate a KML (ground overlay):
dem25m.kml <- GE_SpatialGrid(geoarc)
kmlOverlay(dem25m.kml, "dem25m_c.kml", "dem25m_c.tif", name="Shaded relief in ILWIS")
# export of extracted channels to Google Earth:
proj4string(channels) <- CRS(MGI_Z6)
channels.ll <- spTransform(channels[3], CRS("+proj=longlat +datum=WGS84"))
writeOGR(channels.ll, "channels.kml", "channels", "KML")
# extract the drainage network in GRASS GIS:
library(spgrass6) # version => 0.6-1 ( (71))
# Location of your GRASS installation:
loc <- initGRASS("C:/GRASS", home=tempdir())
# Import the ArcInfo ASCII file to GRASS:
parseGRASS("") # commmand description
execGRASS("", flags="o", parameters=list(input="DEM25m.asc", output="DEM"))
execGRASS("g.region", parameters=list(rast="DEM"))
# extract the drainage network:
execGRASS("r.watershed", flags=c("m", "overwrite"), parameters=list(elevation="DEM", stream="stream", threshold=as.integer(50)))
# thin the raster map so it can be converted to vectors:
execGRASS("r.thin", parameters=list(input="stream", output="streamt"))
# convert to vectors:
execGRASS("", parameters=list(input="streamt", output="streamt", feature="line"))
streamt <- readVECT6("streamt")
######## R script ###############


  1. Bivand, R. 2005. Interfacing GRASS 6 and R. Status and development directions, GRASS Newsletter, 3, 11–16.
  2. Bivand, R., Pebesma, E., Rubio, V., 2008. Applied Spatial Data Analysis with R. Use R Series, p. 400. Springer, Heidelberg, pp. 378.
  3. Brenning, A. 2008. Statistical geocomputing combining R and SAGA: The example of landslide susceptibility analysis with generalized additive models. In: J. Böhner, T. Blaschke & L. Montanarella (eds.), SAGA - Seconds Out (= Hamburger Beiträge zur Physischen Geographie und Landschaftsökologie, 19), 23-32.
  4. Conrad, O. 2007. SAGA — Entwurf, Funktionsumfang und Anwendung eines Systems fur Automatisierte Geowissenschaftliche Analysen, Ph.D. thesis, University of Gottingen, Gottingen.
  5. Grohmann, C.H. 2004. Morphometric analysis in Geographic Information Systems: applications of free software GRASS and R. Computers & Geosciences, 30 (9-10):1055-1067.
  6. Hengl, T., 2009. A Practical Guide to Geostatistical Mapping. 2nd Ed, University of Amsterdam, 291 p.
  7. Neteler, M. and Mitasova, H. 2008. Open Source GIS: A GRASS GIS Approach, Springer, New York, 3rd edn.


dataset script


Latest Posts

ANADEM: A Digital Terrain Model for South America

There is a new paper (open access) describing a Machine Learning-based DTM for South America:

Laipelt L., Andrade B.C., Collischonn W., Teixeira A.A., Paiva R.C.D., Ruhoff A., 2024. ANADEM: A Digital Terrain Model for South America. Remote Sensing 16(13):2321.

GRASS GIS 8.4.0RC1 release

The GRASS GIS 8.4.0RC1 release provides more than 515 improvements and fixes with respect to the release 8.3.2.

Check the full announcement at

Please support in testing this release candidate.

Best BiCubic Method to Compute the Planimetric Misregistration between Images with Sub-Pixel Accuracy: Application to Digital Elevation Models

There is a new paper (open access) describing a novel method to estimate sub-pixel planimetric displacements between two DEMs:

Riazanoff, S.; Corseaux, A.; Albinet, C.; Strobl, P.A.; López-Vázquez, C.; Guth, P.L.; Tadono, T. Best BiCubic Method to Compute the Planimetric Misregistration between Images with Sub-Pixel Accuracy: Application to Digital Elevation Models. ISPRS Int. J. Geo-Inf. 13, 96.