Skip to content

Commit de7261f

Browse files
ConnectedSystemsConnectedSystems
committed
Updated calibration examples
1 parent 57df9aa commit de7261f

2 files changed

Lines changed: 126 additions & 105 deletions

File tree

‎examples/calibration/alt_calibration.jl‎

Lines changed: 9 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -32,12 +32,15 @@ function calibrate(sn, v_id, climate, calib_data)
3232

3333
# Get node parameters (default values and bounds)
3434
p_names, x0, param_bounds = param_info(this_node; with_level=false)
35-
opt = bbsetup(opt_func; SearchRange=param_bounds,
36-
Method=:adaptive_de_rand_1_bin_radiuslimited,
37-
MaxTime=2400.0, # time in seconds to spend
38-
TraceInterval=30.0,
39-
PopulationSize=75,
40-
)
35+
opt = bbsetup(
36+
opt_func;
37+
parameters=x0,
38+
SearchRange=param_bounds,
39+
Method=:adaptive_de_rand_1_bin_radiuslimited,
40+
MaxTime=2400.0, # time in seconds to spend
41+
TraceInterval=30.0,
42+
PopulationSize=75,
43+
)
4144

4245
res = bboptimize(opt)
4346

‎examples/calibration/calibration.jl‎

Lines changed: 117 additions & 99 deletions
Original file line numberDiff line numberDiff line change
@@ -1,107 +1,125 @@
1-
# Import common packages and functions
2-
include("_obj_func_definition.jl")
3-
4-
5-
"""Example calibration function.
6-
7-
Illustrate model calibration using the BlackBoxOptim package.
81
"""
9-
function calibrate(sn, v_id, climate, calib_data)
10-
11-
# Fitness of model is dependent on next node.
12-
ins = inlets(sn, v_id)
13-
next_node_id = outlets(sn, v_id)[1]
14-
15-
# Recurse through and calibrate all nodes upstream
16-
if !isempty(ins)
17-
for nid in ins
18-
calibrate(sn, nid, climate, calib_data)
19-
end
20-
end
21-
22-
this_node = sn[v_id]
23-
24-
# Create new optimization function (see definition inside `_obj_func_definition.jl`)
25-
opt_func = x -> obj_func(x, climate, sn, v_id, next_node_id, calib_data)
26-
27-
# Get node parameters (default values and bounds)
28-
p_names, x0, param_bounds = param_info(this_node; with_level=false)
29-
opt = bbsetup(opt_func; SearchRange=param_bounds,
30-
Method=:adaptive_de_rand_1_bin_radiuslimited,
31-
MaxTime=2400.0, # time in seconds to spend
32-
TraceInterval=30.0,
33-
PopulationSize=75,
34-
)
35-
36-
res = bboptimize(opt)
37-
38-
bs = best_candidate(res)
39-
@info "Calibrated $(v_id) ($(this_node.name)), with score: $(best_fitness(res))"
40-
@info "Best Params:" collect(bs)
41-
42-
# Update node with calibrated parameters
43-
update_params!(this_node, bs...)
44-
45-
return res, opt
46-
end
47-
48-
49-
v_id, node = sn["406219"]
50-
@info "Starting calibration..."
51-
res, opt = calibrate(sn, v_id, climate, hist_data)
2+
Example showcasing calibrating and running a streamflow network.
523
53-
# Stream
54-
55-
best_params = best_candidate(res)
4+
Data is prepped with the script `campaspe_data_prep.jl` in the `test/data/campaspe`
5+
directory.
6+
"""
567

57-
@info best_fitness(res)
58-
@info best_params
8+
using OrderedCollections
9+
using Glob
5910

11+
using Statistics
12+
using CSV, DataFrames, YAML
13+
using Streamfall
6014

6115
using Plots
6216

63-
update_params!(node, best_params...)
64-
dam_id, dam_node = sn["406000"]
65-
Streamfall.run_node!(sn, dam_id, climate; extraction=hist_dam_releases)
66-
67-
h_data = hist_data["406000"]
68-
n_data = dam_node.level
69-
70-
nnse_score = Streamfall.NNSE(h_data, n_data)
71-
nse_score = Streamfall.NSE(h_data, n_data)
72-
rmse_score = Streamfall.RMSE(h_data, n_data)
73-
74-
@info "Downstream Dam Level NNSE:" nnse_score
75-
@info "Downstream Dam Level RMSE:" rmse_score
76-
77-
reset!(dam_node)
78-
79-
nse = round(nse_score, digits=4)
80-
rmse = round(rmse_score, digits=4)
81-
82-
plot(h_data,
83-
legend=:bottomleft,
84-
title="Calibrated IHACRES\n(NSE: $(nse); RMSE: $(rmse))",
85-
label="Historic", xlabel="Day", ylabel="Dam Level [mAHD]")
86-
87-
plot!(n_data, label="IHACRES")
88-
89-
savefig("calibration_ts_comparison.png")
90-
91-
# 1:1 Plot
92-
scatter(h_data, n_data, legend=false,
93-
markerstrokewidth=0, markerstrokealpha=0, alpha=0.2)
94-
plot!(h_data, h_data, color=:red, markersize=.1, markerstrokewidth=0,
95-
xlabel="Historic [mAHD]", ylabel="IHACRES [mAHD]", title="Historic vs Modelled")
96-
97-
savefig("calibration_1to1.png")
9817

99-
# NNSE: 0.9643; RMSE: 1.43553
100-
# d: 84.28015146853407
101-
# d2: 2.4224106535469145
102-
# e: 0.8129590022893607
103-
# f: 2.579276454391652
104-
# a: 5.923379062122229
105-
# b: 0.0989925603647026
106-
# storage_coef: 1.8613364808233752 # gw storage factor
107-
# alpha: 0.7279050097363565
18+
sn = load_network("Example Network", "../test/data/campaspe/campaspe_network.yml")
19+
20+
# The Campaspe catchment is represented as a network of eight nodes, including one dam.
21+
# All nodes use the IHACRES_CMD rainfall-runoff model.
22+
plot_network(sn)
23+
24+
# Load climate data - in this case from a CSV file with data for all nodes.
25+
climate_data = CSV.read(
26+
"../test/data/campaspe/climate/climate.csv",
27+
DataFrame;
28+
comment="#"
29+
)
30+
31+
# Indicate which columns are precipitation and evaporation data based on partial identifiers
32+
climate = Climate(climate_data, "_rain", "_evap")
33+
34+
# Historic flows and dam level data
35+
calib_data = CSV.read(
36+
"../test/data/campaspe/gauges/outflow_and_level.csv",
37+
DataFrame;
38+
comment="#"
39+
)
40+
41+
# Historic extractions from the dam
42+
extraction_data = CSV.read("gauges/dam_extraction.csv", DataFrame; comment="#")
43+
44+
# We now have a dataset for calibration (`calib_data`) and a dataset indicating the
45+
# historic dam extractions (`extraction_data`).
46+
# `extraction_data` may also hold water extractions at each "reach".
47+
48+
# Provide a metric to use to fit models against data.
49+
# Note that calibration always assumes minimization, so if the metric does not
50+
# provide this directionality, it must be wrapped to do so.
51+
metric = (y, y_hat) -> 1.0 - Streamfall.NNSE(y, y_hat)
52+
# metric = (y, y_hat) -> 1.0 - Streamfall.NmKGE(y, y_hat)
53+
# metric = (y, y_hat) -> 1.0 - Streamfall.naive_split_metric(
54+
# y, y_hat; n_members=7, metric=Streamfall.NmKGE, comb_method=mean
55+
# )
56+
# metric = Streamfall.RMSE
57+
58+
# Alternatively, individual metrics for each node in a dictionary
59+
# (key-value pairs in the form of name => function).
60+
# Here, Normalized KGE′ is used for all nodes.
61+
# metrics = Dict{String,Function}(
62+
# n.name => (y, y_hat) -> 1.0 - Streamfall.NmKGE(y, y_hat) for n in sn
63+
# )
64+
65+
# Calibrate all gauges in network using Adaptive Differential Evolution with the
66+
# BlackBoxOptim.jl package. Any BlackBoxOptim keyword arguments are passed through.
67+
# The parameter values provided in the network specification will be used as the initial
68+
# guess.
69+
# The default is to spend 5 mins on each node (MaxTime=300), but for this example we run
70+
# calibration for 1mins/node (MaxTime=60).
71+
# If the downstream node represents a dam, the current node is calibrated by fitting the
72+
# outflows such that it reproduces the observed dam levels.
73+
# The `weighting` parameter controls the weighting between nodes for calibration.
74+
# A choice can be made to calibrate against outflows (a weighting of 1) or dam levels
75+
# (a weighting of 0).
76+
# Here, we calibrate to downstream dam levels only (a zero weighting on node outflows)
77+
calibrate!(
78+
sn, climate, calib_data, metric;
79+
extraction=extraction_data, weighting=0.0,
80+
MaxTime=60.0
81+
);
82+
83+
# Could calibrate a specific node, assuming all nodes upstream have already been calibrated
84+
# Set `calibrate_all=true` to calibrate all upstream nodes as well.
85+
# To produce the results shown below, the node upstream from the dam was calibrated an
86+
# additional 2 hours.
87+
# calibrate!(
88+
# sn, 2, climate, calib_data, metric;
89+
# extraction=extraction_data, weighting=0.0, calibrate_all=false,
90+
# MaxTime=7200.0
91+
# );
92+
93+
# Run all nodes in the catchment
94+
run_catchment!(sn, climate; extraction=extraction_data)
95+
96+
# Get performance metrics for dam levels
97+
dam_obs = aligned_dam_levels[:, "Dam Level [mAHD]"]
98+
dam_node = sn[3]
99+
dam_sim = dam_node.level
100+
101+
Streamfall.RMSE(dam_obs[366:end], dam_sim[366:end])
102+
Streamfall.NSE(dam_obs[366:end], dam_sim[366:end])
103+
Streamfall.mKGE(dam_obs[366:end], dam_sim[366:end])
104+
105+
# Plot results (using a 1-year burn-in period)
106+
f = quickplot(dam_obs, dam_sim, climate, "Modelled - 406000", false; burn_in=366)
107+
savefig(f, "example_dam_level.png")
108+
109+
# Save calibrated network to a file
110+
save_network(sn, "example_network_calibrated.yml")
111+
112+
# Illustrating that the re-loaded network reproduces the results as above
113+
sn2 = load_network("Calibrated Example", "example_network_calibrated.yml")
114+
run_catchment!(sn2, climate; extraction=extraction_data)
115+
116+
dam_node = sn2[3]
117+
dam_sim = dam_node.level
118+
119+
rmse_score = Streamfall.RMSE(dam_obs[366:end], dam_sim[366:end])
120+
nse_score = Streamfall.NSE(dam_obs[366:end], dam_sim[366:end])
121+
mKGE_score = Streamfall.mKGE(dam_obs[366:end], dam_sim[366:end])
122+
123+
@info "Scores: " rmse_score nse_score mKGE_score
124+
125+
f2 = quickplot(dam_obs, dam_sim, climate, "Modelled - 406000", false; burn_in=366)

0 commit comments

Comments
 (0)