Merge pull request #38 from ConnectedSystems/txc-caps

Add cap to bias correction offset

Author: ConnectedSystems

README.md

@@ -181,7 +181,7 @@ Area: 1985.73
 │            f │   2.57928 │        0.01 │         3.0 │                             Plant stress threshold, controls at what moisture deficit plants begin to experience stress. │
 │            a │   5.92338 │         0.1 │        10.0 │                                    Quickflow storage coefficient, where higher values lead to faster quickflow response. │
 │            b │ 0.0989926 │       0.001 │         0.1 │                                            Slowflow storage coefficient, lower values lead to slower baseflow recession. │
-│ storage_coef │   1.86134 │     1.0e-10 │        10.0 │                               Groundwater interaction factor, controling how water is exchanged with deeper groundwater. │
+│ storage_coef │   1.86134 │     1.0e-10 │        10.0 │                              Groundwater interaction factor, controlling how water is exchanged with deeper groundwater. │
 │        alpha │  0.727905 │      1.0e-5 │         1.0 │                                                      Effective rainfall scaling factor, partitions rainfall into runoff. │
 └──────────────┴───────────┴─────────────┴─────────────┴──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
 ```

docs/src/assets/ensemble_bias_corrected.png

@@ -6,9 +6,9 @@ as the ensemble constituents. The default ensemble is a normalized weighted sum.
 The usual setup process is shown here, detailed in previous sections of this guide.
 
 ```julia
-using StatsPlots
-using CSV, DataFrames
 using Statistics
+using CSV, DataFrames
+using StatsPlots
 using Streamfall
 
 data_dir = joinpath(
@@ -86,7 +86,7 @@ low flows as with GR4J.
 
 ![](../../assets/ensemble_model_comparison_quickplots.png)
 
-Comparing the temporal cross section:
+Comparing the temporal cross section to get an idea of seasonality:
 
 ```julia
 ihacres_xs = temporal_cross_section(burn_dates, burn_obs, ihacres_node.outflow[burn_in:end]; title="IHACRES", yscale=:log10)
@@ -101,10 +101,11 @@ A reduction in the median error can be seen with extreme errors reduced somewhat
 ![](../../assets/ensemble_xsection.png)
 
 The median error can then be applied to modelled streamflow (on a month-day basis) as a
-form of bias correction.
+form of bias correction. Here, the correction factor is capped to -80% and +40% of predicted
+outflows.
 
 ```julia
-q_star = Streamfall.apply_temporal_correction(ensemble, climate, Qo[:, "410730"])
+q_star = Streamfall.apply_temporal_correction(ensemble, climate, Qo[:, "410730"]; low_cap=0.8, high_cap=0.4)
 
 bc_ensemble_qp = quickplot(burn_obs, q_star[burn_in:end], climate; label="Bias Corrected Ensemble", log=true)
 
@@ -116,16 +117,17 @@ bias_corrected_xs = temporal_cross_section(
     yscale=:log10
 )
 
-plot(bc_ensemble_qp, bias_corrected_xs; layout=(2,1), size=(800, 800))
+ens_qp = plot(bc_ensemble_qp, bias_corrected_xs; layout=(2,1), size=(800, 800))
 ```
 
-While the median error has increased, its variance has reduced significantly. At the same
-time, performance at the 75 and 95% CI remain steady relative to the original weighted
-ensemble results.
+It can be seen here that low flows are better represented, with a commensurate decrease
+in median error (and its variance). At the same time, performance at the 75 and 95% CI 
+remain steady relative to the original weighted ensemble results.
 
 ![](../../assets/ensemble_bias_corrected.png)
 
 This ensemble approach may be improved further by:
 
 - Using a rolling window to smooth ensemble predictions
 - Defining a custom objective function to target specific conditions
+- Using more advanced ensemble approaches other than the simple weighted mean approach

docs/src/assets/ensemble_model_comparison_quickplots.png

@@ -85,7 +85,7 @@ Area: 130.0
 │            f │   0.8 │        0.01 │         3.0 │                             Plant stress threshold, controls at what moisture deficit plants begin to experience stress. │
 │            a │   0.9 │         0.1 │        10.0 │                                    Quickflow storage coefficient, where higher values lead to faster quickflow response. │
 │            b │   0.1 │       0.001 │         0.1 │                                            Slowflow storage coefficient, lower values lead to slower baseflow recession. │
-│ storage_coef │   2.9 │     1.0e-10 │        10.0 │                               Groundwater interaction factor, controling how water is exchanged with deeper groundwater. │
+│ storage_coef │   2.9 │     1.0e-10 │        10.0 │                              Groundwater interaction factor, controlling how water is exchanged with deeper groundwater. │
 │        alpha │  0.95 │      1.0e-5 │         1.0 │                                                      Effective rainfall scaling factor, partitions rainfall into runoff. │
 └──────────────┴───────┴─────────────┴─────────────┴──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
 ```

docs/src/assets/ensemble_xsection.png

@@ -1,4 +1,4 @@
-import ..Analysis: TemporalCrossSection
+import ..Analysis: TemporalCrossSection, offsets
 
 Base.@kwdef mutable struct WeightedEnsembleNode{N<:NetworkNode,P,A<:Real} <: EnsembleNode
     name::String
@@ -177,20 +177,31 @@ end
     apply_temporal_correction(
         ensemble::WeightedEnsembleNode,
         climate::Climate,
-        obs::Vector{T}
+        obs::Vector{T};
+        low_cap::Float64=1.0,
+        high_cap::Float64=1.0
     ) where {T<:Real}
 
 Correct for model bias using median error.
+By default, the correction is capped to ± 100% of model prediction modified by
+`low_cap` and `high_cap`
 """
 function apply_temporal_correction(
     ensemble::WeightedEnsembleNode,
     climate::Climate,
-    obs::Vector{T}
-) where {T<:Real}
+    obs::Vector{Float64};
+    low_cap::Float64=1.0,
+    high_cap::Float64=1.0
+)
+    low_cap = clamp(low_cap, 0.0, 1.0)
+    high_cap = clamp(high_cap, 0.0, 1.0)
+
     dates = timesteps(climate)
     tcs = TemporalCrossSection(dates, obs, ensemble.outflow)
 
-    return max.(ensemble.outflow .+ Streamfall.Analysis.offsets(tcs), 0.0)
+    x = min.(ensemble.outflow .+ offsets(tcs), ensemble.outflow .* (1.0 + high_cap))
+
+    return max.(x, ensemble.outflow .* (1.0 - low_cap))
 end
 
 """

docs/src/examples/ensembles/weighted_ensembles.md

@@ -20,7 +20,7 @@ Base.@kwdef mutable struct IHACRESBilinearNode{P,A<:AbstractFloat} <: IHACRESNod
     a::P = Param(0.9, bounds=(0.1, 10.0), desc="Quickflow storage coefficient, where higher values lead to faster quickflow response.")  # quickflow storage coefficient == (1/tau_q)
     b::P = Param(0.1, bounds=(1e-3, 0.1), desc="Slowflow storage coefficient, lower values lead to slower baseflow recession.")  # slowflow storage coefficent == (1/tau_s)
 
-    storage_coef::P = Param(2.9, bounds=(1e-10, 10.0), desc="Groundwater interaction factor, controling how water is exchanged with deeper groundwater.")
+    storage_coef::P = Param(2.9, bounds=(1e-10, 10.0), desc="Groundwater interaction factor, controlling how water is exchanged with deeper groundwater.")
     alpha::P = Param(0.1, bounds=(1e-5, 1 - 1 / 10^9), desc="Effective rainfall scaling factor, partitions rainfall into runoff.")
 
     # const level_params::Array{P, 1} = [