Calibrating HBV hydrological model run locally outside container forced with ERA5 forcing data#
In this notebook we will demonstrate how to calibrate the HBV model and works as an example of how to calibrate models in general.
We will use an extention to eWaterCycle: eWaterCycle-DA with DA for Data Assimilation. This package, developed by former MSc student David Haasnoot, adds functionality to deal with ensembles of models in eWaterCycle.
We do now run into a bit of a chicken and egg problem:
calibration of a model needs to be done before running the actual model experiment.
it is better to first demonstrate how to run a model before calibrating. But this requires calibration.
Therefore, too understand how models are run, please have a look at the step 3 notebook first before reading on.
# 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
import json
import os
# Niceties
from rich import print
from tqdm import tqdm
# General eWaterCycle
import ewatercycle
import ewatercycle.models
import ewatercycle.forcing
# We need the ewatercycle_DA package. If that is not available on your machine,
# uncomment the line below to install it
!pip install ewatercycle-da
Requirement already satisfied: ewatercycle-da in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (0.0.7)
Requirement already satisfied: ewatercycle in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle-da) (2.4.0)
Requirement already satisfied: bmipy in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (2.0.1)
Requirement already satisfied: cartopy in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (0.25.0)
Requirement already satisfied: cftime>=1.6.2 in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (1.6.4.post1)
Requirement already satisfied: fiona in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (1.10.1)
Requirement already satisfied: fsspec in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (2024.10.0)
Requirement already satisfied: grpc4bmi>=0.4.0 in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (0.5.0)
Requirement already satisfied: hydrostats in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (0.78)
Requirement already satisfied: matplotlib>=3.5.0 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (3.8.4)
Requirement already satisfied: numpy in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (1.26.4)
Requirement already satisfied: pandas in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (2.2.3)
Requirement already satisfied: pydantic>=2 in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (2.11.9)
Requirement already satisfied: pyoos in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (0.8.4)
Requirement already satisfied: pyopenssl>=24.0.0 in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (25.3.0)
Requirement already satisfied: python-dateutil in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (2.9.0.post0)
Requirement already satisfied: remotebmi in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (0.1.0)
Requirement already satisfied: ruamel-yaml in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (0.18.15)
Requirement already satisfied: shapely in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (2.1.1)
Requirement already satisfied: xarray in /home/mmelotto/.local/lib/python3.12/site-packages (from ewatercycle->ewatercycle-da) (2025.9.0)
Requirement already satisfied: grpcio in /home/mmelotto/.local/lib/python3.12/site-packages (from grpc4bmi>=0.4.0->ewatercycle->ewatercycle-da) (1.75.0)
Requirement already satisfied: grpcio-reflection in /home/mmelotto/.local/lib/python3.12/site-packages (from grpc4bmi>=0.4.0->ewatercycle->ewatercycle-da) (1.62.3)
Requirement already satisfied: grpcio-status in /home/mmelotto/.local/lib/python3.12/site-packages (from grpc4bmi>=0.4.0->ewatercycle->ewatercycle-da) (1.62.3)
Requirement already satisfied: googleapis-common-protos>=1.5.5 in /home/mmelotto/.local/lib/python3.12/site-packages (from grpc4bmi>=0.4.0->ewatercycle->ewatercycle-da) (1.70.0)
Requirement already satisfied: protobuf<5,>=4 in /home/mmelotto/.local/lib/python3.12/site-packages (from grpc4bmi>=0.4.0->ewatercycle->ewatercycle-da) (4.25.8)
Requirement already satisfied: docker in /home/mmelotto/.local/lib/python3.12/site-packages (from grpc4bmi>=0.4.0->ewatercycle->ewatercycle-da) (7.1.0)
Requirement already satisfied: packaging in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from grpc4bmi>=0.4.0->ewatercycle->ewatercycle-da) (24.2)
Requirement already satisfied: typeguard in /home/mmelotto/.local/lib/python3.12/site-packages (from grpc4bmi>=0.4.0->ewatercycle->ewatercycle-da) (4.4.4)
Requirement already satisfied: contourpy>=1.0.1 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from matplotlib>=3.5.0->ewatercycle->ewatercycle-da) (1.3.1)
Requirement already satisfied: cycler>=0.10 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from matplotlib>=3.5.0->ewatercycle->ewatercycle-da) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from matplotlib>=3.5.0->ewatercycle->ewatercycle-da) (4.55.0)
Requirement already satisfied: kiwisolver>=1.3.1 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from matplotlib>=3.5.0->ewatercycle->ewatercycle-da) (1.4.7)
Requirement already satisfied: pillow>=8 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from matplotlib>=3.5.0->ewatercycle->ewatercycle-da) (11.0.0)
Requirement already satisfied: pyparsing>=2.3.1 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from matplotlib>=3.5.0->ewatercycle->ewatercycle-da) (3.2.0)
Requirement already satisfied: annotated-types>=0.6.0 in /home/mmelotto/.local/lib/python3.12/site-packages (from pydantic>=2->ewatercycle->ewatercycle-da) (0.7.0)
Requirement already satisfied: pydantic-core==2.33.2 in /home/mmelotto/.local/lib/python3.12/site-packages (from pydantic>=2->ewatercycle->ewatercycle-da) (2.33.2)
Requirement already satisfied: typing-extensions>=4.12.2 in /home/mmelotto/.local/lib/python3.12/site-packages (from pydantic>=2->ewatercycle->ewatercycle-da) (4.15.0)
Requirement already satisfied: typing-inspection>=0.4.0 in /home/mmelotto/.local/lib/python3.12/site-packages (from pydantic>=2->ewatercycle->ewatercycle-da) (0.4.1)
Requirement already satisfied: cryptography<47,>=45.0.7 in /home/mmelotto/.local/lib/python3.12/site-packages (from pyopenssl>=24.0.0->ewatercycle->ewatercycle-da) (46.0.1)
Requirement already satisfied: cffi>=2.0.0 in /home/mmelotto/.local/lib/python3.12/site-packages (from cryptography<47,>=45.0.7->pyopenssl>=24.0.0->ewatercycle->ewatercycle-da) (2.0.0)
Requirement already satisfied: pycparser in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from cffi>=2.0.0->cryptography<47,>=45.0.7->pyopenssl>=24.0.0->ewatercycle->ewatercycle-da) (2.22)
Requirement already satisfied: six>=1.5 in /home/mmelotto/.local/lib/python3.12/site-packages (from python-dateutil->ewatercycle->ewatercycle-da) (1.17.0)
Requirement already satisfied: black in /home/mmelotto/.local/lib/python3.12/site-packages (from bmipy->ewatercycle->ewatercycle-da) (25.1.0)
Requirement already satisfied: click in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from bmipy->ewatercycle->ewatercycle-da) (8.1.7)
Requirement already satisfied: jinja2 in /home/mmelotto/.local/lib/python3.12/site-packages (from bmipy->ewatercycle->ewatercycle-da) (3.1.6)
Requirement already satisfied: mypy-extensions>=0.4.3 in /home/mmelotto/.local/lib/python3.12/site-packages (from black->bmipy->ewatercycle->ewatercycle-da) (1.1.0)
Requirement already satisfied: pathspec>=0.9.0 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from black->bmipy->ewatercycle->ewatercycle-da) (0.12.1)
Requirement already satisfied: platformdirs>=2 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from black->bmipy->ewatercycle->ewatercycle-da) (4.3.6)
Requirement already satisfied: pyshp>=2.3 in /home/mmelotto/.local/lib/python3.12/site-packages (from cartopy->ewatercycle->ewatercycle-da) (2.3.1)
Requirement already satisfied: pyproj>=3.3.1 in /home/mmelotto/.local/lib/python3.12/site-packages (from cartopy->ewatercycle->ewatercycle-da) (3.7.2)
Requirement already satisfied: certifi in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from pyproj>=3.3.1->cartopy->ewatercycle->ewatercycle-da) (2024.8.30)
Requirement already satisfied: requests>=2.26.0 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from docker->grpc4bmi>=0.4.0->ewatercycle->ewatercycle-da) (2.32.3)
Requirement already satisfied: urllib3>=1.26.0 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from docker->grpc4bmi>=0.4.0->ewatercycle->ewatercycle-da) (2.2.3)
Requirement already satisfied: charset-normalizer<4,>=2 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from requests>=2.26.0->docker->grpc4bmi>=0.4.0->ewatercycle->ewatercycle-da) (3.4.0)
Requirement already satisfied: idna<4,>=2.5 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from requests>=2.26.0->docker->grpc4bmi>=0.4.0->ewatercycle->ewatercycle-da) (3.10)
Requirement already satisfied: attrs>=19.2.0 in /home/mmelotto/.local/lib/python3.12/site-packages (from fiona->ewatercycle->ewatercycle-da) (25.3.0)
Requirement already satisfied: click-plugins>=1.0 in /home/mmelotto/.local/lib/python3.12/site-packages (from fiona->ewatercycle->ewatercycle-da) (1.1.1.2)
Requirement already satisfied: cligj>=0.5 in /home/mmelotto/.local/lib/python3.12/site-packages (from fiona->ewatercycle->ewatercycle-da) (0.7.2)
Requirement already satisfied: numba in /home/mmelotto/.local/lib/python3.12/site-packages (from hydrostats->ewatercycle->ewatercycle-da) (0.61.2)
Requirement already satisfied: scipy in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from hydrostats->ewatercycle->ewatercycle-da) (1.16.3)
Requirement already satisfied: HydroErr in /home/mmelotto/.local/lib/python3.12/site-packages (from hydrostats->ewatercycle->ewatercycle-da) (1.24)
Requirement already satisfied: MarkupSafe>=2.0 in /home/mmelotto/.local/lib/python3.12/site-packages (from jinja2->bmipy->ewatercycle->ewatercycle-da) (3.0.3)
Requirement already satisfied: llvmlite<0.45,>=0.44.0dev0 in /home/mmelotto/.local/lib/python3.12/site-packages (from numba->hydrostats->ewatercycle->ewatercycle-da) (0.44.0)
Requirement already satisfied: pytz>=2020.1 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from pandas->ewatercycle->ewatercycle-da) (2024.1)
Requirement already satisfied: tzdata>=2022.7 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from pandas->ewatercycle->ewatercycle-da) (2024.2)
Requirement already satisfied: paegan>=1.1 in /home/mmelotto/.local/lib/python3.12/site-packages (from pyoos->ewatercycle->ewatercycle-da) (1.1.4)
Requirement already satisfied: OWSLib>=0.8.3 in /home/mmelotto/.local/lib/python3.12/site-packages (from pyoos->ewatercycle->ewatercycle-da) (0.34.1)
Requirement already satisfied: beautifulsoup4>=4.2.1 in /home/mmelotto/.local/lib/python3.12/site-packages (from pyoos->ewatercycle->ewatercycle-da) (4.14.2)
Requirement already satisfied: lxml>=3.2.0 in /home/mmelotto/.conda/envs/ewatercycle_snellius/lib/python3.12/site-packages (from pyoos->ewatercycle->ewatercycle-da) (5.3.0)
Requirement already satisfied: soupsieve>1.2 in /home/mmelotto/.local/lib/python3.12/site-packages (from beautifulsoup4>=4.2.1->pyoos->ewatercycle->ewatercycle-da) (2.8)
Requirement already satisfied: pyyaml in /home/mmelotto/.local/lib/python3.12/site-packages (from OWSLib>=0.8.3->pyoos->ewatercycle->ewatercycle-da) (6.0.3)
Requirement already satisfied: netCDF4>=1.0.2 in /home/mmelotto/.local/lib/python3.12/site-packages (from paegan>=1.1->pyoos->ewatercycle->ewatercycle-da) (1.7.2)
Requirement already satisfied: connexion in /home/mmelotto/.local/lib/python3.12/site-packages (from remotebmi->ewatercycle->ewatercycle-da) (3.2.0)
Requirement already satisfied: httpx>=0.27.2 in /home/mmelotto/.local/lib/python3.12/site-packages (from remotebmi->ewatercycle->ewatercycle-da) (0.28.1)
Requirement already satisfied: uvicorn in /home/mmelotto/.local/lib/python3.12/site-packages (from remotebmi->ewatercycle->ewatercycle-da) (0.35.0)
Requirement already satisfied: anyio in /home/mmelotto/.local/lib/python3.12/site-packages (from httpx>=0.27.2->remotebmi->ewatercycle->ewatercycle-da) (4.11.0)
Requirement already satisfied: httpcore==1.* in /home/mmelotto/.local/lib/python3.12/site-packages (from httpx>=0.27.2->remotebmi->ewatercycle->ewatercycle-da) (1.0.9)
Requirement already satisfied: h11>=0.16 in /home/mmelotto/.local/lib/python3.12/site-packages (from httpcore==1.*->httpx>=0.27.2->remotebmi->ewatercycle->ewatercycle-da) (0.16.0)
Requirement already satisfied: sniffio>=1.1 in /home/mmelotto/.local/lib/python3.12/site-packages (from anyio->httpx>=0.27.2->remotebmi->ewatercycle->ewatercycle-da) (1.3.1)
Requirement already satisfied: asgiref>=3.4 in /home/mmelotto/.local/lib/python3.12/site-packages (from connexion->remotebmi->ewatercycle->ewatercycle-da) (3.9.1)
Requirement already satisfied: inflection>=0.3.1 in /home/mmelotto/.local/lib/python3.12/site-packages (from connexion->remotebmi->ewatercycle->ewatercycle-da) (0.5.1)
Requirement already satisfied: jsonschema>=4.17.3 in /home/mmelotto/.local/lib/python3.12/site-packages (from connexion->remotebmi->ewatercycle->ewatercycle-da) (4.25.1)
Requirement already satisfied: python-multipart>=0.0.15 in /home/mmelotto/.local/lib/python3.12/site-packages (from connexion->remotebmi->ewatercycle->ewatercycle-da) (0.0.20)
Requirement already satisfied: starlette>=0.35 in /home/mmelotto/.local/lib/python3.12/site-packages (from connexion->remotebmi->ewatercycle->ewatercycle-da) (0.48.0)
Requirement already satisfied: werkzeug>=2.2.1 in /home/mmelotto/.local/lib/python3.12/site-packages (from connexion->remotebmi->ewatercycle->ewatercycle-da) (3.1.3)
Requirement already satisfied: jsonschema-specifications>=2023.03.6 in /home/mmelotto/.local/lib/python3.12/site-packages (from jsonschema>=4.17.3->connexion->remotebmi->ewatercycle->ewatercycle-da) (2025.9.1)
Requirement already satisfied: referencing>=0.28.4 in /home/mmelotto/.local/lib/python3.12/site-packages (from jsonschema>=4.17.3->connexion->remotebmi->ewatercycle->ewatercycle-da) (0.36.2)
Requirement already satisfied: rpds-py>=0.7.1 in /home/mmelotto/.local/lib/python3.12/site-packages (from jsonschema>=4.17.3->connexion->remotebmi->ewatercycle->ewatercycle-da) (0.27.1)
Requirement already satisfied: ruamel.yaml.clib>=0.2.7 in /home/mmelotto/.local/lib/python3.12/site-packages (from ruamel-yaml->ewatercycle->ewatercycle-da) (0.2.12)
[notice] A new release of pip is available: 25.2 -> 25.3
[notice] To update, run: pip install --upgrade pip
# eWaterCycle Data assimilation package
from ewatercycle_DA import DA
# Parameters
region_id = None
settings_path = "settings.json"
# Parameters
region_id = "camelsgb_23008"
settings_path = "regions/camelsgb_23008/settings.json"
# Load settings
# Read from the JSON file
with open(settings_path, "r") as json_file:
settings = json.load(json_file)
display(settings)
{'caravan_id': 'camelsgb_23008',
'calibration_start_date': '1994-08-01T00:00:00Z',
'calibration_end_date': '2004-07-31T00:00:00Z',
'validation_start_date': '2004-08-01T00:00:00Z',
'validation_end_date': '2014-07-31T00:00:00Z',
'future_start_date': '2029-08-01T00:00:00Z',
'future_end_date': '2049-08-31T00:00:00Z',
'CMIP_info': {'dataset': ['MPI-ESM1-2-HR'],
'ensembles': ['r1i1p1f1'],
'experiments': ['historical', 'ssp126', 'ssp245', 'ssp370', 'ssp585'],
'project': 'CMIP6',
'frequency': 'day',
'grid': 'gn',
'variables': ['pr', 'tas', 'rsds']},
'base_path': '/gpfs/scratch1/shared/mmelotto/ewatercycleClimateImpact/HBV',
'path_caravan': '/gpfs/scratch1/shared/mmelotto/ewatercycleClimateImpact/HBV/forcing_data/camelsgb_23008/caravan',
'path_ERA5': '/gpfs/scratch1/shared/mmelotto/ewatercycleClimateImpact/HBV/forcing_data/camelsgb_23008/ERA5',
'path_CMIP6': '/gpfs/scratch1/shared/mmelotto/ewatercycleClimateImpact/HBV/forcing_data/camelsgb_23008/CMIP6',
'path_output': '/gpfs/scratch1/shared/mmelotto/ewatercycleClimateImpact/HBV/output_data/camelsgb_23008',
'path_shape': '/gpfs/scratch1/shared/mmelotto/ewatercycleClimateImpact/HBV/forcing_data/camelsgb_23008/caravan/camelsgb_23008.shp',
'downloads': '/gpfs/scratch1/shared/mmelotto/ewatercycleClimateImpact/HBV/downloads/camelsgb_23008'}
# We check if we already ran this test
file_parameters_path = settings['path_output'] + "/" + settings['caravan_id'] + "_params_MC.csv"
need_to_run = True
if os.path.exists(file_parameters_path):
display(f"File already exists: {file_parameters_path}")
need_to_run = False
# sys.exit() # stops notebook execution
'File already exists: /gpfs/scratch1/shared/mmelotto/ewatercycleClimateImpact/HBV/output_data/camelsgb_23008/camelsgb_23008_params_MC.csv'
Pre-generated observations of discharge from caravan#
Here we re-load the disharge data we generated in this notebook.
# Load the caravan forcing object
caravan_data_object = ewatercycle.forcing.sources['CaravanForcing'].load(directory=settings['path_caravan'])
display(caravan_data_object)
CaravanForcing(start_time='1994-08-01T00:00:00Z', end_time='2014-07-31T00:00:00Z', directory=PosixPath('/gpfs/scratch1/shared/mmelotto/ewatercycleClimateImpact/HBV/forcing_data/camelsgb_23008/caravan'), shape=PosixPath('/gpfs/scratch1/shared/mmelotto/ewatercycleClimateImpact/HBV/forcing_data/camelsgb_23008/caravan/camelsgb_23008.shp'), filenames={'evspsblpot': 'camelsgb_23008_1994-08-01_2014-07-31_evspsblpot.nc', 'tasmin': 'camelsgb_23008_1994-08-01_2014-07-31_tasmin.nc', 'tasmax': 'camelsgb_23008_1994-08-01_2014-07-31_tasmax.nc', 'tas': 'camelsgb_23008_1994-08-01_2014-07-31_tas.nc', 'pr': 'camelsgb_23008_1994-08-01_2014-07-31_pr.nc', 'Q': 'camelsgb_23008_1994-08-01_2014-07-31_Q.nc'})
Pre-generated ERA5 forcing data for HBV model#
Here we load the ERA5 data we generated in this notebook
load_location = Path(settings['path_ERA5']) / "work" / "diagnostic" / "script"
ERA5_forcing_object = ewatercycle.forcing.sources["LumpedMakkinkForcing"].load(directory=load_location)
display(ERA5_forcing_object)
LumpedMakkinkForcing(start_time='1994-08-01T00:00:00Z', end_time='2014-07-31T00:00:00Z', directory=PosixPath('/gpfs/scratch1/shared/mmelotto/ewatercycleClimateImpact/HBV/forcing_data/camelsgb_23008/ERA5/work/diagnostic/script'), shape=PosixPath('/gpfs/scratch1/shared/mmelotto/ewatercycleClimateImpact/HBV/forcing_data/camelsgb_23008/ERA5/work/diagnostic/script/camelsgb_23008.shp'), filenames={'pr': 'OBS6_ERA5_reanaly_1_day_pr_1994-2014.nc', 'tas': 'OBS6_ERA5_reanaly_1_day_tas_1994-2014.nc', 'rsds': 'OBS6_ERA5_reanaly_1_day_rsds_1994-2014.nc', 'evspsblpot': 'Derived_Makkink_evspsblpot.nc'})
Calibration basics and objective function#
In model calibration, we are looking for a set of parameters such that when the model is run with that set of parameters, we get the best model output. What “best” means differs per application or research question. In general, we like some model outputs to be as close as possible to observations. For this purpose we create an objective function that takes the model output of interest and observations as inputs and calculates some score that shows goodness of fit. Here we use a RMS difference function:
def calibrationObjective(modelOutput, observation, start_calibration, end_calibration):
'''A function that takes in two dataFrames, interpolates the model output to the
observations and calculates the average absolute difference between the two. '''
# Combine the two in one dataFrame
hydro_data = pd.concat([modelOutput.reindex(observation.index, method = 'ffill'), observation], axis=1,
keys=['model','observation'])
# Only select the calibration period
hydro_data = hydro_data[hydro_data.index > pd.to_datetime(pd.Timestamp(start_calibration).date())]
hydro_data = hydro_data[hydro_data.index < pd.to_datetime(pd.Timestamp(end_calibration).date())]
# Calculate RMS difference
squareDiff = (hydro_data['model'] - hydro_data['observation'])**2
rootMeanSquareDiff = np.sqrt(np.mean(squareDiff))
return rootMeanSquareDiff
Create an ensemble of models#
Instead of single model, we create an ensemble of models. In our case each ensemblemember is a HBV model that will get its own parameters. After running the entire ensemble we will apply the calibration objective function to determine the best set of parameters.
# Set the number of ensemble members. In Data Assimilation "ensemble member" and "particle" is used interchangeably
# Based on which school of DA you come from :-)
n_particles = 1000
# Create an array with parameter values.
# First set minimum and maximum values on the parameters
p_min_initial = np.array([0, 0.2, 40, .5, .001, 1, .01, .0001, 0.01])
p_max_initial = np.array([8, 1, 800, 4, .3, 10, .1, .01, 10.0])
# Create an empty array to store the parameter sets
parameters = np.zeros([len(p_min_initial), n_particles])
# Fill the array with random values between the minimum and maximum
for param in range(len(p_min_initial)):
parameters[param,:] = np.random.uniform(p_min_initial[param],p_max_initial[param],n_particles)
# Print parameter names and values for first ensemble member
param_names = ["Imax", "Ce", "Sumax", "Beta", "Pmax", "Tlag", "Kf", "Ks", "FM"]
display(list(zip(param_names, np.round(parameters[:,0], decimals=3))))
[('Imax', 7.76),
('Ce', 0.66),
('Sumax', 162.464),
('Beta', 3.859),
('Pmax', 0.087),
('Tlag', 6.682),
('Kf', 0.021),
('Ks', 0.005),
('FM', 7.71)]
# Set initial state values
# Si, Su, Sf, Ss, Sp
s_0 = np.array([0, 100, 0, 5, 0])
# Each ensemble member gets their own parameters
# which are set during the initialize phase.
# Here we make a list of 'arguments' to pass to the model
# during initialize.
setup_kwargs_lst = []
for index in range(n_particles):
setup_kwargs_lst.append({'parameters': parameters[:,index],
})
if need_to_run:
# Create the ensemble object
ensemble = DA.Ensemble(N=n_particles)
ensemble.setup()
In the initialize step below we specify which model we will be using and pass the list of parameters. For other purposes (multimodel comparisons) here we could also provide lists of different models each with their own forcing and other setup arguments.
if need_to_run:
# This initializes the models for all ensemble members.
ensemble.initialize(model_name=["HBVLocal"]*n_particles,
forcing=[ERA5_forcing_object]*n_particles,
setup_kwargs=setup_kwargs_lst)
if need_to_run:
# We appoint one of the ensemble members the role "reference model".
# In this use case, this is pure for timekeeping as shown in the next cell
ref_model = ensemble.ensemble_list[0].model
Models run with one command.#
All models can be run with the ensemble.update() command.
if need_to_run:
n_timesteps = int((ref_model.end_time - ref_model.start_time) / ref_model.time_step)
time = []
lst_Q = []
for i in tqdm(range(n_timesteps)):
time.append(pd.Timestamp(ref_model.time_as_datetime.date()))
ensemble.update()
lst_Q.append(ensemble.get_value("Q").flatten())
if need_to_run:
ensemble.finalize()
Find best parameter set#
By calculating the objective function for each model output, we can search the combination of parameters with the lowest objective function.
if need_to_run:
# Create a pandas dataframe to hold all the model outputs
Q_m_arr = np.array(lst_Q).T
df_ensemble = pd.DataFrame(data=Q_m_arr[:,:len(time)].T,index=time,columns=[f'particle {n}' for n in range(n_particles)])
if need_to_run:
# Create a dataframe for the observations
ds_observation = xr.open_mfdataset([caravan_data_object['Q']]).to_pandas()
if need_to_run:
objective_values_calibration = []
for i in tqdm(range(n_particles)):
objective_values_calibration.append(calibrationObjective(df_ensemble.iloc[:,i],ds_observation["Q"],
settings['calibration_start_date'],
settings['calibration_end_date']))
if need_to_run:
# Make some plot of the spread of the objective functions for the different parameters
xFigNr = 2
yFigNr = 4
fig, axs = plt.subplots(xFigNr, yFigNr,figsize = (15,15))
for xFig in range(xFigNr):
for yFig in range(yFigNr):
paramCounter = xFig*yFigNr + yFig
axs[xFig,yFig].plot(parameters[paramCounter,:],objective_values_calibration,'.')
axs[xFig,yFig].set_title(param_names[paramCounter])
if need_to_run:
# Let's also print the minimal values:
parameters_minimum_index = np.argmin(np.array(objective_values_calibration))
parameters_minimum = parameters[:,parameters_minimum_index]
display(list(zip(param_names, np.round(parameters_minimum, decimals=3))))
if need_to_run:
# Make a plot of the model output of the minimum value
ds_observation["Q"].plot()
ax = df_ensemble.iloc[:,parameters_minimum_index].plot(lw=2.5)
plt.legend()
if need_to_run:
df_best = df_ensemble.iloc[:,parameters_minimum_index]
if need_to_run:
df_best.index = df_best.index.tz_localize("UTC")
df_select = df_best.tz_convert("UTC")[settings['validation_start_date']:settings['validation_end_date']]
if need_to_run:
# Make a plot of the model output of the minimum value
ds_observation["Q"].plot()
ax = df_select.plot(lw=2.5)
plt.legend()
plt.xlim(settings['validation_start_date'],settings['validation_end_date'])
Save results#
We want to save these results to file to be able to load them in other studies
if need_to_run:
# Again make sure folder exists
Path(settings["path_output"]).mkdir(exist_ok=True, parents=True)
# Save to csv file
np.savetxt(settings['path_output'] + "/" + settings['caravan_id'] + "_params_MC.csv", parameters_minimum, delimiter=",")