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:
Post a Comment