Digital terrain analysis with RHydro - Revision history

Wbuytaert at 09:12, 28 August 2009

2009-08-28T09:12:56Z

Wbuytaert at 17:48, 13 February 2009

2009-02-13T17:48:39Z

Wbuytaert at 21:38, 26 August 2008

2008-08-26T21:38:56Z

New page

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]]

← Older revision		Revision as of 09:12, 28 August 2009
Line 1:		Line 1:
−	~~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):~~	+	The digital terrain analysis tools are now merged with the topmodel package. A revised guide can be found at [http://wiki.r-project.org/rwiki/doku.php?id=guides:tutorials:hydrological_data_analysis:dta\|the R wiki]
−
−	~~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]~~]

@@ Line 1: / Line 1: @@
-This is an example using the DEM supplied with the package. Load the package and the DEM:
+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(RHydro)
    data(huagrahuma.dem)
@@ Line 54: / Line 55: @@
    topidxmap <- topidx$atb * catchment
-Eliminate the NA values from the topidxmap:
+And make topographic index classes which can be used with topmodel.
-  topidxarray <- as.vector(topidxmap)
+   classes <- make.classes(topidxmap,16)
    plot(classes)
 The object ''classes'' can be used as input for [[topmodel]]