14 April 2014

[R]Hydromad on Cikapundung-update on monthly analysis

After discussing my codes on the previous post about using Hydromad to Cikapundung dataset, It turned out that the codes were for daily analysis. Where as to convert it to monthly analysis, we should set the warm up and objective function.

We add:

  • hydromad.options(warmup=12) to set montly analysis
  • objective=~hmadstat("r.squared") to use Nash-Sutcliffe Efficiency/NSE as the objective function for our calibration.

Then the whole code will look like this,

# open hydromad
library(hydromad)
setwd("~/R practice-cikapundung")

# load data (flow, rain, temp)
Flow <- read.csv("flowlembang2.csv")
head(Flow)

##       Date Flow
## 1 1/1/2007 2.57
## 2 1/2/2007 2.57
## 3 1/3/2007 2.50
## 4 1/4/2007 2.44
## 5 1/5/2007 2.41
## 6 1/6/2007 2.39

Rain <- read.csv("rainlembang.csv")
head(Rain)

##       Date  Rain
## 1 1/1/2007 2.000
## 2 1/2/2007 2.000
## 3 1/3/2007 2.125
## 4 1/4/2007 2.250
## 5 1/5/2007 2.375
## 6 1/6/2007 2.500

Temp <- read.csv("templembang.csv")
head(Temp)

##       Date MaxT
## 1 1/1/2007 18.9
## 2 1/2/2007 18.4
## 3 1/3/2007 18.3
## 4 1/4/2007 17.9
## 5 1/5/2007 19.6
## 6 1/6/2007 20.2


# Convert the date column
Flow$Date <- as.Date(Flow[, 1], "%m/%d/%Y")
Rain$Date <- as.Date(Rain[, 1], "%m/%d/%Y")
Temp$Date <- as.Date(Temp[, 1], "%m/%d/%Y")

# use package zoo
library(zoo)
tsQ <- zoo(Flow$Flow, order.by = Flow$Date, frequency = 1)
tsP <- zoo(Rain$Rain, order.by = Rain$Date, frequency = 1)
tsT <- zoo(Temp$MaxT, order.by = Temp$Date, frequency = 1)

# merge data
Cikapundung <- merge(P = tsP, Q = tsQ, E = tsT, all = F)
head(Cikapundung)

##                P    Q    E
## 2007-01-01 2.000 2.57 18.9
## 2007-01-02 2.000 2.57 18.4
## 2007-01-03 2.125 2.50 18.3
## 2007-01-04 2.250 2.44 17.9
## 2007-01-05 2.375 2.41 19.6
## 2007-01-06 2.500 2.39 20.2

monthlyPQE <- aggregate(Cikapundung, as.yearmon, mean)

# model spec define data period name 'ts.cal', 'ts.val', 'ts.later' split
# the data into different sections
ts.cal <- window(Cikapundung, start = "2007-01-01", end = "2008-12-31")
ts.val <- window(Cikapundung, start = "2008-01-01", end = "2009-12-31")
ts.later <- window(Cikapundung, start = "2009-01-01", end = "2010-12-31")

# using monthly data with Armax and Expuh routing armax
hydromad.options(warmup = 12)
ckpModel.armax <- hydromad(ts.cal, sma = "cmd", routing = "armax", rfit = list("sriv", 
    order = c(n = 2, m = 1)))
ckpFit.armax.NSE <- fitByOptim(ckpModel.armax, objective = ~hmadstat("r.squared")(Q, 
    X), samples = 100, method = "PORT")

## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: did not converge after 12 iterations
## Warning: false convergence (8)

summary(ckpFit.armax.NSE)

## 
## Call:
## hydromad(DATA = ts.cal, sma = "cmd", routing = "armax", rfit = list("sriv", 
##     order = c(n = 2, m = 1)), f = 0.4026, e = 0.8829, d = 146)
## 
## Time steps: 719 (3 missing).
## Runoff ratio (Q/P): (3.92 / 9.34) = 0.419
## rel bias: -0.136
## r squared: -0.35
## r sq sqrt: -0.878
## r sq log: -0.878
## 
## For definitions see ?hydromad.stats

coef(ckpFit.armax.NSE)

##         f         e         d     shape       a_1       a_2       b_0 
##   0.40260   0.88290 146.00000   0.00000   1.98295  -0.98306   0.01065 
##       b_1     delay 
##  -0.01054   2.00000

print(ckpFit.armax.NSE)

## 
## Hydromad model with "cmd" SMA and "armax" routing:
## Start = 2007-01-01, End = 2008-12-31
## 
## SMA Parameters:
##       f        e        d    shape  
##   0.403    0.883  146.000    0.000  
## Routing Parameters:
##     a_1      a_2      b_0      b_1    delay  
##  1.9829  -0.9831   0.0107  -0.0105   2.0000  
## TF Structure: S + Q (two stores in parallel)
##     Poles:0.9915+0.0063i, 0.9915-0.0063i 
## 
## Fit: ($fit.result)
## fitByOptim(MODEL = ckpModel.armax, objective = ~hmadstat("r.squared")(Q, 
##     X), method = "PORT", samples = 100)
##      135 function evaluations in 11.37 seconds 
## 
## Routing fit info:  list(converged = FALSE, iteration = 12) 
## 
## Message: false convergence (8)

xyplot(ckpFit.armax.NSE, with.P = TRUE, type = c("l", "g"))

plot of chunk qplot 


## expuh
hydromad.options(warmup = 12)
ckpModel.expuh <- hydromad(ts.cal, sma = "cwi", routing = "expuh", tau_s = c(5, 
    100), tau_q = c(0, 5), v_s = c(0, 1))
ckpFit.expuh.NSE <- fitByOptim(ckpModel.expuh, objective = ~hmadstat("r.squared")(Q, 
    X), samples = 100, method = "PORT")
summary(ckpFit.expuh.NSE)

## 
## Call:
## hydromad(DATA = ts.cal, tau_s = 100, tau_q = 2.9763, v_s = 0.892532, 
##     sma = "cwi", routing = "expuh", tw = 6.96977, f = 8, scale = 0.00547393)
## 
## Time steps: 719 (1 missing).
## Runoff ratio (Q/P): (3.92 / 9.32) = 0.421
## rel bias: 1.98e-17
## r squared: -0.303
## r sq sqrt: -0.534
## r sq log: -0.545
## 
## For definitions see ?hydromad.stats

coef(ckpFit.expuh.NSE)

##        tw         f     scale         l         p     t_ref     tau_s 
## 6.970e+00 8.000e+00 5.474e-03 0.000e+00 1.000e+00 2.000e+01 1.000e+02 
##     tau_q       v_s 
## 2.976e+00 8.925e-01

print(ckpFit.expuh.NSE)

## 
## Hydromad model with "cwi" SMA and "expuh" routing:
## Start = 2007-01-01, End = 2008-12-31
## 
## SMA Parameters:
##       tw         f     scale         l         p     t_ref  
##  6.96977   8.00000   0.00547   0.00000   1.00000  20.00000  
## Routing Parameters:
##   tau_s    tau_q      v_s  
## 100.000    2.976    0.893  
## TF Structure: S + Q (two stores in parallel)
##     Poles:0.7146, 0.99 
## 
## Fit: ($fit.result)
## fitByOptim(MODEL = ckpModel.expuh, objective = ~hmadstat("r.squared")(Q, 
##     X), method = "PORT", samples = 100)
##      221 function evaluations in 11.87 seconds

xyplot(ckpFit.expuh.NSE, with.P = TRUE, type = c("l", "g"))

plot of chunk qplot

You might also get the ## Warning: did not converge after 12 iterations. Please refer to the following discussion on hydromad user group. The fitByOptim command is the one that iterate the model equation with the best fit parameters. It saves us time for not have to re-input all the parameters from the print(xxxx), coef(xxx), or summary(xxx).

If we look at the plot, we would agree that both of the fitting (pink line) don't capture many parts of the observed data (blue line). However it is normal because we only had 4 years of data in total and only use one year (2007-2008) period as calibration data.

Additional reference:

Good luck {@dasaptaerwin}

Posted by at
Labels: , , ,

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)