Calibrating the MOHYSE hydrological model using OSTRICH on the Raven server

Here we use birdy’s WPS client to calibrate the MOHYSE hydrological model on the server and analyze the calibrated parameter set and hydrograph.

from birdy import WPSClient

from example_data import TESTDATA
import datetime as dt
from urllib.request import urlretrieve
import xarray as xr
import numpy as np
from matplotlib import pyplot as plt
import os

# Set environment variable WPS_URL to "http://localhost:9099" to run on the default local server
url = os.environ.get("WPS_URL", "https://pavics.ouranos.ca/twitcher/ows/proxy/raven/wps")
wps = WPSClient(url)

# The model parameter boundaries. Can either be a string of comma separated values, a list, an array or a named tuple.
low_p = '0.01, 0.01, 0.01, -5.00, 0.01, 0.01, 0.01, 0.01'
high_p = '20.0, 1.0, 20.0, 5.0, 0.5, 1.0, 1.0, 1.0'
low_h = '0.01, 0.01'
high_h = '15.0, 15.0'

# Forcing files
ts=TESTDATA['ostrich-mohyse-nc-ts']

# OSTRICH configuration parameters
config = dict(
    algorithm='DDS',
    max_iterations=10,
    area=4250.6,
    elevation=843.0,
    latitude=54.4848,
    longitude=-123.3659,
    lowerbounds=low_p,
    upperbounds=high_p,
    start_date=dt.datetime(1954, 1, 1),
    duration=208,
    hruslowerbounds=low_h,
    hrusupperbounds=high_h,
    # Comment out the random seed to show different results!
    random_seed=6.67408*10**-11,
    suppress_output=False,
    )

# Let's call Ostrich with the timeseries, calibration parameters and other configuration parameters
resp = wps.ostrich_mohyse(ts=str(ts), **config)

# And get the response
# With `asobj` set to False, only the reference to the output is returned in the response.
# Setting `asobj` to True will retrieve the actual files and copy the locally.
[calibration, hydrograph, storage, solution, diagnostics, calibparams, rv2] = resp.get(asobj=True)

Since we requested output objects, we can simply access the output objects. The dianostics is just a CSV file: # Here there is a problem because we only get the 8 calibrated parameters but not the 2 “hru” calibrated parameters, which means we can’t run the code afterwards with the calibrated HRU parameters.

print(calibparams)

7.721801, 0.8551484, 17.74571, 1.627677, 0.0770245, 0.94096, 0.6941596, 0.820787

The hydrograph and storage outputs are netCDF files storing the time series. These files are opened by default using xarray, which provides convenient and powerful time series analysis and plotting tools.

hydrograph.q_sim

<xarray.DataArray 'q_sim' (time: 209, nbasins: 1)>
array([[ 0.      ],
       [ 0.      ],
       [ 0.      ],
       ...,
       [17.875367],
       [17.30535 ],
       [16.747391]])
Coordinates:
  * time        (time) datetime64[ns] 1954-01-01 1954-01-02 ... 1954-07-28
    basin_name  (nbasins) object ...
Dimensions without coordinates: nbasins
Attributes:
    units:      m**3 s**-1
    long_name:  Simulated outflows

from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()

hydrograph.q_sim.plot()

[<matplotlib.lines.Line2D at 0x7fe70d6b0ac8>]

print("Max: ", hydrograph.q_sim.max())
print("Mean: ", hydrograph.q_sim.mean())
print("Monthly means: ", hydrograph.q_sim.groupby('time.month').mean(dim='time'))

Max:  <xarray.DataArray 'q_sim' ()>
array(61.20512771)
Mean:  <xarray.DataArray 'q_sim' ()>
array(26.91468183)
Monthly means:  <xarray.DataArray 'q_sim' (month: 7, nbasins: 1)>
array([[ 0.        ],
       [ 0.26801909],
       [11.08674251],
       [40.6703697 ],
       [59.69526693],
       [48.9501015 ],
       [26.24312637]])
Coordinates:
    basin_name  (nbasins) object ...
  * month       (month) int64 1 2 3 4 5 6 7
Dimensions without coordinates: nbasins

Now, let’s do a test: let’s run the model again using the same parameters and check to see that the NSE is the same.

First, lets extract the diagnostics from the optimized run, giving the calibration NSE.

diagnostics

'observed data series,filename,DIAG_NASH_SUTCLIFFE,DIAG_RMSE,\nHYDROGRAPH,/tmp/pywps_process_0ohz5jwm/Salmon-River-Near-Prince-George_meteo_daily.nc,0.382681,40.7086,\n'

Run the MOHYSE WPS:

# Model configuration parameters
config = dict(area=4250.6,
    elevation=843.0,
    latitude=54.4848,
    longitude=-123.3659,
    start_date=dt.datetime(1954, 1, 1),
    duration=208,
             )
# Let's call the model with the timeseries, model parameters and other configuration parameters
resp = wps.raven_mohyse(ts=str(ts), params=calibparams, **config)
[hydrograph, storage, solution, diagnostics2, rv2] = resp.get(asobj=True)

Now lets check to see if the diagnostics are the same.

diagnostics2

'observed data series,filename,DIAG_NASH_SUTCLIFFE,DIAG_RMSE,\nHYDROGRAPH,/tmp/pywps_process_x3jsi9es/Salmon-River-Near-Prince-George_meteo_daily.nc,-1.81387,86.9127,\n'