Tutorial on the use of topmodel in R

This is a reconstruction of a set of pages that we originally wrote on the R wiki, until the site disappeared without a trace. Sorry for the "under construction" nature of the documentation.

## install and load the required packages:

install.packages("topmodel")
install.packages("Hmisc")
library(topmodel)
library(Hmisc)

############ PART 1: running the rainfall-runoff model ##############

## Load the example dataset from the Huagrahuma catchment
## and attach it to the search path

data(huagrahuma)
attach(huagrahuma)

## Initial exploration of the data:

str(huagrahuma)
topidx
parameters
rain

plot(rain, type="h")

## run the model and visualise the outcome:

Qsim <- topmodel(parameters, topidx, delay, rain, ET0)
plot(Qsim, type="l", col="red")
points(Qobs)

## Evaluate the model with a performance metric

NSeff(Qobs, Qsim)

############ PART 2: Sensitivity analysis ##############

## let's try first to vary only one parameter
## The function runif() samples randomly from a uniform distribution

parameters["m"] <- runif(1, min = 0, max = 0.1)
parameters["m"]

## Run the model and evaluate with the Nash – Sutcliffe efficiency metric:

Qsim <- topmodel(parameters, topidx, delay, rain, ET0)
NSeff(Qobs, Qsim)

## What value do you get? Do you think this is a good simulation?
## Verify by plotting:

plot(Qsim, type="l", col="red")
points(Qobs)

## Now sample all parameters at random. We take a sample size of 100

n <- 100

qs0 <- runif(n, min = 0.0001, max = 0.00025)
lnTe <- runif(n, min = -2, max = 3)
m <- runif(n, min = 0, max = 0.1)
Sr0 <- runif(n, min = 0, max = 0.2)
Srmax <- runif(n, min = 0, max = 0.1)
td <- runif(n, min = 0, max = 3)
vch <- runif(n, min = 100, max = 2500)
vr <- runif(n, min = 100, max = 2500)
k0 <- runif(n, min = 0, max = 10)
CD <- runif(n, min = 0, max = 5)
dt <- 0.25

parameters <- cbind(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,CD,dt)

## run the model and evaluate with the Nash – Sutcliffe efficiency metric:
## Note: the function accepts a table of parameter sets
## (one parameter set per row)

NS <- topmodel(parameters, topidx, delay, rain, ET0, Qobs = Qobs)
max(NS)

## visualisation of the sensitivity using dotty plots:

plot(lnTe, NS, ylim = c(0,1))

############ PART 3: GLUE uncertainty analysis ##############

## choose a behavioural threshold and remove the “bad” parameter sets:

parameters <- parameters[NS > 0.3,]
NS <- NS[NS > 0.3]

## generate predictions for the behavioural parameter sets:

Qsim <- topmodel(parameters,topidx,delay,rain,ET0)

## (have a look at the predictions for the first time step:)

hist(Qsim[1,])

## construct weights based on the performance measure

weights <- NS - 0.3
weights <- weights / sum(weights)

## make prediction boundaries by weighted quantile calculation
## (we need the Hmisc package for that)

limits <- apply(Qsim, 1, "wtd.quantile", weights = weights,
probs = c(0.05,0.95), normwt=T)

plot(limits[2,], type="l")
points(limits[1,], type="l")
points(Qobs, col="red")

## how many measurements fall outside of the prediction limits?

outside <- (Qobs > limits[2,]) | (Qobs < limits[1,])
summary(outside)

## width of the prediction boundaries

mean(limits[2,] - limits[1,]) / mean(Qobs, na.rm=T)