Source code for epygram.fields.H1DField

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) Météo France (2014-)
# This software is governed by the CeCILL-C license under French law.
# http://www.cecill.info
"""
Contains the class that handle a Horizontal 2D field.
"""

from __future__ import print_function, absolute_import, unicode_literals, division

import numpy
import datetime

import footprints
from bronx.graphics.axes import set_figax, set_nice_time_axis

from epygram import config, util, epygramError
from epygram.base import FieldSet
from epygram.geometries import H1DGeometry
from .D3Field import D3Field

epylog = footprints.loggers.getLogger(__name__)


[docs]class H1DField(D3Field): """ Horizontal 1-Dimensions field class. A field is defined by its identifier 'fid', its data, its geometry (gridpoint and optionally spectral), and its validity. The natural being of a field is gridpoint, so that: a field always has a gridpoint geometry, but it has a spectral geometry only in case it is spectral. """ _collector = ('field',) _footprint = dict( attr=dict( structure=dict( values=set(['H1D'])), geometry=dict( type=H1DGeometry), ) ) ############## # ABOUT DATA # ##############
[docs] def getlevel(self, level=None, k=None): """Returns self. Useless but for compatibility reasons.""" if k is None and level is None: raise epygramError("You must give k or level.") if k is not None and level is not None: raise epygramError("You cannot give, at the same time, k and level") if level is not None: if level not in self.geometry.vcoordinate.levels: raise epygramError("The requested level does not exist.") return self
def _check_operands(self, other): """ Internal method to check compatibility of terms in operations on fields. """ if isinstance(other, self.__class__): if self.spectral != other.spectral: raise epygramError("cannot operate a spectral field with a" + " non-spectral field.") if self.geometry.rectangular_grid: if self.geometry.dimensions['X'] != other.geometry.dimensions['X'] or\ self.geometry.dimensions['Y'] != other.geometry.dimensions['Y']: raise epygramError("operations on fields cannot be done if" + " fields do not share their gridpoint" + " dimensions.") else: if self.geometry.dimensions != other.geometry.dimensions: raise epygramError("operations on fields cannot be done if" + " fields do not share their gridpoint" + " dimensions.") if self.spectral_geometry != other.spectral_geometry: raise epygramError("operations on fields cannot be done if" + " fields do not share their spectral" + " geometry.") else: super(D3Field, self)._check_operands(other) ################### # PRE-APPLICATIVE # ################### # (but useful and rather standard) ! # [so that, subject to continuation through updated versions, # including suggestions/developments by users...]
[docs] def plotfield(self, *args, **kwargs): """ Interface method to methods plotprofiles() and plotverticalhovmoller(), depending on time dimension (or not). Cf. these functions for arguments. """ if len(self.validity) == 1: return self.plottransects(*args, **kwargs) else: return self.plothorizontalhovmoller(*args, **kwargs)
[docs] def plottransects(self, *args, **kwargs): """Cf. eponymous function of module for arguments.""" return plottransects(self, *args, **kwargs)
[docs] def plothorizontalhovmoller(self, *args, **kwargs): """Cf. eponymous function of module for arguments.""" return plothorizontalhovmoller(self, *args, **kwargs)
[docs] def plotanimation(self, *args, **kwargs): """Cf. eponymous function of module for arguments.""" return plotanimation(self, *args, **kwargs)
# FUNCTIONS # #############
[docs]def plothorizontalhovmoller(transect, over=(None, None), fidkey=None, title=None, logscale=False, zoom=None, colorbar='vertical', graphicmode='colorshades', minmax=None, levelsnumber=21, center_cmap_on_0=False, colormap='jet', minmax_in_title=True, contourcolor='k', contourwidth=1, contourlabel=True, datefmt=None, showgrid=True, x_is='distance', figsize=None, rcparams=None): """ Makes a simple vertical Hovmöller plot of the field. :param transect: being a :class:`epygram.fields.H1DField` :param over: any existing figure and/or ax to be used for the plot, given as a tuple (fig, ax), with None for missing objects. *fig* is the frame of the matplotlib figure, containing eventually several subplots (axes); *ax* is the matplotlib axes on which the drawing is done. When given (!= None), these objects must be coherent, i.e. ax being one of the fig axes. :param fidkey: type of fid for entitling the plot with *fid[fidkey]*, if title is *None*; if *None*, labels with raw fid. :param title: title for the plot. :param logscale: to set Y logarithmic scale :param zoom: a dict containing optional limits to zoom on the plot. \n Syntax: e.g. {'ymax':500, ...}. :param colorbar: if *False*, hide colorbar the plot; else, befines the colorbar orientation, among ('horizontal', 'vertical'). Defaults to 'vertical'. :param graphicmode: among ('colorshades', 'contourlines'). :param minmax: defines the min and max values for the plot colorbar. \n Syntax: [min, max]. [0.0, max] also works. Default is min/max of the field. :param levelsnumber: number of levels for contours and colorbar. :param center_cmap_on_0: aligns the colormap center on the value 0. :param colormap: name of the **matplotlib** colormap to use. :param minmax_in_title: if True and minmax != None, adds min and max values in title :param contourcolor: color or colormap to be used for 'contourlines' graphicmode. It can be either a legal html color name, or a colormap name. :param contourwidth: width of contours for 'contourlines' graphicmode. :param contourlabel: displays labels on contours. :param datefmt: date format to use, e.g. "%Y-%m-%d %H:%M:%S %Z" :param showgrid: True/False to show grid or not :param x_is: abscissa to be among: - 'distance': distance from first point of transect - 'lon': longitude of points - 'lat': latitude of points :param figsize: figure sizes in inches, e.g. (5, 8.5). Default figsize is config.plotsizes. :param rcparams: list of (*args, **kwargs) to be passed to pyplot.rc() defaults to [(('font',), dict(family='serif')),] Warning: requires **matplotlib**. """ import matplotlib.pyplot as plt import matplotlib.dates as mdates from mpl_toolkits.axes_grid1 import make_axes_locatable if rcparams is None: rcparams = [(('font',), dict(family='serif')),] for args, kwargs in rcparams: plt.rc(*args, **kwargs) if figsize is None: figsize = config.plotsizes # User colormaps if colormap not in plt.colormaps(): util.load_cmap(colormap) # Figure, ax fig, ax = set_figax(*over, figsize=figsize) # coords y = numpy.zeros((len(transect.validity), transect.geometry.dimensions['X'])) validities = {transect.validity[i].get():i for i in range(len(transect.validity))} yaxis_label = 'Validity' if len(validities) == 1 and len(transect.validity) != 1: yaxis_label = 'Basis' validities = {transect.validity[i].getbasis():i for i in len(transect.validity)} epoch = datetime.datetime(1970, 1, 1) for i in range(len(transect.validity)): d = transect.validity[i].get() if yaxis_label == 'Validity' else transect.validity[i].getbasis() timedelta = d - epoch p = (timedelta.microseconds + (timedelta.seconds + timedelta.days * 24 * 3600) * 10 ** 6) / 1e6 y[i, :] = mdates.epoch2num(p) x = numpy.zeros((len(transect.validity), transect.geometry.dimensions['X'])) lonlat = zip(*transect.geometry.get_lonlat_grid()) p0 = lonlat[0] plast = p0 distance = 0 for i in range(transect.geometry.dimensions['X']): if x_is == 'distance': p = lonlat[i] distance += transect.geometry.distance((plast[0], plast[1]), (p[0], p[1])) x[:, i] = distance plast = p elif x_is == 'lon': x[:, i] = lonlat[i][0] elif x_is == 'lat': x[:, i] = lonlat[i][1] plast = lonlat[-1] data = transect.getdata() # min/max m = data.min() M = data.max() if minmax is not None: if minmax_in_title: minmax_in_title = '(min: ' + \ '{: .{precision}{type}}'.format(m, type='E', precision=3) + \ ' // max: ' + \ '{: .{precision}{type}}'.format(M, type='E', precision=3) + ')' try: m = float(minmax[0]) except Exception: m = data.min() try: M = float(minmax[1]) except Exception: M = data.max() else: minmax_in_title = '' if abs(m - M) > config.epsilon: levels = numpy.linspace(m, M, levelsnumber) vmin = vmax = None if center_cmap_on_0: vmax = max(abs(m), M) vmin = -vmax else: raise epygramError("cannot plot uniform field.") L = int((levelsnumber - 1) // 15) + 1 hlevels = [levels[l] for l in range(len(levels) - L / 3) if l % L == 0] + [levels[-1]] # plot if logscale: ax.set_yscale('log') ax.grid() if graphicmode == 'colorshades': pf = ax.contourf(x, y, data, levels, cmap=colormap, vmin=vmin, vmax=vmax) if colorbar: position = 'right' if colorbar == 'vertical' else 'bottom' cax = make_axes_locatable(ax).append_axes(position, size="5%", pad=0.1) cb = plt.colorbar(pf, orientation=colorbar, ticks=hlevels, cax=cax) if minmax_in_title != '': cb.set_label(minmax_in_title) elif graphicmode == 'contourlines': pf = ax.contour(x, y, data, levels=levels, colors=contourcolor, linewidths=contourwidth) if contourlabel: ax.clabel(pf, colors=contourcolor) # time set_nice_time_axis(ax, 'y', showgrid=showgrid, datefmt=datefmt) # decoration if x_is == 'distance': ax.set_xlabel('Distance from left-end point (m).') elif x_is == 'lon': ax.set_xlabel(u'Longitude (\u00B0).') elif x_is == 'lat': ax.set_xlabel(u'Latitude (\u00B0).') ax.set_ylabel(yaxis_label) if zoom is not None: ykw = {} xkw = {} for pair in (('bottom', 'ymin'), ('top', 'ymax')): try: ykw[pair[0]] = zoom[pair[1]] except Exception: pass for pair in (('left', 'xmin'), ('right', 'xmax')): try: xkw[pair[0]] = zoom[pair[1]] except Exception: pass ax.set_ylim(**ykw) ax.set_xlim(**xkw) if title is None: if fidkey is None: fid = transect.fid[sorted(transect.fid.keys())[0]] else: fid = transect.fid[fidkey] title = u'Horizontal Hovmöller of ' + str(fid) + ' between \n' + \ '<- (' + str(p0[0]) + ', ' + \ str(p0[1]) + ')' + ' and ' + \ '(' + str(plast[0]) + ', ' + \ str(plast[1]) + ') ->' ax.set_title(title) return (fig, ax)
[docs]def plottransects(transects, over=(None, None), labels=None, fidkey=None, unit='SI', title=None, logscale=False, zoom=None, x_is='distance', figsize=None, rcparams=None): """ To plot a series of transects. Returns a tuple of :mod:`matplotlib` (*Figure*, *ax*). :param transects: being a :class:`epygram.base.FieldSet` of :class:`epygram.fields.H1DField`, or a single :class:`epygram.fields.H1DField`. \n All transects are supposed to have the same unit, and the same horizontal coordinate. :param over: any existing figure and/or ax to be used for the plot, given as a tuple (fig, ax), with None for missing objects. *fig* is the frame of the matplotlib figure, containing eventually several subplots (axes); *ax* is the matplotlib axes on which the drawing is done. When given (!= None), these objects must be coherent, e.g. ax being one of the fig axes. :param labels: a list of labels for the profiles (same length and same order). :param fidkey: key of fid for labelling the curve with *fid[fidkey]*; if *None*, labels with raw fid. :param unit: label for X coordinate. :param title: title for the plot. :param logscale: to set Y logarithmic scale :param zoom: a dict containing optional limits to zoom on the plot. \n Syntax: e.g. {'ymax':500, ...}. :param x_is: abscissa to be among: - 'distance': distance from first point of transect - 'lon': longitude of points - 'lat': latitude of points :param figsize: figure sizes in inches, e.g. (5, 8.5). Default figsize is config.plotsizes. :param rcparams: list of (*args, **kwargs) to be passed to pyplot.rc() defaults to [(('font',), dict(family='serif')), (('figure',), dict(autolayout=True))] """ import matplotlib.pyplot as plt if rcparams is None: rcparams = [(('font',), dict(family='serif')), (('figure',), dict(autolayout=True))] for args, kwargs in rcparams: plt.rc(*args, **kwargs) plt.rc('figure', autolayout=True) if figsize is None: figsize = config.plotsizes colors = ['red', 'blue', 'green', 'orange', 'magenta', 'darkolivegreen', 'yellow', 'salmon', 'black'] linestyles = ['-', '--', '-.', ':'] if not isinstance(transects, FieldSet): p = transects.deepcopy() transects = FieldSet() transects.append(p) t0 = transects[0] x = numpy.zeros((t0.geometry.dimensions['X'],)) lonlat = zip(*t0.geometry.get_lonlat_grid()) p0 = lonlat[0] plast = p0 distance = 0 for i in range(t0.geometry.dimensions['X']): if x_is == 'distance': p = lonlat[i] distance += t0.geometry.distance((plast[0], plast[1]), (p[0], p[1])) x[i] = distance plast = p elif x_is == 'lon': x[i] = lonlat[i][0] elif x_is == 'lat': x[i] = lonlat[i][1] plast = lonlat[-1] # Figure fig, ax = set_figax(*over, figsize=figsize) if logscale: ax.set_yscale('log') i = 0 for t in transects: if len(t.validity) != 1: raise epygramError("plottransects can handle only profiles with one validity.") if labels is not None: label = labels[i] else: if fidkey is not None: label = t.fid.get(fidkey, t.fid) else: label = str(t.fid) data = t.getdata() plot_kwargs = {} if len(transects) > 1: plot_kwargs['color'] = colors[i % len(colors)] plot_kwargs['linestyle'] = linestyles[i // len(colors)] ax.plot(x, data.flatten(), label=label, **plot_kwargs) i += 1 # Decoration if zoom is not None: if 'ymin' in zoom.keys(): ax.set_ylim(bottom=zoom['ymin']) if 'ymax' in zoom.keys(): ax.set_ylim(top=zoom['ymax']) if 'xmin' in zoom.keys(): ax.set_xlim(left=zoom['xmin']) if 'xmax' in zoom.keys(): ax.set_xlim(right=zoom['xmax']) if title is None: title = u'Transect between \n' + \ '<- (' + str(p0[0]) + ', ' + \ str(p0[1]) + ')' + ' and ' + \ '(' + str(plast[0]) + ', ' + \ str(plast[1]) + ') ->' ax.set_title(title) legend = ax.legend(loc='upper right', shadow=True) for label in legend.get_texts(): label.set_fontsize('medium') ax.set_ylabel(r'$' + unit + '$') if x_is == 'distance': ax.set_xlabel('Distance from left-end point (m).') elif x_is == 'lon': ax.set_xlabel(u'Longitude (\u00B0).') elif x_is == 'lat': ax.set_xlabel(u'Latitude (\u00B0).') ax.grid() return (fig, ax)
[docs]def plotanimation(transect, title='__auto__', repeat=False, interval=1000, **kwargs): """ To plot a time-dependent transect as an animation. Returns a :class:`matplotlib.animation.FuncAnimation`. :param transect: being a :class:`epygram.fields.H1DField`. :param title: title for the plot. '__auto__' (default) will print the current validity of the time frame. :param repeat: to repeat animation :param interval: number of milliseconds between two validities Other kwargs passed to plottransects(). """ import matplotlib.animation as animation if len(transect.validity) == 1: raise epygramError("plotanimation can handle only transect with several validities.") if title is not None: if title == '__auto__': title_prefix = '' else: title_prefix = title title = title_prefix + '\n' + transect.validity[0].get().isoformat(sep=' ') else: title_prefix = None transect0 = transect.getvalidity(0) mindata = transect.getdata().min() maxdata = transect.getdata().max() mindata -= (maxdata - mindata) / 10 maxdata += (maxdata - mindata) / 10 zoom = kwargs.get('zoom') zoom = util.ifNone_emptydict(zoom) if 'ymax' not in zoom.keys(): zoom.update(ymax=maxdata) if 'ymin' not in zoom.keys(): zoom.update(ymin=mindata) kwargs['zoom'] = zoom fig, ax = plottransects(transect0, title=title, **kwargs) if kwargs.get('colorbar_over') is None: kwargs['colorbar_over'] = fig.axes[-1] # the last being created, in plotfield() kwargs['over'] = (fig, ax) def update(i, ax, myself, transecti, title_prefix, kwargs): print("updating") if i < len(myself.validity): ax.clear() transecti = myself.getvalidity(i) if title_prefix is not None: title = title_prefix + '\n' + transecti.validity.get().isoformat(sep=b' ') transecti.plotfield(title=title, **kwargs) anim = animation.FuncAnimation(fig, update, fargs=[ax, transect, transect0, title_prefix, kwargs], frames=range(len(transect.validity) + 1), # AM: don't really understand why but needed for the last frame to be shown interval=interval, repeat=repeat) return anim