Simhyd documentation

R/simhyd.R

@@ -1,30 +1,74 @@
 ## hydromad: Hydrological Modelling and Analysis of Data
-## Rewrite
-## willem 20151014
+## SimHyd
 ## Coded based on diagram and description in:
 # Chiew et al 2009 WATER RESOURCES RESEARCH, VOL. 45, W10414, doi:10.1029/2008WR007338, 2009
 
-## SimHyd model
-
-
 #' SimHyd model
 #'
-#' Description here. Requires documentation.
+#' Coded based on diagram and description in:
+#' Chiew et al 2009 WATER RESOURCES RESEARCH, 
+#' VOL. 45, W10414, doi:10.1029/2008WR007338, 2009
+#'
+#' @name simhyd
+#' @aliases simhyd.sim simhyd.ranges
+#'
+#' @param DATA time-series-like object with columns P (precipitation, mm) and E
+#' (potential evapo-transpiration, mm).
+#' @param U effective rainfall series.
+#' @param INSC Interception storage capacity (mm)
+#' @param COEFF maximum infiltration loss (mm)
+#' @param SQ Infiltration loss exponent
+#' @param SMSC Soil moisture store capacity (mm)
+#' @param SUB Constant of proportionality in interflow equation
+#' @param CRAK Constant of proportionality in groundwater recharge equation
+#' @param K Baseflow linear recession parameter
+#' @param etmult mutliplier to convert daily maximum temperature to estimates
+#'  of potential evaporation
+#' @param GWt0 Initial state of the groundwater store (mm)
+#' @param SMSt0 Initial state of soil moisture store (mm)
+#' @param return_state to return the series U, S (storage) and ET
+#' (evapotranspiration).
+#' @author Willem Vervoort \email{willemvervoort@@gmail.com} and Joseph Guillaume
+#' \email{josephguillaume@@gmail.com}
+#' @seealso \code{\link{hydromad}(sma = "simhyd", routing = "simhydrouting")} to
+#' work with models as objects (recommended).
+#' @references Chiew, F. H. S., Teng, J., Vaze, J., Post, D. A., Perraud, J. M.,
+#' Kirono, D. G. C., and Viney, N. R. (2009), 
+#' Estimating climate change impact on runoff across southeast Australia: 
+#' Method, results, and implications of the modeling method, 
+#' \emph{Water Resour. Res.}, 45, W10414, \url{doi:10.1029/2008WR007338}.
+#'
+#' @keywords models
+#' @examples
+#'
+#' ## view default parameter ranges:
+#' str(hydromad.getOption("simhyd"))
 #'
-#' @aliases simhyd.sim simhyd.ranges simhydrouting.sim simhydrouting.ranges
+#' data(HydroTestData)
+#' mod0 <- hydromad(HydroTestData, sma = "simhyd", routing = "simhydrouting")
+#' mod0
 #'
-#' @param DATA Placeholder
-#' @param INSC Placeholder
-#' @param COEFF Placeholder
-#' @param SQ Placeholder
-#' @param SMSC Placeholder
-#' @param SUB Placeholder
-#' @param CRAK Placeholder
-#' @param K Placeholder
-#' @param etmult Placeholder
-#' @param GWt0 Placeholder
-#' @param SMSt0 Placeholder
-#' @param return_state Placeholder
+#' data(Cotter)
+#' x <- Cotter[1:1000]
+#'
+#' # Specify simhyd model
+#' mod0 <- hydromad(x, sma = "simhyd", routing = "simhydrouting")
+#' # specify parameter ranges for a few parameters
+#' mod0 <- update(mod0,  COEFF=c(0,400), SQ=c(0.1,5), 
+#' K=c(0,1))
+#' # Allow etmult to vary, because we're using temperature data instead of PET.
+#' mod0 <- update(mod0, etmult = c(0.05, 1.5))
+#' 
+#' mod0
+#'
+#' ## now try to fit the free parameters
+#' set.seed(10)
+#' fit1 <- fitByOptim(mod0)
+#'
+#' fit1
+#' summary(fit1)
+#' xyplot(fit1)
+
 
 #' @useDynLib hydromad _hydromad_simhyd_sim
 #' @useDynLib hydromad
@@ -148,50 +192,6 @@ simhyd.sim <-
     }
   }
 
-# Routing based on Muskinghum
-#' @export
-simhydrouting.sim <- function(U, DELAY = 1, X_m = 0.2,
-                              epsilon = hydromad.getOption("sim.epsilon"),
-                              return_components = FALSE) {
-  X <- rep(0, length(U))
-  inAttr <- attributes(U)
-  U <- as.ts(U)
-  bad <- is.na(U)
-  U[bad] <- 0
-  if (2 * DELAY * X_m < 1 & 2 * DELAY * (1 - X_m) > 1) {
-    # Muskingum components
-    C0 <- (-DELAY * X_m + 0.5) / (DELAY * (1 - X_m) + 0.5)
-    # print(C0)
-    C1 <- (DELAY * X_m + 0.5) / (DELAY * (1 - X_m) + 0.5)
-    # print(C1)
-    C2 <- (DELAY * (1 - X_m) - 0.5) / (DELAY * (1 - X_m) + 0.5)
-    # print(C2)
-  } else {
-    C0 <- 0
-    C1 <- 1
-    C2 <- 0
-    # print("model parameters adjusted")
-  }
-
-  if (C0 + C1 + C2 != 1) {
-    C0 <- 0
-    C1 <- 1
-    C2 <- 0
-    # print("model parameters adjusted again")
-  }
-  # print(C0+C1+C2)
-  # if (round(C0+C1+C2)!=1)  C0 <- 0; C1 <- 1; C2 <- 0
-
-  X[1] <- U[1]
-  for (t in 1:(length(U) - 1)) {
-    X[t + 1] <- C0 * U[t + 1] + C1 * U[t] + C2 * X[t]
-    # print(X[t+1])
-  }
-  X[abs(X) < epsilon] <- 0
-  X[bad] <- NA
-  attributes(X) <- inAttr
-  X
-}
 
 #' @export
 simhyd.ranges <- function() {
@@ -207,10 +207,3 @@ simhyd.ranges <- function() {
   )
 }
 
-#' @export
-simhydrouting.ranges <- function() {
-  list(
-    DELAY = c(0.1, 5),
-    X_m = c(0.01, 0.5)
-  )
-}

R/simhydrouting.R

@@ -0,0 +1,94 @@
+## hydromad: Hydrological Modelling and Analysis of Data
+## SimHydrouting
+## Coded based on diagram and description in:
+# Chiew et al 2009 WATER RESOURCES RESEARCH, VOL. 45, W10414, doi:10.1029/2008WR007338, 2009
+
+#' Routing based on Muskingum
+#' 
+#' Coded based on diagram and description in:
+#' Chiew et al 2009 WATER RESOURCES RESEARCH, 
+#' VOL. 45, W10414, doi:10.1029/2008WR007338, 2009
+
+
+#' @name simhydrouting
+#' @aliases simhydrouting.sim simhydrouting.ranges
+#'
+#' @param U effective rainfall series from a SMA object
+#' @param DELAY delay parameter in the Muskingum routing (Days)
+#' @param X storage weight parameter in Muskingum routing (-)
+#' @param return_components included here to align with hydromad structure
+#' @author Willem Vervoort \email{willemvervoort@@gmail.com} and Joseph Guillaume
+#' \email{josephguillaume@@gmail.com}
+#' @seealso \code{\link{hydromad}(sma = "simhyd", routing = "simhydrouting")} to
+#' work with models as objects (recommended).
+#' @references Chiew, F. H. S., Teng, J., Vaze, J., Post, D. A., Perraud, J. M.,
+#' Kirono, D. G. C., and Viney, N. R. (2009), 
+#' Estimating climate change impact on runoff across southeast Australia: 
+#' Method, results, and implications of the modeling method, 
+#' \emph{Water Resour. Res.}, 45, W10414, \url{doi:10.1029/2008WR007338}.
+#'
+#' @keywords models
+#' @examples
+#'
+#' ## view default parameter ranges:
+#' str(hydromad.getOption("simhydrouting"))
+#' 
+#'
+#' data(HydroTestData)
+#' mod0 <- hydromad(HydroTestData, sma = "simhyd", routing = "simhydrouting")
+#' mod0
+#'
+
+
+#' @export
+simhydrouting.sim <- function(U, DELAY = 1, X_m = 0.2,
+                              epsilon = hydromad.getOption("sim.epsilon"),
+                              return_components = FALSE) {
+  X <- rep(0, length(U))
+  inAttr <- attributes(U)
+  U <- as.ts(U)
+  bad <- is.na(U)
+  U[bad] <- 0
+  if (2 * DELAY * X_m < 1 & 2 * DELAY * (1 - X_m) > 1) {
+    # Muskingum components
+    C0 <- (-DELAY * X_m + 0.5) / (DELAY * (1 - X_m) + 0.5)
+    # print(C0)
+    C1 <- (DELAY * X_m + 0.5) / (DELAY * (1 - X_m) + 0.5)
+    # print(C1)
+    C2 <- (DELAY * (1 - X_m) - 0.5) / (DELAY * (1 - X_m) + 0.5)
+    # print(C2)
+  } else {
+    C0 <- 0
+    C1 <- 1
+    C2 <- 0
+    # print("model parameters adjusted")
+  }
+
+  if (C0 + C1 + C2 != 1) {
+    C0 <- 0
+    C1 <- 1
+    C2 <- 0
+    # print("model parameters adjusted again")
+  }
+  # print(C0+C1+C2)
+  # if (round(C0+C1+C2)!=1)  C0 <- 0; C1 <- 1; C2 <- 0
+
+  X[1] <- U[1]
+  for (t in 1:(length(U) - 1)) {
+    X[t + 1] <- C0 * U[t + 1] + C1 * U[t] + C2 * X[t]
+    # print(X[t+1])
+  }
+  X[abs(X) < epsilon] <- 0
+  X[bad] <- NA
+  attributes(X) <- inAttr
+  X
+}
+
+#' @export
+simhydrouting.ranges <- function() {
+  list(
+    DELAY = c(0.1, 5),
+    X_m = c(0.01, 0.5)
+  )
+}
+