Synthetic Hydrograph Comparisons

hydrographs.py

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

x = np.arange(0, np.pi * 12, 0.025)

q_insitu = 2 * np.sin(0.5 * x) + 4
q_constant = np.ones_like(x) * 4
q_biased_high = 2 * np.sin(0.5 * x) + 4 + 1.5
q_biased_low = 2 * np.sin(0.5 * x) + 4 - 1.5
q_inverse = 2 * np.sin(0.5 * x + np.pi) + 4
q_bad_timing = 2 * np.sin(0.5 * x - 1) + 4
q_higher_frequency = 3 * np.sin(2 * x) + 4
q_nonstationary_mean = 2 * np.sin(5 * x) + 4 + np.sin(x)

df = pd.DataFrame({
    'q_insitu': q_insitu,
    'q_constant': q_constant,
    'q_biased_high': q_biased_high,
    'q_biased_low': q_biased_low,
    'q_inverse': q_inverse,
    'q_bad_timing': q_bad_timing,
    'q_higher_frequency': q_higher_frequency,
    'q_nonstationary_mean': q_nonstationary_mean,
})


def kge2012(simulated_array, observed_array):
    # Means
    sim_mean = np.mean(simulated_array)
    obs_mean = np.mean(observed_array)

    # Standard Deviations
    sim_sigma = np.std(simulated_array)
    obs_sigma = np.std(observed_array)

    # Pearson R
    top_pr = np.sum((observed_array - obs_mean) * (simulated_array - sim_mean))
    bot1_pr = np.sqrt(np.sum((observed_array - obs_mean) ** 2))
    bot2_pr = np.sqrt(np.sum((simulated_array - sim_mean) ** 2))
    pr = top_pr / (bot1_pr * bot2_pr)

    # Ratio between mean of simulated and observed data
    beta = sim_mean / obs_mean

    # CV is the coefficient of variation (standard deviation / mean)
    sim_cv = sim_sigma / sim_mean
    obs_cv = obs_sigma / obs_mean

    # Variability Ratio, or the ratio of simulated CV to observed CV
    gam = sim_cv / obs_cv

    if obs_mean == 0 or obs_sigma == 0 or sim_mean == 0:
        return np.nan

    kge = 1 - np.sqrt((pr - 1) ** 2 + (gam - 1) ** 2 + (beta - 1) ** 2)
    return kge


for qcol in df.columns.drop('q_insitu'):
    fig, ax = plt.subplots(tight_layout=True, figsize=(6, 4))
    (
        df
        [['q_insitu', qcol]]
        .rename(columns={'q_insitu': 'in-situ', qcol: 'model'})
        .plot(figsize=(7, 4), ax=ax, legend=True, color=['k', 'r'])
    )

    # calculate mean error, rmse, R, kling gupta, and nash sutcliffe
    me = np.mean(df['q_insitu'] - df[qcol])
    mse = np.mean((df['q_insitu'] - df[qcol]) ** 2)
    rmse = np.sqrt(np.mean((df['q_insitu'] - df[qcol]) ** 2))
    r = np.corrcoef(df['q_insitu'], df[qcol])[0, 1]
    kg = kge2012(df[qcol], df['q_insitu'])
    ns = 1 - (np.sum((df['q_insitu'] - df[qcol]) ** 2) / np.sum((df['q_insitu'] - np.mean(df['q_insitu'])) ** 2))
    volume_insitu = np.sum(df['q_insitu'])
    volume_qcol = np.sum(df[qcol])
    std_insitu = np.std(df['q_insitu'])
    std_qcol = np.std(df[qcol])
    if np.isnan(r):
        r = 0

    ax.set_ylabel('Discharge (m3/s)')
    ax.set_xlabel('Time')

    ax.set_ylim(0, 8)
    ax.set_xticks([])

    fig.savefig(f'./hydrograph_comparison_nometrics_{qcol}.png', dpi=500)
    ax.set_title(
        f'ME: {round(me, 2)} MSE: {round(mse, 2)} RMSE: {round(rmse, 2)} R: {round(r, 2)} KG: {round(kg, 2)} NS: {round(ns, 2)}' +
        f'\nVolume Insitu: {round(volume_insitu, 2)} Volume Model: {round(volume_qcol, 2)}'
        f'\nStd Insitu: {round(std_insitu, 2)} Std Model: {round(std_qcol, 2)}'
    )
    fig.savefig(f'./hydrograph_comparison_withmetrics_{qcol}.png', dpi=500)

    plt.close(fig)