C. Plotting Code Samples#

C.1. hymod.ipynb#

The following are the plotting functions as described in the hymod.ipynb Jupyter notebook tutorial.

The following are the necessary package imports to run these functions:

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

from matplotlib.lines import Line2D

C.1.1. plot_observed_vs_simulated_streamflow()#

def plot_observed_vs_simulated_streamflow(df, hymod_dict, figsize=[12, 6]):
    """Plot observed versus simulated streamflow.

    :param df:              Dataframe of hymod input data including columns for precip, potential evapotranspiration,
                            and streamflow

    :param hymod_dict:      A dictionary of hymod outputs
    :type hymod_dict:       dict

    :param figsize:         Matplotlib figure size
    :type figsize:          list

    """

    # set plot style
    plt.style.use('seaborn-white')

    # set up figure
    fig, ax = plt.subplots(figsize=figsize)

    # plot observed streamflow
    ax.plot(range(0, len(df['Strmflw'])), df['Strmflw'], color='pink')

    # plot simulated streamflow
    ax.plot(range(0, len(df['Strmflw'])), hymod_dict['Q'], color='black')

    # set axis labels
    ax.set_ylabel('Streamflow($m^3/s$)')
    ax.set_xlabel('Days')

    # set plot title
    plt.title('Observed vs. Simulated Streamflow')

    return ax

C.1.2. plot_observed_vs_sensitivity_streamflow()#

def plot_observed_vs_sensitivity_streamflow(df_obs, df_sim, figsize=[10, 4]):
    """Plot observed streamflow versus simulations generated from sensitivity analysis.

    :param df_obs:          Dataframe of mean monthly hymod input data including columns for precip,
                            potential evapotranspiration, and streamflow

    :param df_sim:          Dataframe of mean monthly simulation data from sensitivity analysis

    :param figsize:         Matplotlib figure size
    :type figsize:          list

    """

    month_list = range(len(df_sim))

    # set up figure
    fig, ax = plt.subplots(figsize=figsize)

    # set labels
    ax.set_xlabel('Days')
    ax.set_ylabel('Flow Discharge (m^3/s)')

    # plots all simulated streamflow cases under different sample sets
    for i in df_sim.columns:
        plt.plot(month_list, df_sim[i], color="pink", alpha=0.2)

    # plot observed streamflow
    plt.plot(month_list, df_obs['Strmflw'], color="black")

    plt.title('Observed vs. Sensitivity Analysis Outputs')

    return ax

C.1.3. plot_monthly_heatmap()#

def plot_monthly_heatmap(arr_sim, df_obs, title='', figsize=[14, 6]):
    """Plot a sensitivity metric overlain by observed flow.

    :param arr_sim:         Numpy array of simulated metrics

    :param df_obs:          Dataframe of mean monthly observed data from sensitivity analysis

    :param title:           Title of plot
    :type title:            str

    :param figsize:         Matplotlib figure size
    :type figsize:          list

    """

    # set up figure
    fig, ax = plt.subplots(figsize=figsize)

    # plot heatmap
    sns.heatmap(arr_sim,
                ax=ax,
                yticklabels=['Kq', 'Ks', 'Alp', 'Huz', 'B'],
                cmap=sns.color_palette("ch:s=-.2,r=.6"))

    # setup overlay axis
    ax2 = ax.twinx()

    # plot line
    ax2.plot(np.arange(0.5, 12.5), df_obs['Strmflw'], color='slateblue')

    # plot points on line
    ax2.plot(np.arange(0.5, 12.5), df_obs['Strmflw'], color='slateblue', marker='o')

    # set axis limits and labels
    ax.set_ylim(0, 5)
    ax.set_xlim(0, 12)
    ax.set_xticklabels(['jan', 'feb', 'mar', 'apr', 'may', 'jun', 'jul', 'aug', 'sep', 'oct', 'nov', 'dec'])
    ax2.set_ylabel('Flow Discharge($m^3/s$)')

    plt.title(title)

    plt.show()

    return ax, ax2

C.1.4. plot_annual_heatmap()#

def plot_annual_heatmap(arr_sim, df_obs, title='', figsize=[14,5]):
    """Plot a sensitivity metric overlain by observed flow..

    :param arr_sim:         Numpy array of simulated metrics

    :param df_obs:          Dataframe of mean monthly observed data from sensitivity analysis

    :param title:           Title of plot
    :type title:            str

    :param figsize:         Matplotlib figure size
    :type figsize:          list

    """

    # set up figure
    fig, ax = plt.subplots(figsize=figsize)

    # plot heatmap
    sns.heatmap(arr_sim, ax=ax, cmap=sns.color_palette("YlOrBr"))

    # setup overlay axis
    ax2 = ax.twinx()

    # plot line
    ax2.plot(np.arange(0.5, 10.5), df_obs['Strmflw'], color='slateblue')

    # plot points on line
    ax2.plot(np.arange(0.5, 10.5), df_obs['Strmflw'], color='slateblue', marker='o')

    # set up axis lables and limits
    ax.set_ylim(0, 5)
    ax.set_xlim(0, 10)
    ax.set_yticklabels(['Kq', 'Ks', 'Alp', 'Huz', 'B'])
    ax.set_xticklabels(range(2000, 2010))
    ax2.set_ylabel('Flow Discharge($m^3/s$)')

    plt.title(title)

    return ax, ax2

C.1.5. plot_varying_heatmap()#

def plot_varying_heatmap(arr_sim, df_obs, title='', figsize=[14,5]):
    """Plot a sensitivity metric overlain by observed flow..

    :param arr_sim:         Numpy array of simulated metrics

    :param df_obs:          Dataframe of mean monthly observed data from sensitivity analysis

    :param title:           Title of plot
    :type title:            str

    :param figsize:         Matplotlib figure size
    :type figsize:          list

    """

    # set up figure
    fig, ax = plt.subplots(figsize=figsize)

    # plot heatmap
    sns.heatmap(arr_sim,
                ax=ax,
                yticklabels=['Kq', 'Ks', 'Alp', 'Huz', 'B'],
                cmap=sns.light_palette("seagreen", as_cmap=True))

    n_years = df_obs.shape[0]

    # setup overlay axis
    ax2 = ax.twinx()

    # plot line
    ax2.plot(range(0, n_years), df_obs['Strmflw'], color='slateblue')

    # plot points on line
    ax2.plot(range(0, n_years), df_obs['Strmflw'], color='slateblue', marker='o')

    # set up axis lables and limits
    ax.set_ylim(0, 5)
    ax.set_xlim(-0.5, 119.5)
    ax2.set_ylabel('Flow Discharge')
    ax.set_xlabel('Number of Months')

    plt.title(title)

    return ax, ax2

C.1.6. plot_precalibration_flow()#

def plot_precalibration_flow(df_sim, df_obs, figsize=[10, 4]):
    """Plot flow discharge provided by the ensemble of parameters sets from Pre-Calibration versus the observed
    flow data.

    :param df_sim:          Dataframe of simulated metrics

    :param df_obs:          Dataframe of mean monthly observed data from sensitivity analysis

    :param figsize:         Matplotlib figure size
    :type figsize:          list

    """

    # set up figure
    fig, ax = plt.subplots(figsize=figsize)

    # set axis labels
    ax.set_xlabel('Days')
    ax.set_ylabel('Flow Discharge')

    # plot pre-calibration results
    for i in range(df_sim.shape[1]):
        plt.plot(range(len(df_sim)), df_sim.iloc[:, i],  color="lightgreen", alpha=0.2)

    # plot observed
    plt.plot(range(len(df_sim)), df_obs['Strmflw'],  color="black")

    plt.title('Observed vs. Pre-Calibration Outputs')

    # customize legend
    custom_lines = [Line2D([0], [0],  color="lightgreen", lw=4),
                    Line2D([0], [0], color="black", lw=4)]
    plt.legend(custom_lines, ['Pre-Calibration', 'Observed'])

    return ax

C.1.7. plot_precalibration_glue()#

def plot_precalibration_glue(df_precal, df_glue, df_obs, figsize=[10, 4]):
    """Plot flow discharge provided by the ensemble of parameters sets from Pre-Calibration versus the observed
    flow data.

    :param df_sim:          Dataframe of simulated metrics

    :param df_obs:          Dataframe of mean monthly observed data from sensitivity analysis

    :param figsize:         Matplotlib figure size
    :type figsize:          list

    """

    # set up figure
    fig, ax = plt.subplots(figsize=figsize)

    # set axis labels
    ax.set_xlabel('Days')
    ax.set_ylabel('Flow Discharge')

    # plot pre-calibration results
    for i in range(df_precal.shape[1]):
        plt.plot(range(len(df_precal)), df_precal.iloc[:, i],  color="lightgreen", alpha=0.2)

    # plot glue
    for i in range(df_glue.shape[1]):
        plt.plot(range(len(df_glue)), df_glue.iloc[:, i], color="lightblue", alpha=0.2)

    # plot observed
    plt.plot(range(len(df_precal)), df_obs['Strmflw'],  color="black")

    plt.title('Observed vs. Sensitivity Analysis Outputs across GLUE/Pre-Calibration')

    # customize legend
    custom_lines = [Line2D([0], [0],  color="lightgreen", lw=4),
                    Line2D([0], [0], color="lightblue", lw=4),
                    Line2D([0], [0], color="black", lw=4)]
    plt.legend(custom_lines, ['Pre-Calibration', 'GLUE', 'Observed'])

    return ax

C.2. fishery_dynamics.ipynb#

The following are the plotting functions as described in the fishery_dynamics.ipynb Jupyter notebook tutorial.

The following are the necessary package imports to run these functions:

import numpy as np
import matplotlib.pyplot as plt

from matplotlib import patheffects as pe

C.2.1. plot_objective_performance()#

def plot_objective_performance(objective_performance, profit_solution, robust_solution, figsize=(18, 9)):
    """Plot the identified solutions with regards to their objective performance
    in a parallel axis plot

    :param objective_performance:               Objective performance array
    :param profit_solution:                     Profitable solutions array
    :param robust_solution:                     Robust solutions array
    :param figsize:                             Figure size
    :type figsize:                              tuple

    """

    # create the figure object
    fig = plt.figure(figsize=figsize)

    # set up subplot axis object
    ax = fig.add_subplot(1, 1, 1)

    # labels where constraint is always 0
    objs_labels = ['Net present\nvalue (NPV)',
                   'Prey population deficit',
                   'Longest duration\nof low harvest',
                   'Worst harvest instance',
                   'Variance of harvest',
                   'Duration of predator\npopulation collapse']

    # normalization across objectives
    mins = objective_performance.min(axis=0)
    maxs = objective_performance.max(axis=0)
    norm_reference = objective_performance.copy()

    for i in range(5):
        mm = objective_performance[:, i].min()
        mx = objective_performance[:, i].max()
        if mm != mx:
            norm_reference[:, i] = (objective_performance[:, i] - mm) / (mx - mm)
        else:
            norm_reference[:, i] = 1

    # colormap from matplotlib
    cmap = plt.cm.get_cmap("Blues")

    # plot all solutions
    for i in range(len(norm_reference[:, 0])):
        ys = np.append(norm_reference[i, :], 1.0)
        xs = range(len(ys))
        ax.plot(xs, ys, c=cmap(ys[0]), linewidth=2)

    # to highlight robust solutions
    ys = np.append(norm_reference[profit_solution, :], 1.0)  # Most profitable
    xs = range(len(ys))
    l1 = ax.plot(xs[0:6],
                 ys[0:6],
                 c=cmap(ys[0]),
                 linewidth=3,
                 label='Most robust in NPV',
                 path_effects=[pe.Stroke(linewidth=6, foreground='darkgoldenrod'), pe.Normal()])

    ys = np.append(norm_reference[robust_solution, :], 1.0)  # Most robust in all criteria
    xs = range(len(ys))
    l2 = ax.plot(xs[0:6],
                 ys[0:6],
                 c=cmap(ys[0]),
                 linewidth=3,
                 label='Most robust across criteria',
                 path_effects=[pe.Stroke(linewidth=6, foreground='gold'), pe.Normal()])

    # build colorbar
    sm = plt.cm.ScalarMappable(cmap=cmap)
    sm.set_array([objective_performance[:, 0].min(), objective_performance[:, 0].max()])
    cbar = fig.colorbar(sm)
    cbar.ax.set_ylabel("\nNet present value (NPV)")

    # tick values
    minvalues = ["{0:.3f}".format(mins[0]),
                 "{0:.3f}".format(-mins[1]),
                 str(-mins[2]),
                 "{0:.3f}".format(-mins[3]),
                 "{0:.2f}".format(-mins[4]),
                 str(0)]

    maxvalues = ["{0:.2f}".format(maxs[0]),
                 "{0:.3f}".format(-maxs[1]),
                 str(-maxs[2]),
                 "{0:.2f}".format(maxs[3]),
                 "{0:.2f}".format(-maxs[4]),
                 str(0)]

    ax.set_ylabel("Preference ->", size=12)
    ax.set_yticks([])
    ax.set_xticks([0, 1, 2, 3, 4, 5])
    ax.set_xticklabels([minvalues[i] + '\n' + objs_labels[i] for i in range(len(objs_labels))])

    # make a twin axis for toplabels
    ax1 = ax.twiny()
    ax1.set_yticks([])
    ax1.set_xticks([0, 1, 2, 3, 4, 5])
    ax1.set_xticklabels([maxvalues[i] for i in range(len(maxs) + 1)])

    return ax, ax1

C.2.2. plot_factor_performance()#

def plot_factor_performance(param_values, collapse_days, b, m, a):
    """Visualize the performance of our policies in three-dimensional
    parametric space.

    :param param_values:                Saltelli sample array
    :param collapse_days:               Simulation array
    :param b:                           b parameter boundary interval
    :param m:                           m parameter boundary interval
    :param a:                           a parameter boundary interval

    """

    # set colormap
    cmap = plt.cm.get_cmap("RdBu_r")

    # build figure object
    fig = plt.figure(figsize=plt.figaspect(0.5), dpi=600, constrained_layout=True)

    # set up scalable colormap
    sm = plt.cm.ScalarMappable(cmap=cmap)

    # set up subplot for profit maximizing policy
    ax1 = fig.add_subplot(1, 2, 1, projection='3d')

    # add point data for profit plot
    sows = ax1.scatter(param_values[:,1],
                       param_values[:,6],
                       param_values[:,0],
                       c=collapse_days[:,0],
                       cmap=cmap,
                       s=0.5)

    # add surface data for boundary separating successful and failed states of the world
    pts_ineq = ax1.plot_surface(b, m, a, color='black', alpha=0.25, zorder=1)

    # add reference point to plot
    pt_ref = ax1.scatter(0.5, 0.7, 0.005, c='black', s=50, zorder=0)

    # set up plot aesthetics and labels
    ax1.set_xlabel("b")
    ax1.set_ylabel("m")
    ax1.set_zlabel("a")
    ax1.set_zlim([0.0, 2.0])
    ax1.set_xlim([0.0, 1.0])
    ax1.set_ylim([0.0, 1.5])
    ax1.xaxis.set_view_interval(0,  0.5)
    ax1.set_facecolor('white')
    ax1.view_init(12, -17)
    ax1.set_title('Profit maximizing policy')

    # set up subplot for robust policy
    ax2 = fig.add_subplot(1, 2, 2, projection='3d')

    # add point data for robust plot
    sows = ax2.scatter(param_values[:,1],
                       param_values[:,6],
                       param_values[:,0],
                       c=collapse_days[:,1],
                       cmap=cmap,
                       s=0.5)

    # add surface data for boundary separating successful and failed states of the world
    pts_ineq = ax2.plot_surface(b, m, a, color='black', alpha=0.25, zorder=1)

    # add reference point to plot
    pt_ref = ax2.scatter(0.5, 0.7, 0.005, c='black', s=50, zorder=0)

    # set up plot aesthetics and labels
    ax2.set_xlabel("b")
    ax2.set_ylabel("m")
    ax2.set_zlabel("a")
    ax2.set_zlim([0.0, 2.0])
    ax2.set_xlim([0.0, 1.0])
    ax2.set_ylim([0.0, 1.5])
    ax2.xaxis.set_view_interval(0, 0.5)
    ax2.set_facecolor('white')
    ax2.view_init(12, -17)
    ax2.set_title('Robust policy')

    # set up colorbar
    sm.set_array([collapse_days.min(), collapse_days.max()])
    cbar = fig.colorbar(sm)
    cbar.set_label('Days with predator collapse')

    return ax1, ax2