.. _TutorialCalibration: ================================ Running a calibration with Shyft ================================ :ref:`setup-shyft` This notebook is guiding through the simulation process of a catchment. The following steps are described: 1. **Loading required python modules and setting path to Shyft installation** 2. **Configuration of a Shyft calibration** 3. **Running a Shyft calibration** 4. **Inspecting the calibration results** 1. Loading required python modules and setting path to Shyft installation ------------------------------------------------------------------------- Shyft requires a number of different modules to be loaded as part of the package. Below, we describe the required steps for loading the modules, and note that some steps are only required for the use of the tutorial. .. code-block:: python # Pure python modules and jupyter functionality # first you should import the third-party python modules which you'll use later on # the first line enables that figures are shown inline, directly in the notebook %matplotlib inline import os import datetime as dt import pandas as pd from os import path import sys from matplotlib import pyplot as plt The Shyft Environment ^^^^^^^^^^^^^^^^^^^^^ Set shyft_data_path: .. code-block:: python # try to auto-configure the path, -will work in all cases where doc and data # are checked out at same level from shyft.hydrology import shyftdata_dir shyft_data_path = os.path.abspath("/workspaces/shyft-data/") if path.exists(shyft_data_path) and 'SHYFT_DATA' not in os.environ: os.environ['SHYFT_DATA']=shyft_data_path print(shyft_data_path) print(shyftdata_dir) .. code-block:: python # importing the shyft modules needed for running a calibration from shyft.hydrology.repository.default_state_repository import DefaultStateRepository from shyft.hydrology.orchestration.configuration.yaml_configs import YAMLCalibConfig, YAMLSimConfig from shyft.hydrology.orchestration.simulators.config_simulator import ConfigCalibrator, ConfigSimulator 2. Configuration of a Shyft calibration --------------------------------------- .. code-block:: python # conduct a configured simulation first. config_file_path = os.path.join(shyft_data_path, "neanidelv/yaml_config/neanidelva_simulation.yaml") cfg = YAMLSimConfig(config_file_path, "neanidelva") simulator = ConfigSimulator(cfg) # run the model, and we'll just pull the `api.model` from the `simulator` simulator.run() state = simulator.region_model.state Now that we have the initial state, we'll run the calibration (this is not a strictly required step, but we use it later) .. code-block:: python # set up configuration using *.yaml configuration files config_file_path = os.path.join(shyft_data_path,"neanidelv/yaml_config/neanidelva_calibration.yaml") # here is the *.yaml file cfg = YAMLCalibConfig(config_file_path, "neanidelva") .. code-block:: python # initialize an instance of the orchestration's ConfigCalcalibrator class, which has all the functionality needed # to run a calibration using the above initiated configuration calib = ConfigCalibrator(cfg) n_cells = calib.region_model.size() state_repos = DefaultStateRepository(calib.region_model) # Notice that this repository needs the real model # so that it's able to generate a precise # default state-with-id vector for this # specific model # # 3. Running a Shyft calibration ------------------------------ Calibration may take up to 10 minutes. .. code-block:: python # once the calibrator is set up, all you need to do is running the calibration... # the calibrated parameters are stored in a model.yaml. results = calib.calibrate(cfg.sim_config.time_axis, state_repos.get_state(0), cfg.optimization_method['name'], cfg.optimization_method['params']) 4. Inspecting the calibration results ------------------------------------- First the Nash-Suttcliffe-efficiency of the calibrated simulation is computed to see the quality of the calibration. Then the calibrated model parameters are accessed and printed out. .. code-block:: python # Get NSE of calibrated run: result_params = [] for i in range(results.size()): result_params.append(results.get(i)) print("Final NSE =", 1-calib.optimizer.calculate_goal_function(result_params)) # Final NSE = 0.8521010139040602 .. code-block:: python # Check out the calibrated parameters. diff = 1.0E-3 print("{0:30s} {1:10s}".format("PARAM-NAME", "CALIB-VALUE")) for i in range(results.size()): print("{0:30s} {1:10f}".format(results.get_name(i), results.get(i))) Output: .. code-block:: none PARAM-NAME CALIB-VALUE kirchner.c1 -2.865322 kirchner.c2 0.641488 kirchner.c3 -0.097761 ae.ae_scale_factor 1.500000 gs.tx -0.144348 gs.wind_scale 1.841114 gs.max_water 0.100000 gs.wind_const 1.000000 gs.fast_albedo_decay_rate 13.656055 gs.slow_albedo_decay_rate 27.828786 gs.surface_magnitude 30.000000 gs.max_albedo 0.900000 gs.min_albedo 0.600000 gs.snowfall_reset_depth 5.000000 gs.snow_cv 0.400000 gs.glacier_albedo 0.400000 p_corr.scale_factor 1.000000 gs.snow_cv_forest_factor 0.000000 gs.snow_cv_altitude_factor 0.000000 pt.albedo 0.200000 pt.alpha 1.260000 gs.initial_bare_ground_fraction 0.040000 gs.winter_end_day_of_year 100.000000 gs.calculate_iso_pot_energy 0.000000 gm.dtf 5.000000 routing.velocity 0.000000 routing.alpha 0.900000 routing.beta 3.000000 gs.n_winter_days 100.000000 gm.direct_response 0.000000 msp.reservoir_direct_response_fraction 1.000000 Plotting simulated and observed discharge ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ We are now plotting the simulated and observed discharge timeseries over the course of the melt period. .. code-block:: python # get the target vector and discharge statistics from the configured calibrator target_obs = calib.tv[0] disch_sim = calib.region_model.statistics.discharge(target_obs.catchment_indexes) disch_obs = target_obs.ts.values ts_timestamps_sim = [dt.datetime.utcfromtimestamp(p) for p in disch_sim.time_axis.time_points] ts_timestamps_sim = ts_timestamps_sim[:-1] ts_timestamps_obs = [dt.datetime.utcfromtimestamp(p) for p in target_obs.ts.time_axis.time_points] ts_timestamps_obs = ts_timestamps_obs[:-1] .. code-block:: python # plot up the results fig, ax = plt.subplots(1, figsize=(15,10)) ax.plot(ts_timestamps_sim, disch_sim.values, lw=2, label = "sim") ax.plot(ts_timestamps_obs, disch_obs, lw=2, ls='--', label = "obs") ax.set_title("observed and simulated discharge") ax.legend() ax.set_ylabel("discharge [m3 s-1]") .. image:: plotting_simulated_vs_observed.png 5. Changing parameters on-the-fly --------------------------------- Instead of changing model parameters in the yaml-configs, reload the configuration and re-run the model, we can also just change the parameters "on-the-fly", and rerun the model. This makes it easy to investigate the influence of certain model parameters on the simulation results. 5a. Snow or rain? The parameter gs.tx sets the threshold temperature at which precipitation is treated as snow fall. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ In the following we'll investigate the impact of manually manipulating this parameter. .. code-block:: python parameters = calib.region_model.get_region_parameter() # fetching parameters from the simulator object print(u"Calibrated rain/snow threshold temp: {} C".format(parameters.gs.tx)) # print current value of hs.tx # output: Calibrated rain/snow threshold temp: -0.14434796435811803 C In the following, we first set the hs.tx parameter to a higher, and then to a lower value compared to the value the calibration results suggest. We re-run the simulation, respicetively, and plot the results. .. code-block:: python calib.optimizer.calculate_goal_function(result_params) # reset the parameters to the values of the calibration parameters.gs.tx = 4.0 # setting a higher value for tx s_init = state.extract_state([]) # type(state) # s0=state_repos.get_state(0) # s0.state_vector # state.apply_state(s0, []) calib.run(state=s_init) # rerun the model, with new parameter disch_sim_p_high = calib.region_model.statistics.discharge(target_obs.catchment_indexes) # fetch discharge ts parameters.gs.tx = -4.0 # setting a higher value for tx calib.run(state=s_init) # rerun the model, with new parameter disch_sim_p_low = calib.region_model.statistics.discharge(target_obs.catchment_indexes) # fetch discharge ts fig, ax = plt.subplots(1, figsize=(15,10)) ax.plot(ts_timestamps_sim, disch_sim.values, lw=2, label = "calib") ax.plot(ts_timestamps_sim, disch_sim_p_high.values, lw=2, label = "high") ax.plot(ts_timestamps_sim, disch_sim_p_low.values, lw=2, label = "low") ax.plot(ts_timestamps_obs, disch_obs, lw=2, ls='--', label = "obs") ax.set_title("investigating parameter gs.tx") ax.legend() ax.set_ylabel("discharge [m3 s-1]") .. image:: snow_or_rain.png .. code-block:: python s_init = state.extract_state([]) # reset the max water parameter parameters.gs.max_water = 1.0 # setting a higher value for tx calib.run(state=s_init) # rerun the model, with new parameter disch_sim_p_high = calib.region_model.statistics.discharge(target_obs.catchment_indexes) # fetch discharge ts parameters.gs.max_water = .001 # setting a higher value for tx calib.run(state=s_init) # rerun the model, with new parameter disch_sim_p_low = calib.region_model.statistics.discharge(target_obs.catchment_indexes) # fetch discharge ts # plot the results fig, ax = plt.subplots(1, figsize=(15,10)) ax.plot(ts_timestamps_sim, disch_sim.values, lw=2, label = "calib") ax.plot(ts_timestamps_sim, disch_sim_p_high.values, lw=2, label = "high") ax.plot(ts_timestamps_sim, disch_sim_p_low.values, lw=2, label = "low") ax.plot(ts_timestamps_obs, disch_obs, lw=2, ls='--', label = "obs") ax.set_title("investigating parameter gs.max_water") ax.legend() ax.set_ylabel("discharge [m3 s-1]") .. image:: snow_or_rain_2.png