v1.4.2.0 NEW: added calculate_Imax fun to computethe maximum capacity of the...

v1.4.2.0 NEW: added calculate_Imax fun to computethe maximum capacity of the GR5H interception store #13

DESCRIPTION

 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.4.1.23
+Version: 1.4.2.0
 Date: 2019-12-10
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),

NAMESPACE

@@ -28,6 +28,7 @@ export(ErrorCrit_KGE)
 export(ErrorCrit_KGE2)
 export(ErrorCrit_NSE)
 export(ErrorCrit_RMSE)
+export(calculate_Imax)
 export(PE_Oudin)
 export(PEdaily_Oudin)
 export(RunModel)

NEWS.md

@@ -2,12 +2,14 @@
 
 
 
-### 1.4.1.23 Release Notes (2019-12-10)
+### 1.4.2.0 Release Notes (2019-12-10)
 
 
 #### New features
 
 - <code>plot.Outputsmodel()</code> now allows to draw actual evapotranspiration when <code>which = "ActuEvap"</code> or <code>which = "All"</code> (overlaid to potential evapotranspiration if already drawn).
+- Added <code>RunModel_GR5H()</code> function to run the hourly model GR5H. This model prestents an interception store
+- Added <code>calculate_Imax()</code> which allows to estimate the maximum capacity of the GR5H interception store.
 
 
 #### Bug fixes

R/calculate_Imax.R

+
+calculate_Imax <- function(InputsModel, 
+                           IndPeriod_Run, 
+                           tested_values = seq(0.1, 3, 0.1)) {
+  
+  ##_____Arguments_check_____________________________________________________________________
+  if (!inherits(InputsModel, "InputsModel")) {
+    stop("'InputsModel' must be of class 'InputsModel'")
+  }  
+  if (!inherits(InputsModel, "hourly")) {
+    stop("'InputsModel' must be of class 'hourly'")
+  }  
+  if (!(class(tested_values) == "numeric")){
+    stop("'tested_values' must be 'numeric'")
+  }
+  
+  ##check_IndPeriod_Run
+  if (!is.vector(IndPeriod_Run)) {
+    stop("'IndPeriod_Run' must be a vector of numeric values")
+  }
+  if (!is.numeric(IndPeriod_Run)) {
+    stop("'IndPeriod_Run' must be a vector of numeric values")
+  }
+  if (!identical(as.integer(IndPeriod_Run), as.integer(seq(from = IndPeriod_Run[1], to = tail(IndPeriod_Run, 1), by = 1)))) {
+    stop("'IndPeriod_Run' must be a continuous sequence of integers")
+  }
+  if (storage.mode(IndPeriod_Run) != "integer") {
+    stop("'IndPeriod_Run' should be of type integer")
+  }
+  
+  ##aggregate data at the daily time step  
+  TabSeries <- data.frame(DatesR = InputsModel$DatesR[IndPeriod_Run], 
+                          Precip = InputsModel$Precip[IndPeriod_Run], 
+                          PotEvap = InputsModel$PotEvap[IndPeriod_Run])
+  daily_data <- SeriesAggreg(TabSeries, "hourly", "daily", 
+                             ConvertFun = c("sum", "sum"))
+  
+  ##calculate total interception of daily GR models on the period
+  cum_daily <- sum(pmin(daily_data$Precip, daily_data$PotEvap))
+  
+  ##calculate total interception of the GR5H interception store on the period
+  ## and compute difference with daily values
+  differences <- array(NA, c(length(tested_values)))
+  for (Imax in tested_values) {
+    C0 <- 0
+    cum_hourly <- 0
+    for (i in IndPeriod_Run) {
+      
+      Ec  <- min(InputsModel$PotEvap[i], InputsModel$Precip[i]+C0)
+      Pth <- max(0, InputsModel$Precip[i]-(Imax-C0)-Ec)
+      C0  <- C0 + InputsModel$Precip[i] - Ec - Pth
+      cum_hourly <- cum_hourly + Ec
+      
+    }
+    differences[which(Imax == tested_values)] <- abs(cum_hourly - cum_daily)
+  }
+  
+  ##return the Imax value that minimises the difference
+  return(tested_values[which.min(differences)])
+  
+}