Appendix C: Calibration code#
In this appendix is the code used to calibrate the model.
Startup#
# General python
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
import numpy as np
from pathlib import Path
import pandas as pd
import matplotlib.pyplot as plt
import xarray as xr
# Niceties
from rich import print
from ipywidgets import IntProgress
from IPython.display import display
# General eWaterCycle
import ewatercycle
import ewatercycle.models
import ewatercycle.forcing
Defining variables & pathways#
# Choosing time period & basin si
experiment_start_date = "2000-01-01T00:00:00Z"
experiment_end_date = "2005-12-31T00:00:00Z"
calibration_start = "1999-01-01T00:00:00Z"
calibration_end = "2019-12-31T00:00:00Z"
# validation_start = ""
# validation_end = ""
# Create pathways for ERA 5 forcings
forcing_path_ERA5 = Path.home() / "BEP-maxime" / "book" / "thesis_projects" / "BSc" / "2026_Q4_MaximedeBekker_CEG" / "Workyard" / "forcings" / "ERA5" / "ERA5-1999-2019"
discharge_file = Path.home() / "BEP-maxime" / "book" / "thesis_projects" / "BSc" / "2026_Q4_MaximedeBekker_CEG" / "Workyard" / "07DA001_discharge_daily_withoutmissing.csv"
shape_file = Path.home() / "BEP-maxime" / "book" / "thesis_projects" / "BSc" / "2026_Q4_MaximedeBekker_CEG" / "Workyard" / "Shapefiles" / "07DA001_basin.shp"
basin_size = 132572
hbv_config = Path.home() / "BEP-maxime" / "book" / "thesis_projects" / "BSc" / "2026_Q4_MaximedeBekker_CEG" / "Workyard" / "hbv_config"
hbv_config.mkdir(parents=True, exist_ok=True)
Load & prepare discharge data#
# Load CSV discharge 07DA001
q_obs = pd.read_csv(discharge_file, skiprows=1)
q_obs = q_obs[["Date", "Value"]].copy()
q_obs["Date"] = pd.to_datetime(q_obs["Date"])
q_obs = q_obs.rename(columns={"Value": "discharge_m3s"})
# Filter q_obs to start & end date
# start_date = pd.to_datetime(experiment_start_date.replace("Z", ""))
# end_date = pd.to_datetime(experiment_end_date.replace("Z", ""))
# q_obs = q_obs[(q_obs["Date"] >= start_date) & (q_obs["Date"] <= end_date)]
# observed_output = pd.Series(data=q_obs["discharge_m3s"].to_numpy(), name="Observed discharge", index=q_obs["Date"])
Generate or Load ERA5 data#
# Generate ERA5 data
# ERA5_forcing = ewatercycle.forcing.sources['LumpedMakkinkForcing'].generate(
# dataset="ERA5",
# start_time=experiment_start_date,
# end_time=experiment_end_date,
# shape=shape_file,
# )
# Load ERA5 data
ERA5_forcing = ewatercycle.forcing.sources["LumpedMakkinkForcing"].load(directory=forcing_path_ERA5)
print(ERA5_forcing)
LumpedMakkinkForcing( start_time='1999-01-01T00:00:00Z', end_time='2019-12-31T00:00:00Z', directory=PosixPath('/home/maxime/BEP-maxime/Workyard/forcings/ERA5-1999-2019'), shape=PosixPath('/home/maxime/BEP-maxime/Workyard/Shapefiles/07DA001_basin.shp'), filenames={ 'evspsblpot': 'combined_ERA5_1999_2019_evspsblpot.nc', 'pr': 'combined_ERA5_1999_2019_pr.nc', 'rsds': 'combined_ERA5_1999_2019_rsds.nc', 'tas': 'combined_ERA5_1999_2019_tas.nc' } )
Calibration Methods#
Method 1: NSE#
def NSE(sim, obs):
sim = sim.to_numpy()
obs = obs.to_numpy()
return 1 - np.sum((obs - sim) ** 2) / np.sum((obs - np.mean(obs)) ** 2)
Method 2: Log_NSE#
def log_NSE(sim, obs):
# Convert for easy calculations
sim = sim.to_numpy()
obs = obs.to_numpy()
# Avoid log(0)
eps = 0.00000000001
log_sim = np.log(sim + eps)
log_obs = np.log(obs + eps)
return 1 - np.sum((log_obs - log_sim) ** 2) / np.sum((log_obs - np.mean(log_obs)) ** 2)
Define Parameters & Storages and Times#
# Define initial storages
# Si, Su, Sf, Ss, Sp
s_0 = np.array([0, 100, 0, 5, 0])
# Define parameters (range)
parameter_ranges = {
"Imax": (5.5, 8), # Maximum interception storage
"Ce": (0.35, 0.6), # Evaporation correction factor
"Sumax": (150, 220), # Maximum soil moisture storage
"Beta": (1.6, 2), # Soil runoff parameter
"Pmax": (0.1, 0.6), # Maximum percolation rate
"Tlag": (5.2, 7), # Time lag
"Kf": (0.01, 0.1), # Fast reservoir recession coefficient
"Ks": (0.001, 0.1), # Slow reservoir recession coefficient
"FM": (0.3, 1.4), # Snowmelt factor
}
# Convert to arrays
parameter_names = list(parameter_ranges.keys())
p_min = np.array([parameter_ranges[name][0] for name in parameter_names])
p_max = np.array([parameter_ranges[name][1] for name in parameter_names])
Create samples#
### Create samples
N = 10
par_samples = np.random.rand(N, len(parameter_ranges)) # Generates random values between 0-1
parameter_sets = p_min + par_samples * (p_max - p_min) # Fits value to parameter range
print(parameter_sets[:][0])
[7.85285914e+00 3.89990216e-01 1.76127959e+02 1.92058774e+00 5.56024109e-01 6.99810960e+00 3.01314564e-02 6.63774925e-02 9.62585104e-01]
Define time period#
# Define time period
calibration_start_date = pd.to_datetime(calibration_start.replace("Z", ""))
calibration_end_date = pd.to_datetime(calibration_end.replace("Z", ""))
# Skip 1 year for calibration
evaluation_start = pd.to_datetime(f"{calibration_start_date.year + 1}-01-01")
# Align q_obs
q_obs = q_obs[(q_obs["Date"] >= calibration_start_date) & (q_obs["Date"] <= calibration_end_date)]
observed_output = pd.Series(data=q_obs["discharge_m3s"].to_numpy(), name="Observed discharge", index=q_obs["Date"])
Calibration start function#
# Create function to run the model
def run_hbv(parameters, initial_storages, forcing):
# Create model object
model = ewatercycle.models.HBV(forcing=forcing)
# Create config file
config_file, _ = model.setup(
parameters=parameters,
initial_storages=initial_storages,
cfg_dir=hbv_config)
# Initialising model
model.initialize(config_file)
# Define & update outputs
Q_m = []
time = []
while model.time < model.end_time:
model.update()
Q_m.append(model.get_value("Q")[0])
time.append(pd.Timestamp(model.time_as_datetime))
model.finalize()
# Convert mm/day to m3/s
model_output_mmday = pd.Series(
data=Q_m,
index=time,
name="Modelled discharge")
model_output_m3s = model_output_mmday * basin_size * 1000 / 86400
return model_output_m3s
Run HBV loop#
results = []
# Progress bar for visualization
f = IntProgress(min=0, max=N)
display(f)
for i in range(N):
print(f"Running parameter set {i+1}/{N}")
simulated = run_hbv(
parameters=parameter_sets[i],
initial_storages=s_0,
forcing=ERA5_forcing)
# Filter data by day only, not by day & time to prevent alignment issues
simulated_daily = simulated
observed_daily = observed_output
simulated_daily.index = pd.to_datetime(simulated_daily.index).tz_localize(None).normalize()
observed_daily.index = pd.to_datetime(observed_daily.index).tz_localize(None).normalize()
# Align simulated and observed data
combined_data = pd.DataFrame(
{"Modelled discharge": simulated_daily,
"Observed discharge": observed_daily}).dropna()
# Skip first year as warm-up
combined_data = combined_data[
(combined_data.index >= evaluation_start) &
(combined_data.index <= calibration_end_date)]
# Exclude ice affected winter period
combined_data = combined_data[
~combined_data.index.month.isin([12, 1, 2, 3, 4])]
# Append results and drop nan
if combined_data.empty:
print(f"No overlapping data for parameter set {i+1}")
nse = np.nan
log_nse = np.nan
else:
nse = NSE(
combined_data["Modelled discharge"],
combined_data["Observed discharge"])
log_nse = log_NSE(
combined_data["Modelled discharge"],
combined_data["Observed discharge"])
results.append({
"run": i,
"NSE": nse,
"log_NSE": log_nse,
"parameters": parameter_sets[i]})
Running parameter set 1/600
Running parameter set 2/600
Running parameter set 3/600
Running parameter set 4/600
Running parameter set 5/600
Running parameter set 6/600
Running parameter set 7/600
Running parameter set 8/600
Running parameter set 9/600
Running parameter set 10/600
Running parameter set 11/600
Running parameter set 12/600
Running parameter set 13/600
Running parameter set 14/600
Running parameter set 15/600
Running parameter set 16/600
Running parameter set 17/600
Running parameter set 18/600
Running parameter set 19/600
Running parameter set 20/600
Running parameter set 21/600
Running parameter set 22/600
Running parameter set 23/600
Running parameter set 24/600
Running parameter set 25/600
Running parameter set 26/600
Running parameter set 27/600
Running parameter set 28/600
Running parameter set 29/600
Running parameter set 30/600
Running parameter set 31/600
Running parameter set 32/600
Running parameter set 33/600
Running parameter set 34/600
Running parameter set 35/600
Running parameter set 36/600
Running parameter set 37/600
Running parameter set 38/600
Running parameter set 39/600
Running parameter set 40/600
Running parameter set 41/600
Running parameter set 42/600
Running parameter set 43/600
Running parameter set 44/600
Running parameter set 45/600
Running parameter set 46/600
Running parameter set 47/600
Running parameter set 48/600
Running parameter set 49/600
Running parameter set 50/600
Running parameter set 51/600
Running parameter set 52/600
Running parameter set 53/600
Running parameter set 54/600
Running parameter set 55/600
Running parameter set 56/600
Running parameter set 57/600
Running parameter set 58/600
Running parameter set 59/600
Running parameter set 60/600
Running parameter set 61/600
Running parameter set 62/600
Running parameter set 63/600
Running parameter set 64/600
Running parameter set 65/600
Running parameter set 66/600
Running parameter set 67/600
Running parameter set 68/600
Running parameter set 69/600
Running parameter set 70/600
Running parameter set 71/600
Running parameter set 72/600
Running parameter set 73/600
Running parameter set 74/600
Running parameter set 75/600
Running parameter set 76/600
Running parameter set 77/600
Running parameter set 78/600
Running parameter set 79/600
Running parameter set 80/600
Running parameter set 81/600
Running parameter set 82/600
Running parameter set 83/600
Running parameter set 84/600
Running parameter set 85/600
Running parameter set 86/600
Running parameter set 87/600
Running parameter set 88/600
Running parameter set 89/600
Running parameter set 90/600
Running parameter set 91/600
Running parameter set 92/600
Running parameter set 93/600
Running parameter set 94/600
Running parameter set 95/600
Running parameter set 96/600
Running parameter set 97/600
Running parameter set 98/600
Running parameter set 99/600
Running parameter set 100/600
Running parameter set 101/600
Running parameter set 102/600
Running parameter set 103/600
Running parameter set 104/600
Running parameter set 105/600
Running parameter set 106/600
Running parameter set 107/600
Running parameter set 108/600
Running parameter set 109/600
Running parameter set 110/600
Running parameter set 111/600
Running parameter set 112/600
Running parameter set 113/600
Running parameter set 114/600
Running parameter set 115/600
Running parameter set 116/600
Running parameter set 117/600
Running parameter set 118/600
Running parameter set 119/600
Running parameter set 120/600
Running parameter set 121/600
Running parameter set 122/600
Running parameter set 123/600
Running parameter set 124/600
Running parameter set 125/600
Running parameter set 126/600
Running parameter set 127/600
Running parameter set 128/600
Running parameter set 129/600
Running parameter set 130/600
Running parameter set 131/600
Running parameter set 132/600
Running parameter set 133/600
Running parameter set 134/600
Running parameter set 135/600
Running parameter set 136/600
Running parameter set 137/600
Running parameter set 138/600
Running parameter set 139/600
Running parameter set 140/600
Running parameter set 141/600
Running parameter set 142/600
Running parameter set 143/600
Running parameter set 144/600
Running parameter set 145/600
Running parameter set 146/600
Running parameter set 147/600
Running parameter set 148/600
Running parameter set 149/600
Running parameter set 150/600
Running parameter set 151/600
Running parameter set 152/600
Running parameter set 153/600
Running parameter set 154/600
Running parameter set 155/600
Running parameter set 156/600
Running parameter set 157/600
Running parameter set 158/600
Running parameter set 159/600
Running parameter set 160/600
Running parameter set 161/600
Running parameter set 162/600
Running parameter set 163/600
Running parameter set 164/600
Running parameter set 165/600
Running parameter set 166/600
Running parameter set 167/600
Running parameter set 168/600
Running parameter set 169/600
Running parameter set 170/600
Running parameter set 171/600
Running parameter set 172/600
Running parameter set 173/600
Running parameter set 174/600
Running parameter set 175/600
Running parameter set 176/600
Running parameter set 177/600
Running parameter set 178/600
Running parameter set 179/600
Running parameter set 180/600
Running parameter set 181/600
Running parameter set 182/600
Running parameter set 183/600
Running parameter set 184/600
Running parameter set 185/600
Running parameter set 186/600
Running parameter set 187/600
Running parameter set 188/600
Running parameter set 189/600
Running parameter set 190/600
Running parameter set 191/600
Running parameter set 192/600
Running parameter set 193/600
Running parameter set 194/600
Running parameter set 195/600
Running parameter set 196/600
Running parameter set 197/600
Running parameter set 198/600
Running parameter set 199/600
Running parameter set 200/600
Running parameter set 201/600
Running parameter set 202/600
Running parameter set 203/600
Running parameter set 204/600
Running parameter set 205/600
Running parameter set 206/600
Running parameter set 207/600
Running parameter set 208/600
Running parameter set 209/600
Running parameter set 210/600
Running parameter set 211/600
Running parameter set 212/600
Running parameter set 213/600
Running parameter set 214/600
Running parameter set 215/600
Running parameter set 216/600
Running parameter set 217/600
Running parameter set 218/600
Running parameter set 219/600
Running parameter set 220/600
Running parameter set 221/600
Running parameter set 222/600
Running parameter set 223/600
Running parameter set 224/600
Running parameter set 225/600
Running parameter set 226/600
Running parameter set 227/600
Running parameter set 228/600
Running parameter set 229/600
Running parameter set 230/600
Running parameter set 231/600
Running parameter set 232/600
Running parameter set 233/600
Running parameter set 234/600
Running parameter set 235/600
Running parameter set 236/600
Running parameter set 237/600
Running parameter set 238/600
Running parameter set 239/600
Running parameter set 240/600
Running parameter set 241/600
Running parameter set 242/600
Running parameter set 243/600
Running parameter set 244/600
Running parameter set 245/600
Running parameter set 246/600
Running parameter set 247/600
Running parameter set 248/600
Running parameter set 249/600
Running parameter set 250/600
Running parameter set 251/600
Running parameter set 252/600
Running parameter set 253/600
Running parameter set 254/600
Running parameter set 255/600
Running parameter set 256/600
Running parameter set 257/600
Running parameter set 258/600
Running parameter set 259/600
Running parameter set 260/600
Running parameter set 261/600
Running parameter set 262/600
Running parameter set 263/600
Running parameter set 264/600
Running parameter set 265/600
Running parameter set 266/600
Running parameter set 267/600
Running parameter set 268/600
Running parameter set 269/600
Running parameter set 270/600
Running parameter set 271/600
Running parameter set 272/600
Running parameter set 273/600
Running parameter set 274/600
Running parameter set 275/600
Running parameter set 276/600
Running parameter set 277/600
Running parameter set 278/600
Running parameter set 279/600
Running parameter set 280/600
Running parameter set 281/600
Running parameter set 282/600
Running parameter set 283/600
Running parameter set 284/600
Running parameter set 285/600
Running parameter set 286/600
Running parameter set 287/600
Running parameter set 288/600
Running parameter set 289/600
Running parameter set 290/600
Running parameter set 291/600
Running parameter set 292/600
Running parameter set 293/600
Running parameter set 294/600
Running parameter set 295/600
Running parameter set 296/600
Running parameter set 297/600
Running parameter set 298/600
Running parameter set 299/600
Running parameter set 300/600
Running parameter set 301/600
Running parameter set 302/600
Running parameter set 303/600
Running parameter set 304/600
Running parameter set 305/600
Running parameter set 306/600
Running parameter set 307/600
Running parameter set 308/600
Running parameter set 309/600
Running parameter set 310/600
Running parameter set 311/600
Running parameter set 312/600
Running parameter set 313/600
Running parameter set 314/600
Running parameter set 315/600
Running parameter set 316/600
Running parameter set 317/600
Running parameter set 318/600
Running parameter set 319/600
Running parameter set 320/600
Running parameter set 321/600
Running parameter set 322/600
Running parameter set 323/600
Running parameter set 324/600
Running parameter set 325/600
Running parameter set 326/600
Running parameter set 327/600
Running parameter set 328/600
Running parameter set 329/600
Running parameter set 330/600
Running parameter set 331/600
Running parameter set 332/600
Running parameter set 333/600
Running parameter set 334/600
Running parameter set 335/600
Running parameter set 336/600
Running parameter set 337/600
Running parameter set 338/600
Running parameter set 339/600
Running parameter set 340/600
Running parameter set 341/600
Running parameter set 342/600
Running parameter set 343/600
Running parameter set 344/600
Running parameter set 345/600
Running parameter set 346/600
Running parameter set 347/600
Running parameter set 348/600
Running parameter set 349/600
Running parameter set 350/600
Running parameter set 351/600
Running parameter set 352/600
Running parameter set 353/600
Running parameter set 354/600
Running parameter set 355/600
Running parameter set 356/600
Running parameter set 357/600
Running parameter set 358/600
Running parameter set 359/600
Running parameter set 360/600
Running parameter set 361/600
Running parameter set 362/600
Running parameter set 363/600
Running parameter set 364/600
Running parameter set 365/600
Running parameter set 366/600
Running parameter set 367/600
Running parameter set 368/600
Running parameter set 369/600
Running parameter set 370/600
Running parameter set 371/600
Running parameter set 372/600
Running parameter set 373/600
Running parameter set 374/600
Running parameter set 375/600
Running parameter set 376/600
Running parameter set 377/600
Running parameter set 378/600
Running parameter set 379/600
Running parameter set 380/600
Running parameter set 381/600
Running parameter set 382/600
Running parameter set 383/600
Running parameter set 384/600
Running parameter set 385/600
Running parameter set 386/600
Running parameter set 387/600
Running parameter set 388/600
Running parameter set 389/600
Running parameter set 390/600
Running parameter set 391/600
Running parameter set 392/600
Running parameter set 393/600
Running parameter set 394/600
Running parameter set 395/600
Running parameter set 396/600
Running parameter set 397/600
Running parameter set 398/600
Running parameter set 399/600
Running parameter set 400/600
Running parameter set 401/600
Running parameter set 402/600
Running parameter set 403/600
Running parameter set 404/600
Running parameter set 405/600
Running parameter set 406/600
Running parameter set 407/600
Running parameter set 408/600
Running parameter set 409/600
Running parameter set 410/600
Running parameter set 411/600
Running parameter set 412/600
Running parameter set 413/600
Running parameter set 414/600
Running parameter set 415/600
Running parameter set 416/600
Running parameter set 417/600
Running parameter set 418/600
Running parameter set 419/600
Running parameter set 420/600
Running parameter set 421/600
Running parameter set 422/600
Running parameter set 423/600
Running parameter set 424/600
Running parameter set 425/600
Running parameter set 426/600
Running parameter set 427/600
Running parameter set 428/600
Running parameter set 429/600
Running parameter set 430/600
Running parameter set 431/600
Running parameter set 432/600
Running parameter set 433/600
Running parameter set 434/600
Running parameter set 435/600
Running parameter set 436/600
Running parameter set 437/600
Running parameter set 438/600
Running parameter set 439/600
Running parameter set 440/600
Running parameter set 441/600
Running parameter set 442/600
Running parameter set 443/600
Running parameter set 444/600
Running parameter set 445/600
Running parameter set 446/600
Running parameter set 447/600
Running parameter set 448/600
Running parameter set 449/600
Running parameter set 450/600
Running parameter set 451/600
Running parameter set 452/600
Running parameter set 453/600
Running parameter set 454/600
Running parameter set 455/600
Running parameter set 456/600
Running parameter set 457/600
Running parameter set 458/600
Running parameter set 459/600
Running parameter set 460/600
Running parameter set 461/600
Running parameter set 462/600
Running parameter set 463/600
Running parameter set 464/600
Running parameter set 465/600
Running parameter set 466/600
Running parameter set 467/600
Running parameter set 468/600
Running parameter set 469/600
Running parameter set 470/600
Running parameter set 471/600
Running parameter set 472/600
Running parameter set 473/600
Running parameter set 474/600
Running parameter set 475/600
Running parameter set 476/600
Running parameter set 477/600
Running parameter set 478/600
Running parameter set 479/600
Running parameter set 480/600
Running parameter set 481/600
Running parameter set 482/600
Running parameter set 483/600
Running parameter set 484/600
Running parameter set 485/600
Running parameter set 486/600
Running parameter set 487/600
Running parameter set 488/600
Running parameter set 489/600
Running parameter set 490/600
Running parameter set 491/600
Running parameter set 492/600
Running parameter set 493/600
Running parameter set 494/600
Running parameter set 495/600
Running parameter set 496/600
Running parameter set 497/600
Running parameter set 498/600
Running parameter set 499/600
Running parameter set 500/600
Running parameter set 501/600
Running parameter set 502/600
Running parameter set 503/600
Running parameter set 504/600
Running parameter set 505/600
Running parameter set 506/600
Running parameter set 507/600
Running parameter set 508/600
Running parameter set 509/600
Running parameter set 510/600
Running parameter set 511/600
Running parameter set 512/600
Running parameter set 513/600
Running parameter set 514/600
Running parameter set 515/600
Running parameter set 516/600
Running parameter set 517/600
Running parameter set 518/600
Running parameter set 519/600
Running parameter set 520/600
Running parameter set 521/600
Running parameter set 522/600
Running parameter set 523/600
Running parameter set 524/600
Running parameter set 525/600
Running parameter set 526/600
Running parameter set 527/600
Running parameter set 528/600
Running parameter set 529/600
Running parameter set 530/600
Running parameter set 531/600
Running parameter set 532/600
Running parameter set 533/600
Running parameter set 534/600
Running parameter set 535/600
Running parameter set 536/600
Running parameter set 537/600
Running parameter set 538/600
Running parameter set 539/600
Running parameter set 540/600
Running parameter set 541/600
Running parameter set 542/600
Running parameter set 543/600
Running parameter set 544/600
Running parameter set 545/600
Running parameter set 546/600
Running parameter set 547/600
Running parameter set 548/600
Running parameter set 549/600
Running parameter set 550/600
Running parameter set 551/600
Running parameter set 552/600
Running parameter set 553/600
Running parameter set 554/600
Running parameter set 555/600
Running parameter set 556/600
Running parameter set 557/600
Running parameter set 558/600
Running parameter set 559/600
Running parameter set 560/600
Running parameter set 561/600
Running parameter set 562/600
Running parameter set 563/600
Running parameter set 564/600
Running parameter set 565/600
Running parameter set 566/600
Running parameter set 567/600
Running parameter set 568/600
Running parameter set 569/600
Running parameter set 570/600
Running parameter set 571/600
Running parameter set 572/600
Running parameter set 573/600
Running parameter set 574/600
Running parameter set 575/600
Running parameter set 576/600
Running parameter set 577/600
Running parameter set 578/600
Running parameter set 579/600
Running parameter set 580/600
Running parameter set 581/600
Running parameter set 582/600
Running parameter set 583/600
Running parameter set 584/600
Running parameter set 585/600
Running parameter set 586/600
Running parameter set 587/600
Running parameter set 588/600
Running parameter set 589/600
Running parameter set 590/600
Running parameter set 591/600
Running parameter set 592/600
Running parameter set 593/600
Running parameter set 594/600
Running parameter set 595/600
Running parameter set 596/600
Running parameter set 597/600
Running parameter set 598/600
Running parameter set 599/600
Running parameter set 600/600
Display best results#
results = pd.DataFrame(results)
print(results[["run", "NSE", "log_NSE"]])
# # Best NSE
# best_run_nse = results["NSE"].idxmax()
# best_result_nse = results.loc[best_run_nse]
# best_parameters_nse = best_result_nse["parameters"]
# best_nse = best_result_nse["NSE"]
# print(f'Best run: {best_run_nse} with NSE: {best_nse}, with parameters:')
# print(list(zip(parameter_names, best_parameters_nse)))
# Best log_NSE
best_run_log_nse = results["log_NSE"].idxmax()
best_result_log_nse = results.loc[best_run_log_nse]
best_parameters_log_nse = best_result_log_nse["parameters"]
best_log_nse = best_result_log_nse["log_NSE"]
print(f'Best run: {best_run_log_nse} with log_NSE: {best_log_nse}, with parameters:')
print(list(zip(parameter_names, best_parameters_log_nse)))
run NSE log_NSE 0 0 0.288376 0.337195 1 1 0.279102 0.480497 2 2 0.384565 0.462463 3 3 -0.085379 0.265060 4 4 0.118825 0.395402 .. ... ... ... 595 595 0.378628 0.445890 596 596 0.190509 0.402412 597 597 0.034397 -0.105571 598 598 0.174568 0.090100 599 599 0.379411 0.461890 [600 rows x 3 columns]
Best run: 232 with NSE: 0.5042075946097913, with parameters:
[ ('Imax', 6.488362330365561), ('Ce', 0.4845571293512798), ('Sumax', 218.83336737778694), ('Beta', 1.7699130278847444), ('Pmax', 0.1358146729030282), ('Tlag', 5.44235281882723), ('Kf', 0.04585512368848791), ('Ks', 0.008840391568068273), ('FM', 0.8983450417631202) ]
Best run: 451 with log_NSE: 0.5985790718618934, with parameters:
[ ('Imax', 6.9919295474161505), ('Ce', 0.44128061023639964), ('Sumax', 183.49045936834898), ('Beta', 1.62907119850059), ('Pmax', 0.5726586861631707), ('Tlag', 5.2154764587256395), ('Kf', 0.06546844116401723), ('Ks', 0.0025930490563726475), ('FM', 1.1427426360454793) ]
# Get ensemble
best_5_runs = results["log_NSE"].nlargest(5).index
best_5_parameters = []
best_5_LNSE = []
for i in range(len(best_5_runs)):
run = best_5_runs[i]
result = results.loc[run]
parameters = result["parameters"]
lnse = result["log_NSE"]
best_5_parameters.append(parameters)
best_5_LNSE.append(log_nse)
print(f"Rank {i+1}: run {run} with LNSE: {lnse}")
print(list(zip(parameter_names, parameters)))
print()
Rank 1: run 451 with LNSE: 0.5985790718618934
[ ('Imax', 6.9919295474161505), ('Ce', 0.44128061023639964), ('Sumax', 183.49045936834898), ('Beta', 1.62907119850059), ('Pmax', 0.5726586861631707), ('Tlag', 5.2154764587256395), ('Kf', 0.06546844116401723), ('Ks', 0.0025930490563726475), ('FM', 1.1427426360454793) ]
Rank 2: run 193 with LNSE: 0.5871245105302814
[ ('Imax', 5.882658193683639), ('Ce', 0.47241275652770587), ('Sumax', 194.30386240959916), ('Beta', 1.781621276884348), ('Pmax', 0.4608178778064105), ('Tlag', 5.215364575839981), ('Kf', 0.06419932503784231), ('Ks', 0.015979806886031615), ('FM', 1.0919722363345776) ]
Rank 3: run 71 with LNSE: 0.5802931655594039
[ ('Imax', 6.279358279165555), ('Ce', 0.45738099679628363), ('Sumax', 218.06826364855027), ('Beta', 1.8389403651683596), ('Pmax', 0.49366891001367486), ('Tlag', 6.037707687593043), ('Kf', 0.07007118899924838), ('Ks', 0.024409700594446738), ('FM', 1.2556224532809368) ]
Rank 4: run 152 with LNSE: 0.5697643546613578
[ ('Imax', 7.301551796430246), ('Ce', 0.49391712556044987), ('Sumax', 211.00486173599964), ('Beta', 1.9612109256650778), ('Pmax', 0.46785716423961), ('Tlag', 6.467897643950326), ('Kf', 0.055226399867357515), ('Ks', 0.03429434734237983), ('FM', 1.191697747701756) ]
Rank 5: run 232 with LNSE: 0.5666729995545778
[ ('Imax', 6.488362330365561), ('Ce', 0.4845571293512798), ('Sumax', 218.83336737778694), ('Beta', 1.7699130278847444), ('Pmax', 0.1358146729030282), ('Tlag', 5.44235281882723), ('Kf', 0.04585512368848791), ('Ks', 0.008840391568068273), ('FM', 0.8983450417631202) ]
Plot Observed vs Modelled#
# Overall Plot
q_critical = 500
simulated_calibrated_nse = run_hbv(
parameters=best_parameters_nse,
initial_storages=s_0,
forcing=ERA5_forcing)
simulated_calibrated_log_nse = run_hbv(
parameters=best_parameters_log_nse,
initial_storages=s_0,
forcing=ERA5_forcing)
# Align dates
simulated_plot_nse = simulated_calibrated_nse
simulated_plot_log_nse = simulated_calibrated_log_nse
observed_plot = observed_output
simulated_plot_nse.index = pd.to_datetime(simulated_plot_nse.index).tz_localize(None).normalize()
simulated_plot_log_nse.index = pd.to_datetime(simulated_plot_log_nse.index).tz_localize(None).normalize()
observed_plot.index = pd.to_datetime(observed_plot.index).tz_localize(None).normalize()
# Combine modelled and observed data
plot_data = pd.DataFrame({
"Modelled discharge NSE": simulated_plot_nse,
"Modelled discharge log_NSE": simulated_plot_log_nse,
"Observed discharge": observed_plot}).dropna()
# Calculate NSE & log NSE for plotted data
plot_nse = NSE(
plot_data["Modelled discharge NSE"],
plot_data["Observed discharge"])
plot_log_nse = log_NSE(
plot_data["Modelled discharge log_NSE"],
plot_data["Observed discharge"])
# Plot
plt.figure(figsize=(14, 6))
plt.plot(plot_data.index, plot_data["Observed discharge"], label="Observed discharge")
# plt.plot(plot_data.index, plot_data["Modelled discharge NSE"], label="Modelled discharge NSE")
plt.plot(plot_data.index, plot_data["Modelled discharge log_NSE"], label="Modelled discharge log_NSE")
plt.xlabel("Date")
plt.ylabel("Discharge (m³/s)")
plt.title(f"Observed vs modelled discharge at 07DA001, NSE = {best_nse:.3f}, log NSE = {best_log_nse:.3f}")
plt.axhline(y=q_critical, linestyle=":", label=f'Critical discharge ({q_critical} m³/s', color='black')
plt.legend()
plt.grid(True)
plt.show()
# Plot for 1 year
selected_year = 2015
plot_data_year = plot_data[plot_data.index.year == selected_year]
plt.figure(figsize=(14, 6))
plt.plot(plot_data_year.index,plot_data_year["Observed discharge"],label="Observed discharge")
# plt.plot(plot_data_year.index,plot_data_year["Modelled discharge NSE"],label="Modelled discharge NSE")
plt.plot(plot_data_year.index,plot_data_year["Modelled discharge log_NSE"],label="Modelled discharge log_NSE", color="green")
plt.axhline(y=q_critical,linestyle=":",color="black",label=f"Critical discharge ({q_critical} m³/s)")
plt.xlabel("Date")
plt.ylabel("Discharge (m³/s)")
plt.title(f"Observed vs modelled discharge at 07DA001 in {selected_year}")
plt.legend()
plt.grid(True)
plt.show()
