06 - Calibration of a Raven hydrological model

Calibration of a Raven model

In this notebook, we show how to calibrate a Raven model using the GR4J-CN predefined structure. Users can refer to the documentation for the parameterization of other hydrological model structures.

Let’s start by importing the packages that will do the work.

import datetime as dt

import spotpy

from ravenpy.config import commands as rc
from ravenpy.config.emulators import GR4JCN
from ravenpy.utilities.calibration import SpotSetup

Preparing the model to be calibrated on a given watershed

Our test watershed from the last notebook is selected for this test. It can be replaced with any desired watershed.

from ravenpy.utilities.testdata import get_file

# We get the netCDF for testing on a server. You can replace the getfile method by a string containing the path to your own netCDF
nc_file = get_file(
    "raven-gr4j-cemaneige/Salmon-River-Near-Prince-George_meteo_daily.nc"
)

# Display the dataset that we will be using
print(nc_file)

The process is very similar to setting up a hydrological model. We first need to create the model with its configuration. We must provide the same information as before, except for the model parameters since those need to be calibrated.

# Here, we need to give the name of your different dataset in order to match with Raven models.
alt_names = {
    "RAINFALL": "rain",
    "SNOWFALL": "snow",
    "TEMP_MIN": "tmin",
    "TEMP_MAX": "tmax",
    "PET": "pet",
    "HYDROGRAPH": "qobs",
}

# The HRU of your watershed
hru = dict(area=4250.6, elevation=843.0, latitude=54.4848, longitude=-123.3659)

# You can decide the evaluation metrics that will be used to calibrate the parameters of your model. You need atleast
# 1 evaluation metric, but you can do any combinaison of evaluations from this list :
#
# NASH_SUTCLIFFE,
# LOG_NASH,
# RMSE,
# PCT_BIAS,
# ABSERR,
# ABSMAX,
# PDIFF,
# TMVOL,
# RCOEFF,
# NSC,
# KLING_GUPTA
eval_metrics = ("NASH_SUTCLIFFE",)


# We need to create the desired model with its parameters the same way as in the Notebook 04_Emulating_hydrological_models.
model_config = GR4JCN(
    ObservationData=[rc.ObservationData.from_nc(nc_file, alt_names="qobs")],
    Gauge=[
        rc.Gauge.from_nc(
            nc_file,
            alt_names=alt_names,
            data_kwds={"ALL": {"elevation": hru["elevation"]}},
        )
    ],
    HRUs=[hru],
    StartDate=dt.datetime(1990, 1, 1),
    EndDate=dt.datetime(1999, 12, 31),
    RunName="test",
    EvaluationMetrics=eval_metrics,  # We add this code to tell Raven which objective function we want to pass.
    SuppressOutput=True,
)

Spotpy Calibration

Once you’ve created your model, you need to create a SpotSetup, which will be used to calibrate your model.

# In order to calibrate your model, you need to give the lower and higher bounds of the model. In this case, we are passing
# the boundaries for a GR4JCN, but it's important to change them, if you are using another model. Note that the list of these
# boundaries for each model is at the end of this notebook.
low_params = (0.01, -15.0, 10.0, 0.0, 1.0, 0.0)
high_params = (2.5, 10.0, 700.0, 7.0, 30.0, 1.0)


spot_setup = SpotSetup(
    config=model_config,
    low=low_params,
    high=high_params,
)

Now that the model is setup and configured and that SpotSetup object exists, we need to create a sampler from spotpy module which will optimize the hydrological model paramaters. You can see that we are using the DDS algorithm to optimize the parameters:

[Tolson, B.A. and Shoemaker, C.A., 2007. Dynamically dimensioned search algorithm for computationally efficient watershed model calibration. Water Resources Research, 43(1)].

If you want to use another algorithm, please refer to the Spotpy documentation here : https://spotpy.readthedocs.io/

Finally, we run the sampler by the amount of desired repetitions.

# Number of total model evaluations in the calibration. This value should be over 500 for real optimisation,
# and upwards of 10000 evaluations for models with many parameters. This will take a LONG period of time so
# be sure of all the configuration above before executing with a high number of model evaluations.
model_evaluations = 10

# Set up the spotpy sampler with the method, the setup configuration, a run name and other options. Please refer to
# the spotpy documentation for more options. We recommend sticking to this format for efficiency of most applications.
sampler = spotpy.algorithms.dds(
    spot_setup, dbname="RAVEN_model_run", dbformat="ram", save_sim=False
)

# Launch the actual optimization. Multiple trials can be launched, where the entire process is repeated and
# the best overall value from all trials is returned.
sampler.sample(model_evaluations, trials=1)

Analysing the calibration results

The best parameters as well as the objective functions can be analyzed.

# Get all the values of each iteration
results = sampler.getdata()

print("The best Nash-Sutcliffe value is : ")

# Get the raw resutlts directly in an array
bestindex, bestobjfun = spotpy.analyser.get_maxlikeindex(
    results
)  # Want to get the MAX NSE (change for min for RMSE)
best_model_run = list(
    results[bestindex][0]
)  # Get the parameter set returning the best NSE
optimized_parameters = best_model_run[
    1:-1
]  # Remove the NSE value (position 0) and the ID at the last position to get the actual parameter set.

print("\nThe best parameters are : ")
# Display the parameter set ready to use in a future run:
print(optimized_parameters)

Next steps

In the next notebooks, we will apply the model to specific use-cases, including making and using hotstart files for forecasting, performing hindcasting and forecasting, applying data assimilation and evaluating the impacts of climate change on the hydrology of a watershed. In the meantime, you can explore calibration with any of the emulated models below with the provided low and high bounds. You can also provide your own for specific cases.

List of Model-Boundaries

GR4J-CN :

low = (0.01, -15.0, 10.0, 0.0, 1.0, 0.0),
high = (2.5, 10.0, 700.0, 7.0, 30.0, 1.0)

HMETS :

low = (0.3, 0.01, 0.5, 0.15, 0.0, 0.0, -2.0, 0.01, 0.0, 0.01, 0.005, -5.0, 0.0, 0.0, 0.0, 0.0, 0.00001, 0.0, 0.00001, 0.0, 0.0),
high = (20.0, 5.0, 13.0, 1.5, 20.0, 20.0, 3.0, 0.2, 0.1, 0.3, 0.1, 2.0, 5.0, 1.0, 3.0, 1.0, 0.02, 0.1, 0.01, 0.5, 2.0)

Mohyse :

low = (0.01, 0.01, 0.01, -5.00, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01),
high = (20.0, 1.0, 20.0, 5.0, 0.5, 1.0, 1.0, 1.0, 15.0, 15.0)

HBV-EC :

low = (-3.0, 0.0, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0, 0.01, 0.05, 0.01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.01, 0.0, 0.05, 0.8, 0.8),
high = (3.0, 8.0, 8.0, 0.1, 1.0, 1.0, 7.0, 100.0, 1.0, 0.1, 6.0, 5.0, 5.0, 0.2, 1.0, 30.0, 3.0, 2.0, 1.0, 1.5, 1.5)

CanadianShield :

low = (0.01, 0.01, 0.01, 0.0, 0.0, 0.05, 0.0, -5.0, -5.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.01, 0.005, -3.0, 0.5, 5.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.8, 0.0),
high = (0.5, 2.0, 3.0, 3.0, 0.05, 0.45, 7.0, -1.0, -1.0, 2.0, 2.0, 100.0, 100.0, 100.0, 0.4, 0.1, 0.3, 0.1, 3.0, 4.0, 500.0, 5.0, 5.0, 1.0, 8.0, 20.0, 1.5, 0.2, 0.2, 10.0, 10.0, 1.2, 1.2, 1.0)

HYPR :

low = (-1.0, -3.0, 0.0, 0.3, -1.3, -2.0, 0.0, 0.1, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.01, 0.0, 0.0, 1.5, 0.0, 0.0, 0.8),
high = (1.0, 3.0, 0.8, 1.0, 0.3, 0.0, 30.0, 0.8, 2.0, 100.0, 0.5, 5.0, 1.0, 1000.0, 6.0, 7.0, 8.0, 3.0, 5.0, 5.0, 1.2)

SACSMA :

low = (-3.0, -1.52287874, -0.69897, 0.025, 0.01, 0.075, 0.015, 0.04, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3, 0.01, 0.8, 0.8),
high = (-1.82390874, -0.69897, -0.30102999, 0.125, 0.075, 0.3, 0.3, 0.6, 0.5, 3.0, 80.0, 0.8, 0.05, 0.2, 0.1, 0.4, 8.0, 20.0, 5.0, 1.2, 1.2)

Blended :

low = (0.0, 0.1, 0.5, -5.0, 0.0, 0.5, 5.0, 0.0, 0.0, 0.0, -5.0, 0.5, 0.0, 0.01, 0.005, -5.0, 0.0, 0.0, 0.0, 0.3, 0.01, 0.5, 0.15, 1.5, 0.0, -1.0, 0.01, 0.00001, 0.0, 0.0, -3.0, 0.5, 0.8, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
high = (1.0, 3.0, 3.0, -1.0, 100.0, 2.0, 10.0, 3.0, 0.05, 0.45, -2.0, 2.0, 0.1, 0.3, 0.1, 2.0, 1.0, 5.0, 0.4, 20.0, 5.0, 13.0, 1.5, 3.0, 5.0, 1.0, 0.2, 0.02, 0.5, 2.0, 3.0, 4.0, 1.2, 1.2, 0.02, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)