FUSE-RHydro tutorial 3

FUSE_RHydro03.R

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