Home
Overview
Monitoring
Projects
CV
Publications
R code
Links

Miscellaneous functions

NSeff()

Convenient implementation of the Nash and Sutcliffe efficiency calculation

Eff <- NSeff(Qobs,Qsim)

topidx()

Calculation of the topographic index starting from a DEM matrix (for a more complete example see below)

topographic_index <- topidx(DEM, resolution=25)



Data import and preprocessing for topmodel()

In this example, all input files should be plain text (.txt) files. (Each decent spreadsheet converts easily to text files with "save as"). All inputs are given in meters per time unit

Rainfall

In this case we have the rainfall data from a tipping bucket gauge. This means that data are stored as the moments of the tipping of the bucket. The resolution of the gauges is 0.2 mm, so each entry represents 0.2 mm of rainfall. An input file, say rain.txt could look like this:

24/05/2001 00:52:46
24/05/2001 01:07:37
24/05/2001 01:22:51
24/05/2001 01:27:00
25/05/2001 01:05:17
...

The code we use to import these data into R is:

rain <- read.table("rain.txt")
rain <- paste(rain[,1],rain[,2])
rain <- strptime(as.character(rain), format = "%d/%m/%Y %H:%M:%S")
rain <- as.POSIXct(rain, tz="GMT")

NOTE: be sure to set the time zone correctly. Importing data from a region without summer time on a computer configured for a time zone with summer time is problematic. You can change the time zone in R with:

Sys.putenv(TZ="GMT")

Then we transform the data to rainfall per time unit, say, per 15 min:

rain <- rain+15*60
rain <- cut(c(trunc(rain[1],"hour"),rain), "15 min")
rain <- tapply(rep(0.2,length(rain)),rain,"sum")
rain[1] <- rain[1]-0.2

Note the first line of this code block: we shift all the entries 15 min. to the future. This makes sure that the rainfall associated with each time step is the rainfall that occurred over the past time interval (and not the coming interval).

Discharge

Discharge is easier to import and does not require much preprocessing. The same goes for rainfall which is already a time series. For instance, if the input file is:

04/09/2001 15:30 0.0492
04/09/2001 16:00 0.0495
04/09/2001 16:30 0.0503
04/09/2001 17:00 0.0504
04/09/2001 17:30 0.0500
...

Import these data:

input <- read.table("discharge.txt")
time <- paste(input[,1],input[,2])
time <- strptime(as.character(time), format = "%d/%m/%Y %H:%M")
time <- as.POSIXct(time, tz="GMT")
level <- input[,3]
discharge <- data.frame(time,level)
rm(input,time,level)

Strictly, we do not need the time information in the discharge and rain objects, but it is handy to have them around, for instance to check that both series are synchronous.

Evapotranspiration

Same procedure as discharge

The topographic index

Some GIS programs can calculate the topographic index, but the package provides a function too. If you use the R function, you should delineate your catchment in your preferred GIS application and export a DEM of your catchment as a text file (ascii). Pixels outside the catchment area should be given a distinct value that can be set to NA in R.

For now, the GIS functions of this package are quite limited. The DEM has to be imported as a matrix, which can then be processed by topidx(). Take for instance this minimalistic DEM, saved in a test file called "DEM.txt". Values outside the catchment are given the value -9999 (this can be any other value):

-9999 -9999 828.9 835.6 -9999
818.3 826.0 830.7 834.5 836.0
817.1 824.0 825.2 833.3 836.9
816.5 820.0 824.1 330.8 -9999
810.7 815.6 822.2 -9999 -9999

This file can be imported and processed in R with:

DEM <- read.table("DEM.txt")
DEM <- as.matrix(DEM)
# remove the values outside the catchment:
DEM[DEM==-9999] <- NA
# plot the DEM:
image(DEM)
topindex <- topidx(DEM, resolution=25)
topidxclasses <- topidxclasses(topindex,16)