Digital terrain analysis with RHydro

This is an example using the DEM supplied with the package. Load the package and the DEM (we also need the package topmodel for the make.classes() function):

 library(topmodel)
 library(RHydro)
 data(huagrahuma.dem)

Fill depression and look at the difference between the original and filled map:

 dem.filled <- sinkfind(huagrahuma.dem, cellsize=25, degree=0.1)
 library(lattice)
 levelplot(dem.filled)
 
 difference <- dem.filled - huagrahuma.dem
 max(difference)
 min(difference)
 levelplot(difference)

Topographic index calculation:

 topidx <- atb(dem.filled, cellsize=25)
 class(topidx)
 str(topidx)
 levelplot(topidx$atb)
 levelplot(topidx$area)

If you want, you can compare with the analysis on the non-filled DEM that contains obvious artefacts

 old <- atb(huagrahuma.dem, cellsize=25)
 levelplot(old$area)
 levelplot(topidx$area - old$area)

Have a look at the potential outlet:

 outlet(topidx$area,c(28,8),2)

We know from manual catchment delineation (on a topographic map) that the area of the catchment is around 2.4 km2, so the following cells are likely to contain the outlet: c(28,8); c(29,8); c(29,9). Just try them out with the subcatch() function and take the one you like.

 catchment <- subcatch(dem.filled, c(29,8)) 
 image(catchment)
 catchment[catchment == 0] <- NA

 flowlength <- flowlength(dem.filled, c(29,8))

Headwater cells are identified on the basis of a topographic index or an accumulation area threshold. This is an “intuitive” choice:

 rivers <- river(dem.filled,topidx$atb,topidx$area,cellsize=25,thatb=12.35,tharea=10000)
 image(rivers)

Play around with the values of thatb and tharea and look at the impact…

Mask the rivers and topographic index maps outside the catchment:

 rivers <- rivers * catchment
 topidxmap <- topidx$atb * catchment

And make topographic index classes which can be used with topmodel.

 classes <- make.classes(topidxmap,16)
 plot(classes)

The object classes can be used as input for topmodel