02 April 2014

[R] Following Hydromad Tutorial


The following is the code, following Hydromad Tutorial by Felix Andrews.

---
# R script following Hydromad Tutorial by Felix Andrews
# link: http://hydromad.catchment.org/downloads/tutorial.pdf
# ------------------------------

# open hydromad
library(hydromad)

# open data
data(Cotter)

# checking data
xyplot(Cotter, screens=c("streamflow (mm/day)", 
                         "areal rain (mm/day)", 
                         "temperature (deg.c)"), xlab=NULL)
xyplot(window(Cotter, start="1974-01-01", 
              end="1975-01-01"))
monthlyPQE <- aggregate="" as.yearmon="" font="" mean="" otter="">
xyplot(monthlyPQE, screens=c("streamflow (mm/day)", 
                             "areal rain (mm/day)", "temperature (deg.c)"), xlab=NULL)
ok <- 1:2="" complete.cases="" font="" otter="">
with(Cotter, sum(Q[ok])/sum(P[ok]))
head(Cotter)
tail(Cotter)


# calculate rolling cross-correlation
# between rainfall and streamflow rises
# from 1980 to 1990 data period
x <- font="" otter="" rollccf="">
xyplot(x, xlim=extendrange(as.Date(c("1980-01-01", 
                                     "1990-01-01"))))

# estimateDelay(x)
delay <- estimatedelay="" font="" otter="">
delay

# model spec
# define data period name "data70s", "data80s", "data90s"
data70s <- font="" nbsp="" otter="" start="1970-01-01" window="">
                  end="1979-12-31")
data80s <- font="" nbsp="" otter="" start="1980-01-01" window="">
                  end="1989-12-31")
data90s <- font="" nbsp="" otter="" start="1990-01-01" window="">
                  end="1999-12-31")

# define routing and parameters
# using: data90s periode, "expuh" routing
# cwi (catchment wetness index) as 
# sma (soil moisture accounting)
# tau_s=5-100 (slow flow component), 
# tau_q=0-5 (fast flow comp), 
# v_s=0-1 (fraction volume in slow comp)
cotterModel <- data90s="" hydromad="" routing="expuh" sma="cwi" tau_q="c(0,5)," tau_s="c(5,100)," v_s="c(0,1))</font">
print(cotterModel)

# saving result to txt file
out<-capture .output="" cottermodel="" font="" print="">
cat(out,file="cotterModel1.txt",sep="\n",append=TRUE)

# model calibration using"fitByOptim"
cotterModel <- cottermodel="" font="" nbsp="" rfit="list(" routing="armax" sriv="" update="">
               order=c(n=2, m=1)))
cotterFit <- cottermodel="" fitbyoptim="" font="" method="PORT" samples="100,">

# try to simulate other period
# remember data90s is the calib period, 
# then we try to simulate data70s and data80s 
sim70s <- cotterfit="" newdata="data70s)</font" update="">
sim80s <- cotterfit="" newdata="data80s)</font" update="">
simAll <- cotterfit="" newdata="Cotter)</font" update="">

# run verification to whole dataset (70s,80s), 
# excluding data90s
dataVerif <- cotter="" dataverif="" font="" frame="" from="" make="">
dataVerif$Q[time(data90s)] <- and="" data90ss="" excluding="" font="" na="" nbsp="" q="" verify="" we="">
simVerif <- cotterfit="" dataverif="" font="" in="" newdata="dataVerif)" run="" save="" simverif="" update="">

# grouping all models using "runlist" func
allModels <- calibration="cotterFit," font="" runlist="" sim70s="" sim80s="" simverif="">
summary(allModels) # summary of model performance
print(allModels)

# plotting cotterFit with defined period
xyplot(cotterFit, with.P=TRUE, xlim=as.Date(c("1994-01-01", "1997-01-01")))

# plotting allModels
xyplot(allModels[2:3], scales=list(y=list(log=TRUE)))
summary(simAll, breaks="5 years")

# plotting performance over time
twoYrStats <- breaks="2 years" font="" simall="" summary="">
statSeries <- c="" font="" nbsp="" r.squared="" twoyrstats="">
               "r.sq.sqrt", "rel.bias", "runoff")]
statSeries[, 1:2] <- 0="" 1:2="" each="" font="" in="" max="" pick="" pmax="" segment="" statseries="" values="" will="">
c(xyplot(statSeries, type="s", lwd=2, 
         ylab="statistic", xlab=NULL), 
         'observed streamflow'=xyplot(observed(simAll)),
         layout=c(1,5), x.same=TRUE)+layer_(panel.refline(h=0, 
         v=time(statSeries)))

# to plot the flow duration curve for 
# modelled vs observed data
qqmath(cotterFit, scales=list(y=list(log=TRUE)), 
       type=c("l", "g"))
qqmath(allModels, type = c("l", "g"), scales = list(y = list(log = TRUE)),
       xlab = "Standard normal variate", ylab = "Flow (mm/day)",
       f.value = ppoints(100), tails.n = 50, as.table = TRUE)

# inverse fitting method
ihSpec <- data90s="" f="1," font="" hydromad="" routing="armax" sma="cwi" tw="10,">
osumm <- font="" ihspec="" nbsp="" rfit="sriv" trymodelorders="" update="">
          n=0:3, m=0:3, delay=0) # fixed delay and variable n and m

summary(osumm)

No comments: