Digital terrain analysis with RHydro
This is an example using the DEM supplied with the package. Load the package and the DEM:
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
Eliminate the NA values from the topidxmap:
topidxarray <- as.vector(topidxmap) topidxarray <- topidxarray[!is.na(topidxarray)]
And make topographic index classes which can be used with topmodel:
classes <- make.classes(topidxarray,16) plot(classes)
The object classes can be used as input for topmodel