Create a gist now

Instantly share code, notes, and snippets.

Embed
What would you like to do?
FUSE-RHydro tutorial 3
if(!require(RHydro)) install.packages("RHydro",repos="http://R-Forge.R-project.org")
library(RHydro)
temp <- read.csv("dummyData.csv")
DATA <- zooreg(temp[,2:4], order.by=temp[,1])
myDELTIM <- 1
# Step A: define the parameter ranges + mid range
mids <- c(60, 230, 342, 426)
require(tgp)
DefaultRanges <- data.frame(rbind(rferr_add = c(0,0),
rferr_mlt = c(1,1),
maxwatr_1 = c(25,500),
maxwatr_2 = c(50,5000),
fracten = c(0.05,0.95),
frchzne = c(0.05,0.95),
fprimqb = c(0.05,0.95),
rtfrac1 = c(0.05,0.95),
percrte = c(0.01,1000),
percexp = c(1,20),
sacpmlt = c(1,250),
sacpexp = c(1,5),
percfrac = c(0.05,0.95),
iflwrte = c(0.01,1000),
baserte = c(0.001,1000),
qb_powr = c(1,10),
qb_prms = c(0.001,0.25),
qbrate_2a = c(0.001,0.25),
qbrate_2b = c(0.001,0.25),
sareamax = c(0.05,0.95),
axv_bexp = c(0.001,3),
loglamb = c(5,10),
tishape = c(2,5),
timedelay = c(0.01,5) )
)
names(DefaultRanges) <- c("Min","Max")
nRuns <- 100
parameters <- lhs( nRuns, as.matrix(DefaultRanges) )
parameters <- data.frame(parameters)
names(parameters) <- c("rferr_add","rferr_mlt","maxwatr_1",
"maxwatr_2","fracten","frchzne",
"fprimqb","rtfrac1","percrte",
"percexp","sacpmlt","sacpexp",
"percfrac","iflwrte","baserte",
"qb_powr","qb_prms","qbrate_2a",
"qbrate_2b","sareamax","axv_bexp",
"loglamb","tishape","timedelay")
# Step B: run a multi-model calibration (use the Nash-Sutcliffe efficiency as objective function)
require(qualV)
indices <- rep(NA,4*nRuns)
discharges <- matrix(NA,ncol=4*nRuns,nrow=dim(DATA)[1])
kCounter <- 0
for (m in 1:4){
myMID <- mids[m]
for (pid in 1:nRuns){
kCounter <- kCounter + 1
ParameterSet <- as.list(parameters[pid,])
# Run FUSE Soil Moisture Accounting module
Qinst <- fusesma.sim(DATA,
mid=myMID,
modlist,
deltim=myDELTIM,
states=FALSE, fluxes=FALSE, fracstate0=0.25,
ParameterSet$rferr_add,
ParameterSet$rferr_mlt,
ParameterSet$frchzne,
ParameterSet$fracten,
ParameterSet$maxwatr_1,
ParameterSet$percfrac,
ParameterSet$fprimqb,
ParameterSet$qbrate_2a,
ParameterSet$qbrate_2b,
ParameterSet$qb_prms,
ParameterSet$maxwatr_2,
ParameterSet$baserte,
ParameterSet$rtfrac1,
ParameterSet$percrte,
ParameterSet$percexp,
ParameterSet$sacpmlt,
ParameterSet$sacpexp,
ParameterSet$iflwrte,
ParameterSet$axv_bexp,
ParameterSet$sareamax,
ParameterSet$loglamb,
ParameterSet$tishape,
ParameterSet$qb_powr)
# Run FUSE Routing module
Qrout <- fuserouting.sim(Qinst, mid=myMID,
modlist=modlist,
timedelay=ParameterSet$timedelay,
deltim=myDELTIM)
indices[kCounter] <- EF(DATA$Q,Qrout)
discharges[,kCounter] <- Qrout
}
}
# Step C: compare results
bestRun <- which(indices == max(indices))
bestModel <- function(runNumber){
if (runNumber<(nRuns+1)) myBestModel <- "TOPMODEL"
if (runNumber>(nRuns+1) & runNumber<(2*nRuns+1)) myBestModel <- "ARNOXVIC"
if (runNumber>(2*nRuns+1) & runNumber<(3*nRuns+1)) myBestModel <- "PRMS"
if (runNumber>(3*nRuns+1) & runNumber<(4*nRuns+1)) myBestModel <- "SACRAMENTO"
return(myBestModel)
}
bestModel(bestRun)
plot(coredata(DATA$Q),type="l",xlab="",ylab="Streamflow [mm/day]", lwd=0.5)
for(pid in 1:(4*nRuns)){
lines(discharges[,pid], col="gray", lwd=3)
}
lines(coredata(DATA$Q),col="black", lwd=1)
lines(discharges[,bestRun],col="red", lwd=1)
# How the best simulation of each model structure compare to each other?
bestRun0060 <- which(indices[1:nRuns] == max(indices[1:nRuns]))
bestRun0230 <- nRuns + which(indices[(nRuns+1):(2*nRuns)] == max(indices[(nRuns+1):(2*nRuns)]))
bestRun0342 <- 2*nRuns + which(indices[(2*nRuns+1):(3*nRuns)] == max(indices[(2*nRuns+1):(3*nRuns)]))
bestRun0426 <- 3*nRuns + which(indices[(3*nRuns+1):(4*nRuns)] == max(indices[(3*nRuns+1):(4*nRuns)]))
plot(coredata(DATA$Q),type="l",xlab="",ylab="Streamflow [mm/day]", lwd=1)
lines(discharges[,bestRun0060], col="green", lwd=1)
lines(discharges[,bestRun0230], col="blue", lwd=1)
lines(discharges[,bestRun0342], col="pink", lwd=1)
lines(discharges[,bestRun0426], col="orange", lwd=1)
legend("top",
c("TOPMODEL", "ARNOXVIC", "PRMS","SACRAMENTO"),
col = c("green", "blue", "pink", "orange"),
lty = c(1, 1, 1, 1))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment