Computing Objective Functions on the Raven server

Here we use birdy’s WPS client to compare simulated and observed discharge using various objective functions on the server.

import os
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

# 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)

First, we run GR4JCN to generate a NetCDF file containing q_obs and q_sim. This file contains q_obs from the initial dataset, plus the simulated discharge q_sim from the model.

# The model parameters. Can either be a string of comma separated values, a list, an array or a named tuple.
params = '0.529, -3.396, 407.29, 1.072, 16.9, 0.947'

# Forcing files
ts=TESTDATA['raven-gr4j-cemaneige-nc-ts']

# Model configuration parameters
config = dict(
    start_date=dt.datetime(2000, 1, 1),
    end_date=dt.datetime(2002, 1, 1),
    area=4250.6,
    elevation=843.0,
    latitude=54.4848,
    longitude=-123.3659,
    )

# Let's call the model
resp = wps.raven_gr4j_cemaneige(ts=str(ts), params = params, **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.
# Here we use false, as we do not want the file, we only want the path to the file.
[hydrograph, storage, solution, diagnostics, rv] = resp.get(asobj=False)

Now we can call the objective function WPS service on the Raven Server and collect the results using ‘asobj=True’ to collect the actual results in a dict format.

resp = wps.objective_function(obs=hydrograph, sim=hydrograph)
results=resp.get(asobj=True)

Let’s see the results!

print(results.metrics)

{'agreementindex': 0.4596628392956652, 'bias': 3.264414783955057, 'correlationcoefficient': 0.20692429176680854, 'covariance': 126.58541412759337, 'decomposed_mse': 1336.777456482738, 'kge': 0.041555657877441976, 'log_p': -86.28634365686393, 'lognashsutcliffe': -1.7249359884422084, 'mae': 24.444255724443455, 'mse': 1336.777456482738, 'nashsutcliffe': -0.037104793191898855, 'pbias': -11.666398787696835, 'rmse': 36.56196734973021, 'rrmse': 1.3066553112711805, 'rsquared': 0.04281766252319531, 'rsr': 1.018383421502873, 'volume_error': -0.11666398787696836}

There are 17 objective functions. We can also call the WPS server with any one of those objective functions:

resp = wps.objective_function(obs=hydrograph, sim=hydrograph, name='agreementindex')
results=resp.get(asobj=True)
print(results.metrics)

{'agreementindex': 0.4596628392956652}

Now that we have numerical values for the goodness-of-fit, let’s perform a visual inspection.

# Rerun the model but this time get the object (asobj=True)
resp2 = wps.raven_gr4j_cemaneige(ts=str(ts), params = params, **config)
[hydrograph, storage, solution, diagnostics, rv] = resp2.get(asobj=False)

# Call the objective function graphing wps while sending the hydrograph as input, which includes 'q_sim' and 'q_obs'
resp3=wps.graph_objective_function_fit(sims=hydrograph)
[graph_objfun_fit, graph_objfun_annual_fit] = resp3.get(asobj=True)

Now that we have run the graphing service, display the plots.

# Graph of full time series
graph_objfun_fit