Add meteogram plots

[jrleeman]

It would be helpful to have the ability to generate a meteogram given a set of station data. Possibly similar to the mesonet style.

[dopplershift]

Yeah, I've had the same thought. I actually have the code from MetPy's attic sitting in my home to keep the original meteogram I wrote handy. I'll put it here for as a starting point for anyone wanting to make a cut (it matches the OK mesonet graphics, at least circa 2009):

from datetime import timedelta
from pytz import UTC
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import NullLocator, NullFormatter
from matplotlib.dates import HourLocator, AutoDateLocator, DateFormatter
from metpy.cbook import iterable, get_title

#Default units for certain variables
default_units = {'temperature':'C', 'dewpoint':'C', 'relative humidity':'%',
    'pressure':'mb', 'wind speed':'m/s', 'solar radiation':'$W/m^2$',
    'rainfall':'mm', 'wind gusts':'m/s'}

def _rescale_yaxis(ax, bounds):
    # Manually tweak the limits here to ignore the low bottom set
    # for fill_between
    XY = np.array([[0]*len(bounds), bounds]).T
    ax.ignore_existing_data_limits = True
    ax.update_datalim(XY, updatex=False, updatey=True)
    ax.autoscale_view()

#TODO: REWRITE AS CLASS
def meteogram(data, fig=None, num_panels=3, time_range=None, layout=None,
    styles=None, limits=None, units=None, tz=UTC):
    '''
    Plots a meteogram (collection of time series) for a data set. This
    is broken down into a series of panels (defaults to 3), each of which
    can plot multiple variables, with sensible defaults, but can also be
    specified using *layout*.

    *data* : numpy record array
        A numpy record array containing time series for individual variables
        in each field.

    *fig* : :class:`matplotlib.figure.Figure` instance or None.
        A matplotlib Figure on which to draw.  If None, a new figure
        will be created.

    *num_panels* : int
        The number of panels to use in the plot.

    *time_range* : sequence, datetime.timedetla, or *None*
        If a sequence, the starting and ending times for the x-axis.  If
        a :class:`datetime.timedelta` object, it represents the time span
        to plot, which will end with the last data point.  It defaults to
        the last 24 hours of data.

    *layout* : dictionary
        A dictionary that maps panel numbers to lists of variables.
        If a panel is not found in the dictionary, a default (up to panel 5)
        will be used.  *None* can be included in the list to denote that
        :func:`pyplot.twinx` should be called, and the remaining variables
        plotted.

    *styles* : dictionary
        A dictionary that maps variable names to dictionary of matplotlib
        style arguments.  Also, the keyword `fill` can be included, to
        indicate that a filled plot should be used.  Any variable not
        specified will use a default (if available).

    *limits* : dictionary
        A dictionary that maps variable names to plot limits.  These limits
        are given by tuples with at least two items, which specify the
        start and end limits.  Either can be *None* which implies using the
        automatically determined value.  Optional third and fourth values
        can be given in the tuple, which is a list of tick values and labels,
        respectively.

    *units* : dictionary
        A dictionary that maps variable names to unit strings for axis labels.

    *tz* : datetime.tzinfo instance
        A :class:`datetime.tzinfo instance specifying the timezone to use
        for plotting the x-axis.  See the docs for :module:`datetime` and
        :module:`pytz` for how to construct and use these objects.  The
        default is UTC.

    Returns : list
        A list of the axes objects that were created.
    '''

    if fig is None:
        fig = plt.figure()

    # Get the time variable
    time = data['datetime']

    # Process time_range.
    major_ticker = AutoDateLocator(tz=tz)
    minor_ticker = NullLocator()
    major_formatter = DateFormatter('%H', tz=tz)
    minor_formatter = NullFormatter()
    if time_range is None:
        time_range = timedelta(hours=24)
        major_ticker = HourLocator(byhour=np.arange(0, 25, 3), tz=tz)
        minor_ticker = HourLocator(tz=tz)

    if not iterable(time_range):
        end = time[-1]
        start = end - time_range
        time_range = (start, end)

    #List of variables in each panel.  None denotes that at that point, twinx
    #should be called and the remaining variables plotted on the other axis
    default_layout = {
        0:['temperature', 'dewpoint'],
        1:['wind gusts', 'wind speed', None, 'wind direction'],
        2:['pressure'],
        3:['rainfall'],
        4:['solar radiation']}

    if layout is not None:
        default_layout.update(layout)
    layout = default_layout

    #Default styles for each variable
    default_styles = {
        'relative humidity':dict(color='#255425', linestyle='--'),
        'dewpoint':dict(facecolor='#265425', edgecolor='None', fill=True),
        'temperature':dict(facecolor='#C14F53', edgecolor='None', fill=True),
        'pressure':dict(facecolor='#895125', edgecolor='None', fill=True),
        'dewpoint':dict(facecolor='#265425', edgecolor='None', fill=True),
        'wind speed':dict(facecolor='#1C2386', edgecolor='None', fill=True),
        'wind gusts':dict(facecolor='#8388FC', edgecolor='None', fill=True),
        'wind direction':dict(markeredgecolor='#A9A64B', marker='D',
            linestyle='', markerfacecolor='None', markeredgewidth=1,
            markersize=3),
        'rainfall':dict(facecolor='#37CD37', edgecolor='None', fill=True),
        'solar radiation':dict(facecolor='#FF8529', edgecolor='None',
            fill=True),
        'windchill':dict(color='#8388FC', linewidth=1.5),
        'heat index':dict(color='#671A5C')}

    if styles is not None:
        default_styles.update(styles)
    styles = default_styles

    #Default data limits
    default_limits = {
        'wind direction':(0, 360, np.arange(0,400,45,),
            ['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW', 'N']),
        'wind speed':(0, None),
        'wind gusts':(0, None),
        'rainfall':(0, None),
        'solar radiation':(0, 1000, np.arange(0,1050,200))
        }
    if limits is not None:
        default_limits.update(limits)
    limits = default_limits

    #Set data units
    def_units = default_units.copy()
    if units is not None:
        def_units.update(units)
    units = def_units

    #Get the station name
    site = data['site'][0]

    #Get strings for the start and end times
    start = time_range[0].astimezone(tz).strftime('%H%M %d %b %Y')
    end = time_range[1].astimezone(tz).strftime('%H%M %d %b %Y %Z')

    axes = []
    for panel in range(num_panels):
        if panel > 0:
            ax = fig.add_subplot(num_panels, 1, panel+1, sharex=ax)
        else:
            ax = fig.add_subplot(num_panels, 1, panel+1)
            ax.set_title('%s\n%s to %s' % (site, start, end))

        panel_twinned = False

        var_min = []
        var_max = []
        for varname in layout[panel]:
            if varname is None:
                _rescale_yaxis(ax, var_min + var_max)
                ax = ax.twinx()
                panel_twinned = True
                var_min = []
                var_max = []
                continue

            # Get the linestyle for this variable
            style = styles.get(varname, dict())

            #Get the variable from the data and plot
            var = data[varname]

            #Set special limits if necessary
            lims = limits.get(varname, (None, None))

            #Store the max and min for auto scaling
            if var.max() is not np.ma.masked:
                var_max.append(var.max())
                var_min.append(var.min())

            if style.pop('fill', False):
                #Plot the filled area.  Need latest Matplotlib for date support
                #with fill_betweeen
                lower = -500 if lims[0] is None else lims[0]
                ax.fill_between(time, lower, var, where=~var.mask, **style)
                #TODO: Matplotlib SVN has support for turning off the
                #automatic scaling of the y-axis.  Can we use that somehow
                #to simplify our code??
                _rescale_yaxis(ax, var_min + var_max)
            else:
                ax.plot(time, var, **style)

            #If then length > 2, then we have ticks and (maybe) labels
            if len(lims) > 2:
                other = lims[2:]
                lims = lims[:2]
                #Separate out the ticks and perhaps labels
                if len(other) == 1:
                    ax.set_yticks(other[0])
                else:
                    ticks,labels = other
                    ax.set_yticks(ticks)
                    ax.set_yticklabels(labels)
            ax.set_ylim(*lims)

            # Set the label to the title-cased nice-version from the
            # field info with units, if given.
            if varname in units and units[varname]:
                unit_str = ' (%s)' % units[varname]
                if '^' in unit_str and '$' not in unit_str:
                    unit_str = '$' + unit_str + '$'
            else:
                unit_str = ''

            descr = get_title(data, varname)
            ax.set_ylabel(descr.title() + unit_str)

        ax.xaxis.set_major_locator(major_ticker)
        ax.xaxis.set_major_formatter(major_formatter)
        ax.xaxis.set_minor_locator(minor_ticker)
        ax.xaxis.set_minor_formatter(minor_formatter)
        if not panel_twinned:
            ax.yaxis.set_ticks_position('both')
            for tick in ax.yaxis.get_major_ticks():
                tick.label2On = True
        axes.append(ax)

    # Set the xlabel as appropriate depending on the timezone
    ax.set_xlabel('Hour (%s)' % time[-1].astimezone(tz).tzname())
    ax.set_xlim(*time_range)

    return axes

[jrleeman]

@mtb-za This is very relevant to our plotting discussion on SWUNG.

[dopplershift]

Wow, I think my sample code has workarounds for matplotlib pre-1.0 behavior.

[ahaberlie]

I'm working on a similar project and will pull request to metpy eventually if no one finishes this before me. Until then feel free to borrow code/inspiration/etc. from the project if anyone is working on this.

https://github.com/ahaberlie/ndfd_meteogram/

Basically trying to remake the NWS website meteograms at this point. Attached is what I've got so far.

[dopplershift]

That looks nice. It's certainly useful to see that, since it gives another use-case for what we want in MetPy.

[dopplershift]

Some ideas from #318 :

1. Use GridSpec
2. Have Meteogram implement __getitem__ as the generic way to access a panel, like meteo[0].plot(...). Just not sure if that should be spelled meteo['thermo'].plot(...)