diff --git a/requirements/dev.yml b/requirements/dev.yml index 5150458..ed8b07d 100644 --- a/requirements/dev.yml +++ b/requirements/dev.yml @@ -3,6 +3,7 @@ channels: - conda-forge dependencies: - matplotlib + - shapely - numpy - scipy - pip diff --git a/requirements/rtd.yml b/requirements/rtd.yml index 8ae66ae..f520371 100644 --- a/requirements/rtd.yml +++ b/requirements/rtd.yml @@ -3,6 +3,7 @@ channels: - conda-forge dependencies: - matplotlib + - shapely - numpy - scipy - sphinx diff --git a/tephi/__init__.py b/tephi/__init__.py index e65ed75..e80b5b9 100644 --- a/tephi/__init__.py +++ b/tephi/__init__.py @@ -1,271 +1,40 @@ -# Copyright Tephi contributors -# -# This file is part of Tephi and is released under the BSD license. -# See LICENSE in the root of the repository for full licensing details. -""" -The tephi module provides tephigram plotting of pressure, temperature and wind -barb data. - -.. warning:: - This is a beta release module and is liable to change. - -""" -from collections import namedtuple -from collections.abc import Iterable -from functools import partial from matplotlib.font_manager import FontProperties import matplotlib.pyplot as plt +from mpl_toolkits.axisartist import Subplot from mpl_toolkits.axisartist.grid_helper_curvelinear import ( GridHelperCurveLinear, ) -from mpl_toolkits.axisartist import Subplot -import numbers import numpy as np import os.path - -from . import isopleths -from . import transforms - +from . import artists, isopleths, transforms __version__ = "0.4.0.dev0" - -# -# Miscellaneous constants. -# -DEFAULT_WIDTH = 700 # in pixels - -ISOBAR_SPEC = [(25, 0.03), (50, 0.10), (100, 0.25), (200, 1.5)] -ISOBAR_LINE = {"color": "blue", "linewidth": 0.5, "clip_on": True} -ISOBAR_TEXT = { - "size": 8, - "color": "blue", - "clip_on": True, - "va": "bottom", - "ha": "right", -} -ISOBAR_FIXED = [50, 1000] - -WET_ADIABAT_SPEC = [(1, 0.05), (2, 0.15), (4, 1.5)] -WET_ADIABAT_LINE = {"color": "orange", "linewidth": 0.5, "clip_on": True} -WET_ADIABAT_TEXT = { - "size": 8, - "color": "orange", - "clip_on": True, - "va": "bottom", - "ha": "left", -} -WET_ADIABAT_FIXED = None - -MIXING_RATIO_SPEC = [(1, 0.05), (2, 0.18), (4, 0.3), (8, 1.5)] -MIXING_RATIO_LINE = {"color": "green", "linewidth": 0.5, "clip_on": True} -MIXING_RATIO_TEXT = { - "size": 8, - "color": "green", - "clip_on": True, - "va": "bottom", - "ha": "right", -} -MIXING_RATIOS = [ - 0.001, - 0.002, - 0.005, - 0.01, - 0.02, - 0.03, - 0.05, - 0.1, - 0.15, - 0.2, - 0.3, - 0.4, - 0.5, - 0.6, - 0.8, - 1.0, - 1.5, - 2.0, - 2.5, - 3.0, - 4.0, - 5.0, - 6.0, - 7.0, - 8.0, - 9.0, - 10.0, - 12.0, - 14.0, - 16.0, - 18.0, - 20.0, - 24.0, - 28.0, - 32.0, - 36.0, - 40.0, - 44.0, - 48.0, - 52.0, - 56.0, - 60.0, - 68.0, - 80.0, -] -MIXING_RATIO_FIXED = None - -MIN_PRESSURE = 50 # mb = hPa -MAX_PRESSURE = 1000 # mb = hPa -MIN_THETA = 0 # degC -MAX_THETA = 250 # degC -MIN_WET_ADIABAT = 1 # degC -MAX_WET_ADIABAT = 60 # degC -MIN_TEMPERATURE = -50 # degC - - RESOURCES_DIR = os.path.join(os.path.abspath(os.path.dirname(__file__)), "etc") DATA_DIR = os.path.join(RESOURCES_DIR, "test_data") -def loadtxt(*filenames, **kwargs): +class _FormatterTheta(object): """ - Load one or more text files of pressure, temperature, wind speed and wind - direction value sets. - - Each line should contain, at minimum, a single pressure value (mb or hPa), - and a single temperature value (degC), but may also contain a dewpoint - value (degC), wind speed (knots) and wind direction value (degrees from - north). - - Note that blank lines and comment lines beginning with a '#' are ignored. - - For example: - - >>> import os.path - >>> import tephi - - >>> winds = os.path.join(tephi.DATA_DIR, 'barbs.txt') - >>> columns = ('pressure', 'dewpoint', 'wind_speed', 'wind_direction') - >>> data = tephi.loadtxt(winds, column_titles=columns) - >>> pressure = data.pressure - >>> dews = data.dewpoint - >>> wind_speed = data.wind_speed - >>> wind_direction = data.wind_direction - - .. seealso:: :func:`numpy.loadtxt`. - - Args: - - * filenames: one or more filenames. - - Kwargs: - - * column_titles: - List of iterables, or None. If specified, should contain one title - string for each column of data per specified file. If all of multiple - files loaded have the same column titles, then only one tuple of column - titles need be specified. - - * delimiter: - The string used to separate values. This is passed directly to - :func:`np.loadtxt`, which defaults to using any whitespace as delimiter - if this keyword is not specified. - - * dtype: - The datatype to cast the data in the text file to. Passed directly to - :func:`np.loadtxt`. - - Returns: - A :func:`collections.namedtuple` instance containing one tuple, named - with the relevant column title if specified, for each column of data - in the text file loaded. If more than one file is loaded, a sequence - of namedtuples is returned. + Dry adiabats potential temperature axis tick formatter. """ - def _repr(nt): - """An improved representation of namedtuples over the default.""" - - typename = nt.__class__.__name__ - fields = nt._fields - n_fields = len(fields) - return_str = "{}(\n".format(typename) - for i, t in enumerate(fields): - gap = " " * 4 - if i == n_fields - 1: - ender = "" - else: - ender = "\n" - return_str += "{}{}={!r}{}".format(gap, t, getattr(nt, t), ender) - return_str += ")" - return return_str - - column_titles = kwargs.pop("column_titles", None) - delimiter = kwargs.pop("delimiter", None) - dtype = kwargs.pop("dtype", "f4") - - if column_titles is not None: - fields = column_titles[0] - if not isinstance(column_titles, str): - if isinstance(fields, Iterable) and not isinstance(fields, str): - # We've an iterable of iterables - multiple titles is True. - multiple_titles = True - if len(column_titles) > len(filenames): - msg = "Received {} files but {} sets of column titles." - raise ValueError( - msg.format(len(column_titles), len(filenames)) - ) - elif isinstance(fields, str): - # We've an iterable of title strings - use for namedtuple. - tephidata = namedtuple("tephidata", column_titles) - multiple_titles = False - else: - # Whatever we've got it isn't iterable, so raise TypeError. - msg = "Expected title to be string, got {!r}." - raise TypeError(msg.format(type(column_titles))) - else: - msg = "Expected column_titles to be iterable, got {!r}." - raise TypeError(msg.format(type(column_titles))) - - else: - tephidata = namedtuple("tephidata", ("pressure", "temperature")) - multiple_titles = False - - data = [] - for ct, arg in enumerate(filenames): - if isinstance(arg, str): - if os.path.isfile(arg): - if multiple_titles: - tephidata = namedtuple("tephidata", column_titles[ct]) - tephidata.__repr__ = _repr - payload = np.loadtxt(arg, dtype=dtype, delimiter=delimiter) - item = tephidata(*payload.T) - data.append(item) - else: - msg = "Item {} is either not a file or does not exist." - raise OSError(msg.format(arg)) - - if len(data) == 1: - data = data[0] - - return data - - -class _FormatterTheta: - """Dry adiabats potential temperature axis tick formatter.""" - def __call__(self, direction, factor, values): - return [r"$\theta={:.1f}$".format(value) for value in values] + return [r"$\theta={}$".format(value) for value in values] -class _FormatterIsotherm: - """Isotherms temperature axis tick formatter.""" +class _FormatterIsotherm(object): + """ + Isotherms temperature axis tick formatter. + + """ def __call__(self, direction, factor, values): - return [r" $T={:.1f}$".format(value) for value in values] + return [r"$T={}$".format(value) for value in values] -class Locator: +class Locator(object): """ Determine the fixed step axis tick locations when called with a tick range. @@ -286,546 +55,167 @@ def __init__(self, step): Args: - * step: the step value for each axis tick. + * step: + The step value for each axis tick. """ self.step = int(step) def __call__(self, start, stop): - """Calculate the axis ticks given the provided tick range.""" + """ + Calculate the axis ticks given the provided tick range. + """ step = self.step - start = (int(start) // step) * step - stop = (int(stop) // step) * step - ticks = np.arange(start, stop + step, step, dtype=int) + start = (int(start) / step) * step + stop = (int(stop) / step) * step + ticks = np.arange(start, stop + step, step) return ticks, len(ticks), 1 -def _refresh_isopleths(axes): - """ - Refresh the plot isobars, wet adiabats and mixing ratios and associated - text labels. - - Args: - - * axes: - Tephigram plotting :class:`matplotlib.axes.AxesSubplot` instance. - - Returns: - Boolean, whether the plot has changed. - - """ - changed = False - - # Determine the current zoom level. - xlim = axes.get_xlim() - delta_xlim = xlim[1] - xlim[0] - ylim = axes.get_ylim() - zoom = delta_xlim / axes.tephigram_original_delta_xlim - - # Determine the display mid-point. - x_point = xlim[0] + delta_xlim * 0.5 - y_point = ylim[0] + (ylim[1] - ylim[0]) * 0.5 - xy = np.array([[x_point, y_point]]) - xy_point = axes.tephigram_inverse.transform(xy)[0] - - for profile in axes.tephigram_profiles: - profile.refresh() - - for isopleth in axes.tephigram_isopleths: - changed = isopleth.refresh(zoom, xy_point) or changed +class TephiAxes(Subplot): + name = "tephigram" - return changed - - -def _handler(event): - """Matplotlib event handler.""" - - for axes in event.canvas.figure.axes: - if hasattr(axes, "tephigram"): - if _refresh_isopleths(axes): - event.canvas.figure.show() - - -class _PlotGroup(dict): - """ - Container for a related group of tephigram isopleths. - - Manages the creation and plotting of all isopleths within the group. - - """ - - def __init__( - self, - axes, - plot_func, - text_kwargs, - step, - zoom, - tags, - fixed=None, - xfocus=None, - ): - self.axes = axes - self.text_kwargs = text_kwargs - self.step = step - self.zoom = zoom - - pairs = [] - for tag in tags: - text = plt.text(0, 0, str(tag), **text_kwargs) - text.set_bbox( - dict( - boxstyle="Round,pad=0.3", - facecolor="white", - edgecolor="white", - alpha=0.5, - clip_on=True, - clip_box=self.axes.bbox, + def __init__(self, *args, **kwargs): + # Validate the subplot arguments. + if len(args) == 0: + args = (1, 1, 1) + elif len(args) == 1 and isinstance(args[0], int): + args = tuple([int(c) for c in str(args[0])]) + if len(args) != 3: + msg = ( + "Integer subplot specification must be a " + "three digit number. Not {}.".format(len(args)) ) - ) - pairs.append((tag, [plot_func(tag), text])) - - dict.__init__(self, pairs) - for line, text in self.values(): - line.set_visible(True) - text.set_visible(True) - self._visible = True - - if fixed is None: - fixed = [] - - if not isinstance(fixed, Iterable): - fixed = [fixed] - - if zoom is None: - self.fixed = set(tags) - else: - self.fixed = set(tags) & set(fixed) - - self.xfocus = xfocus - - def __setitem__(self, tag, item): - emsg = "Cannot add or set an item into the plot group {!r}" - raise ValueError(emsg.format(self.step)) - - def __getitem__(self, tag): - if tag not in self.keys(): - emsg = "Tag item {!r} is not a member of the plot group {!r}" - raise KeyError(emsg.format(tag, self.step)) - return dict.__getitem__(self, tag) - - def refresh(self, zoom, xy_point): - """ - Refresh all isopleths within the plot group. - - Args: - - * zoom: - Zoom level of the current plot, relative to the initial plot. - * xy_point: - The center point of the current point, transformed into - temperature and potential temperature. - - Returns: - Boolean, whether the plot group has changed. - - """ - if self.zoom is None or zoom <= self.zoom: - changed = self._item_on() - else: - changed = self._item_off() - self._refresh_text(xy_point) - return changed - - def _item_on(self, zoom=None): - changed = False - if zoom is None or self.zoom is None or zoom <= self.zoom: - if not self._visible: - for line, text in self.values(): - line.set_visible(True) - text.set_visible(True) - changed = True - self._visible = True - return changed - - def _item_off(self, zoom=None): - changed = False - if self.zoom is not None and (zoom is None or zoom > self.zoom): - if self._visible: - for tag, (line, text) in self.items(): - if tag not in self.fixed: - line.set_visible(False) - text.set_visible(False) - changed = True - self._visible = False - return changed - - def _generate_text(self, tag, xy_point): - line, text = self[tag] - x_data = line.get_xdata() - y_data = line.get_ydata() - - if self.xfocus: - delta = np.power(x_data - xy_point[0], 2) - else: - delta = np.power(x_data - xy_point[0], 2) + np.power( - y_data - xy_point[1], 2 - ) - index = np.argmin(delta) - text.set_position((x_data[index], y_data[index])) - - def _refresh_text(self, xy_point): - if self._visible: - for tag in self: - self._generate_text(tag, xy_point) - elif self.fixed: - for tag in self.fixed: - self._generate_text(tag, xy_point) - - -class _PlotCollection: - """ - Container for tephigram isopleths. - - Manages the creation and plotting of all tephigram isobars, mixing ratio - lines and pseudo saturated wet adiabats. - - """ - - def __init__( - self, - axes, - spec, - stop, - plot_func, - text_kwargs, - fixed=None, - minimum=None, - xfocus=None, - ): - if isinstance(stop, Iterable): - if minimum and minimum > max(stop): - emsg = "Minimum value of {!r} exceeds all other values" - raise ValueError(emsg.format(minimum)) - - items = [ - [step, zoom, set(stop[step - 1 :: step])] - for step, zoom in sorted(spec, reverse=True) - ] + raise ValueError(msg) else: - if minimum and minimum > stop: - emsg = "Minimum value of {!r} exceeds maximum threshold {!r}" - raise ValueError(emsg.format(minimum, stop)) - - items = [ - [step, zoom, set(range(step, stop + step, step))] - for step, zoom in sorted(spec, reverse=True) - ] - - for index, item in enumerate(items): - if minimum: - item[2] = set([value for value in item[2] if value >= minimum]) - - for subitem in items[index + 1 :]: - subitem[2] -= item[2] - - self.groups = { - item[0]: _PlotGroup( - axes, plot_func, text_kwargs, *item, fixed=fixed, xfocus=xfocus - ) - for item in items - if item[2] - } - - if not self.groups: - emsg = "The plot collection failed to generate any plot groups" - raise ValueError(emsg) - - def refresh(self, zoom, xy_point): - """ - Refresh all isopleth groups within the plot collection. - - Args: - - * zoom: - Zoom level of the current plot, relative to the initial plot. - * xy_point: - The center point of the current plot, transformed into - temperature and potential temperature. - - Returns: - Boolean, whether any plot group has changed. - - """ - changed = False - - for group in self.groups.values(): - changed = group.refresh(zoom, xy_point) or changed + msg = "Invalid arguments: " + ", ".join(["{}" for _ in len(args)]) + raise ValueError(msg.format(*args)) - return changed + # Process the kwargs. + figure = kwargs.get("figure") + isotherm_locator = kwargs.get("isotherm_locator") + dry_adiabat_locator = kwargs.get("dry_adiabat_locator") + anchor = None + if "anchor" in kwargs: + anchor = kwargs.pop("anchor") - -class Tephigram: - """ - Generate a tephigram of one or more pressure and temperature data sets. - - """ - - def __init__( - self, - figure=None, - isotherm_locator=None, - dry_adiabat_locator=None, - anchor=None, - ): - """ - Initialise the tephigram transformation and plot axes. - - Kwargs: - - * figure: - An existing :class:`matplotlib.figure.Figure` instance for the - tephigram plot. If a figure is not provided, a new figure will - be created by default. - * isotherm_locator: - A :class:`tephi.Locator` instance or a numeric step size - for the isotherm lines. - * dry_adiabat_locator: - A :class:`tephi.Locator` instance or a numeric step size - for the dry adiabat lines. - * anchor: - A sequence of two (pressure, temperature) pairs specifying the extent - of the tephigram plot in terms of the bottom right-hand corner, and - the top left-hand corner. Pressure data points must be in units of - mb or hPa, and temperature data points must be in units of degC. - - For example: - - .. plot:: - :include-source: - - import matplotlib.pyplot as plt - from numpy import column_stack - import os.path - import tephi - from tephi import Tephigram - - dew_point = os.path.join(tephi.DATA_DIR, 'dews.txt') - dry_bulb = os.path.join(tephi.DATA_DIR, 'temps.txt') - dew_data, temp_data = tephi.loadtxt(dew_point, dry_bulb) - dews = column_stack((dew_data.pressure, dew_data.temperature)) - temps = column_stack((temp_data.pressure, temp_data.temperature)) - tpg = Tephigram() - tpg.plot(dews, label='Dew-point', color='blue', linewidth=2) - tpg.plot(temps, label='Dry-bulb', color='red', linewidth=2) - plt.show() - - """ - if not figure: - # Create a default figure. - self.figure = plt.figure(0, figsize=(9, 9)) - else: - self.figure = figure + # Get the figure. + if figure is None: + figure = plt.gcf() # Configure the locators. - if isotherm_locator and not isinstance(isotherm_locator, Locator): - if not isinstance(isotherm_locator, numbers.Number): - raise ValueError("Invalid isotherm locator") - locator_isotherm = Locator(isotherm_locator) - else: - locator_isotherm = isotherm_locator - - if dry_adiabat_locator and not isinstance( - dry_adiabat_locator, Locator - ): - if not isinstance(dry_adiabat_locator, numbers.Number): - raise ValueError("Invalid dry adiabat locator") - locator_theta = Locator(dry_adiabat_locator) - else: - locator_theta = dry_adiabat_locator - - # Define the tephigram coordinate-system transformation. - self.tephi_transform = transforms.TephiTransform() - ghelper = GridHelperCurveLinear( - self.tephi_transform, + locator_isotherm = isotherm_locator + if locator_isotherm and not isinstance(locator_isotherm, Locator): + if not isinstance(locator_isotherm, int): + raise ValueError("Invalid isotherm locator.") + locator_isotherm = Locator(locator_isotherm) + locator_theta = dry_adiabat_locator + if locator_theta and not isinstance(locator_theta, Locator): + if not isinstance(locator_theta, int): + raise ValueError("Invalid dry adiabat locator.") + + from mpl_toolkits.axisartist.grid_finder import MaxNLocator + + locator_isotherm = MaxNLocator(nbins=20, steps=[10], integer=True) + locator_theta = MaxNLocator(nbins=20, steps=[10], integer=True) + + gridder = GridHelperCurveLinear( + transforms.TephiTransform(), tick_formatter1=_FormatterIsotherm(), grid_locator1=locator_isotherm, tick_formatter2=_FormatterTheta(), grid_locator2=locator_theta, ) - self.axes = Subplot(self.figure, 1, 1, 1, grid_helper=ghelper) - self.transform = self.tephi_transform + self.axes.transData - self.axes.axis["isotherm"] = self.axes.new_floating_axis(1, 0) - self.axes.axis["theta"] = self.axes.new_floating_axis(0, 0) - self.axes.axis["left"].get_helper().nth_coord_ticks = 0 - self.axes.axis["left"].toggle(all=True) - self.axes.axis["bottom"].get_helper().nth_coord_ticks = 1 - self.axes.axis["bottom"].toggle(all=True) - self.axes.axis["top"].get_helper().nth_coord_ticks = 0 - self.axes.axis["top"].toggle(all=False) - self.axes.axis["right"].get_helper().nth_coord_ticks = 1 - self.axes.axis["right"].toggle(all=True) - self.axes.gridlines.set_linestyle("solid") - - self.figure.add_subplot(self.axes) - - # Configure default axes. - axis = self.axes.axis["left"] + super(TephiAxes, self).__init__( + figure, *args, grid_helper=gridder, **kwargs + ) + + # The tephigram cache. + transform = transforms.TephiTransform() + self.transData + self.tephi = dict( + anchor=anchor, + figure=figure.add_subplot(self), + profiles=isopleths.ProfileList(), + transform=transform, + ) + + # Create each axis. + self.axis["isotherm"] = self.new_floating_axis(1, 0) + self.axis["theta"] = self.new_floating_axis(0, 0) + self.axis["left"].get_helper().nth_coord_ticks = 0 + self.axis["left"].toggle(all=True) + self.axis["bottom"].get_helper().nth_coord_ticks = 1 + self.axis["bottom"].toggle(all=True) + self.axis["top"].get_helper().nth_coord_ticks = 0 + self.axis["top"].toggle(all=False) # Turned-off + self.axis["right"].get_helper().nth_coord_ticks = 1 + self.axis["right"].toggle(all=True) + self.gridlines.set_linestyle("solid") + + # Configure each axis. + axis = self.axis["left"] axis.major_ticklabels.set_fontsize(10) axis.major_ticklabels.set_va("baseline") axis.major_ticklabels.set_rotation(135) - axis = self.axes.axis["right"] + axis = self.axis["right"] axis.major_ticklabels.set_fontsize(10) axis.major_ticklabels.set_va("baseline") axis.major_ticklabels.set_rotation(-135) - self.axes.axis["top"].major_ticklabels.set_fontsize(10) - axis = self.axes.axis["bottom"] + self.axis["top"].major_ticklabels.set_fontsize(10) + axis = self.axis["bottom"] axis.major_ticklabels.set_fontsize(10) axis.major_ticklabels.set_ha("left") - axis.major_ticklabels.set_va("top") + axis.major_ticklabels.set_va("bottom") axis.major_ticklabels.set_rotation(-45) # Isotherms: lines of constant temperature (degC). - axis = self.axes.axis["isotherm"] + axis = self.axis["isotherm"] axis.set_axis_direction("right") axis.set_axislabel_direction("-") axis.major_ticklabels.set_rotation(90) - axis.major_ticklabels.set_fontsize(10) + axis.major_ticklabels.set_fontsize(8) axis.major_ticklabels.set_va("bottom") axis.major_ticklabels.set_color("grey") - axis.major_ticklabels.set_visible(False) # turned-off + axis.major_ticklabels.set_visible(False) # Turned-off + axis.major_ticklabels.set_clip_box(self.bbox) # Dry adiabats: lines of constant potential temperature (degC). - axis = self.axes.axis["theta"] + axis = self.axis["theta"] axis.set_axis_direction("right") axis.set_axislabel_direction("+") - axis.major_ticklabels.set_fontsize(10) + axis.major_ticklabels.set_fontsize(8) axis.major_ticklabels.set_va("bottom") axis.major_ticklabels.set_color("grey") - axis.major_ticklabels.set_visible(False) # turned-off + axis.major_ticklabels.set_visible(False) # Turned-off + axis.major_ticklabels.set_clip_box(self.bbox) axis.line.set_linewidth(3) axis.line.set_linestyle("--") # Lock down the aspect ratio. - self.axes.set_aspect(1.0) - self.axes.grid(True) + self.set_aspect("equal") + self.grid(True) # Initialise the text formatter for the navigation status bar. - self.axes.format_coord = self._status_bar - - # Factor in the tephigram transform. - ISOBAR_TEXT["transform"] = self.transform - WET_ADIABAT_TEXT["transform"] = self.transform - MIXING_RATIO_TEXT["transform"] = self.transform - - # Create plot collections for the tephigram isopleths. - func = partial( - isopleths.isobar, - MIN_THETA, - MAX_THETA, - self.axes, - self.transform, - ISOBAR_LINE, - ) - self._isobars = _PlotCollection( - self.axes, - ISOBAR_SPEC, - MAX_PRESSURE, - func, - ISOBAR_TEXT, - fixed=ISOBAR_FIXED, - minimum=MIN_PRESSURE, - ) - - func = partial( - isopleths.wet_adiabat, - MAX_PRESSURE, - MIN_TEMPERATURE, - self.axes, - self.transform, - WET_ADIABAT_LINE, - ) - self._wet_adiabats = _PlotCollection( - self.axes, - WET_ADIABAT_SPEC, - MAX_WET_ADIABAT, - func, - WET_ADIABAT_TEXT, - fixed=WET_ADIABAT_FIXED, - minimum=MIN_WET_ADIABAT, - xfocus=True, - ) - - func = partial( - isopleths.mixing_ratio, - MIN_PRESSURE, - MAX_PRESSURE, - self.axes, - self.transform, - MIXING_RATIO_LINE, - ) - self._mixing_ratios = _PlotCollection( - self.axes, - MIXING_RATIO_SPEC, - MIXING_RATIOS, - func, - MIXING_RATIO_TEXT, - fixed=MIXING_RATIO_FIXED, - ) - - # Initialise for the tephigram plot event handler. - plt.connect("motion_notify_event", _handler) - self.axes.tephigram = True - self.axes.tephigram_original_delta_xlim = DEFAULT_WIDTH - self.original_delta_xlim = DEFAULT_WIDTH - self.axes.tephigram_transform = self.tephi_transform - self.axes.tephigram_inverse = self.tephi_transform.inverted() - self.axes.tephigram_isopleths = [ - self._isobars, - self._wet_adiabats, - self._mixing_ratios, - ] - - # The tephigram profiles. - self._profiles = [] - self.axes.tephigram_profiles = self._profiles + self.format_coord = self._status_bar # Center the plot around the anchor extent. - self._anchor = anchor - if self._anchor is not None: - self._anchor = np.asarray(anchor) - if ( - self._anchor.ndim != 2 - or self._anchor.shape[-1] != 2 - or len(self._anchor) != 2 - ): + if anchor is not None: + anchor = np.asarray(anchor) + if anchor.shape != (2, 2): msg = ( - "Invalid anchor, expecting [(bottom-right-pressure, " - "bottom-right-temperature), (top-left-pressure, " - "top-left-temperature)]" + "Invalid anchor, expecting [(BLHC-T, BLHC-t)," + "(TRHC-T, TRHC-t)]" ) raise ValueError(msg) - ( - (bottom_pressure, bottom_temp), - (top_pressure, top_temp), - ) = self._anchor - - if (bottom_pressure - top_pressure) < 0: - raise ValueError("Invalid anchor pressure range") - if (bottom_temp - top_temp) < 0: - raise ValueError("Invalid anchor temperature range") - - self._anchor = isopleths.Profile(anchor, self.axes) - self._anchor.plot(visible=False) - xlim, ylim = self._calculate_extents() - self.axes.set_xlim(xlim) - self.axes.set_ylim(ylim) + xlim, ylim = transforms.convert_Tt2xy(anchor[:, 0], anchor[:, 1]) + self.set_xlim(xlim) + self.set_ylim(ylim) + self.tephi["anchor"] = xlim, ylim def plot(self, data, **kwargs): """ - Plot the environmental lapse rate profile of the pressure and - temperature data points. + Plot the profile of the pressure and temperature data points. The pressure and temperature data points are transformed into potential temperature and temperature data points before plotting. @@ -839,51 +229,47 @@ def plot(self, data, **kwargs): Args: - * data: (pressure, temperature) pair data points. + * data: + Pressure and temperature pair data points. .. note:: All keyword arguments are passed through to :func:`matplotlib.pyplot.plot`. - For example: - .. plot:: :include-source: import matplotlib.pyplot as plt - from tephi import Tephigram + from tephi import TephiAxes - tpg = Tephigram() + ax = TephiAxes() data = [[1006, 26.4], [924, 20.3], [900, 19.8], [850, 14.5], [800, 12.9], [755, 8.3]] - profile = tpg.plot(data, color='red', linestyle='--', - linewidth=2, marker='o') + profile = ax.plot(data, color='red', linestyle='--', + linewidth=2, marker='o') barbs = [(10, 45, 900), (20, 60, 850), (25, 90, 800)] profile.barbs(barbs) plt.show() - For associating wind barbs with an environmental lapse rate profile, - see :meth:`~tephi.isopleths.Profile.barbs`. + For associating wind barbs with the profile, see + :meth:`~tephi.isopleths.Profile.barbs`. """ - profile = isopleths.Profile(data, self.axes) + profile = isopleths.Profile(self, data) profile.plot(**kwargs) - self._profiles.append(profile) + self.tephi["profiles"].append(profile) # Center the tephigram plot around all the profiles. - if self._anchor is None: + if self.tephi["anchor"] is None: xlim, ylim = self._calculate_extents(xfactor=0.25, yfactor=0.05) - self.axes.set_xlim(xlim) - self.axes.set_ylim(ylim) - - # Refresh the tephigram plot isopleths. - _refresh_isopleths(self.axes) + self.set_xlim(xlim) + self.set_ylim(ylim) # Show the plot legend. if "label" in kwargs: font_properties = FontProperties(size="x-small") plt.legend( - loc="upper left", + loc="upper right", fancybox=True, shadow=True, prop=font_properties, @@ -891,46 +277,109 @@ def plot(self, data, **kwargs): return profile + def add_isobars( + self, + ticks=None, + line=None, + text=None, + min_theta=None, + max_theta=None, + nbins=None, + ): + artist = artists.IsobarArtist( + ticks=ticks, + line=line, + text=text, + min_theta=min_theta, + max_theta=max_theta, + nbins=nbins, + ) + self.add_artist(artist) + + def add_wet_adiabats( + self, + ticks=None, + line=None, + text=None, + min_temperature=None, + max_pressure=None, + nbins=None, + ): + artist = artists.WetAdiabatArtist( + ticks=ticks, + line=line, + text=text, + min_temperature=min_temperature, + max_pressure=max_pressure, + nbins=nbins, + ) + self.add_artist(artist) + + def add_humidity_mixing_ratios( + self, + ticks=None, + line=None, + text=None, + min_pressure=None, + max_pressure=None, + nbins=None, + ): + artist = artists.HumidityMixingRatioArtist( + ticks=ticks, + line=line, + text=text, + min_pressure=min_pressure, + max_pressure=max_pressure, + nbins=nbins, + ) + self.add_artist(artist) + def _status_bar(self, x_point, y_point): - """Generate text for the interactive backend navigation status bar.""" + """ + Generate text for the interactive backend navigation status bar. + """ temperature, theta = transforms.convert_xy2Tt(x_point, y_point) pressure, _ = transforms.convert_Tt2pT(temperature, theta) - xlim = self.axes.get_xlim() - zoom = (xlim[1] - xlim[0]) / self.original_delta_xlim - msg = "T:{:.2f}, theta:{:.2f}, phi:{:.2f} (zoom:{:.3f})" - text = msg.format( - float(temperature), float(theta), float(pressure), zoom - ) - - return text + text = "T={:.2f}\u00b0C, \u03b8={:.2f}\u00b0C, p={:.2f}hPa" + return text.format(float(temperature), float(theta), float(pressure)) def _calculate_extents(self, xfactor=None, yfactor=None): - min_x = min_y = 1e10 - max_x = max_y = -1e-10 - profiles = self._profiles - transform = self.tephi_transform.transform - - if self._anchor is not None: - profiles = [self._anchor] - - for profile in profiles: - temperature = profile.temperature.reshape(-1, 1) - theta = profile.theta.reshape(-1, 1) - xy_points = transform(np.concatenate((temperature, theta), axis=1)) - x_points = xy_points[:, 0] - y_points = xy_points[:, 1] - min_x = np.min([min_x, np.min(x_points)]) - min_y = np.min([min_y, np.min(y_points)]) - max_x = np.max([max_x, np.max(x_points)]) - max_y = np.max([max_y, np.max(y_points)]) - - if xfactor is not None: - delta_x = max_x - min_x - min_x, max_x = min_x - xfactor * delta_x, max_x + xfactor * delta_x - - if yfactor is not None: - delta_y = max_y - min_y - min_y, max_y = min_y - yfactor * delta_y, max_y + yfactor * delta_y - - return ([min_x, max_x], [min_y, max_y]) + min_x = min_y = np.inf + max_x = max_y = -np.inf + + if self.tephi["anchor"] is not None: + xlim, ylim = self.tephi["anchor"] + else: + for profile in self.tephi["profiles"]: + temperature = profile.points.temperature + theta = profile.points.theta + x_points, y_points = transforms.convert_Tt2xy( + temperature, theta + ) + min_x, min_y = ( + np.min([min_x, np.min(x_points)]), + np.min([min_y, np.min(y_points)]), + ) + max_x, max_y = ( + np.max([max_x, np.max(x_points)]), + np.max([max_y, np.max(y_points)]), + ) + + if xfactor is not None: + delta_x = max_x - min_x + min_x, max_x = ( + (min_x - xfactor * delta_x), + (max_x + xfactor * delta_x), + ) + + if yfactor is not None: + delta_y = max_y - min_y + min_y, max_y = ( + (min_y - yfactor * delta_y), + (max_y + yfactor * delta_y), + ) + + xlim, ylim = (min_x, max_x), (min_y, max_y) + + return xlim, ylim diff --git a/tephi/artists.py b/tephi/artists.py new file mode 100644 index 0000000..2759cf2 --- /dev/null +++ b/tephi/artists.py @@ -0,0 +1,307 @@ +import matplotlib.artist +import numpy as np +from scipy.interpolate import interp1d +from shapely.geometry import LineString, Polygon +from shapely.prepared import prep + +from .constants import default +from .isopleths import Isobar, WetAdiabat, HumidityMixingRatio +from .transforms import convert_xy2Tt, convert_Tt2pT + + +class IsoplethArtist(matplotlib.artist.Artist): + def __init__(self): + super(IsoplethArtist, self).__init__() + self._isopleths = None + + def _locator(self, x0, x1, y0, y1): + temperature, theta = convert_xy2Tt([x0, x0, x1, x1], [y0, y1, y1, y0]) + bbox = prep(Polygon(zip(temperature, theta))) + mask = [bbox.intersects(item.geometry) for item in self._isopleths] + mask = np.asarray(mask) + + if self.nbins: + indices = np.where(mask)[0] + if indices.size: + if self.nbins < indices.size: + mask[:] = False + upint = indices.size + self.nbins - 1 + # this is an ugly solution, I'm sure there must be better ones + mask[indices[:: upint // self.nbins + 1]] = True + + return mask + + +class IsobarArtist(IsoplethArtist): + def __init__( + self, + ticks=None, + line=None, + text=None, + min_theta=None, + max_theta=None, + nbins=None, + ): + super(IsobarArtist, self).__init__() + if ticks is None: + ticks = default.get("isobar_ticks") + self.ticks = ticks + self._kwargs = {} + if line is None: + line = default.get("isobar_line") + self._kwargs["line"] = line + if text is None: + text = default.get("isobar_text") + self._kwargs["text"] = text + if min_theta is None: + min_theta = default.get("isobar_min_theta") + self.min_theta = min_theta + if max_theta is None: + max_theta = default.get("isobar_max_theta") + self.max_theta = max_theta + if nbins is None: + nbins = default.get("isobar_nbins") + elif nbins < 2 or isinstance(nbins, str): + nbins = None + self.nbins = nbins + + @matplotlib.artist.allow_rasterization + def draw( + self, renderer, line=None, text=None, min_theta=None, max_theta=None + ): + if not self.get_visible(): + return + axes = self.axes + draw_kwargs = dict(self._kwargs["line"]) + if line is not None: + draw_kwargs.update(line) + text_kwargs = dict(self._kwargs["text"]) + if text is not None: + text_kwargs.update(text) + if min_theta is None: + min_theta = self.min_theta + if max_theta is None: + max_theta = self.max_theta + + if self._isopleths is None: + isobars = [] + for tick in self.ticks: + isobars.append(Isobar(axes, tick, min_theta, max_theta)) + self._isopleths = np.asarray(isobars) + + (x0, x1), (y0, y1) = axes.get_xlim(), axes.get_ylim() + mask = self._locator(x0, x1, y0, y1) + + mx = x0 + axes.viewLim.width * 0.5 + temperature, theta = convert_xy2Tt([mx, mx], [y0, y1]) + text_line = LineString(zip(temperature, theta)) + + temperature, theta = convert_xy2Tt([mx] * 50, np.linspace(y0, y1, 50)) + pressure, _ = convert_Tt2pT(temperature, theta) + func = interp1d(pressure, theta, bounds_error=False) + + for isobar in self._isopleths[mask]: + isobar.draw(renderer, **draw_kwargs) + point = text_line.intersection(isobar.geometry) + if point: + isobar.refresh( + point.x, point.y, renderer=renderer, **text_kwargs + ) + else: + if func(isobar.data) < isobar.extent.theta.lower: + T = isobar.points.temperature[isobar.index.theta.lower] + t = isobar.extent.theta.lower + else: + T = isobar.points.temperature[isobar.index.theta.upper] + t = isobar.extent.theta.upper + isobar.refresh(T, t, renderer=renderer, **text_kwargs) + + +class WetAdiabatArtist(IsoplethArtist): + def __init__( + self, + ticks=None, + line=None, + text=None, + min_temperature=None, + max_pressure=None, + nbins=None, + ): + super(WetAdiabatArtist, self).__init__() + if ticks is None: + ticks = default.get("wet_adiabat_ticks") + self.ticks = sorted(ticks) + self._kwargs = {} + if line is None: + line = default.get("wet_adiabat_line") + self._kwargs["line"] = line + if text is None: + text = default.get("wet_adiabat_text") + self._kwargs["text"] = text + if min_temperature is None: + min_temperature = default.get("wet_adiabat_min_temperature") + self.min_temperature = min_temperature + if max_pressure is None: + max_pressure = default.get("wet_adiabat_max_pressure") + self.max_pressure = max_pressure + if nbins is None: + nbins = default.get("wet_adiabat_nbins") + if nbins < 2 or isinstance(nbins, str): + nbins = None + self.nbins = nbins + + @matplotlib.artist.allow_rasterization + def draw( + self, + renderer, + line=None, + text=None, + min_temperature=None, + max_pressure=None, + ): + if not self.get_visible(): + return + axes = self.axes + draw_kwargs = dict(self._kwargs["line"]) + if line is not None: + draw_kwargs.update(line) + text_kwargs = dict(self._kwargs["text"]) + if text is not None: + text_kwargs.update(text) + if min_temperature is None: + min_temperature = self.min_temperature + if max_pressure is None: + max_pressure = self.max_pressure + + if self._isopleths is None: + adiabats = [] + for tick in self.ticks: + adiabats.append( + WetAdiabat(axes, tick, min_temperature, max_pressure) + ) + self._isopleths = np.asarray(adiabats) + + (x0, x1), (y0, y1) = axes.get_xlim(), axes.get_ylim() + mask = self._locator(x0, x1, y0, y1) + + mx = x0 + axes.viewLim.width * 0.5 + my = y0 + axes.viewLim.height * 0.5 + temperature, theta = convert_xy2Tt([x0, mx, x1], [y0, my, y1]) + text_line = LineString(zip(temperature, theta)) + mT = temperature[1] + snap = None + + for adiabat in self._isopleths[mask]: + adiabat.draw(renderer, **draw_kwargs) + point = text_line.intersection(adiabat.geometry) + if point: + adiabat.refresh( + point.x, point.y, renderer=renderer, **text_kwargs + ) + else: + upper = abs(adiabat.extent.temperature.upper - mT) + lower = abs(adiabat.extent.temperature.lower - mT) + if snap == "upper" or upper < lower: + T = adiabat.extent.temperature.upper + t = adiabat.points.theta[adiabat.index.temperature.upper] + snap = "upper" + else: + T = adiabat.extent.temperature.lower + t = adiabat.points.theta[adiabat.index.temperature.lower] + snap = "lower" + adiabat.refresh(T, t, renderer=renderer, **text_kwargs) + + +class HumidityMixingRatioArtist(IsoplethArtist): + def __init__( + self, + ticks=None, + line=None, + text=None, + min_pressure=None, + max_pressure=None, + nbins=None, + ): + super(HumidityMixingRatioArtist, self).__init__() + if ticks is None: + ticks = default.get("mixing_ratio_ticks") + self.ticks = ticks + self._kwargs = {} + if line is None: + line = default.get("mixing_ratio_line") + self._kwargs["line"] = line + if text is None: + text = default.get("mixing_ratio_text") + self._kwargs["text"] = text + if min_pressure is None: + min_pressure = default.get("mixing_ratio_min_pressure") + self.min_pressure = min_pressure + if max_pressure is None: + max_pressure = default.get("mixing_ratio_max_pressure") + self.max_pressure = max_pressure + if nbins is None: + nbins = default.get("mixing_ratio_nbins") + if nbins < 2 or isinstance(nbins, str): + nbins = None + self.nbins = nbins + + @matplotlib.artist.allow_rasterization + def draw( + self, + renderer, + line=None, + text=None, + min_pressure=None, + max_pressure=None, + ): + if not self.get_visible(): + return + axes = self.axes + draw_kwargs = dict(self._kwargs["line"]) + if line is not None: + draw_kwargs.update(line) + text_kwargs = dict(self._kwargs["text"]) + if text is not None: + text_kwargs.update(text) + if min_pressure is None: + min_pressure = self.min_pressure + if max_pressure is None: + max_pressure = self.max_pressure + + if self._isopleths is None: + ratios = [] + for tick in self.ticks: + ratios.append( + HumidityMixingRatio(axes, tick, min_pressure, max_pressure) + ) + self._isopleths = np.asarray(ratios) + + (x0, x1), (y0, y1) = axes.get_xlim(), axes.get_ylim() + mask = self._locator(x0, x1, y0, y1) + + mx = x0 + axes.viewLim.width * 0.5 + my = y0 + axes.viewLim.height * 0.5 + temperature, theta = convert_xy2Tt([x0, mx, x1], [y1, my, y0]) + text_line = LineString(zip(temperature, theta)) + mt = theta[1] + snap = None + + for ratio in self._isopleths[mask]: + ratio.draw(renderer, **draw_kwargs) + point = text_line.intersection(ratio.geometry) + if point: + ratio.refresh( + point.x, point.y, renderer=renderer, **text_kwargs + ) + else: + upper = abs(ratio.extent.theta.upper - mt) + lower = abs(ratio.extent.theta.lower - mt) + if snap == "upper" or upper < lower: + T = ratio.points.temperature[ratio.index.theta.upper] + t = ratio.extent.theta.upper + snap = "upper" + else: + T = ratio.points.temperature[ratio.index.theta.lower] + t = ratio.extent.theta.lower + snap = "lower" + ratio.refresh(T, t, renderer=renderer, **text_kwargs) diff --git a/tephi/constants.py b/tephi/constants.py new file mode 100644 index 0000000..d162a1f --- /dev/null +++ b/tephi/constants.py @@ -0,0 +1,135 @@ +# Copyright Tephi contributors +# +# This file is part of Tephi and is released under the BSD license. +# See LICENSE in the root of the repository for full licensing details. +"""Tephigram transform and isopleth constants.""" + +# The specific heat capacity of dry air at a constant pressure, +# in units of J kg-1 K-1. +# TBC: This was originally set to 1.01e3 +Cp = 1004.0 + +# Dimensionless ratio: Rd / Cp. +K = 0.286 + +# Conversion offset between degree Celsius and Kelvin. +KELVIN = 273.15 + +# The specific latent heat of vapourisation of water at 0 degC, +# in units of J kg-1. +L = 2.501e6 + +MA = 300.0 + +# The specific gas constant for dry air, in units of J kg-1 K-1. +Rd = 287.0 + +# The specific gas constant for water vapour, in units of J kg-1 K-1. +Rv = 461.0 + +# Dimensionless ratio: Rd / Rv. +E = 0.622 + +# Base surface pressure. +P_BASE = 1000.0 + +# TODO: add in hodograph and mode defaults +default = { + "barbs_gutter": 0.1, + "barbs_length": 7, + "barbs_linewidth": 1.5, + "barbs_zorder": 10, + "isobar_line": dict(color="blue", linewidth=0.5, clip_on=True), + "isobar_min_theta": 0, + "isobar_max_theta": 250, + "isobar_nbins": None, + "isobar_text": dict( + size=8, color="blue", clip_on=True, va="bottom", ha="right" + ), + "isobar_ticks": [ + 1050, + 1000, + 950, + 900, + 850, + 800, + 700, + 600, + 500, + 400, + 300, + 250, + 200, + 150, + 100, + 70, + 50, + 40, + 30, + 20, + 10, + ], + "isopleth_picker": 3, + "isopleth_zorder": 10, + "mixing_ratio_line": dict(color="green", linewidth=0.5, clip_on=True), + "mixing_ratio_text": dict( + size=8, color="green", clip_on=True, va="bottom", ha="right" + ), + "mixing_ratio_min_pressure": 10, + "mixing_ratio_max_pressure": P_BASE, + "mixing_ratio_nbins": 10, + "mixing_ratio_ticks": [ + 0.001, + 0.002, + 0.005, + 0.01, + 0.02, + 0.03, + 0.05, + 0.1, + 0.15, + 0.2, + 0.3, + 0.4, + 0.5, + 0.6, + 0.8, + 1.0, + 1.5, + 2.0, + 2.5, + 3.0, + 4.0, + 5.0, + 6.0, + 7.0, + 8.0, + 9.0, + 10.0, + 12.0, + 14.0, + 16.0, + 18.0, + 20.0, + 24.0, + 28.0, + 32.0, + 36.0, + 40.0, + 44.0, + 48.0, + 52.0, + 56.0, + 60.0, + 68.0, + 80.0, + ], + "wet_adiabat_line": dict(color="orange", linewidth=0.5, clip_on=True), + "wet_adiabat_min_temperature": -50, + "wet_adiabat_max_pressure": P_BASE, + "wet_adiabat_nbins": 10, + "wet_adiabat_text": dict( + size=8, color="orange", clip_on=True, va="top", ha="left" + ), + "wet_adiabat_ticks": range(1, 61), +} diff --git a/tephi/isopleths.py b/tephi/isopleths.py index f6d698e..edeba04 100644 --- a/tephi/isopleths.py +++ b/tephi/isopleths.py @@ -7,21 +7,29 @@ environment profiles and barbs. """ + +from __future__ import absolute_import, division, print_function + +from abc import ABCMeta, abstractmethod +from collections import namedtuple import math +import matplotlib.artist from matplotlib.collections import PathCollection from matplotlib.path import Path import matplotlib.pyplot as plt -import matplotlib.transforms as mtransforms +import matplotlib.transforms as mtrans +from mpl_toolkits.axisartist import Subplot import numpy as np +from shapely.geometry import LineString from scipy.interpolate import interp1d -from ._constants import CONST_CP, CONST_L, CONST_KELVIN, CONST_RD, CONST_RV -from . import transforms +import tephi.constants as constants +from tephi.constants import default +import tephi.transforms as transforms # Wind barb speed (knots) ranges used since 1 January 1955. _BARB_BINS = np.arange(20) * 5 + 3 -_BARB_GUTTER = 0.1 _BARB_DTYPE = np.dtype( dict( names=("speed", "angle", "pressure", "barb"), @@ -29,216 +37,51 @@ ) ) -# -# Reference: http://www-nwp/~hadaa/tephigram/tephi_plot.html -# - - -def mixing_ratio( - min_pressure, max_pressure, axes, transform, kwargs, mixing_ratio_value -): - """ - Generate and plot a humidity mixing ratio line. - - A line of constant saturation mixing ratio with respect to a - plane water surface (g kg-1). - - Args: - - * min_pressure: - Minumum pressure, in mb or hPa, for the mixing ratio line extent. - - * max_pressure: - Maximum pressure, in mb or hPa, for the mixing ratio line extent. - - * axes: - Tephigram plotting :class:`matplotlib.axes.AxesSubplot` instance. - - * transform: - Tephigram plotting transformation - :class:`matplotlib.transforms.CompositeGenericTransform` instance. - - * kwargs: - Keyword arguments for the mixing ratio :class:`matplotlib.lines.Line2D` - instance. - - * mixing_ratio_value: - The mixing ratio value to be plotted. - - Returns: - The mixing ratio :class:`matplotlib.lines.Line2D` instance. - - """ - pressures = np.linspace(min_pressure, max_pressure, 100) - temps = transforms.convert_pw2T(pressures, mixing_ratio_value) - _, thetas = transforms.convert_pT2Tt(pressures, temps) - (line,) = axes.plot(temps, thetas, transform=transform, **kwargs) - - return line - - -def isobar(min_theta, max_theta, axes, transform, kwargs, pressure): - """ - Generate and plot an isobar line. - - A line of constant pressure (mb). - - Args: - - * min_theta: - Minimum potential temperature, in degC, for the isobar extent. - - * max_theta: - Maximum potential temperature, in degC, for the isobar extent. - - * axes: - Tephigram plotting :class:`matplotlib.axes.AxesSubplot` instance. - - * transform: - Tephigram plotting transformation - :class:`matplotlib.transforms.CompositeGenericTransform` instance. - - * kwargs: - Keyword arguments for the isobar :class:`matplotlib.lines.Line2D` - instance. - - * pressure: - The isobar pressure value, in mb or hPa, to be plotted. - - Returns: - The isobar :class:`matplotlib.lines.Line2D` instance. - - """ - steps = 100 - thetas = np.linspace(min_theta, max_theta, steps) - _, temps = transforms.convert_pt2pT([pressure] * steps, thetas) - (line,) = axes.plot(temps, thetas, transform=transform, **kwargs) - - return line - - -def _wet_adiabat_gradient(min_temperature, pressure, temperature, dp): - """ - Calculate the wet adiabat change in pressure and temperature. - - Args: - - * min_temperature: - Minimum potential temperature, in degC, for the wet adiabat line - extent. - - * pressure: - Pressure point value, in mb or hPa, from which to calculate the - gradient difference. - - * temperature: - Potential temperature point value, in degC, from which to calculate - the gradient difference. - - * dp: - The wet adiabat change in pressure, in mb or hPa, from which to - calculate the gradient difference. - - Returns: - The gradient change as a (pressure, potential-temperature) value pair. - - """ - - # TODO: Discover the meaning of the magic numbers. - - kelvin = temperature + CONST_KELVIN - lsbc = (CONST_L / CONST_RV) * ((1.0 / CONST_KELVIN) - (1.0 / kelvin)) - rw = 6.11 * np.exp(lsbc) * (0.622 / pressure) - lrwbt = (CONST_L * rw) / (CONST_RD * kelvin) - nume = ((CONST_RD * kelvin) / (CONST_CP * pressure)) * (1.0 + lrwbt) - deno = 1.0 + (lrwbt * ((0.622 * CONST_L) / (CONST_CP * kelvin))) - gradi = nume / deno - dt = dp * gradi - - if (temperature + dt) < min_temperature: - dt = min_temperature - temperature - dp = dt / gradi - - return dp, dt +# Isopleth defaults. +_DRY_ADIABAT_STEPS = 50 +_HUMIDITY_MIXING_RATIO_STEPS = 50 +_ISOBAR_STEPS = 50 +_ISOTHERM_STEPS = 50 +_SATURATION_ADIABAT_PRESSURE_DELTA = -5.0 +BOUNDS = namedtuple("BOUNDS", "lower upper") +POINTS = namedtuple("POINTS", "temperature theta pressure") -def wet_adiabat( - max_pressure, min_temperature, axes, transform, kwargs, temperature -): - """ - Generate and plot a pseudo saturated wet adiabat line. - A line of constant equivalent potential temperature for saturated - air parcels (degC). - - Args: - - * max_pressure: - Maximum pressure, in mb or hPa, for the wet adiabat line extent. - - * min_temperature: - Minimum potential temperature, in degC, for the wet adiabat line - extent. - - * axes: - Tephigram plotting :class:`matplotlib.axes.AxesSubplot` instance. - - * transform: - Tephigram plotting transformation - :class:`matplotlib.transforms.CompositeGenericTransform` instance. - - * kwargs: - Keyword arguments for the mixing ratio :class:`matplotlib.lines.Line2D` - instance. - - * temperature: - The wet adiabat value, in degC, to be plotted. - - Returns: - The wet adiabat :class:`matplotlib.lines.Line2D` instance. - - """ - temps = [temperature] - pressures = [max_pressure] - dp = -5.0 - - for i in range(200): - dp, dt = _wet_adiabat_gradient( - min_temperature, pressures[i], temps[i], dp +class BarbArtist(matplotlib.artist.Artist): + def __init__(self, barbs, **kwargs): + super(BarbArtist, self).__init__() + self._gutter = kwargs.pop("gutter", default.get("barbs_gutter")) + self._kwargs = dict( + length=default.get("barbs_length"), + zorder=default.get("barbs_zorder", 10), ) - temps.append(temps[i] + dt) - pressures.append(pressures[i] + dp) - - _, thetas = transforms.convert_pT2Tt(pressures, temps) - (line,) = axes.plot(temps, thetas, transform=transform, **kwargs) - - return line - - -class Barbs: - """Generate a wind arrow barb.""" - - def __init__(self, axes): - """ - Create a wind arrow barb for the given axes. - - Args: - - * axes: - A :class:`matplotlib.axes.AxesSubplot` instance. - - """ - self.axes = axes - self.barbs = None - self._gutter = None - self._transform = axes.tephigram_transform + axes.transData - self._kwargs = None - self._custom_kwargs = None - self._custom = dict( + self._kwargs.update(kwargs) + self.set_zorder(self._kwargs["zorder"]) + self._path_kwargs = dict( + color=None, + linewidth=default.get("barbs_linewidth"), + zorder=self._kwargs["zorder"], + ) + alias_by_kwarg = dict( color=["barbcolor", "color", "edgecolor", "facecolor"], linewidth=["lw", "linewidth"], linestyle=["ls", "linestyle"], ) + for kwarg, alias in iter(alias_by_kwarg.items()): + common = set(alias).intersection(kwargs) + if common: + self._path_kwargs[kwarg] = kwargs[sorted(common)[0]] + barbs = np.asarray(barbs) + if barbs.ndim != 2 or barbs.shape[-1] != 3: + msg = ( + "The barbs require to be a sequence of wind speed, " + "wind direction and pressure value triples." + ) + raise ValueError(msg) + self.barbs = np.empty(barbs.shape[0], dtype=_BARB_DTYPE) + for i, barb in enumerate(barbs): + self.barbs[i] = tuple(barb) + (None,) @staticmethod def _uv(magnitude, angle): @@ -281,6 +124,7 @@ def _uv(magnitude, angle): def _make_barb(self, temperature, theta, speed, angle): """Add the barb to the plot at the specified location.""" + transform = self.axes.tephi["transform"] u, v = self._uv(speed, angle) if 0 < speed < _BARB_BINS[0]: # Plot the missing barbless 1-2 knots line. @@ -289,8 +133,9 @@ def _make_barb(self, temperature, theta, speed, angle): pivot = self._kwargs.get("pivot", "tip") offset = pivot_points[pivot] verts = [(0.0, offset), (0.0, length + offset)] - rangle = math.radians(-angle) - verts = mtransforms.Affine2D().rotate(rangle).transform(verts) + verts = ( + mtrans.Affine2D().rotate(math.radians(-angle)).transform(verts) + ) codes = [Path.MOVETO, Path.LINETO] path = Path(verts, codes) size = length**2 / 4 @@ -299,165 +144,296 @@ def _make_barb(self, temperature, theta, speed, angle): [path], (size,), offsets=xy, - transOffset=self._transform, - **self._custom_kwargs, + transOffset=transform, + **self._path_kwargs, ) - barb.set_transform(mtransforms.IdentityTransform()) - self.axes.add_collection(barb) + barb.set_transform(mtrans.IdentityTransform()) else: - barb = plt.barbs( - temperature, - theta, - u, - v, - transform=self._transform, - **self._kwargs, + barb = self.axes.barbs( + temperature, theta, u, v, transform=transform, **self._kwargs ) + collections = list(self.axes.collections).remove(barb) + if collections: + self.axes.collections = tuple(collections) return barb - def refresh(self): - """Refresh the plot with the barbs.""" - if self.barbs is not None: - xlim = self.axes.get_xlim() - ylim = self.axes.get_ylim() - y = np.linspace(*ylim)[::-1] - xdelta = xlim[1] - xlim[0] - x = np.ones(y.size) * (xlim[1] - (xdelta * self._gutter)) - xy = np.column_stack((x, y)) - points = self.axes.tephigram_inverse.transform(xy) - temperature, theta = points[:, 0], points[:, 1] - pressure, _ = transforms.convert_Tt2pT(temperature, theta) - min_pressure, max_pressure = np.min(pressure), np.max(pressure) - func = interp1d(pressure, temperature) - for i, (speed, angle, pressure, barb) in enumerate(self.barbs): - if min_pressure < pressure < max_pressure: - p2T = func(pressure) - temperature, theta = transforms.convert_pT2Tt( - pressure, p2T - ) - if barb is None: - self.barbs[i]["barb"] = self._make_barb( - temperature, theta, speed, angle - ) - else: - barb.set_offsets(np.array([[temperature, theta]])) - barb.set_visible(True) + @matplotlib.artist.allow_rasterization + def draw(self, renderer): + if not self.get_visible(): + return + axes = self.axes + x0, x1 = axes.get_xlim() + y0, y1 = axes.get_ylim() + y = np.linspace(y0, y1)[::-1] + x = np.asarray([x1 - ((x1 - x0) * self._gutter)] * y.size) + temperature, theta = transforms.convert_xy2Tt(x, y) + pressure, _ = transforms.convert_Tt2pT(temperature, theta) + min_pressure, max_pressure = np.min(pressure), np.max(pressure) + func = interp1d(pressure, temperature) + for i, (speed, angle, pressure, barb) in enumerate(self.barbs): + if min_pressure < pressure < max_pressure: + temperature, theta = transforms.convert_pT2Tt( + pressure, func(pressure) + ) + if barb is None: + barb = self._make_barb(temperature, theta, speed, angle) + self.barbs[i]["barb"] = barb else: - if barb is not None: - barb.set_visible(False) + barb.set_offsets(np.array([[temperature, theta]])) + barb.draw(renderer) - def plot(self, barbs, **kwargs): - """ - Plot the sequence of barbs. - Args: +class Isopleth(object): + __metaclass__ = ABCMeta - * barbs: - Sequence of speed, direction and pressure value triples for - each barb. Where speed is measured in units of knots, direction - in units of degrees (clockwise from north), and pressure must - be in units of mb or hPa. + def __init__(self, axes): + self.axes = axes + self._transform = axes.tephi["transform"] + self.points = self._generate_points() + self.geometry = LineString( + np.vstack((self.points.temperature, self.points.theta)).T + ) + self.line = None + self.label = None + self._kwargs = dict(line={}, text={}) + Tmin, Tmax = ( + np.argmin(self.points.temperature), + np.argmax(self.points.temperature), + ) + tmin, tmax = ( + np.argmin(self.points.theta), + np.argmax(self.points.theta), + ) + pmin, pmax = ( + np.argmin(self.points.pressure), + np.argmax(self.points.pressure), + ) + self.index = POINTS( + BOUNDS(Tmin, Tmax), BOUNDS(tmin, tmax), BOUNDS(pmin, pmax) + ) + self.extent = POINTS( + BOUNDS( + self.points.temperature[Tmin], self.points.temperature[Tmax] + ), + BOUNDS(self.points.theta[tmin], self.points.theta[tmax]), + BOUNDS(self.points.pressure[pmin], self.points.pressure[pmax]), + ) - Kwargs: + @abstractmethod + def _generate_points(self): + pass + + def draw(self, renderer, **kwargs): + if self.line is None: + if "zorder" not in kwargs: + kwargs["zorder"] = default.get("isopleth_zorder") + draw_kwargs = dict(self._kwargs["line"]) + draw_kwargs.update(kwargs) + self.line = plt.Line2D( + self.points.temperature, + self.points.theta, + transform=self._transform, + **draw_kwargs, + ) + self.line.set_clip_box(self.axes.bbox) + self.line.draw(renderer) + return self.line + + def plot(self, **kwargs): + """ + Plot the points of the isopleth. - * gutter: - Proportion offset from the right hand side axis to plot the - barbs. Defaults to 0.1 + Kwargs: + See :func:`matplotlib.pyplot.plot`. - Also see :func:`matplotlib.pyplot.barbs` + Returns: + The isopleth :class:`matplotlib.lines.Line2D` """ - self._gutter = kwargs.pop("gutter", _BARB_GUTTER) - # zorder of 4.1 is higher than all MPL defaults, excluding legend. Also - # higher than tephi default for plot-lines. - self._kwargs = dict(length=7, zorder=4.1) - self._kwargs.update(kwargs) - self._custom_kwargs = dict( - color=None, linewidth=1.5, zorder=self._kwargs["zorder"] + if self.line is not None: + if self.line in self.axes.lines: + self.axes.lines.remove(self.line) + if "zorder" not in kwargs: + kwargs["zorder"] = default.get("isopleth_zorder") + if "picker" not in kwargs: + kwargs["picker"] = default.get("isopleth_picker") + plot_kwargs = dict(self._kwargs["line"]) + plot_kwargs.update(kwargs) + (self.line,) = Subplot.plot( + self.axes, + self.points.temperature, + self.points.theta, + transform=self._transform, + **plot_kwargs, ) - for key, values in self._custom.items(): - common = set(values).intersection(kwargs) - if common: - self._custom_kwargs[key] = kwargs[sorted(common)[0]] - if hasattr(barbs, "__next__"): - barbs = list(barbs) - barbs = np.asarray(barbs) - if barbs.ndim != 2 or barbs.shape[-1] != 3: - msg = ( - "The barbs require to be a sequence of wind speed, " - "wind direction and pressure value triples." - ) - raise ValueError(msg) - self.barbs = np.empty(barbs.shape[0], dtype=_BARB_DTYPE) - for i, barb in enumerate(barbs): - self.barbs[i] = tuple(barb) + (None,) - self.refresh() + return self.line + def text(self, temperature, theta, text, **kwargs): + if "zorder" not in kwargs: + kwargs["zorder"] = default.get("isopleth_zorder", 10) + 1 + text_kwargs = dict(self._kwargs["text"]) + text_kwargs.update(kwargs) + if self.label is not None and self.label in self.axes.texts: + self.axes.lines.remove(self.label) + self.label = self.axes.text( + temperature, + theta, + str(text), + transform=self._transform, + **text_kwargs, + ) + self.label.set_bbox( + dict( + boxstyle="Round,pad=0.3", + facecolor="white", + edgecolor="white", + alpha=0.5, + clip_on=True, + clip_box=self.axes.bbox, + ) + ) + return self.label + + def refresh(self, temperature, theta, renderer=None, **kwargs): + if self.label is None: + self.text(temperature, theta, self.data, **kwargs) + if renderer is not None: + try: + self.axes.tests = tuple( + list(self.axes.texts).remove(self.label) + ) + except TypeError: + self.axes.tests = None + else: + self.label.set_position((temperature, theta)) + if renderer is not None: + self.label.draw(renderer) -class Profile: - """Generate an environmental lapse rate profile.""" - def __init__(self, data, axes): - """ - Create an environmental lapse rate profile from the sequence of - pressure and temperature point data. +class DryAdiabat(Isopleth): + def __init__(self, axes, theta, min_pressure, max_pressure): + self.data = theta + self.bounds = BOUNDS(min_pressure, max_pressure) + self._steps = _DRY_ADIABAT_STEPS + super(DryAdiabat, self).__init__(axes) - Args: + def _generate_points(self): + pressure = np.linspace( + self.bounds.lower, self.bounds.upper, self._steps + ) + theta = np.asarray([self.data] * self._steps) + _, temperature = transforms.convert_pt2pT(pressure, theta) + return POINTS(temperature, theta, pressure) - * data: - Sequence of pressure and temperature points defining the - environmental lapse rate. - * axes: - The axes on which to plot the profile. +class HumidityMixingRatio(Isopleth): + def __init__(self, axes, mixing_ratio, min_pressure, max_pressure): + self.data = mixing_ratio + self.bounds = BOUNDS(min_pressure, max_pressure) + self._step = _HUMIDITY_MIXING_RATIO_STEPS + super(HumidityMixingRatio, self).__init__(axes) - """ - if hasattr(data, "__next__"): - data = list(data) - self.data = np.asarray(data) - if self.data.ndim != 2 or self.data.shape[-1] != 2: - msg = ( - "The environment profile data requires to be a sequence " - "of (pressure, temperature) value pairs." - ) - raise ValueError(msg) - self.axes = axes - self._transform = axes.tephigram_transform + axes.transData - self.pressure = self.data[:, 0] - self.temperature = self.data[:, 1] - _, self.theta = transforms.convert_pT2Tt( - self.pressure, self.temperature + def _generate_points(self): + pressure = np.linspace( + self.bounds.lower, self.bounds.upper, self._step ) - self.line = None - self._barbs = Barbs(axes) + temperature = transforms.convert_pw2T(pressure, self.data) + _, theta = transforms.convert_pT2Tt(pressure, temperature) + return POINTS(temperature, theta, pressure) + + +class Isobar(Isopleth): + def __init__(self, axes, pressure, min_theta, max_theta): + self.data = pressure + self.bounds = BOUNDS(min_theta, max_theta) + self._steps = _ISOBAR_STEPS + super(Isobar, self).__init__(axes) + self._kwargs["line"] = default.get("isobar_line") + self._kwargs["text"] = default.get("isobar_text") + + def _generate_points(self): + pressure = np.asarray([self.data] * self._steps) + theta = np.linspace(self.bounds.lower, self.bounds.upper, self._steps) + _, temperature = transforms.convert_pt2pT(pressure, theta) + return POINTS(temperature, theta, pressure) + + +class Isotherm(Isopleth): + def __init__(self, axes, temperature, min_pressure, max_pressure): + self.data = temperature + self.bounds = BOUNDS(min_pressure, max_pressure) + self._steps = _ISOTHERM_STEPS + super(Isotherm, self).__init__(axes) + + def _generate_points(self): + pressure = np.linspace( + self.bounds.lower, self.bounds.upper, self._steps + ) + temperature = np.asarray([self.data] * self._steps) + _, theta = transforms.convert_pT2Tt(pressure, temperature) + return POINTS(temperature, theta, pressure) - def plot(self, **kwargs): + +class Profile(Isopleth): + def __init__(self, axes, data): """ - Plot the environmental lapse rate profile. + Create a profile from the sequence of pressure and temperature points. - Kwargs: + Args: - See :func:`matplotlib.pyplot.plot`. + * axes: + The tephigram axes on which to plot the profile. - Returns: - The profile :class:`matplotlib.lines.Line2D` + * data: + Sequence of pressure and temperature points defining + the profile. """ - if self.line is not None and self.line in self.axes.lines: - self.axes.lines.remove(self.line) - - # zorder of 4 is higher than all MPL defaults, excluding legend. - if "zorder" not in kwargs: - kwargs["zorder"] = 4 + self.data = np.asarray(data) + super(Profile, self).__init__(axes) + self._barbs = None + self._highlight = None + + def has_highlight(self): + return self._highlight is not None + + def highlight(self, state=None): + if state is None: + state = not self.has_highlight() + if state: + if self._highlight is None: + linewidth = self.line.get_linewidth() * 7 + zorder = default.get("isopleth_zorder", 10) - 1 + kwargs = dict( + linewidth=linewidth, + color="grey", + alpha=0.3, + transform=self._transform, + zorder=zorder, + ) + (self._highlight,) = Subplot.plot( + self.axes, + self.points.temperature, + self.points.theta, + **kwargs, + ) + else: + if self._highlight is not None: + self.axes.lines.remove(self._highlight) + self._highlight = None - (self.line,) = self.axes.plot( - self.temperature, self.theta, transform=self._transform, **kwargs - ) - return self.line + def _generate_points(self): + if self.data.ndim != 2 or self.data.shape[-1] != 2: + msg = ( + "The profile data requires to be a sequence " + "of pressure, temperature value pairs." + ) + raise ValueError(msg) - def refresh(self): - """Refresh the plot with the profile and any associated barbs.""" - self._barbs.refresh() + pressure = self.data[:, 0] + temperature = self.data[:, 1] + _, theta = transforms.convert_pT2Tt(pressure, temperature) + return POINTS(temperature, theta, pressure) def barbs(self, barbs, **kwargs): """ @@ -473,10 +449,85 @@ def barbs(self, barbs, **kwargs): Kwargs: + * kwargs: See :func:`matplotlib.pyplot.barbs` """ colors = ["color", "barbcolor", "edgecolor", "facecolor"] if not set(colors).intersection(kwargs): kwargs["color"] = self.line.get_color() - self._barbs.plot(barbs, **kwargs) + self._barbs = BarbArtist(barbs, **kwargs) + self.axes.add_artist(self._barbs) + + def get_barbs(self): + return self._barbs.barbs + + +class WetAdiabat(Isopleth): + def __init__(self, axes, theta_e, min_temperature, max_pressure): + self.data = theta_e + self.bounds = BOUNDS(min_temperature, max_pressure) + self._delta_pressure = _SATURATION_ADIABAT_PRESSURE_DELTA + super(WetAdiabat, self).__init__(axes) + + def _gradient(self, pressure, temperature, dp): + stop = False + + kelvin = temperature + constants.KELVIN + lsbc = (constants.L / constants.Rv) * ( + (1.0 / constants.KELVIN) - (1.0 / kelvin) + ) + rw = 6.11 * np.exp(lsbc) * (constants.E / pressure) + lrwbt = (constants.L * rw) / (constants.Rd * kelvin) + numerator = ((constants.Rd * kelvin) / (constants.Cp * pressure)) * ( + 1.0 + lrwbt + ) + denominator = 1.0 + ( + lrwbt * ((constants.E * constants.L) / (constants.Cp * kelvin)) + ) + grad = numerator / denominator + dt = dp * grad + + if (temperature + dt) < self.bounds.lower: + dt = self.bounds.lower - temperature + dp = dt / grad + stop = True + + return dp, dt, stop + + def _generate_points(self): + temperature = [self.data] + pressure = [self.bounds.upper] + stop = False + dp = self._delta_pressure + + while not stop: + dp, dT, stop = self._gradient(pressure[-1], temperature[-1], dp) + pressure.append(pressure[-1] + dp) + temperature.append(temperature[-1] + dT) + + _, theta = transforms.convert_pT2Tt(pressure, temperature) + return POINTS(temperature, theta, pressure) + + +class ProfileList(list): + def __new__(cls, profiles=None): + profile_list = list.__new__(cls, profiles) + if not all(isinstance(profile, Profile) for profile in profile_list): + msg = "All items in the list must be a Profile instance." + raise TypeError(msg) + return profile_list + + def highlighted(self): + profiles = [profile for profile in self if profile.has_highlight()] + return profiles + + def picker(self, artist): + result = None + for profile in self: + if profile.line == artist: + result = profile + break + if result is None: + raise ValueError("Picker cannot find the profile.") + return result diff --git a/tephi/old__init__.py b/tephi/old__init__.py new file mode 100644 index 0000000..e65ed75 --- /dev/null +++ b/tephi/old__init__.py @@ -0,0 +1,936 @@ +# Copyright Tephi contributors +# +# This file is part of Tephi and is released under the BSD license. +# See LICENSE in the root of the repository for full licensing details. +""" +The tephi module provides tephigram plotting of pressure, temperature and wind +barb data. + +.. warning:: + This is a beta release module and is liable to change. + +""" +from collections import namedtuple +from collections.abc import Iterable +from functools import partial +from matplotlib.font_manager import FontProperties +import matplotlib.pyplot as plt +from mpl_toolkits.axisartist.grid_helper_curvelinear import ( + GridHelperCurveLinear, +) +from mpl_toolkits.axisartist import Subplot +import numbers +import numpy as np +import os.path + +from . import isopleths +from . import transforms + + +__version__ = "0.4.0.dev0" + + +# +# Miscellaneous constants. +# +DEFAULT_WIDTH = 700 # in pixels + +ISOBAR_SPEC = [(25, 0.03), (50, 0.10), (100, 0.25), (200, 1.5)] +ISOBAR_LINE = {"color": "blue", "linewidth": 0.5, "clip_on": True} +ISOBAR_TEXT = { + "size": 8, + "color": "blue", + "clip_on": True, + "va": "bottom", + "ha": "right", +} +ISOBAR_FIXED = [50, 1000] + +WET_ADIABAT_SPEC = [(1, 0.05), (2, 0.15), (4, 1.5)] +WET_ADIABAT_LINE = {"color": "orange", "linewidth": 0.5, "clip_on": True} +WET_ADIABAT_TEXT = { + "size": 8, + "color": "orange", + "clip_on": True, + "va": "bottom", + "ha": "left", +} +WET_ADIABAT_FIXED = None + +MIXING_RATIO_SPEC = [(1, 0.05), (2, 0.18), (4, 0.3), (8, 1.5)] +MIXING_RATIO_LINE = {"color": "green", "linewidth": 0.5, "clip_on": True} +MIXING_RATIO_TEXT = { + "size": 8, + "color": "green", + "clip_on": True, + "va": "bottom", + "ha": "right", +} +MIXING_RATIOS = [ + 0.001, + 0.002, + 0.005, + 0.01, + 0.02, + 0.03, + 0.05, + 0.1, + 0.15, + 0.2, + 0.3, + 0.4, + 0.5, + 0.6, + 0.8, + 1.0, + 1.5, + 2.0, + 2.5, + 3.0, + 4.0, + 5.0, + 6.0, + 7.0, + 8.0, + 9.0, + 10.0, + 12.0, + 14.0, + 16.0, + 18.0, + 20.0, + 24.0, + 28.0, + 32.0, + 36.0, + 40.0, + 44.0, + 48.0, + 52.0, + 56.0, + 60.0, + 68.0, + 80.0, +] +MIXING_RATIO_FIXED = None + +MIN_PRESSURE = 50 # mb = hPa +MAX_PRESSURE = 1000 # mb = hPa +MIN_THETA = 0 # degC +MAX_THETA = 250 # degC +MIN_WET_ADIABAT = 1 # degC +MAX_WET_ADIABAT = 60 # degC +MIN_TEMPERATURE = -50 # degC + + +RESOURCES_DIR = os.path.join(os.path.abspath(os.path.dirname(__file__)), "etc") +DATA_DIR = os.path.join(RESOURCES_DIR, "test_data") + + +def loadtxt(*filenames, **kwargs): + """ + Load one or more text files of pressure, temperature, wind speed and wind + direction value sets. + + Each line should contain, at minimum, a single pressure value (mb or hPa), + and a single temperature value (degC), but may also contain a dewpoint + value (degC), wind speed (knots) and wind direction value (degrees from + north). + + Note that blank lines and comment lines beginning with a '#' are ignored. + + For example: + + >>> import os.path + >>> import tephi + + >>> winds = os.path.join(tephi.DATA_DIR, 'barbs.txt') + >>> columns = ('pressure', 'dewpoint', 'wind_speed', 'wind_direction') + >>> data = tephi.loadtxt(winds, column_titles=columns) + >>> pressure = data.pressure + >>> dews = data.dewpoint + >>> wind_speed = data.wind_speed + >>> wind_direction = data.wind_direction + + .. seealso:: :func:`numpy.loadtxt`. + + Args: + + * filenames: one or more filenames. + + Kwargs: + + * column_titles: + List of iterables, or None. If specified, should contain one title + string for each column of data per specified file. If all of multiple + files loaded have the same column titles, then only one tuple of column + titles need be specified. + + * delimiter: + The string used to separate values. This is passed directly to + :func:`np.loadtxt`, which defaults to using any whitespace as delimiter + if this keyword is not specified. + + * dtype: + The datatype to cast the data in the text file to. Passed directly to + :func:`np.loadtxt`. + + Returns: + A :func:`collections.namedtuple` instance containing one tuple, named + with the relevant column title if specified, for each column of data + in the text file loaded. If more than one file is loaded, a sequence + of namedtuples is returned. + + """ + + def _repr(nt): + """An improved representation of namedtuples over the default.""" + + typename = nt.__class__.__name__ + fields = nt._fields + n_fields = len(fields) + return_str = "{}(\n".format(typename) + for i, t in enumerate(fields): + gap = " " * 4 + if i == n_fields - 1: + ender = "" + else: + ender = "\n" + return_str += "{}{}={!r}{}".format(gap, t, getattr(nt, t), ender) + return_str += ")" + return return_str + + column_titles = kwargs.pop("column_titles", None) + delimiter = kwargs.pop("delimiter", None) + dtype = kwargs.pop("dtype", "f4") + + if column_titles is not None: + fields = column_titles[0] + if not isinstance(column_titles, str): + if isinstance(fields, Iterable) and not isinstance(fields, str): + # We've an iterable of iterables - multiple titles is True. + multiple_titles = True + if len(column_titles) > len(filenames): + msg = "Received {} files but {} sets of column titles." + raise ValueError( + msg.format(len(column_titles), len(filenames)) + ) + elif isinstance(fields, str): + # We've an iterable of title strings - use for namedtuple. + tephidata = namedtuple("tephidata", column_titles) + multiple_titles = False + else: + # Whatever we've got it isn't iterable, so raise TypeError. + msg = "Expected title to be string, got {!r}." + raise TypeError(msg.format(type(column_titles))) + else: + msg = "Expected column_titles to be iterable, got {!r}." + raise TypeError(msg.format(type(column_titles))) + + else: + tephidata = namedtuple("tephidata", ("pressure", "temperature")) + multiple_titles = False + + data = [] + for ct, arg in enumerate(filenames): + if isinstance(arg, str): + if os.path.isfile(arg): + if multiple_titles: + tephidata = namedtuple("tephidata", column_titles[ct]) + tephidata.__repr__ = _repr + payload = np.loadtxt(arg, dtype=dtype, delimiter=delimiter) + item = tephidata(*payload.T) + data.append(item) + else: + msg = "Item {} is either not a file or does not exist." + raise OSError(msg.format(arg)) + + if len(data) == 1: + data = data[0] + + return data + + +class _FormatterTheta: + """Dry adiabats potential temperature axis tick formatter.""" + + def __call__(self, direction, factor, values): + return [r"$\theta={:.1f}$".format(value) for value in values] + + +class _FormatterIsotherm: + """Isotherms temperature axis tick formatter.""" + + def __call__(self, direction, factor, values): + return [r" $T={:.1f}$".format(value) for value in values] + + +class Locator: + """ + Determine the fixed step axis tick locations when called with a tick range. + + """ + + def __init__(self, step): + """ + Set the fixed step value for the axis tick locations. + + Generate tick location specification when called with a tick range. + + For example: + + >>> from tephi import Locator + >>> locator = Locator(10) + >>> locator(-45, 23) + (array([-50, -40, -30, -20, -10, 0, 10, 20]), 8, 1) + + Args: + + * step: the step value for each axis tick. + + """ + self.step = int(step) + + def __call__(self, start, stop): + """Calculate the axis ticks given the provided tick range.""" + + step = self.step + start = (int(start) // step) * step + stop = (int(stop) // step) * step + ticks = np.arange(start, stop + step, step, dtype=int) + return ticks, len(ticks), 1 + + +def _refresh_isopleths(axes): + """ + Refresh the plot isobars, wet adiabats and mixing ratios and associated + text labels. + + Args: + + * axes: + Tephigram plotting :class:`matplotlib.axes.AxesSubplot` instance. + + Returns: + Boolean, whether the plot has changed. + + """ + changed = False + + # Determine the current zoom level. + xlim = axes.get_xlim() + delta_xlim = xlim[1] - xlim[0] + ylim = axes.get_ylim() + zoom = delta_xlim / axes.tephigram_original_delta_xlim + + # Determine the display mid-point. + x_point = xlim[0] + delta_xlim * 0.5 + y_point = ylim[0] + (ylim[1] - ylim[0]) * 0.5 + xy = np.array([[x_point, y_point]]) + xy_point = axes.tephigram_inverse.transform(xy)[0] + + for profile in axes.tephigram_profiles: + profile.refresh() + + for isopleth in axes.tephigram_isopleths: + changed = isopleth.refresh(zoom, xy_point) or changed + + return changed + + +def _handler(event): + """Matplotlib event handler.""" + + for axes in event.canvas.figure.axes: + if hasattr(axes, "tephigram"): + if _refresh_isopleths(axes): + event.canvas.figure.show() + + +class _PlotGroup(dict): + """ + Container for a related group of tephigram isopleths. + + Manages the creation and plotting of all isopleths within the group. + + """ + + def __init__( + self, + axes, + plot_func, + text_kwargs, + step, + zoom, + tags, + fixed=None, + xfocus=None, + ): + self.axes = axes + self.text_kwargs = text_kwargs + self.step = step + self.zoom = zoom + + pairs = [] + for tag in tags: + text = plt.text(0, 0, str(tag), **text_kwargs) + text.set_bbox( + dict( + boxstyle="Round,pad=0.3", + facecolor="white", + edgecolor="white", + alpha=0.5, + clip_on=True, + clip_box=self.axes.bbox, + ) + ) + pairs.append((tag, [plot_func(tag), text])) + + dict.__init__(self, pairs) + for line, text in self.values(): + line.set_visible(True) + text.set_visible(True) + self._visible = True + + if fixed is None: + fixed = [] + + if not isinstance(fixed, Iterable): + fixed = [fixed] + + if zoom is None: + self.fixed = set(tags) + else: + self.fixed = set(tags) & set(fixed) + + self.xfocus = xfocus + + def __setitem__(self, tag, item): + emsg = "Cannot add or set an item into the plot group {!r}" + raise ValueError(emsg.format(self.step)) + + def __getitem__(self, tag): + if tag not in self.keys(): + emsg = "Tag item {!r} is not a member of the plot group {!r}" + raise KeyError(emsg.format(tag, self.step)) + return dict.__getitem__(self, tag) + + def refresh(self, zoom, xy_point): + """ + Refresh all isopleths within the plot group. + + Args: + + * zoom: + Zoom level of the current plot, relative to the initial plot. + * xy_point: + The center point of the current point, transformed into + temperature and potential temperature. + + Returns: + Boolean, whether the plot group has changed. + + """ + if self.zoom is None or zoom <= self.zoom: + changed = self._item_on() + else: + changed = self._item_off() + self._refresh_text(xy_point) + return changed + + def _item_on(self, zoom=None): + changed = False + if zoom is None or self.zoom is None or zoom <= self.zoom: + if not self._visible: + for line, text in self.values(): + line.set_visible(True) + text.set_visible(True) + changed = True + self._visible = True + return changed + + def _item_off(self, zoom=None): + changed = False + if self.zoom is not None and (zoom is None or zoom > self.zoom): + if self._visible: + for tag, (line, text) in self.items(): + if tag not in self.fixed: + line.set_visible(False) + text.set_visible(False) + changed = True + self._visible = False + return changed + + def _generate_text(self, tag, xy_point): + line, text = self[tag] + x_data = line.get_xdata() + y_data = line.get_ydata() + + if self.xfocus: + delta = np.power(x_data - xy_point[0], 2) + else: + delta = np.power(x_data - xy_point[0], 2) + np.power( + y_data - xy_point[1], 2 + ) + index = np.argmin(delta) + text.set_position((x_data[index], y_data[index])) + + def _refresh_text(self, xy_point): + if self._visible: + for tag in self: + self._generate_text(tag, xy_point) + elif self.fixed: + for tag in self.fixed: + self._generate_text(tag, xy_point) + + +class _PlotCollection: + """ + Container for tephigram isopleths. + + Manages the creation and plotting of all tephigram isobars, mixing ratio + lines and pseudo saturated wet adiabats. + + """ + + def __init__( + self, + axes, + spec, + stop, + plot_func, + text_kwargs, + fixed=None, + minimum=None, + xfocus=None, + ): + if isinstance(stop, Iterable): + if minimum and minimum > max(stop): + emsg = "Minimum value of {!r} exceeds all other values" + raise ValueError(emsg.format(minimum)) + + items = [ + [step, zoom, set(stop[step - 1 :: step])] + for step, zoom in sorted(spec, reverse=True) + ] + else: + if minimum and minimum > stop: + emsg = "Minimum value of {!r} exceeds maximum threshold {!r}" + raise ValueError(emsg.format(minimum, stop)) + + items = [ + [step, zoom, set(range(step, stop + step, step))] + for step, zoom in sorted(spec, reverse=True) + ] + + for index, item in enumerate(items): + if minimum: + item[2] = set([value for value in item[2] if value >= minimum]) + + for subitem in items[index + 1 :]: + subitem[2] -= item[2] + + self.groups = { + item[0]: _PlotGroup( + axes, plot_func, text_kwargs, *item, fixed=fixed, xfocus=xfocus + ) + for item in items + if item[2] + } + + if not self.groups: + emsg = "The plot collection failed to generate any plot groups" + raise ValueError(emsg) + + def refresh(self, zoom, xy_point): + """ + Refresh all isopleth groups within the plot collection. + + Args: + + * zoom: + Zoom level of the current plot, relative to the initial plot. + * xy_point: + The center point of the current plot, transformed into + temperature and potential temperature. + + Returns: + Boolean, whether any plot group has changed. + + """ + changed = False + + for group in self.groups.values(): + changed = group.refresh(zoom, xy_point) or changed + + return changed + + +class Tephigram: + """ + Generate a tephigram of one or more pressure and temperature data sets. + + """ + + def __init__( + self, + figure=None, + isotherm_locator=None, + dry_adiabat_locator=None, + anchor=None, + ): + """ + Initialise the tephigram transformation and plot axes. + + Kwargs: + + * figure: + An existing :class:`matplotlib.figure.Figure` instance for the + tephigram plot. If a figure is not provided, a new figure will + be created by default. + * isotherm_locator: + A :class:`tephi.Locator` instance or a numeric step size + for the isotherm lines. + * dry_adiabat_locator: + A :class:`tephi.Locator` instance or a numeric step size + for the dry adiabat lines. + * anchor: + A sequence of two (pressure, temperature) pairs specifying the extent + of the tephigram plot in terms of the bottom right-hand corner, and + the top left-hand corner. Pressure data points must be in units of + mb or hPa, and temperature data points must be in units of degC. + + For example: + + .. plot:: + :include-source: + + import matplotlib.pyplot as plt + from numpy import column_stack + import os.path + import tephi + from tephi import Tephigram + + dew_point = os.path.join(tephi.DATA_DIR, 'dews.txt') + dry_bulb = os.path.join(tephi.DATA_DIR, 'temps.txt') + dew_data, temp_data = tephi.loadtxt(dew_point, dry_bulb) + dews = column_stack((dew_data.pressure, dew_data.temperature)) + temps = column_stack((temp_data.pressure, temp_data.temperature)) + tpg = Tephigram() + tpg.plot(dews, label='Dew-point', color='blue', linewidth=2) + tpg.plot(temps, label='Dry-bulb', color='red', linewidth=2) + plt.show() + + """ + if not figure: + # Create a default figure. + self.figure = plt.figure(0, figsize=(9, 9)) + else: + self.figure = figure + + # Configure the locators. + if isotherm_locator and not isinstance(isotherm_locator, Locator): + if not isinstance(isotherm_locator, numbers.Number): + raise ValueError("Invalid isotherm locator") + locator_isotherm = Locator(isotherm_locator) + else: + locator_isotherm = isotherm_locator + + if dry_adiabat_locator and not isinstance( + dry_adiabat_locator, Locator + ): + if not isinstance(dry_adiabat_locator, numbers.Number): + raise ValueError("Invalid dry adiabat locator") + locator_theta = Locator(dry_adiabat_locator) + else: + locator_theta = dry_adiabat_locator + + # Define the tephigram coordinate-system transformation. + self.tephi_transform = transforms.TephiTransform() + ghelper = GridHelperCurveLinear( + self.tephi_transform, + tick_formatter1=_FormatterIsotherm(), + grid_locator1=locator_isotherm, + tick_formatter2=_FormatterTheta(), + grid_locator2=locator_theta, + ) + self.axes = Subplot(self.figure, 1, 1, 1, grid_helper=ghelper) + self.transform = self.tephi_transform + self.axes.transData + self.axes.axis["isotherm"] = self.axes.new_floating_axis(1, 0) + self.axes.axis["theta"] = self.axes.new_floating_axis(0, 0) + self.axes.axis["left"].get_helper().nth_coord_ticks = 0 + self.axes.axis["left"].toggle(all=True) + self.axes.axis["bottom"].get_helper().nth_coord_ticks = 1 + self.axes.axis["bottom"].toggle(all=True) + self.axes.axis["top"].get_helper().nth_coord_ticks = 0 + self.axes.axis["top"].toggle(all=False) + self.axes.axis["right"].get_helper().nth_coord_ticks = 1 + self.axes.axis["right"].toggle(all=True) + self.axes.gridlines.set_linestyle("solid") + + self.figure.add_subplot(self.axes) + + # Configure default axes. + axis = self.axes.axis["left"] + axis.major_ticklabels.set_fontsize(10) + axis.major_ticklabels.set_va("baseline") + axis.major_ticklabels.set_rotation(135) + axis = self.axes.axis["right"] + axis.major_ticklabels.set_fontsize(10) + axis.major_ticklabels.set_va("baseline") + axis.major_ticklabels.set_rotation(-135) + self.axes.axis["top"].major_ticklabels.set_fontsize(10) + axis = self.axes.axis["bottom"] + axis.major_ticklabels.set_fontsize(10) + axis.major_ticklabels.set_ha("left") + axis.major_ticklabels.set_va("top") + axis.major_ticklabels.set_rotation(-45) + + # Isotherms: lines of constant temperature (degC). + axis = self.axes.axis["isotherm"] + axis.set_axis_direction("right") + axis.set_axislabel_direction("-") + axis.major_ticklabels.set_rotation(90) + axis.major_ticklabels.set_fontsize(10) + axis.major_ticklabels.set_va("bottom") + axis.major_ticklabels.set_color("grey") + axis.major_ticklabels.set_visible(False) # turned-off + + # Dry adiabats: lines of constant potential temperature (degC). + axis = self.axes.axis["theta"] + axis.set_axis_direction("right") + axis.set_axislabel_direction("+") + axis.major_ticklabels.set_fontsize(10) + axis.major_ticklabels.set_va("bottom") + axis.major_ticklabels.set_color("grey") + axis.major_ticklabels.set_visible(False) # turned-off + axis.line.set_linewidth(3) + axis.line.set_linestyle("--") + + # Lock down the aspect ratio. + self.axes.set_aspect(1.0) + self.axes.grid(True) + + # Initialise the text formatter for the navigation status bar. + self.axes.format_coord = self._status_bar + + # Factor in the tephigram transform. + ISOBAR_TEXT["transform"] = self.transform + WET_ADIABAT_TEXT["transform"] = self.transform + MIXING_RATIO_TEXT["transform"] = self.transform + + # Create plot collections for the tephigram isopleths. + func = partial( + isopleths.isobar, + MIN_THETA, + MAX_THETA, + self.axes, + self.transform, + ISOBAR_LINE, + ) + self._isobars = _PlotCollection( + self.axes, + ISOBAR_SPEC, + MAX_PRESSURE, + func, + ISOBAR_TEXT, + fixed=ISOBAR_FIXED, + minimum=MIN_PRESSURE, + ) + + func = partial( + isopleths.wet_adiabat, + MAX_PRESSURE, + MIN_TEMPERATURE, + self.axes, + self.transform, + WET_ADIABAT_LINE, + ) + self._wet_adiabats = _PlotCollection( + self.axes, + WET_ADIABAT_SPEC, + MAX_WET_ADIABAT, + func, + WET_ADIABAT_TEXT, + fixed=WET_ADIABAT_FIXED, + minimum=MIN_WET_ADIABAT, + xfocus=True, + ) + + func = partial( + isopleths.mixing_ratio, + MIN_PRESSURE, + MAX_PRESSURE, + self.axes, + self.transform, + MIXING_RATIO_LINE, + ) + self._mixing_ratios = _PlotCollection( + self.axes, + MIXING_RATIO_SPEC, + MIXING_RATIOS, + func, + MIXING_RATIO_TEXT, + fixed=MIXING_RATIO_FIXED, + ) + + # Initialise for the tephigram plot event handler. + plt.connect("motion_notify_event", _handler) + self.axes.tephigram = True + self.axes.tephigram_original_delta_xlim = DEFAULT_WIDTH + self.original_delta_xlim = DEFAULT_WIDTH + self.axes.tephigram_transform = self.tephi_transform + self.axes.tephigram_inverse = self.tephi_transform.inverted() + self.axes.tephigram_isopleths = [ + self._isobars, + self._wet_adiabats, + self._mixing_ratios, + ] + + # The tephigram profiles. + self._profiles = [] + self.axes.tephigram_profiles = self._profiles + + # Center the plot around the anchor extent. + self._anchor = anchor + if self._anchor is not None: + self._anchor = np.asarray(anchor) + if ( + self._anchor.ndim != 2 + or self._anchor.shape[-1] != 2 + or len(self._anchor) != 2 + ): + msg = ( + "Invalid anchor, expecting [(bottom-right-pressure, " + "bottom-right-temperature), (top-left-pressure, " + "top-left-temperature)]" + ) + raise ValueError(msg) + ( + (bottom_pressure, bottom_temp), + (top_pressure, top_temp), + ) = self._anchor + + if (bottom_pressure - top_pressure) < 0: + raise ValueError("Invalid anchor pressure range") + if (bottom_temp - top_temp) < 0: + raise ValueError("Invalid anchor temperature range") + + self._anchor = isopleths.Profile(anchor, self.axes) + self._anchor.plot(visible=False) + xlim, ylim = self._calculate_extents() + self.axes.set_xlim(xlim) + self.axes.set_ylim(ylim) + + def plot(self, data, **kwargs): + """ + Plot the environmental lapse rate profile of the pressure and + temperature data points. + + The pressure and temperature data points are transformed into + potential temperature and temperature data points before plotting. + + By default, the tephigram will automatically center the plot around + all profiles. + + .. warning:: + Pressure data points must be in units of mb or hPa, and temperature + data points must be in units of degC. + + Args: + + * data: (pressure, temperature) pair data points. + + .. note:: + All keyword arguments are passed through to + :func:`matplotlib.pyplot.plot`. + + For example: + + .. plot:: + :include-source: + + import matplotlib.pyplot as plt + from tephi import Tephigram + + tpg = Tephigram() + data = [[1006, 26.4], [924, 20.3], [900, 19.8], + [850, 14.5], [800, 12.9], [755, 8.3]] + profile = tpg.plot(data, color='red', linestyle='--', + linewidth=2, marker='o') + barbs = [(10, 45, 900), (20, 60, 850), (25, 90, 800)] + profile.barbs(barbs) + plt.show() + + For associating wind barbs with an environmental lapse rate profile, + see :meth:`~tephi.isopleths.Profile.barbs`. + + """ + profile = isopleths.Profile(data, self.axes) + profile.plot(**kwargs) + self._profiles.append(profile) + + # Center the tephigram plot around all the profiles. + if self._anchor is None: + xlim, ylim = self._calculate_extents(xfactor=0.25, yfactor=0.05) + self.axes.set_xlim(xlim) + self.axes.set_ylim(ylim) + + # Refresh the tephigram plot isopleths. + _refresh_isopleths(self.axes) + + # Show the plot legend. + if "label" in kwargs: + font_properties = FontProperties(size="x-small") + plt.legend( + loc="upper left", + fancybox=True, + shadow=True, + prop=font_properties, + ) + + return profile + + def _status_bar(self, x_point, y_point): + """Generate text for the interactive backend navigation status bar.""" + + temperature, theta = transforms.convert_xy2Tt(x_point, y_point) + pressure, _ = transforms.convert_Tt2pT(temperature, theta) + xlim = self.axes.get_xlim() + zoom = (xlim[1] - xlim[0]) / self.original_delta_xlim + msg = "T:{:.2f}, theta:{:.2f}, phi:{:.2f} (zoom:{:.3f})" + text = msg.format( + float(temperature), float(theta), float(pressure), zoom + ) + + return text + + def _calculate_extents(self, xfactor=None, yfactor=None): + min_x = min_y = 1e10 + max_x = max_y = -1e-10 + profiles = self._profiles + transform = self.tephi_transform.transform + + if self._anchor is not None: + profiles = [self._anchor] + + for profile in profiles: + temperature = profile.temperature.reshape(-1, 1) + theta = profile.theta.reshape(-1, 1) + xy_points = transform(np.concatenate((temperature, theta), axis=1)) + x_points = xy_points[:, 0] + y_points = xy_points[:, 1] + min_x = np.min([min_x, np.min(x_points)]) + min_y = np.min([min_y, np.min(y_points)]) + max_x = np.max([max_x, np.max(x_points)]) + max_y = np.max([max_y, np.max(y_points)]) + + if xfactor is not None: + delta_x = max_x - min_x + min_x, max_x = min_x - xfactor * delta_x, max_x + xfactor * delta_x + + if yfactor is not None: + delta_y = max_y - min_y + min_y, max_y = min_y - yfactor * delta_y, max_y + yfactor * delta_y + + return ([min_x, max_x], [min_y, max_y]) diff --git a/tephi/transforms.py b/tephi/transforms.py index a904cf8..7183e66 100644 --- a/tephi/transforms.py +++ b/tephi/transforms.py @@ -6,15 +6,11 @@ Tephigram transform support. """ + from matplotlib.transforms import Transform import numpy as np -from ._constants import CONST_K, CONST_KELVIN, CONST_L, CONST_MA, CONST_RV - - -# -# Reference: http://www-nwp/~hadaa/tephigram/tephi_plot.html -# +import tephi.constants as constants def convert_Tt2pT(temperature, theta): @@ -37,12 +33,12 @@ def convert_Tt2pT(temperature, theta): temperature, theta = np.asarray(temperature), np.asarray(theta) # Convert temperature and theta from degC to kelvin. - kelvin = temperature + CONST_KELVIN - theta = theta + CONST_KELVIN + kelvin = temperature + constants.KELVIN + theta = theta + constants.KELVIN # Calculate the associated pressure given the temperature and # potential temperature. - pressure = 1000.0 * np.power(kelvin / theta, 1 / CONST_K) + pressure = constants.P_BASE * np.power(kelvin / theta, 1 / constants.K) return pressure, temperature @@ -67,13 +63,13 @@ def convert_pT2Tt(pressure, temperature): pressure, temperature = np.asarray(pressure), np.asarray(temperature) # Convert temperature from degC to kelvin. - kelvin = temperature + CONST_KELVIN + kelvin = temperature + constants.KELVIN # Calculate the potential temperature given the pressure and temperature. - theta = kelvin * ((1000.0 / pressure) ** CONST_K) + theta = kelvin * ((constants.P_BASE / pressure) ** constants.K) # Convert potential temperature from kelvin to degC. - return temperature, theta - CONST_KELVIN + return temperature, theta - constants.KELVIN def convert_pt2pT(pressure, theta): @@ -95,14 +91,14 @@ def convert_pt2pT(pressure, theta): pressure, theta = np.asarray(pressure), np.asarray(theta) # Convert potential temperature from degC to kelvin. - theta = theta + CONST_KELVIN + theta = theta + constants.KELVIN - # Calculate the temperature given the pressure and - # potential temperature. - kelvin = theta * (pressure**CONST_K) / (1000.0**CONST_K) + # Calculate the temperature given the pressure and potential temperature. + denom = constants.P_BASE**constants.K + kelvin = theta * (pressure**constants.K) / denom # Convert temperature from kelvin to degC. - return pressure, kelvin - CONST_KELVIN + return pressure, kelvin - constants.KELVIN def convert_Tt2xy(temperature, theta): @@ -125,13 +121,13 @@ def convert_Tt2xy(temperature, theta): temperature, theta = np.asarray(temperature), np.asarray(theta) # Convert potential temperature from degC to kelvin. - theta = theta + CONST_KELVIN + theta = theta + constants.KELVIN theta = np.clip(theta, 1, 1e10) phi = np.log(theta) - x_data = phi * CONST_MA + temperature - y_data = phi * CONST_MA - temperature + x_data = phi * constants.MA + temperature + y_data = phi * constants.MA - temperature return x_data, y_data @@ -155,10 +151,10 @@ def convert_xy2Tt(x_data, y_data): """ x_data, y_data = np.asarray(x_data), np.asarray(y_data) - phi = (x_data + y_data) / (2 * CONST_MA) + phi = (x_data + y_data) / (2 * constants.MA) temperature = (x_data - y_data) / 2.0 - theta = np.exp(phi) - CONST_KELVIN + theta = np.exp(phi) - constants.KELVIN return temperature, theta @@ -173,21 +169,22 @@ def convert_pw2T(pressure, mixing_ratio): Pressure in mb in hPa. * mixing_ratio: - Dimensionless mixing ratios. + Mixing ratio in g kg-1. Returns: Temperature in degC. """ - pressure = np.array(pressure) + pressure = np.asarray(pressure) # Calculate the dew-point. - vapp = pressure * (8.0 / 5.0) * (mixing_ratio / 1000.0) + vapp = pressure * (8.0 / 5.0) * (mixing_ratio / constants.P_BASE) temp = 1.0 / ( - (1.0 / CONST_KELVIN) - ((CONST_RV / CONST_L) * np.log(vapp / 6.11)) + (1.0 / constants.KELVIN) + - ((constants.Rv / constants.L) * np.log(vapp / 6.11)) ) - return temp - CONST_KELVIN + return temp - constants.KELVIN class TephiTransform(Transform): @@ -214,7 +211,7 @@ def transform_non_affine(self, values): """ return np.concatenate( - convert_Tt2xy(values[:, 0:1], values[:, 1:2]), axis=1 + convert_Tt2xy(values[:, 0:1], values[:, 1:2]), axis=-1 ) def inverted(self): @@ -247,7 +244,7 @@ def transform_non_affine(self, values): """ return np.concatenate( - convert_xy2Tt(values[:, 0:1], values[:, 1:2]), axis=1 + convert_xy2Tt(values[:, 0:1], values[:, 1:2]), axis=-1 ) def inverted(self):