Blog

Geomorphometry in R + SAGA + ILWIS + GRASS

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;
library(maptools)
library(rgdal)
library(RSAGA)
# ! first download and install SAGA GIS [http://www.saga-gis.org], ILWIS GIS [https://52north.org/download/Ilwis/52n-Ilwis-v3-05-02.zip] and GRASS GIS [http://grass.itc.it] to your machine!
rsaga.env(path="C:/Progra~1/saga_vc")
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 [http://spatial-analyst.net/wiki/index.php?title=MGI_/_Balkans_coordinate_systems]
setwd("C:/tmp")
# obtain the data:
download.file("http://geomorphometry.org/system/files/BaranjaHill.zip", destfile=paste(getwd(), "BaranjaHill.zip", sep="/"))
fname <- zip.file.extract(file="DEM25m.asc", zipname="BaranjaHill.zip")
file.copy(fname, "./DEM25m.asc", overwrite=TRUE)
list.files(pattern="asc")
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("http://spatial-analyst.net/CRS/gk_6.csy", 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:
rsaga.esri.to.sgrd(in.grids="dem25m.asc", out.sgrds="dem25m.sgrd", in.path=getwd())
# First, fill the spurious sinks:
rsaga.get.modules("ta_preprocessor")
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.modules("ta_channels")
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:
rsaga.sgrd.to.esri(in.sgrds="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
gridparameters(geoarc)
gridparameters(dem25m)
# 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 (http://spatial.nhh.no/R/Devel/spgrass6_0.6-1.zip (71))
# Location of your GRASS installation:
loc <- initGRASS("C:/GRASS", home=tempdir())
loc
# Import the ArcInfo ASCII file to GRASS:
parseGRASS("r.in.gdal") # commmand description
execGRASS("r.in.gdal", flags="o", parameters=list(input="DEM25m.asc", output="DEM"))
execGRASS("g.region", parameters=list(rast="DEM"))
gmeta6()
# 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("r.to.vect", parameters=list(input="streamt", output="streamt", feature="line"))
streamt <- readVECT6("streamt")
plot(streamt)
######## R script ###############

MORE READING:

  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.

Attachment:

SAGA_ILWIS_GRASS_0.zip

Book - Geomorphometry: Concepts, Software, Applications

A title in the Developments in Soil Science series volume 33. Read this book using the AmazonOnlineReader!

Table of content:

  • Foreword
  • **Part 1 **
  • 1. Geomorphometry: A Brief Guide
    1. Mathematical and Digital Models of the Land Surface
    1. DEM Production Methods and Sources
    1. Preparation of DEMs for Geomorphometric Analysis
    1. Geostatistical Simulation and Error Propagation in Geomorphometry
    1. Basic Land-Surface Parameters
    1. Land-Surface Parameters and Objects in Hydrology
    1. Land-Surface Parameters Specific to Topo-Climatology
    1. Landforms and Landform Elements in Geomorphometry
  • **Part 2 **
    1. Overview of Software Packages Used in Geomorphometry
    1. Geomorphometry in ESRI Packages
    1. Geomorphometry in SAGA
    1. Geomorphometry in ILWIS
    1. Geomorphometry in LandSerf
    1. Geomorphometry in MicroDEM
    1. Geomorphometry in TAS GIS
    1. Geomorphometry in GRASS GIS
    1. Geomorphometry in RiverTools
  • **Part 3 **
    1. Geomorphometry - A Key to Landscape Mapping and Modelling
    1. Soil Mapping Applications
    1. Vegetation Mapping Applications
    1. Applications in Geomorphology
    1. Modelling Mass Movements and Landslide Susceptibility
    1. Automated Predictive Mapping of Ecological Entities
    1. Geomorphometry and Spatial Hydrologic Modelling
    1. Applications in Meteorology
    1. Applications in Precision Agriculture
    1. The Future of Geomorphometry
  • Bibliography
  • Index
  • Colour Plate Section

How to cite?

The book can be cited as:

Hengl, T., Reuter, H.I. (eds) 2008. Geomorphometry: Concepts, Software, Applications. Developments in Soil Science, vol. 33, Elsevier, 772 pp.

To refer to a specific chapter, please write e.g.:

Pike, R.J., Evans, I.S., Hengl, T., 2008. Geomorphometry: a Brief Guide. In: Hengl, T. and Reuter, H.I. (Eds), Geomorphometry: Geomorphometry: Concepts, Software, Applications. Developments in Soil Science, vol. 33, Elsevier, 1-28 pp.

*Note that each specific chapter of the book has an unique DOI.

Specific chapter of the book can be obtained separately from the Elsevier’s ScienceDirect service.


How to obtain the data set?

The Baranja Hill data set used in this book can be obtained here. It comprises various data layers ranging from a TOPO DEM, SRTM DEM, original point measurements of heights, contour lines from the 1:25k and 1:5k scale topo-maps and various other thematic layers.


Attachment:

Pike_2008_Geomorphometry_ch1.pdf

Geomorphometry_bookflyer_Hengl.pdf

Latest Posts

Cover Design Contest for the Upcoming Book on Geomorphometry

Dear geomorphometry community,

We are pleased to invite submissions for a cover design contest for the second edition of the Geomorphometry book, to be published in 2026.

The submissions will be gathered in a poll, and the entire community will be able to vote for their favorite design.

If your design is selected, you will receive the appropriate credits, but would need to provide the necessary permissions to use the image.

You can submit your design by email before October 17th. Please ensure that the image is of at least 300 dpi resolution.

Get designing!

The editors,
Hannes Reuter
Carlos Grohmann
Vincent Lecours

Coffee Talk - Recent Research Progress in Geomorphometry in China

Recent Research Progress in Geomorphometry in China

Dr. Li-Yang Xiong
Nanjing Normal University, China

October 1st , 2025
8:00 MDT (UTC -6), 10:00 EDT (UTC -4), 11:00 BRT (UTC - 3), 15:00 BST (UTC +1), 16:00 CEST (UTC +2), 17:00 EEST (UTC +3), 22:00 CST (UTC +8)

Recording available in our YouTube channel

Bio: Dr. Li-Yang Xiong is a professor at the School of Geographical Science, Nanjing Normal University (NNU), China. He is currently responsible for managing NNU’s research in Digital Terrain Model and Digital Terrain Analysis. His main research interests include AI based terrain modelling, loess terrain feature characterization, landform evolution modeling, paleotopography reconstruction and geomorphological process mining. His recent work involves deep learning-based DEM reconstruction, geomorphology-oriented digital terrain analysis, and value-added digital terrain applications for geoscience. He also serves as Associate Editor for the journal Earth Surface Processes and Landforms and as an Editorial Board Member for International Journal of Geographical Information Science.

Abstract: In this talk, I will present some recent research achievements related to terrain modelling theory, terrain analysis method and terrain application in China. This terrain modeling theory focused on how we understand terrain knowledge and integrate it into AI methods for terrain reconstruction. In term of the terrain analysis method, the mathematical vector operation we believe should be highlighted in the research of Geomorphometry, which is suitable for multi-source data structure by considering the directional property of terrain parameters. Actually, this directional property should be made a full consideration for process- oriented geographical modeling and simulation. Lastly, I will show some terrain applications towards different typical geographical areas in China as well as global scale application.

PHD position in Italy

Dear colleagues,

I’m grateful if you can circulate information on this PhD opportunity in Italy. The potential candidates can contact me (strevisani@iuav.it) for further information. Here the main elements of the position:

Research topics: Predicting and supporting benthic and pelagic biodiversity through geomorphometry and machine learning

Link to the call (Italian and English): https://www.unipa.it/didattica/dottorati/dottorato-xli/bando-di-accesso-ciclo-41/

Position code [BIODIV.OGS]

Research headquarters OGS Trieste and University of Palermo

Funded by OGS - Istituto Nazionale di Oceanografia e di Geofisica Sperimentale

Key dates: Deadline: 7th August 2025 - 14:59 (Italian time)