Digital terrain analysis with RHydro

From SourceWiki
Revision as of 21:38, 26 August 2008 by Wbuytaert (talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

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