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