r.topmodel

1   Data requirements

• Creating a topographic index statistics file
  • Input
    • topidx: Topographic index map
    • ntopidxclasses: Number of topographic index classes
  • Output
    • outtopidxstats: Output topographic index statistics file, which can be used for the topidxstats parameter later.
• Running TOPMODEL
  • Input
    • parameters: Model parameters file
    • topidxstats: Topographic index statistics file
    • input: Model input file
    • timestep: OPTIONAL for generating output for this specific time step
    • topidxclass: OPTIONAL for generating output for this specific topographic index class
  • Output
    • output: Model output file
• Creating a topographic index statistics file and running TOPMODEL using this statistics file
  • Input
    • topidx: Topographic index map
    • ntopidxclasses: Number of topographic index classes
    • parameters: Model parameters file
    • topidxstats: Topographic index statistics file. Give the same value as outtopidxstats to use this statistics file as input in the same run.
    • input: Model input file
    • timestep: OPTIONAL for generating output for this specific time step
    • topidxclass: OPTIONAL for generating output for this specific topographic index class
  • Output
    • outtopidxstats: Output topographic index statistics file, which will be used for the topidxstats parameter internally.
    • output: Model output file

2   Parameters file

The nch parameter (the number of distance increments) is automatically determined from pairs of d and Ad_r to avoid a mistake that the modeler forgets to update it after adding or removing them.

# Subcatchment name
Subcatchment 1

################################################################################
# A [m^2]: Total subcatchment area
3.31697E+07

################################################################################
# qs0 [m/h]: Initial subsurface flow per unit area
#		"The first streamflow input is assumed to represent
#		 only the subsurface flow contribution in the watershed."
#								- Liaw (1988)
0.000075

# lnTe [ln(m^2/h)]: Areal average of ln(T0)
4.

# m [m]: Scaling parameter
0.0125

# Sr0 [m]: Initial root zone storage deficit
0.0025

# Srmax [m]: Maximum root zone storage deficit
0.041

# td [h]: Unsaturated zone time delay per unit storage deficit if greater than 0
#  OR
# -alpha: Effective vertical hydraulic gradient if not greater than 0.
#
# For example, -10 means alpha=10.
60.

# vch [m/h]: Main channel routing velocity
20000.

# vr [m/h]: Internal subcatchment routing velocity
10000.

################################################################################
# infex: Calculate infiltration excess if not zero (integer)
0

# K0 [m/h]: Surface hydraulic conductivity
2.

# psi [m]: Wetting front suction
0.1

# dtheta: Water content change across the wetting front
0.1

################################################################################
# d [m]: Distance from the catchment outlet
#		The first value should be the mainstream distance from
#		the subcatchment outlet to the catchment outlet.
# Ad_r:  Cumulative area ratio of subcatchment (0.0 to 1.0)
#		The first and last values should be 0 and 1, respectively.

#   d  Ad_r
    0   0.0
 1000   0.2
 2000   0.4
 3000   0.6
 4000   0.8
 5000   1.0

3   Input file

For unit consistency over all input parameters, lengths are in meters and times are in hours except for rainfall and potential evapotranspiration time series, which are in time steps (dt).

The input file contains weather data like rainfall and potential evapotranspiration. Its format is like the following (a comment line starts with #):

# dt [h]: Time step
24

################################################################################
# R [m/dt]:  Rainfall
# Ep [m/dt]: Potential evapotranspiration

# R             Ep
0.000033        0.000000
0.000053        0.011938
0.004821        0.000000
.
.
.

4   Creating input data

It is very important to set the region and resolution correctly to avoid known issues with r.topidx. The following command will allow you to use the region and resolution of the elevation map (dem in this document) for further analysis.

g.region rast=dem

Create a drainage map (drain). Note that you don’t need to fill sinks in the elevation map because r.watershed uses a least-cost algorithm.

r.watershed elevation=dem drainage=drain

You will need a watershed boundary map (basin) to define the study area. You can find the coordinates of the watershed outlet (coord) by using the GUI of r.water.outlet interactively.

r.water.outlet input=drain output=basin coord=-109050.647255,1141527.27768

Now, you have to find the longest flow path (lfp vector) using this method.

Once you got lfp, use the following command to create outlet points (outlets) for subbasins along the longest flow path (LFP). This example creates subbasin outlets at every 10% of LFP. You can change this distance as you like. First, find the LFP category.

v.db.select -fc map=lfp columns=cat | head -1

Let’s say 123 is the category of LFP reported above. Create outlet points.

v.segment input=lfp output=outlets<<EOF
P 1 123 -10%
P 2 123 -20%
P 3 123 -30%
P 4 123 -40%
P 5 123 -50%
P 6 123 -60%
P 7 123 -70%
P 8 123 -80%
P 9 123 -90%
EOF

Now, you can create a subbasin vector at each outlet. Print the coordinates of all the outlets using the following command.

v.to.db -p map=outlets option=coor

Take x, y coordinates above and create subbasins at those outlets.

r.water.outlet input=drain output=subbasin1 coord=x1,y1 
r.water.outlet input=drain output=subbasin2 coord=x2,y2 
...

Or, you can use the following shell script to automate the above two steps.

#!/bin/sh
i=1
for coord in `v.to.db -p map=outlets option=coor | sed '/cat/d; s/^.*|\(.*\)|\(.*\)|.*$/\1,\2/'`
do
        r.water.outlet input=drain output=subbasin$i coord=$coord --o
        i=`expr $i + 1`
done

Pairs of distance and area, d and Ad_r pairs in the parameters file, can be obtained as follows. The following shell script assumes that the entire study area is located at the end point of LFP (i.e., the first line in output is 0 0). This script divides LFP into 10 pieces (“0 0” + “$i -le 9”), but you can change 9 to any number you want.

#!/bin/sh
echo "0 0"
total_area=`r.stats -an --q input=basin | sed 's/^[^ ]* //'`
total_length=`v.to.db -p --q map=lfp option=length | sed '/cat/d; s/.*|//; q'`
i=1
while [ $i -le 9 ]
do
        area=`r.stats -an --q input=subbasin$i | sed 's/^[^ ]* //'`
        perl -e 'printf "%f %f\n", '"0.1*$i*$total_length, ($total_area-$area)/$total_area"
        i=`expr $i + 1`
done
echo "$total_length 1"

Now, you have a list of distance and area pairs that define the study area. Create a parameters file (params.txt) and an input file (input.txt) based on the data requirements section in this document. Use the watershed boundary map (basin) as a mask to avoid unnecessary calculations outside the study area.

r.mask raster=basin

Create a topographic index map (topidx) using r.topidx.

r.topidx input=dem output=topidx

Create a topographic index statistics file (topidxstats.txt) by dividing the topographic index map into a number of classes (ntopidxclasses 30).

r.topmodel -p topidx=topidx ntopidxclasses=30 outtopidxstats=topidxstats.txt

Finally, you can run TOPMODEL and create an output file (output.txt).

r.topmodel param=params.txt topidxstats=topidxstats.txt input=input.txt output=output.txt

5   Calibrating model parameters using R

The following R script can read and write the parameters file:

write_rtopmodel_params <- function(file="params.txt",
				   qs0=NULL, lnTe=NULL, m=NULL, Sr0=NULL,
				   Srmax=NULL, td=NULL, vch=NULL, vr=NULL,
				   infex=NULL, K0=NULL, psi=NULL, dtheta=NULL){
	lines <- readLines(file)
	params <- read_rtopmodel_params_from_lines(lines)

	if(is.null(qs0) && is.null(lnTe) && is.null(m) && is.null(Sr0) &&
	   is.null(Srmax) && is.null(td) && is.null(vch) && is.null(vr) &&
	   is.null(infex) && is.null(K0) && is.null(psi) && is.null(dtheta))
		params$params

	if(!is.null(qs0)){
		params$params$qs0 <- qs0
		lines[params$linenos[1]] <- as.character(qs0)
	}
	if(!is.null(lnTe)){
		params$params$lnTe <- lnTe
		lines[params$linenos[2]] <- as.character(lnTe)
	}
	if(!is.null(m)){
		params$params$m <- m
		lines[params$linenos[3]] <- as.character(m)
	}
	if(!is.null(Sr0)){
		params$params$Sr0 <- Sr0
		lines[params$linenos[4]] <- as.character(Sr0)
	}
	if(!is.null(Srmax)){
		params$params$Srmax <- Srmax
		lines[params$linenos[5]] <- as.character(Srmax)
	}
	if(!is.null(td)){
		params$params$td <- td
		lines[params$linenos[6]] <- as.character(td)
	}
	if(!is.null(vch)){
		params$params$vch <- vch
		lines[params$linenos[7]] <- as.character(vch)
	}
	if(!is.null(vr)){
		params$params$vr <- vr
		lines[params$linenos[8]] <- as.character(vr)
	}
	if(!is.null(infex)){
		params$params$infex <- infex
		lines[params$linenos[9]] <- as.character(infex)
	}
	if(!is.null(K0)){
		params$params$K0 <- K0
		lines[params$linenos[10]] <- as.character(K0)
	}
	if(!is.null(psi)){
		params$params$psi <- psi
		lines[params$linenos[11]] <- as.character(psi)
	}
	if(!is.null(dtheta)){
		params$params$dtheta <- dtheta
		lines[params$linenos[12]] <- as.character(dtheta)
	}
	writeLines(lines, file)

	params$params
}

read_rtopmodel_params <- function(file="params.txt"){
	lines <- readLines(file)
	params <- read_rtopmodel_params_from_lines(lines)

	params$params
}

read_rtopmodel_params_from_lines <- function(lines){
	linenos <- c()
	params <- list()
	l <- 0
	for(i in 1:length(lines)){
		line <- lines[i]
		if(line == "" || substr(line, 1, 1) == "#")
			next
		l <- l + 1
		if(l == 1)
			next

		x <- as.numeric(strsplit(gsub("^[ \t]+|[ \t]+$", "", line),
					 "[ \t]+")[[1]])
		if(l >= 3 && l <= 14 && length(x) >= 1 && !is.na(x[1])){
			linenos[l-2] <- i
			if(l == 3)
				params$qs0 <- x[1]
			else if(l == 4)
				params$lnTe <- x[1]
			else if(l == 5)
				params$m <- x[1]
			else if(l == 6)
				params$Sr0 <- x[1]
			else if(l == 7)
				params$Srmax <- x[1]
			else if(l == 8)
				params$td <- x[1]
			else if(l == 9)
				params$vch <- x[1]
			else if(l == 10)
				params$vr <- x[1]
			else if(l == 11)
				params$infex <- x[1]
			else if(l == 12)
				params$K0 <- x[1]
			else if(l == 13)
				params$psi <- x[1]
			else
				params$dtheta <- x[1]
		}
	}

	list(linenos=linenos, params=params)
}

read_rtopmodel_output <- function(file="output.txt", name="Qt"){
	column <- ifelse(name=="Qt", 12,
		  ifelse(name=="qt", 23,
		  ifelse(name=="qo", 34,
		  ifelse(name=="qs", 45,
		  ifelse(name=="qv", 56,
		  ifelse(name=="S_mean", 67,
		  ifelse(name=="f", 78,
		  ifelse(name=="fex", 89, 12))))))))
	lines <- readLines(file)
	start <- 0
	output <- c()
	for(i in 1:length(lines)){
		line <- lines[i]
		if(substr(line, 1, 10) == "  timestep"){
			start <- i
			next
		}
		if(!start)
			next
		output[i-start] <- as.numeric(substr(line, column, column+9))
	}

	output
}

You can use any optimization algorithm for calibrating the model parameters. The following R scripts are meant to be used with ispso:

par.name <- c("qs0", "lnTe", "m", "Sr0", "Srmax", "td", "vr", "K0", "psi", "dtheta")
par.dim <- length(par.name);
par.min <- c(0, -7, 0.001, 0, 0.005, 0.001, 50, 0.0001, 0.01, 0.01)
par.max <- c(50, 10, 0.25, 0.01, 0.08, 40, 2000, 0.2, 0.5, 0.6)

run_rtopmodel <- function(x){
	cmd <- function(...) system(sprintf(...))

	if(length(x) != par.dim)
		stop()

	x[x<0] <- 0
	x[x>1] <- 1

	sink("x.txt", append=TRUE)
	cat(x, "\n")
	sink()

	parval <- par.min+(par.max-par.min)*x

	sink("parval.txt", append=TRUE)
	cat(parval, "\n")
	sink()

	write_text <- "write_rtopmodel_params(\"params.txt\""
	for(i in 1:par.dim)
		write_text <- paste(write_text, ",", par.name[i], "=", parval[i], sep="")
	write_text <- paste(write_text, ")", sep="")
	eval(parse(text=write_text))
	cmd("GRASS_VERBOSE=-1 r.topmodel --o param=params.txt topidxstats=topidxstats.txt input=input.txt output=output.txt")
}

source("ispso.R")
source("read_write_rtopmodel.R")
source("run_rtopmodel.R")

obs <- read.table("obs.txt")[[1]]

cmd <- function(...) system(sprintf(...))

ns <- c()
run <- 0

set_parameters <- function(s){
	########################################################################
	# Plot method
	#s$.plot_method <- "density"
	#s$.plot_method <- "movement"
	#s$.plot_method <- sprintf("%s,species", s$.plot_method)
	s$.plot_method <- ""
	s$.plot_delay <- 0
	########################################################################

	# Swarm size
	#s$S <- 10 + floor(2*sqrt(s$D))
	s$S <- 20

	# Maximum particle velocity
	s$vmax <- (s$xmax-s$xmin)*0.1

	# Maximum initial particle velocity
	s$vmax0 <- diagonal(s)*0.001

	# Stopping criteria: Stop if the number of exclusions per particle
	# since the last minimum is greater than exclusion_factor * max sol
	# iter / average sol iter. The more difficult the problem is (i.e.,
	# high max sol iter / average sol iter), the more iterations the
	# algoritm requires to stop.
	s$exclusion_factor <- 3
	# Maximum iteration
	s$maxiter <- 1000
	# Small positive number close to 0
	s$xeps <- 0.001
	s$feps <- 0.0001

	# Search radius for preys: One particle has two memories (i.e., x and
	# pbest).  When two particles collide with each other within prey, one
	# particle takes more desirable x and pbest from the two particles'
	# memories, and the other particle is replaced with a quasi-random
	# particle using scrambled Sobol' sequences (PREY).
	s$rprey <- diagonal(s)*0.0001

	# Nesting criteria for global and local optima using particles' ages
	# (NEST_BY_AGE).
	s$age <- 10

	# Speciation radius: Li (2004) recommends 0.05*L<=rspecies<=0.1*L.
	s$rspecies <- diagonal(s)*0.2

	# Nesting radius
	s$rnest <- diagonal(s)*0.01

	invisible(s)
}

s <- list()
s$f <- function(x){
	run <<- run + 1
	run_rtopmodel(x)

	sim <- read_rtopmodel_output(name="Qt")

	ns[run] <<- 1-sum((sim-obs)^2)/sum((mean(obs)-obs)^2)
	sink("ns.txt", append=TRUE)
	cat(ns[run], "\n")
	sink()
	printf("%d: %f\n", run, ns[run])

	return(-ns[run])
}
s$D <- par.dim
s$xmin <- rep(0, s$D)
s$xmax <- rep(1, s$D)
s <- set_parameters(s)

ret <- ispso(s)

bestx <- ret$pop[which(ret$pop[,"f"]==min(ret$pop[,"f"])),1:s$D]
bestpar <- par.min+(par.max-par.min)*bestx
names(bestpar) <- par.name

print("Best Parameters")
print(bestpar)

You should run doit.R within a GRASS session by sourcing the file.

> source("doit.R")

6   Comparison of data files with the original TOPMODEL

Title of the simulation
inputs.dat
subcat.dat
params.dat
topmod.out

NSTEP DT
R PE QOBS
.
.
.

NSC IMAP IOUT
SUBCAT
NAC AREA
AC ST
.
.
.
NCH
ACH D . . .
MAPFILE

SUBCAT
SZM T0 TD CHV RV SRMAX Q0 SR0 INFEX XK0 HF DTH

7   Example

r.topmodel_ex.zip

This example will give you a better idea of how to prepare input data. The zip file contains r.topmodel_ex/tmod9502 and r.topmodel_ex/r.topmodel directories. Basically, the two models are the same except that the former is for TMOD9502 and the latter is for r.topmodel.

8   Manual pages

• r.topidx man page
• r.topmodel man page

9   References

Beven, K., Lamb, R., Quinn, P., Romanowicz, R., Freer, J., 1995. TOPMODEL. In: Singh, V.P. (Ed.), Computer Models of Watershed Hydrology. Water Resources Publications, pp. 627-668.

Buytaert, W., 2009. TOPMODEL. https://source.ggy.bris.ac.uk/wiki/Topmodel, accessed on November 9, 2015.

Cho, H., 2000. GIS Hydrological Modeling System by Using Programming Interface of GRASS. Master’s Thesis, Department of Civil Engineering, Kyungpook National University, Korea.

Cho, H., Kim, D., Lee, K., 2014. Efficient Uncertainty Analysis of TOPMODEL Using Particle Swarm Optimization. Journal of the Korean Water Resources Association 47 (3), 285-295. doi:10.3741/JKWRA.2014.47.3.285.

Cho, H., Park, J., Kim, D., 2019. Evaluation of Four GLUE Likelihood Measures and Behavior of Large Parameter Samples in ISPSO-GLUE for TOPMODEL. Water 11 (3), 447. doi:10.3390/w11030447.

Conrad, O., 2003. SAGA-GIS Module Library Documentation (v2.1.3) Module TOPMODEL. http://www.saga-gis.org/saga_module_doc/2.1.3/sim_hydrology_2.html, accessed on November 9, 2015.

Liaw, S.C., 1988. Streamflow Simulation Using a Physically Based Hydrologic Model in Humid Forested Watersheds. Dissertation, Colorado State University, CO. p163.