|
| 1 | +Streamfall supports ensemble modeling using individual instances of rainfall-runoff models |
| 2 | +as the ensemble constituents. The default ensemble is a normalized weighted sum. |
| 3 | + |
| 4 | +The usual setup process is shown here, detailed in previous sections of this guide. |
| 5 | + |
| 6 | +```julia |
| 7 | +using YAML, DataFrames, CSV, Plots |
| 8 | +using Statistics |
| 9 | +using Streamfall |
| 10 | + |
| 11 | + |
| 12 | +climate = Climate("../test/data/campaspe/climate/climate.csv", "_rain", "_evap") |
| 13 | + |
| 14 | +# Historic flows and dam level data |
| 15 | +obs_data = CSV.read( |
| 16 | + "../test/data/cotter/climate/CAMELS-AUS_410730.csv", |
| 17 | + DataFrame; |
| 18 | + comment="#" |
| 19 | +) |
| 20 | + |
| 21 | +Qo = extract_flow(obs_data, "410730") |
| 22 | +climate = extract_climate(obs_data) |
| 23 | +``` |
| 24 | + |
| 25 | +Specific node representations can then be created, each representing the same sub-catchment. |
| 26 | +A stream network is not considered for this demonstration. |
| 27 | + |
| 28 | +```julia |
| 29 | + |
| 30 | +# Create one instance each of IHACRES_CMD and GR4J |
| 31 | +ihacres_node = create_node(IHACRESBilinearNode, "410730", 129.2) |
| 32 | +gr4j_node = create_node(GR4JNode, "410730", 129.2) |
| 33 | + |
| 34 | +# Create a weighted ensemble with equal weights |
| 35 | +# The default behavior is to combine component predictions with a normalized weighted sum. |
| 36 | +ensemble = create_node(WeightedEnsembleNode, [ihacres_node, gr4j_node], [0.5, 0.5]) |
| 37 | +``` |
| 38 | + |
| 39 | +In previous examples, the `calibrate!()` method was used to calibrate nodes to available |
| 40 | +data. For WeightedEnsembleNodes, `calibrate!()` optimizes just the weights to allow for |
| 41 | +pre-calibrated instances to be provided for use in an ensemble. This may be useful if |
| 42 | +it is desired for the component models to be calibrated according to different criteria and |
| 43 | +objective functions. |
| 44 | + |
| 45 | +A separate `calibrate_instances!()` method is available to calibrate individual component |
| 46 | +models (and the weights used). |
| 47 | + |
| 48 | +```julia |
| 49 | +# The default behaviour for WeightedEnsembleNodes are to calibrate just the weights. |
| 50 | +# res, opt = calibrate!(ensemble, climate, Qo, (obs, sim) -> 1.0 .- Streamfall.NmKGE(obs, sim); MaxTime=180) |
| 51 | + |
| 52 | +# Here, the component models are uncalibrated so we calibrate these and the weights. |
| 53 | +res, opt = calibrate_instances!(ensemble, climate, Qo, (obs, sim) -> 1.0 .- Streamfall.NmKGE(obs, sim); MaxTime=180) |
| 54 | +``` |
| 55 | + |
| 56 | +Running the calibrated models directly, the amount of improvement to model performance can |
| 57 | +be assessed. Here, a 1-year burn in period is used. |
| 58 | + |
| 59 | + |
| 60 | +```julia |
| 61 | +run_node!(ihacres_node, climate) |
| 62 | +run_node!(gr4j_node, climate) |
| 63 | +run_node!(ensemble, climate) |
| 64 | + |
| 65 | +burn_in = 365 |
| 66 | +burn_dates = timesteps(climate)[burn_in:end] |
| 67 | +burn_obs = Qo[burn_in:end, "410730"] |
| 68 | + |
| 69 | +ihacres_qp = quickplot(burn_obs, ihacres_node.outflow[burn_in:end], climate, "IHACRES", true) |
| 70 | +gr4j_qp = quickplot(burn_obs, gr4j_node.outflow[burn_in:end], climate, "GR4J", true) |
| 71 | +ensemble_qp = quickplot(burn_obs, ensemble.outflow[burn_in:end], climate, "Weighted Ensemble", true) |
| 72 | + |
| 73 | +plot(ihacres_qp, gr4j_qp, ensemble_qp; layout=(3, 1), size=(800, 1200)) |
| 74 | +``` |
| 75 | + |
| 76 | +Below a small improvement to model performance based on the modified Kling-Gupta Efficiency |
| 77 | +score can be seen. Comparing the Q-Q plots, IHACRES had a tendency to underestimate low |
| 78 | +flows and high flows, whereas GR4J had a tendency to overestimate low flows. |
| 79 | + |
| 80 | +The weighted ensemble combined characteristics of both, with a tendency to overestimate |
| 81 | +low flows as with GR4J. |
| 82 | + |
| 83 | + |
| 84 | + |
| 85 | +Comparing the temporal cross section: |
| 86 | + |
| 87 | +```julia |
| 88 | +ihacres_xs = temporal_cross_section(burn_dates, burn_obs, ihacres_node.outflow[burn_in:end]; title="IHACRES") |
| 89 | +gr4j_xs = temporal_cross_section(burn_dates, burn_obs, gr4j_node.outflow[burn_in:end]; title="GR4J") |
| 90 | +ensemble_xs = temporal_cross_section(burn_dates, burn_obs, ensemble.outflow[burn_in:end]; title="Weighted Ensemble (IHACRES-GR4J)") |
| 91 | + |
| 92 | +plot(ihacres_xs, gr4j_xs, ensemble_xs; layout=(3, 1), size=(800, 1200)) |
| 93 | +``` |
| 94 | + |
| 95 | +A reduction in the median error can be seen with extreme errors reduced somewhat (according to |
| 96 | +the 95% CI). |
| 97 | + |
| 98 | + |
| 99 | + |
| 100 | +The median error can then be applied to modelled streamflow (on a month-day basis) as a |
| 101 | +form of bias correction. |
| 102 | + |
| 103 | +```julia |
| 104 | +q_star = Streamfall.apply_temporal_correction(ensemble, climate, Qo[:, "410730"]) |
| 105 | + |
| 106 | +bc_ensemble_qp = quickplot(burn_obs, q_star[burn_in:end], climate, "Bias Corrected Ensemble") |
| 107 | + |
| 108 | +bias_corrected_xs = temporal_cross_section( |
| 109 | + burn_dates, |
| 110 | + burn_obs, |
| 111 | + q_star[burn_in:end]; |
| 112 | + title="Bias Corrected Ensemble" |
| 113 | +) |
| 114 | + |
| 115 | +plot(bc_ensemble_qp, bias_corrected_xs; layout=(2,1), size=(800, 800)) |
| 116 | +``` |
| 117 | + |
| 118 | +While the median error has increased, its variance has reduced significantly. At the same |
| 119 | +time, performance at the 75 and 95% CI remain steady relative to the original weighted |
| 120 | +ensemble results. |
| 121 | + |
| 122 | + |
0 commit comments