Notes from plot_thresholds_all() code review on 2015-03-10

gistfile1.py

from collections import OrderedDict

# Here's a possible data structure to hold constant values for the sites
# that we report ssh for. The outer dict has 2 keys: DFO and NOAA. Each of
# those keys has an OrderedDict as its value. OrderedDicts are defined by
# lists of key-value 2-tuples. Here the keys of the DFO OrderedDict are the
# names of the DFO sites, and the values are a dict of constant values.
# So, iterating over SITES['DFO'].items() will return the constants for
# Point Atkinson, Campbell River, and Victoria, in that order. The Campbell
# River extreme_ssh is accessed by SITES['DFO']['Campbell River']['extreme_ssh'].
SITES = {
    'DFO': OrderedDict([
        ('Point Atkinson', {
            'lat': -123.1,
            'lon': 49.4,
            'stn_number': 7795,
            'extreme_ssh': 5.61,
        }),
        ('Campbell River', {
            'lat': -125.1,
            'lon': 50.1,
            'stn_number': 8074,
            'extreme_ssh': 5.35,
        }),
        ('Victoria', {
            'lat': -124,
            'lon': 48.43,
            'stn_number': 7120,
            'extreme_ssh': 3.76,
        }),
    ]),
    'NOAA': OrderedDict([
    ]),
}


# This is the function to plot the small square map that is used to show
# station locations. It is used in 3 figure functions with slightly different
# x/y limits in 1 instance, and a different title in that same instance.
# So, we can use default argument values to make 2 of the calls of this function
# somewhat succinct. In a slight departure from our discussion I chose to
# make the map title a required arg so that the calls will be a little more
# self-documenting. The function is marked as private with a leading
# underscore on its name because it is not intended to be useful outside of
# the figures module.
def _plot_stations_map(
    ax, grid_B, PNW_coastline, title, xlim=(-125.4, -122.2), ylim=(48, 50.3)
):
    """The docstring that a less lazy person would write...
    """
    plot_map(ax, grid_B, PNW_coastline, coastline='full', fill=1, domain=0)
    ax0.set_xlim(xlim)
    ax0.set_ylim(ylim)
    ax0.set_title(title, **title_font)
    axis_colors(ax0, 'gray')
# In plot_thresholds_all() we'll call this function like:
#     _plot_stations_map(
#         ax0, grid_B, title='Degree of Flood Risk',
#         xlim=(-125.4, -122.2), ylim=(48, 50.3))
# and in Sandheads_winds() and compare_water_levels() the calls will looke
# like:
#     _plot_stations_map(ax0, grid_B, PNW_coastline, title='Station Locations')


# We had a sideline discussion about converting arrow objects to datetime
# objects that matplotlib can work with. A list comprehension is a nice
# Pythonic way of doing that. If arrows is a list of arrow objects:
#     [a.datetime for a in arrows]
# will give a list of datetime objects that have tzinfo attributes, and
#     [a.naive for a in arrows]
# will give a list of datetime objects without tzinfo atrributes.


# Comments marked with ## in the function are explanatory

def plot_thresholds_all(
    grid_T, grid_B, model_path, PNW_coastline, PST=1, MSL=1,
    figsize=(20, 15.5),
):
    """Plots sea surface height over one day with respect to warning
    thresholds.

    This function applies only to Point Atkinson, Campbell River, and Victoria.
    There are three different warning thresholds.
    The locations of stations are colored depending on the threshold in
    which they fall: green, yellow, red.

    :arg grid_T: Hourly tracer results dataset from NEMO.
    :type grid_T: :class:`netCDF4.Dataset`

    :arg grid_B: Bathymetry dataset for the Salish Sea NEMO model.
    :type grid_B: :class:`netCDF4.Dataset`

    :arg model_path: The directory where the model wind files are stored.
    :type model_path: string

    :arg PST: Specifies if plot should be presented in PST.
              1 = plot in PST, 0 = plot in UTC.
    :type PST: 0 or 1

    :arg MSL: Specifies if the plot should be centred about mean sea level.
              1=centre about MSL, 0=centre about 0.
    :type MSL: 0 or 1

    :arg figsize:  Figure size (width, height) in inches.
    :type figsize: 2-tuple

    :returns: matplotlib figure object instance (fig).
    """
    ## plt.figure() accepts a facecolor arg
    fig = plt.figure(figsize=figsize, facecolor='#2B3E50')
    # fig.patch.set_facecolor('#2B3E50') - not needed
    gs = gridspec.GridSpec(3, 2, width_ratios=[1.5, 1])
    gs.update(wspace=0.13, hspace=0.2)

    # Map of region
    ax0 = fig.add_subplot(gs[:, 1])
    _plot_stations_map(
        ax0, grid_B, title='Degree of Flood Risk',
        xlim=(-125.4, -122.2), ylim=(48, 50.3))
    ## Next 6 lines are replaced by _plot_stations_map() call above
    # plot_map(ax0, grid_B, PNW_coastline, 'full', 1, 0)
    # ax0.set_xlim(-125.4, -122.2)
    # ax0.set_ylim(48, 50.3)
    # ax0.set_title('Degree of Flood Risk', **title_font)
    # axis_colors(ax0, 'gray')

    ## Provided by the SITES module constant
    # Stations information
    # lats, lons = station_coords() - not needed

    # Bathymetry - redundant comment?
    bathy, X, Y = tidetools.get_bathy_data(grid_B)
    ## Moved time range and timezone indicator out of loop
    # Time range - redundant comment?
    t_orig, t_final, t = get_model_time_variables(grid_T)
    ## Cleaner way of choosing timezone indicator
    tzone = '[PST]' if PST else '[UTC]'
    ## Next line replaced by above
    # tzone = PST * '[PST]' + abs((PST - 1)) * '[UTC]'
    ## Calculate the time shift in the same way
    t_shift = time_shift if PST else 0

    ## Next 3 lines are provided by the SITES module constant
    # m = np.arange(3)
    # names = ['Point Atkinson', 'Campbell River', 'Victoria']
    # extreme_sshs = [5.61, 5.35, 3.76]

    # for M, name, extreme_ssh in zip(m, names, extreme_sshs): - replaced
    for M, (name, site_values) in enumerate(SITES['DFO'].items()):
        # Get sea surface height
        j, i = tidetools.find_closest_model_point(
            site_values['lon'], site_values['lat'], X, Y, bathy,
            allow_land=False)
        ## Next 2 lines replaced by above
        # j, i = tidetools.find_closest_model_point(
        #     lons[name], lats[name], X, Y, bathy, allow_land=False)
        ## Calculate (possibly corrected) ssh at site all in 1 place
        ssh_loc = grid_T.variables['sossheig'][:, j, i]
        ssh_loc += MSL_DATUMS[name] if MSL else 0
        ## Next 2 lines replaced by above
        # ssh = grid_T.variables['sossheig']
        # ssh_loc = ssh[:, j, i]

        # Sea surface height plot
        ax = fig.add_subplot(gs[M, 0])
        # Plot tides, corrected model and original model
        if name == 'Point Atkinson':
            plot_PA_observations(ax, PST)
        ttide = plot_tides(ax, name, PST, MSL)
        ## This might need to change since ssh_loc now may include MSL_DATUMS
        ssh_corr = plot_corrected_model(
            ax, t, ssh_loc, ttide, PST, MSL, MSL_DATUMS[name])
        ## Simplified thanks to t_shift & ssh_loc calculations above
        ax.plot(
            t + t_shift, ssh_loc,
            '--', c=model_c, linewidth=1, label='Model')
        ## Replaced by above call
        # ax.plot(
        #     t + PST * time_shift, ssh_loc + MSL_DATUMS[name] * MSL,
        #     '--', c=model_c, linewidth=1, label='Model')

        # Axis
        ax.set_xlim(t_orig + PST * time_shift, t_final + PST * time_shift)
        ax.set_ylim([-1, 6])
        ax.set_title(
            'Hourly Sea Surface Height at {name}: {t_orig:%d-%b-%Y}'
            .format(name=name, t_orig=t_orig),
            **title_font)
        ax.set_xlabel('Time {}'.format(tzone), **axis_font)
        ax.set_ylabel('Water Level above Chart Datum (m)', **axis_font)
        ax.grid()
        axis_colors(ax, 'gray')
        ax.xaxis.set_major_formatter(hfmt)
        fig.autofmt_xdate()

        # Map
        bbox_args = dict(boxstyle='square', facecolor='white', alpha=0.3)
        ax0.annotate('Point Atkinson', (-123.1, 49.4),
                     fontsize=15, color='black', bbox=bbox_args)
        ax0.annotate('Campbell River', (-125.1, 50.1),
                     fontsize=15, color='black', bbox=bbox_args)
        ax0.annotate('Victoria', (-124, 48.43),
                     fontsize=15, color='black', bbox=bbox_args)

        # Define thresholds in sea surface height plots
        # [max_tides, mid_tides, extreme_ssh] = plot_threshold_map(
        max_tides, mid_tides, _ = plot_threshold_map(
            ax0, ttide, ssh_corr, 'D', 10, 1.0, name)

        # Plot thresholds in sea surface height plots
        ax.axhline(
            y=max_tides,
            color='Gold', linewidth=2, ls='solid', label='Maximum tides')
        ax.axhline(
            y=mid_tides,
            color='Red', linewidth=2, ls='solid', label='Extreme water')
        ax.axhline(
            y=extreme_ssh,
            color='DarkRed', linewidth=2, ls='solid',
            label='Historical maximum')

        # Legend
        if M == 0:
            legend = ax.legend(
                bbox_to_anchor=(1.285, 1), loc=2, borderaxespad=0.,
                prop={'size': 15}, title=r'Legend')
            legend.get_title().set_fontsize('20')

    # Citation
    ax0.text(0.03, -0.45,
             'Tidal predictions calculated with t_tide: '
             'http://www.eos.ubc.ca/~rich/#T_Tide \n'
             'Observed water levels from Fisheries and Oceans, Canada \n'
             'via Scott Tinis at stormsurgebc.ca',
             horizontalalignment='left',
             verticalalignment='top',
             transform=ax0.transAxes, color='white')

    return fig