"""
Climate contour/vector plots using cf-python, matplotlib and cartopy.
Andy Heaps NCAS-CMS September 2018
"""
import cf
import numpy as np
import subprocess
from scipy import interpolate
import time
import matplotlib
from copy import deepcopy
import os
import sys
import matplotlib.pyplot as plot
from matplotlib.path import Path
import matplotlib.patches as mpatches
from matplotlib.collections import PolyCollection
from distutils.version import StrictVersion
import cartopy
import cartopy.crs as ccrs
import cartopy.util as cartopy_util
import cartopy.feature as cfeature
class pvars(object):
def __init__(self, **kwargs):
'''Initialize a new Pvars instance'''
for attr, value in kwargs.iteritems():
setattr(self, attr, value)
def __str__(self):
'''x.__str__() <==> str(x)'''
out = ['%s = %s' % (a, repr(v))]
for a, v in self.__dict__.iteritems():
return '\n'.join(out)
cf_version_min = '1.0.1'
cf_errstr = '\n cf-python > ' + cf_version_min + \
' needs to be installed to use cf-plot \n'
try:
import cf
if StrictVersion(cf.__version__) < StrictVersion(cf_version_min):
raise Warning(cf_errstr)
except ImportError:
raise Warning(cf_errstr)
# Check for a display and use the Agg backing store if none is present
# This is for batch mode processing
try:
disp = os.environ["DISPLAY"]
except:
matplotlib.use('Agg')
# Code to check if the ImageMagick display command is available
def which(program):
def is_exe(fpath):
return os.path.exists(fpath) and os.access(fpath, os.X_OK)
def ext_candidates(fpath):
yield fpath
for ext in os.environ.get("PATHEXT", "").split(os.pathsep):
yield fpath + ext
for path in os.environ["PATH"].split(os.pathsep):
exe_file = os.path.join(path, program)
for candidate in ext_candidates(exe_file):
if is_exe(candidate):
return candidate
return None
# Default colour scales
# cscale1 is a differential data scale - blue to red
cscale1 = ['#0a3278', '#0f4ba5', '#1e6ec8', '#3ca0f0', '#50b4fa', '#82d2ff',
'#a0f0ff', '#c8faff', '#e6ffff', '#fffadc', '#ffe878', '#ffc03c',
'#ffa000', '#ff6000', '#ff3200', '#e11400', '#c00000', '#a50000']
# viridis is a continuous data scale - blue, green, yellow
viridis = ['#440154', '#440255', '#440357', '#450558', '#45065a', '#45085b',
'#46095c', '#460b5e', '#460c5f', '#460e61', '#470f62', '#471163',
'#471265', '#471466', '#471567', '#471669', '#47186a', '#48196b',
'#481a6c', '#481c6e', '#481d6f', '#481e70', '#482071', '#482172',
'#482273', '#482374', '#472575', '#472676', '#472777', '#472878',
'#472a79', '#472b7a', '#472c7b', '#462d7c', '#462f7c', '#46307d',
'#46317e', '#45327f', '#45347f', '#453580', '#453681', '#443781',
'#443982', '#433a83', '#433b83', '#433c84', '#423d84', '#423e85',
'#424085', '#414186', '#414286', '#404387', '#404487', '#3f4587',
'#3f4788', '#3e4888', '#3e4989', '#3d4a89', '#3d4b89', '#3d4c89',
'#3c4d8a', '#3c4e8a', '#3b508a', '#3b518a', '#3a528b', '#3a538b',
'#39548b', '#39558b', '#38568b', '#38578c', '#37588c', '#37598c',
'#365a8c', '#365b8c', '#355c8c', '#355d8c', '#345e8d', '#345f8d',
'#33608d', '#33618d', '#32628d', '#32638d', '#31648d', '#31658d',
'#31668d', '#30678d', '#30688d', '#2f698d', '#2f6a8d', '#2e6b8e',
'#2e6c8e', '#2e6d8e', '#2d6e8e', '#2d6f8e', '#2c708e', '#2c718e',
'#2c728e', '#2b738e', '#2b748e', '#2a758e', '#2a768e', '#2a778e',
'#29788e', '#29798e', '#287a8e', '#287a8e', '#287b8e', '#277c8e',
'#277d8e', '#277e8e', '#267f8e', '#26808e', '#26818e', '#25828e',
'#25838d', '#24848d', '#24858d', '#24868d', '#23878d', '#23888d',
'#23898d', '#22898d', '#228a8d', '#228b8d', '#218c8d', '#218d8c',
'#218e8c', '#208f8c', '#20908c', '#20918c', '#1f928c', '#1f938b',
'#1f948b', '#1f958b', '#1f968b', '#1e978a', '#1e988a', '#1e998a',
'#1e998a', '#1e9a89', '#1e9b89', '#1e9c89', '#1e9d88', '#1e9e88',
'#1e9f88', '#1ea087', '#1fa187', '#1fa286', '#1fa386', '#20a485',
'#20a585', '#21a685', '#21a784', '#22a784', '#23a883', '#23a982',
'#24aa82', '#25ab81', '#26ac81', '#27ad80', '#28ae7f', '#29af7f',
'#2ab07e', '#2bb17d', '#2cb17d', '#2eb27c', '#2fb37b', '#30b47a',
'#32b57a', '#33b679', '#35b778', '#36b877', '#38b976', '#39b976',
'#3bba75', '#3dbb74', '#3ebc73', '#40bd72', '#42be71', '#44be70',
'#45bf6f', '#47c06e', '#49c16d', '#4bc26c', '#4dc26b', '#4fc369',
'#51c468', '#53c567', '#55c666', '#57c665', '#59c764', '#5bc862',
'#5ec961', '#60c960', '#62ca5f', '#64cb5d', '#67cc5c', '#69cc5b',
'#6bcd59', '#6dce58', '#70ce56', '#72cf55', '#74d054', '#77d052',
'#79d151', '#7cd24f', '#7ed24e', '#81d34c', '#83d34b', '#86d449',
'#88d547', '#8bd546', '#8dd644', '#90d643', '#92d741', '#95d73f',
'#97d83e', '#9ad83c', '#9dd93a', '#9fd938', '#a2da37', '#a5da35',
'#a7db33', '#aadb32', '#addc30', '#afdc2e', '#b2dd2c', '#b5dd2b',
'#b7dd29', '#bade27', '#bdde26', '#bfdf24', '#c2df22', '#c5df21',
'#c7e01f', '#cae01e', '#cde01d', '#cfe11c', '#d2e11b', '#d4e11a',
'#d7e219', '#dae218', '#dce218', '#dfe318', '#e1e318', '#e4e318',
'#e7e419', '#e9e419', '#ece41a', '#eee51b', '#f1e51c', '#f3e51e',
'#f6e61f', '#f8e621', '#fae622', '#fde724']
# Read in defaults if they exist and overlay
# for contour options of fill, blockfill and lines
global_fill=True
global_lines=True
global_blockfill=False
global_degsym=False
defaults_file=os.path.expanduser("~")+'/.cfplot_defaults'
if os.path.exists(defaults_file):
with open(defaults_file) as file:
for line in file:
vals=line.split(' ')
com, val = vals
v=False
if val.strip() == 'True':
v=True
if com == 'blockfill':
global_blockfill=v
if com == 'lines':
global_lines=v
if com == 'fill':
global_fill=v
if com == 'degsym':
global_degsym=v
#####################################
# plotvars - global plotting variables
#####################################
plotvars = pvars(lonmin=-180, lonmax=180, latmin=-90, latmax=90, proj='cyl',
resolution='110m', plot_type=1, boundinglat=0, lon_0=0,
levels=None,
levels_min=None, levels_max=None, levels_step=None,
norm=None, levels_extend='both', xmin=None,
xmax=None, ymin=None, ymax=None, xlog=None, ylog=None,
rows=1, columns=1, file=None, orientation='landscape',
user_mapset=0, user_gset=0, cscale_flag=0, user_levs=0,
user_plot=0, master_plot=None, plot=None, cs=cscale1,
cs_user='cscale1', mymap=None, xticks=None, yticks=None,
xticklabels=None, yticklabels=None, xstep=None, ystep=None,
xlabel=None, ylabel=None, title=None, title_fontsize=15,
axis_label_fontsize=11, text_fontsize=11,
text_fontweight='normal', axis_label_fontweight='normal',
colorbar_fontsize=11, colorbar_fontweight='normal',
title_fontweight='normal', continent_thickness=None,
continent_color=None, continent_linestyle=None,
pos=1, viewer='display',
tspace_year=None, tspace_month=None, tspace_day=None,
tspace_hour=None, xtick_label_rotation=0,
xtick_label_align='center', ytick_label_rotation=0,
ytick_label_align='right', legend_text_size=11,
legend_text_weight='normal', tight=False,
cs_uniform=True, master_title=None,
master_title_location=[0.5, 0.95], master_title_fontsize=30,
master_title_fontweight='normal', dpi=None,
plot_xmin=None, plot_xmax=None, plot_ymin=None,
plot_ymax=None, land_color=None, ocean_color=None,
lake_color=None, twinx=False, twiny=False,
rotated_grid_thickness=2.0, rotated_grid_spacing=10,
rotated_deg_spacing=0.75, rotated_continents=True,
rotated_grid=True, rotated_labels=True,
legend_frame=True, legend_frame_edge_color='k',
legend_frame_face_color=None, degsym=global_degsym,
axis_width=None, grid=True, grid_spacing=1,
grid_lons=None, grid_lats=None,
grid_colour='grey', grid_linestyle='--',
grid_thickness=1.0)
# Check for iPython notebook inline
# and set the viewer to None if found
is_inline = 'inline' in matplotlib.get_backend()
if is_inline:
plotvars.viewer=None
[docs]def con(f=None, x=None, y=None, fill=global_fill, lines=global_lines, line_labels=True,
title=None, colorbar_title=None, colorbar=True,
colorbar_label_skip=None, ptype=0, negative_linestyle='solid',
blockfill=global_blockfill, zero_thick=False, colorbar_shrink=None,
colorbar_orientation=None, colorbar_position=None, xlog=False,
ylog=False, axes=True, xaxis=True, yaxis=True, xticks=None,
xticklabels=None, yticks=None, yticklabels=None, xlabel=None,
ylabel=None, colors='k', swap_axes=False, verbose=None,
linewidths=None, alpha=1.0, colorbar_text_up_down=False,
colorbar_text_down_up=False, colorbar_drawedges=True,
linestyles=None, zorder=None):
"""
| con is the interface to contouring in cf-plot. The minimum use is con(f)
| where f is a 2 dimensional array. If a cf field is passed then an
| appropriate plot will be produced i.e. a longitude-latitude or
| latitude-height plot for example. If a 2d numeric array is passed then
| the optional arrays x and y can be used to describe the x and y data
| point locations.
|
| f - array to contour
| x - x locations of data in f (optional)
| y - y locations of data in f (optional)
| fill=True - colour fill contours
| lines=True - draw contour lines and labels
| line_labels=True - label contour lines
| title=title - title for the plot
| ptype=0 - plot type - not needed for cf fields.
| 0 = no specific plot type,
| 1 = longitude-latitude,
| 2 = latitude - height,
| 3 = longitude - height,
| 4 = latitude - time,
| 5 = longitude - time
| 6 = rotated pole
| negative_linestyle='solid' - set to one of 'solid', 'dashed'
| zero_thick=False - add a thick zero contour line. Set to 3 for example.
| blockfill=False - set to True for a blockfill plot
| colorbar_title=colbar_title - title for the colour bar
| colorbar=1 - add a colour bar if a filled contour plot
| colorbar_label_skip=None - skip colour bar labels. Set to 2 to skip
| every other label.
| colorbar_orientation=None - options are 'horizontal' and 'vertical'
| The default for most plots is horizontal but
| for polar stereographic plots this is vertical.
| colorbar_shrink=None - value to shrink the colorbar by. If the colorbar
| exceeds the plot area then values of 1.0, 0.55
| or 0.5m ay help it better fit the plot area.
| colorbar_position=None - position of colorbar
| [xmin, ymin, x_extent,y_extent] in normalised
| coordinates. Use when a common colorbar
| is required for a set of plots. A typical set
| of values would be [0.1, 0.05, 0.8, 0.02]
| colorbar_text_up_down=False - if True horizontal colour bar labels alternate
| above (start) and below the colour bar
| colorbar_text_down_up=False - if True horizontal colour bar labels alternate
| below (start) and above the colour bar
| colorbar_drawedges=True - Draw internal divisions in the colorbar
| colors='k' - contour line colors - takes one or many values.
| xlog=False - logarithmic x axis
| ylog=False - logarithmic y axis
| axes=True - plot x and y axes
| xaxis=True - plot xaxis
| yaxis=True - plot y axis
| xticks=None - xtick positions
| xticklabels=None - xtick labels
| yticks=None - y tick positions
| yticklabels=None - ytick labels
| xlabel=None - label for x axis
| ylabel=None - label for y axis
| swap_axes=False - swap plotted axes - only valid for X, Y, Z vs T plots
| verbose=None - change to 1 to get a verbose listing of what con
| is doing
| linewidths=None - contour linewidths. Either a single number for all
| lines or array of widths
| linestyles=None - takes 'solid', 'dashed', 'dashdot' or 'dotted'
| alpha=1.0 - transparency setting. The default is no transparency.
| colorbar_text_up_down=False - on a horizontal colorbar alternate the
| labels top and bottom starting in the up position
| colorbar_text_down_up=False - on a horizontal colorbar alternate the
| labels bottom and top starting in the bottom position
| colorbar_drawedges=True - draw internal delimeter lines in the colorbar
| zorder=None - order of drawing
:Returns:
None
"""
# Turn off divide warning in contour routine which is a numpy issue
old_settings = np.seterr(all='ignore')
np.seterr(divide='ignore')
# Set potential user axis labels
user_xlabel = xlabel
user_ylabel = ylabel
# Extract required data for contouring
# If a cf-python field
if isinstance(f, cf.Field):
# Extract data
if verbose:
print 'con - calling cf_data_assign'
field, x, y, ptype, colorbar_title, xlabel, ylabel, xpole, ypole =\
cf_data_assign(f, colorbar_title, verbose=verbose)
if user_xlabel is not None:
xlabel = user_xlabel
if user_ylabel is not None:
ylabel = user_ylabel
elif isinstance(f, cf.FieldList):
raise TypeError("Can't plot a field list")
else:
if verbose:
print 'con - using user assigned data'
field = f # field data passed in as f
if x is None:
x = np.arange(np.shape(field)[1])
if y is None:
y = np.arange(np.shape(field)[0])
check_data(field, x, y)
xlabel = ''
ylabel = ''
# Set contour line styles
matplotlib.rcParams['contour.negative_linestyle'] = negative_linestyle
# Set contour lines off on block plots
if blockfill:
fill = False
field_orig = field
x_orig = x
y_orig = y
# Check number of colours and levels match if user has modified the
# number of colours
if plotvars.cscale_flag == 2:
if plotvars.user_levs == 0:
errstr = "\n\nBlockfill error - need to match number of "
errstr += "colours and contour intervals\n"
raise TypeError(errstr)
ncols=np.size(plotvars.cs)
nintervals=np.size(plotvars.levels)-1
if plotvars.levels_extend == 'min':
nintervals += 1
if plotvars.levels_extend == 'max':
nintervals += 1
if plotvars.levels_extend == 'both':
nintervals += 2
if ncols != nintervals:
errstr = "\n\nBlockfill error - need to match number of "
errstr += "colours and contour intervals\n"
errstr += "Don't forget to take account of the colorbar "
errstr += "extensions"
raise TypeError(errstr)
# Turn off colorbar if fill is turned off
if not fill and not blockfill:
colorbar = False
# Revert to default colour scale if cscale_flag flag is set
if plotvars.cscale_flag == 0:
plotvars.cs = cscale1
# Set the orientation of the colorbar
if plotvars.plot_type == 1:
if plotvars.proj == 'npstere' or plotvars.proj == 'spstere':
if colorbar_orientation is None:
colorbar_orientation = 'vertical'
if colorbar_orientation is None:
colorbar_orientation = 'horizontal'
# Store original map resolution
resolution_orig=plotvars.resolution
proj = plotvars.proj
# Set size of color bar if not specified
if colorbar_shrink is None:
colorbar_shrink = 1.0
if plotvars.plot_type == 1:
scale = (plotvars.lonmax - plotvars.lonmin) / \
(plotvars.latmax - plotvars.latmin)
if scale <= 1:
colorbar_shrink = 0.55
if plotvars.orientation == 'landscape':
if (proj == 'cyl' and colorbar_orientation == 'vertical'):
colorbar_shrink = 0.5
if (proj == 'cyl' and colorbar_orientation == 'horizontal'):
colorbar_shrink = 1.0
if plotvars.orientation == 'portrait':
if (proj == 'cyl' and colorbar_orientation == 'vertical'):
colorbar_shrink = 0.25
if (proj == 'cyl' and colorbar_orientation == 'horizontal'):
colorbar_shrink = 1.0
if plotvars.proj == 'npstere' or plotvars.proj == 'spstere':
if plotvars.orientation == 'landscape':
if colorbar_orientation == 'horizontal':
colorbar_shrink = 1.0
if colorbar_orientation == 'vertical':
colorbar_shrink = 1.0
if plotvars.orientation == 'portrait':
if colorbar_orientation == 'horizontal':
colorbar_shrink = 1.0
if colorbar_orientation == 'vertical':
colorbar_shrink = 1.0
# Set plot type if user specified
if (ptype is not None):
plotvars.plot_type = ptype
# Get contour levels
includes_zero = 0
if plotvars.user_levs == 1:
# User defined
if verbose:
print 'con - using user defined contour levels'
clevs = plotvars.levels
mult = 0
fmult = 1
else:
if plotvars.levels_step is None:
# Automatic levels
if verbose:
print 'con - generating automatic contour levels'
dmin = np.nanmin(field)
dmax = np.nanmax(field)
clevs, mult = gvals(dmin=dmin, dmax=dmax, tight=0)
fmult = 10**-mult
else:
# Use step to generate the levels
step = plotvars.levels_step
dmin = np.nanmin(field)
dmax = np.nanmax(field)
if isinstance(step, (int, long)):
dmin = int(np.nanmin(field))
dmax = int(np.nanmax(field))
fmult = 1
mult = 0
clevs = []
if dmin < 0:
clevs = ((np.arange(-1*dmin/step+1)*-step)[::-1])
if dmax > 0:
if np.size(clevs) > 0:
clevs = np.concatenate((clevs[:-1], np.arange(dmax/step+1)*step))
else:
clevs = np.arange(dmax/step+1)*step
# Strip out any outlying values
pts = np.where(np.logical_and(clevs >= dmin, clevs <= dmax))
clevs = clevs[pts]
# Throw an error if less than two levels
if np.size(clevs) < 2:
errstr = "\n\ncfp.con error - need more than two levels "
errstr += "to make a contour plot\n"
raise TypeError(errstr)
# Set the colour scale
# Nothing defined
if plotvars.cscale_flag == 0:
includes_zero = 0
col_zero = 0
for cval in clevs:
if includes_zero == 0:
col_zero = col_zero + 1
if cval == 0:
includes_zero = 1
if includes_zero == 1:
cs_below = col_zero
cs_above = np.size(clevs) - col_zero + 1
if plotvars.levels_extend == 'max' or plotvars.levels_extend == 'neither':
cs_below = cs_below - 1
if plotvars.levels_extend == 'min' or plotvars.levels_extend == 'neither':
cs_above = cs_above - 1
uniform = True
if plotvars.cs_uniform is False:
uniform = False
cscale('scale1', below=cs_below, above=cs_above, uniform=uniform)
else:
ncols=ncols=np.size(clevs)+1
if plotvars.levels_extend == 'min' or plotvars.levels_extend == 'max':
ncols=ncols-1
if plotvars.levels_extend == 'neither':
ncols=ncols-2
cscale('viridis', ncols=ncols)
plotvars.cscale_flag = 0
# User selected colour map but no mods so fit to levels
if plotvars.cscale_flag == 1:
ncols=ncols=np.size(clevs)+1
if plotvars.levels_extend == 'min' or plotvars.levels_extend == 'max':
ncols=ncols-1
if plotvars.levels_extend == 'neither':
ncols=ncols-2
cscale(plotvars.cs_user, ncols=ncols)
plotvars.cscale_flag = 1
# Set colorbar labels
# Set a sensible label spacing if the user hasn't already done so
if colorbar_label_skip is None:
if colorbar_orientation == 'horizontal':
nchars = 0
for lev in clevs:
nchars = nchars + len(str(lev))
colorbar_label_skip = nchars / 80 + 1
if plotvars.columns > 1:
colorbar_label_skip = nchars * (plotvars.columns) / 80
else:
colorbar_label_skip = 1
if colorbar_label_skip > 1:
if includes_zero:
# include zero in the colorbar labels
zero_pos = [i for i, mylev in enumerate(clevs) if mylev == 0][0]
colorbar_labels = clevs[zero_pos]
i = zero_pos + colorbar_label_skip
while i <= len(clevs) - 1:
colorbar_labels = np.append(colorbar_labels, clevs[i])
i = i + colorbar_label_skip
i = zero_pos - colorbar_label_skip
if i >= 0:
while i >= 0:
colorbar_labels = np.append(clevs[i], colorbar_labels)
i = i - colorbar_label_skip
else:
colorbar_labels = clevs[0]
i = colorbar_label_skip
while i <= len(clevs) - 1:
colorbar_labels = np.append(colorbar_labels, clevs[i])
i = i + colorbar_label_skip
else:
colorbar_labels = clevs
# Make a list of strings of the colorbar levels for later labeling
clabels=[]
for i in colorbar_labels:
clabels.append(str(i))
# Add mult to colorbar_title if used
if (colorbar_title is None):
colorbar_title = ''
if (mult != 0):
colorbar_title = colorbar_title + ' *10$^{' + str(mult) + '}$'
# Catch null titles
if title is None:
title = ''
if plotvars.title is not None:
title = plotvars.title
# Set plot variables
title_fontsize = plotvars.title_fontsize
text_fontsize = plotvars.text_fontsize
colorbar_fontsize = plotvars.colorbar_fontsize
axis_label_fontsize = plotvars.axis_label_fontsize
continent_thickness = plotvars.continent_thickness
continent_color = plotvars.continent_color
continent_linestyle = plotvars.continent_linestyle
land_color = plotvars.land_color
ocean_color = plotvars.ocean_color
lake_color = plotvars.lake_color
text_fontweight = plotvars.text_fontweight
colorbar_fontweight = plotvars.colorbar_fontweight
title_fontweight = plotvars.title_fontweight
axis_label_fontweight = plotvars.axis_label_fontweight
if continent_thickness is None:
continent_thickness = 1.5
if continent_color is None:
continent_color = 'k'
if continent_linestyle is None:
continent_linestyle = 'solid'
cb_orient = colorbar_orientation
# Retrieve any user defined axis labels
if xlabel == '' and plotvars.xlabel is not None:
xlabel = plotvars.xlabel
if ylabel == '' and plotvars.ylabel is not None:
ylabel = plotvars.ylabel
if xticks is None and plotvars.xticks is not None:
xticks = plotvars.xticks
if plotvars.xticklabels is not None:
xticklabels = plotvars.xticklabels
else:
xticklabels = map(str, xticks)
if yticks is None and plotvars.yticks is not None:
yticks = plotvars.yticks
if plotvars.yticklabels is not None:
yticklabels = plotvars.yticklabels
else:
yticklabels = map(str, yticks)
##########
# Map plot
##########
if ptype == 1:
if verbose:
print 'con - making a map plot'
# Open a new plot if necessary
if plotvars.user_plot == 0:
gopen(user_plot=0)
# Set up mapping
lonrange = np.nanmax(x) - np.nanmin(x)
latrange = np.nanmax(y) - np.nanmin(y)
# Reset mapping
if plotvars.user_mapset == 0:
plotvars.lonmin = -180
plotvars.lonmax = 180
plotvars.latmin = -90
plotvars.latmax = 90
if (lonrange > 350 and latrange > 170) or plotvars.user_mapset == 1:
set_map()
else:
mapset(lonmin=np.nanmin(x), lonmax=np.nanmax(x),
latmin=np.nanmin(y), latmax=np.nanmax(y),
user_mapset=0, resolution=resolution_orig)
set_map()
mymap = plotvars.mymap
user_mapset = plotvars.user_mapset
lonrange = np.nanmax(x) - np.nanmin(x)
if lonrange > 350 and np.ndim(y) == 1:
# Add cyclic information if missing.
if lonrange < 360:
field, x = cartopy_util.add_cyclic_point(field, x)
lonrange = np.nanmax(x) - np.nanmin(x)
if x[-1] - x[0] == 360.0:
x[-1] = x[-1] + 0.001 # **cartopy line drawing fix
# Shift grid if needed
if plotvars.lonmin < np.nanmin(x):
x = x - 360
if plotvars.lonmin > np.nanmax(x):
x = x + 360
# Flip latitudes and field if latitudes are in descending order
if np.ndim(y) == 1:
if y[0] > y[-1]:
y = y[::-1]
field = np.flipud(field)
# Plotting a sub-area of the grid produces stray contour labels
# in polar plots. Subsample the grid to remove this problem
if plotvars.proj == 'npstere':
myypos = find_pos_in_array(vals=y, val=plotvars.boundinglat)
if myypos != -1:
y = y[myypos:]
field = field[myypos:, :]
if plotvars.proj == 'spstere':
myypos = find_pos_in_array(
vals=y, val=plotvars.boundinglat, above=1)
if myypos != -1:
y = y[0:myypos + 1]
field = field[0:myypos + 1, :]
# Create the meshgrid if required
if (np.ndim(field) == 1 and np.ndim(x) == 1 and np.ndim(y) == 1):
lons = x
lats = y
if (np.ndim(field) == 2 and np.ndim(x) == 2 and np.ndim(y) == 2):
lons = x
lats = y
if (np.ndim(field) == 2 and np.ndim(x) == 1 and np.ndim(y) == 1):
lons, lats = x, y
# Set the plot limits
if lonrange > 350:
gset(
xmin=plotvars.lonmin,
xmax=plotvars.lonmax,
ymin=plotvars.latmin,
ymax=plotvars.latmax,
user_gset=0)
else:
if user_mapset == 1:
gset(
xmin=plotvars.lonmin,
xmax=plotvars.lonmax,
ymin=plotvars.latmin,
ymax=plotvars.latmax,
user_gset=0)
else:
gset(
xmin=np.nanmin(lons),
xmax=np.nanmax(lons),
ymin=np.nanmin(lats),
ymax=np.nanmax(lats),
user_gset=0)
# Filled contours
if fill or blockfill:
if verbose:
print 'con - adding filled contours'
# Get colour scale for use in contouring
# If colour bar extensions are enabled then the colour map goes
# from 1 to ncols-2. The colours for the colour bar extensions
# are then changed on the colourbar and plot after the plot is made
cscale_ncols = np.size(plotvars.cs)
colmap = cscale_get_map()
cmap = matplotlib.colors.ListedColormap(colmap)
if (plotvars.levels_extend ==
'min' or plotvars.levels_extend == 'both'):
cmap.set_under(plotvars.cs[0])
if (plotvars.levels_extend ==
'max' or plotvars.levels_extend == 'both'):
cmap.set_over(plotvars.cs[-1])
# filled colour contours
cfill = mymap.contourf(lons, lats, field * fmult, clevs,
extend=plotvars.levels_extend,
cmap=cmap, norm=plotvars.norm,
alpha=alpha, transform=ccrs.PlateCarree(),
zorder=zorder)
# Block fill
if blockfill:
if verbose:
print 'con - adding blockfill'
if isinstance(f, cf.Field):
if f.ref('transverse_mercator'):
# Special case for transverse mercator
bfill(
f=f,
clevs=clevs,
alpha=alpha, zorder=zorder)
else:
if getattr(f.coord('lon'), 'hasbounds', False):
xpts = np.squeeze(f.coord('lon').bounds.array[:, 0])
ypts = np.squeeze(f.coord('lat').bounds.array[:, 0])
# Add last longitude point
xpts = np.append(xpts, f.coord('lon').bounds.array[-1, 1])
# Add last latitude point
ypts = np.append(ypts, f.coord('lat').bounds.array[-1, 1])
bfill(
f=field_orig *
fmult,
x=xpts,
y=ypts,
clevs=clevs,
lonlat=1,
bound=1,
alpha=alpha,
zorder=zorder)
else:
bfill(
f=field_orig *
fmult,
x=x_orig,
y=y_orig,
clevs=clevs,
lonlat=1,
bound=0,
alpha=alpha,
zorder=zorder)
else:
bfill(
f=field_orig *
fmult,
x=x_orig,
y=y_orig,
clevs=clevs,
lonlat=1,
bound=0,
alpha=alpha,
zorder=zorder)
# Contour lines and labels
if lines:
if verbose:
print 'con - adding contour lines and labels'
cs = mymap.contour(lons, lats, field * fmult, clevs, colors=colors,
linewidths=linewidths,
linestyles=linestyles, alpha=alpha,
transform=ccrs.PlateCarree(), zorder=zorder)
if line_labels:
nd = ndecs(clevs)
fmt = '%d'
if nd != 0:
fmt = '%1.' + str(nd) + 'f'
plotvars.plot.clabel(cs, fmt=fmt, colors=colors,
fontsize=text_fontsize,
fontweight=text_fontweight, zorder=1000)
# Thick zero contour line
if zero_thick:
cs = mymap.contour(lons, lats, field * fmult, [-1e-32, 0],
colors=colors, linewidths=zero_thick,
linestyles=linestyles, alpha=alpha,
transform=ccrs.PlateCarree(), zorder=zorder)
# axes
plot_map_axes(axes=axes, xaxis=xaxis, yaxis=yaxis,
xticks=xticks, xticklabels=xticklabels,
yticks=yticks, yticklabels=yticklabels,
user_xlabel=user_xlabel, user_ylabel=user_ylabel,
verbose=verbose)
# Coastlines
feature = cfeature.NaturalEarthFeature(
name='land', category='physical',
scale=plotvars.resolution,
facecolor='none')
mymap.add_feature(feature, edgecolor=continent_color,
linewidth=continent_thickness,
linestyle=continent_linestyle, zorder=zorder)
if ocean_color is not None:
mymap.add_feature(cfeature.OCEAN, edgecolor='face', facecolor=ocean_color,
zorder=999)
if land_color is not None:
mymap.add_feature(cfeature.LAND, edgecolor='face', facecolor=land_color,
zorder=999)
if lake_color is not None:
mymap.add_feature(cfeature.LAKES, edgecolor='face', facecolor=lake_color,
zorder=999)
# Title
if title != '':
map_title(title)
# Color bar
if colorbar:
cfplot_colorbar(cfill=cfill, colorbar_labels=colorbar_labels,
colorbar_orientation=colorbar_orientation,
colorbar_position=colorbar_position,
colorbar_shrink=colorbar_shrink,
cb_orient=cb_orient,colorbar_title=colorbar_title,
text_fontsize=text_fontsize,
text_fontweight=text_fontweight,
colorbar_text_up_down=colorbar_text_up_down,
colorbar_text_down_up=colorbar_text_down_up,
colorbar_drawedges=colorbar_drawedges,
verbose=verbose)
# Reset plot limits if not a user plot
if plotvars.user_gset == 0:
gset()
########################################
# Latitude, longitude or time vs Z plots
########################################
if ptype == 2 or ptype == 3 or ptype == 7:
if verbose:
if ptype == 2:
print 'con - making a latitude-pressure plot'
if ptype == 3:
print 'con - making a longitude-pressure plot'
if ptype == 7:
print 'con - making a time-pressure plot'
# Work out which way is up
positive = None
if isinstance(f, cf.Field):
if hasattr(f.item('Z'), 'positive'):
positive = f.item('Z').positive
else:
errstr = '\ncf-plot - data error \n'
errstr += 'data needs a vertical coordinate direction'
errstr += ' as required in CF data conventions'
errstr += '\nMaking a contour plot assuming positive is up\n\n'
errstr += 'If this is incorrect the data needs to be modified to \n'
errstr += 'include a correct value for the direction attribute\n\n'
print errstr
positive='up'
else:
positive = 'down'
if 'theta' in ylabel.split(' '):
positive = 'up'
if 'height' in ylabel.split(' '):
positive = 'up'
if plotvars.user_plot == 0:
gopen(user_plot=0)
# Use gset parameter of ylog if user has set this
if plotvars.ylog is True or plotvars.ylog == 1:
ylog = True
# Set plot limits
user_gset = plotvars.user_gset
if user_gset == 0:
# Program selected data plot limits
xmin = np.nanmin(x)
if xmin < -80 and xmin >= -90:
xmin = -90
xmax = np.nanmax(x)
if xmax > 80 and xmax <= 90:
xmax = 90
if positive == 'down':
ymin = np.nanmax(y)
ymax = np.nanmin(y)
if ymax < 10:
ymax = 0
else:
ymin = np.nanmin(y)
ymax = np.nanmax(y)
else:
# Use user specified plot limits
xmin = plotvars.xmin
xmax = plotvars.xmax
ymin = plotvars.ymin
ymax = plotvars.ymax
xstep = 10
if (xmin == -90 and xmax == 90):
xstep = 30
ystep = 100
myrange = abs(ymax - ymin)
if myrange < 1:
ystep = abs(ymax - ymin)/10.
if abs(ymax - ymin) > 1:
ystep = 1
if abs(ymax - ymin) > 10:
ystep = 10
if abs(ymax - ymin) > 100:
ystep = 100
if abs(ymax - ymin) > 1000:
ystep = 200
if abs(ymax - ymin) > 2000:
ystep = 500
if abs(ymax - ymin) > 5000:
ystep = 1000
if abs(ymax - ymin) > 15000:
ystep = 5000
# Work out ticks and tick labels
if ylog is False or ylog == 0:
heightticks = gvals(dmin=min(ymin, ymax),
dmax=max(ymin, ymax),
tight=1, mystep=ystep, mod=0)[0]
if myrange < 1 and myrange > 0.1:
heightticks=np.arange(10)/10.0
else:
heightticks = []
for tick in 1000, 100, 10, 1:
if tick >= min(ymin, ymax) and tick <= max(ymin, ymax):
heightticks.append(tick)
heightlabels = heightticks
if axes:
if xaxis:
if xticks is not None:
if xticklabels is None:
xticklabels = xticks
else:
xticks = [100000000]
xticklabels = xticks
xlabel = ''
if yaxis:
if yticks is not None:
heightticks = yticks
heightlabels = yticks
if yticklabels is not None:
heightlabels = yticklabels
else:
heightticks = [100000000]
ylabel = ''
else:
xticks = [100000000]
xticklabels = xticks
heightticks = [100000000]
heightlabels = heightticks
xlabel = ''
ylabel = ''
if yticks is None:
yticks = heightticks
yticklabels = heightlabels
if ptype == 7:
if isinstance(f, cf.Field):
if plotvars.user_gset == 0:
tmin = f.item('T').dtarray[0]
tmax = f.item('T').dtarray[-1]
else:
# Use user set values if present
tmin = plotvars.xmin
tmax = plotvars.xmax
ref_time = f.item('T').units
ref_calendar = f.item('T').calendar
ref_time_origin = str(f.item('T').Units.reftime)
time_units = cf.Units(ref_time, ref_calendar)
t = cf.Data(cf.dt(tmin), units=time_units)
xmin = t.array
t = cf.Data(cf.dt(tmax), units=time_units)
xmax = t.array
if xticks is None and xaxis:
if ptype == 2:
xticks, xticklabels = mapaxis(
min=xmin, max=xmax, type=2) # lat-pressure
if ptype == 3:
xticks, xticklabels = mapaxis(
min=xmin, max=xmax, type=1) # lon-pressure
if ptype == 7:
# time-pressure
if isinstance(f, cf.Field):
# Change plotvars.xmin and plotvars.xmax from a date string
# to a number
ref_time = f.item('T').units
ref_calendar = f.item('T').calendar
ref_time_origin = str(f.item('T').Units.reftime)
time_units = cf.Units(ref_time, ref_calendar)
t = cf.Data(cf.dt(tmin), units=time_units)
xmin = t.array
t = cf.Data(cf.dt(tmax), units=time_units)
xmax = t.array
taxis = cf.Data(
[cf.dt(tmin), cf.dt(tmax)], units=time_units)
time_ticks, time_labels, tlabel = timeaxis(taxis)
# Use user supplied labels if present
if user_xlabel is None:
xlabel = tlabel
if xticks is None:
xticks = time_ticks
if xticklabels is None:
xticklabels = time_labels
else:
errstr = '\nNot a CF field\nPlease use ptype=0 and '
errstr = srrstr + 'specify axis labels manually\n'
raise Warning(errstr)
# Set plot limits
if ylog is False or ylog == 0:
gset(
xmin=xmin,
xmax=xmax,
ymin=ymin,
ymax=ymax,
user_gset=user_gset)
else:
if ymax == 0:
ymax = 1 # Avoid zero in a log plot
gset(
xmin=xmin,
xmax=xmax,
ymin=ymin,
ymax=ymax,
ylog=1,
user_gset=user_gset)
# Label axes
axes_plot(xticks=xticks, xticklabels=xticklabels, yticks=heightticks,
yticklabels=heightlabels, xlabel=xlabel, ylabel=ylabel)
# Get colour scale for use in contouring
# If colour bar extensions are enabled then the colour map goes
# from 1 to ncols-2. The colours for the colour bar extensions are
# then changed on the colourbar and plot after the plot is made
cscale_ncols = np.size(plotvars.cs)
colmap = cscale_get_map()
# Filled contours
if fill or blockfill:
colmap = cscale_get_map()
cmap = matplotlib.colors.ListedColormap(colmap)
if (plotvars.levels_extend ==
'min' or plotvars.levels_extend == 'both'):
cmap.set_under(plotvars.cs[0])
if (plotvars.levels_extend ==
'max' or plotvars.levels_extend == 'both'):
cmap.set_over(plotvars.cs[-1])
cfill = plotvars.plot.contourf(x, y, field * fmult, clevs,
extend=plotvars.levels_extend,
cmap=cmap,
norm=plotvars.norm, alpha=alpha,
zorder=zorder)
# Block fill
if blockfill:
if isinstance(f, cf.Field):
hasbounds = True
if ptype == 2:
if getattr(f.coord('Y'), 'hasbounds', False):
hasbounds = False
xpts = np.squeeze(f.coord('Y').bounds.array)[:, 0]
xpts = np.append(xpts, f.coord('Y').bounds.array[-1, 1])
if ptype == 3:
if getattr(f.coord('X'), 'hasbounds', False):
hasbounds = False
xpts = np.squeeze(f.coord('X').bounds.array)[:, 0]
xpts = np.append(xpts, f.coord('X').bounds.array[-1, 1])
if ptype == 7:
if getattr(f.coord('T'), 'hasbounds', False):
hasbounds = False
xpts = np.squeeze(f.coord('T').bounds.array)[:, 0]
xpts = np.append(xpts, f.coord('T').bounds.array[-1, 1])
if getattr(f.coord('Z'), 'hasbounds', False):
hasbounds = False
ypts = np.squeeze(f.coord('Z').bounds.array)[:, 0]
ypts = np.append(ypts, f.coord('Z').bounds.array[-1, 1])
if hasbounds is False:
bfill(
f=field_orig *
fmult,
x=xpts,
y=ypts,
clevs=clevs,
lonlat=0,
bound=1,
alpha=alpha,
zorder=zorder)
else:
bfill(
f=field_orig *
fmult,
x=x_orig,
y=y_orig,
clevs=clevs,
lonlat=0,
bound=0,
alpha=alpha,
zorder=zorder)
else:
bfill(
f=field_orig *
fmult,
x=x_orig,
y=y_orig,
clevs=clevs,
lonlat=0,
bound=0,
alpha=alpha,
zorder=zorder)
# Contour lines and labels
if lines:
cs = plotvars.plot.contour(
x, y, field * fmult, clevs, colors=colors,
linewidths=linewidths, linestyles=linestyles, zorder=zorder)
if line_labels:
nd = ndecs(clevs)
fmt = '%d'
if nd != 0:
fmt = '%1.' + str(nd) + 'f'
plotvars.plot.clabel(
cs,
fmt=fmt,
colors=colors,
fontsize=text_fontsize,
fontweight=text_fontweight)
# Thick zero contour line
if zero_thick:
cs = plotvars.plot.contour(x, y, field * fmult,
[-1e-32, 0], colors=colors,
linewidths=zero_thick,
linestyles=linestyles, alpha=alpha,
zorder=zorder)
# Color bar
if colorbar:
cfplot_colorbar(cfill=cfill, colorbar_labels=colorbar_labels,
colorbar_orientation=colorbar_orientation,
colorbar_position=colorbar_position,
colorbar_shrink=colorbar_shrink,
cb_orient=cb_orient,colorbar_title=colorbar_title,
text_fontsize=text_fontsize,
text_fontweight=text_fontweight,
colorbar_text_up_down=colorbar_text_up_down,
colorbar_text_down_up=colorbar_text_down_up,
colorbar_drawedges=colorbar_drawedges,
verbose=verbose)
# Title
plotvars.plot.set_title(title, y=1.03, fontsize=title_fontsize,
fontweight=title_fontweight)
# Reset plot limits to those supplied by the user
if user_gset == 1 and ptype == 7:
gset(xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax,
user_gset=user_gset)
# reset plot limits if not a user plot
if plotvars.user_gset == 0:
gset()
#################
# Hovmuller plots
#################
if (ptype == 4 or ptype == 5):
if verbose:
print 'con - making a Hovmuller plot'
yplotlabel = 'Time'
if ptype == 4:
xplotlabel = 'Longitude'
if ptype == 5:
xplotlabel = 'Latitude'
user_gset = plotvars.user_gset
# Time strings set to None initially
tmin = None
tmax = None
# Set plot limits
if all(val is not None for val in [
plotvars.xmin, plotvars.xmax, plotvars.ymin, plotvars.ymax]):
# Store time strings for later use
tmin = plotvars.ymin
tmax = plotvars.ymax
# Check data has CF attributes needed
check_units = check_units = True
check_calendar = True
check_Units_reftime = True
if hasattr(f.item('T'), 'units') is False:
check_units = False
if hasattr(f.item('T'), 'calendar') is False:
check_calendar = False
if hasattr(f.item('T'), 'Units'):
if not hasattr(f.item('T').Units, 'reftime'):
check_Units_reftime = False
else:
check_Units_reftime = False
if False in [check_units, check_calendar, check_Units_reftime]:
print '\nThe required CF time information to make the plot'
print 'is not available please fix the following before'
print 'trying to plot again'
if check_units is False:
print 'Time axis missing: units'
if check_calendar is False:
print 'Time axis missing: calendar'
if check_Units_reftime is False:
print 'Time axis missing: Units.reftime'
return
# Change from date string in ymin and ymax to date as a float
ref_time = f.item('T').units
ref_calendar = f.item('T').calendar
ref_time_origin = str(f.item('T').Units.reftime)
time_units = cf.Units(ref_time, ref_calendar)
t = cf.Data(cf.dt(plotvars.ymin), units=time_units)
ymin = t.array
t = cf.Data(cf.dt(plotvars.ymax), units=time_units)
ymax = t.array
xmin = plotvars.xmin
xmax = plotvars.xmax
else:
xmin = np.nanmin(x)
xmax = np.nanmax(x)
ymin = np.nanmin(y)
ymax = np.nanmax(y)
# Extract axis labels
time_ticks, time_labels, ylabel = timeaxis(f.item('T'))
if ptype == 4:
lonlatticks, lonlatlabels = mapaxis(min=xmin, max=xmax, type=1)
if ptype == 5:
lonlatticks, lonlatlabels = mapaxis(min=xmin, max=xmax, type=2)
if axes:
if xaxis:
if xticks is not None:
lonlatticks = xticks
lonlatlabels = xticks
if xticklabels is not None:
lonlatlabels = xticklabels
else:
lonlatticks = [100000000]
xlabel = ''
if yaxis:
if yticks is not None:
timeticks = yticks
timelabels = yticks
if yticklabels is not None:
timelabels = yticklabels
else:
timeticks = [100000000]
ylabel = ''
else:
latticks = [100000000]
timeticks = [100000000]
xplotlabel = ''
yplotlabel = ''
if user_xlabel is not None:
xplotlabel = user_xlabel
if user_ylabel is not None:
yplotlabel = user_ylabel
# Use the automatically generated labels if none are supplied
if ylabel is None:
yplotlabel = time_axis_label
if np.size(time_ticks) > 0:
timeticks = time_ticks
if np.size(time_labels) > 0:
timelabels = time_labels
# Swap axes if requested
if swap_axes:
x, y = y, x
field = np.flipud(np.rot90(field))
xmin, ymin = ymin, xmin
xmax, ymax = ymax, xmax
xplotlabel, yplotlabel = yplotlabel, xplotlabel
lonlatticks, timeticks = timeticks, lonlatticks
lonlatlabels, timelabels = timelabels, lonlatlabels
# Set plot limits
if plotvars.user_plot == 0:
gopen(user_plot=0)
gset(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, user_gset=user_gset)
# Revert to time strings if set
if all(val is not None for val in [tmin, tmax]):
plotvars.ymin = tmin
plotvars.ymax = tmax
# Set and label axes
axes_plot(xticks=lonlatticks, xticklabels=lonlatlabels,
yticks=timeticks, yticklabels=timelabels,
xlabel=xplotlabel, ylabel=yplotlabel)
# Get colour scale for use in contouring
# If colour bar extensions are enabled then the colour map goes
# from 1 to ncols-2. The colours for the colour bar extensions are
# then changed on the colourbar and plot after the plot is made
cscale_ncols = np.size(plotvars.cs)
colmap = cscale_get_map()
# Filled contours
if fill or blockfill:
colmap = cscale_get_map()
cmap = matplotlib.colors.ListedColormap(colmap)
if (plotvars.levels_extend ==
'min' or plotvars.levels_extend == 'both'):
cmap.set_under(plotvars.cs[0])
if (plotvars.levels_extend ==
'max' or plotvars.levels_extend == 'both'):
cmap.set_over(plotvars.cs[-1])
cfill = plotvars.plot.contourf(x, y, field * fmult, clevs,
extend=plotvars.levels_extend,
cmap=cmap,
norm=plotvars.norm, alpha=alpha,
zorder=zorder)
# Block fill
if blockfill:
if isinstance(f, cf.Field):
if f.coord('lon').hasbounds:
if ptype == 4:
xpts = np.squeeze(f.coord('X').bounds.array)[:, 0]
xpts = np.append(xpts, f.coord('X').bounds.array[-1, 1])
if ptype == 5:
xpts = np.squeeze(f.coord('Y').bounds.array)[:, 0]
xpts = np.append(xpts, f.coord('Y').bounds.array[-1, 1])
ypts = np.squeeze(f.coord('T').bounds.array)[:, 0]
ypts = np.append(ypts, f.coord('T').bounds.array[-1, 1])
if swap_axes:
xpts, ypts = ypts, xpts
field_orig = np.flipud(np.rot90(field_orig))
bfill(
f=field_orig *
fmult,
x=xpts,
y=ypts,
clevs=clevs,
lonlat=0,
bound=1,
alpha=alpha,
zorder=zorder)
else:
if swap_axes:
x_orig, y_orig = y_orig, x_orig
field_orig = np.flipud(np.rot90(field_orig))
bfill(
f=field_orig *
fmult,
x=x_orig,
y=y_orig,
clevs=clevs,
lonlat=0,
bound=0,
alpha=alpha,
zorder=zorder)
else:
if swap_axes:
x_orig, y_orig = y_orig, x_orig
field_orig = np.flipud(np.rot90(field_orig))
bfill(
f=field_orig *
fmult,
x=x_orig,
y=y_orig,
clevs=clevs,
lonlat=0,
bound=0,
alpha=alpha,
zorder=zorder)
# Contour lines and labels
if lines:
cs = plotvars.plot.contour(
x, y, field * fmult, clevs, colors=colors,
linewidths=linewidths, linestyles=linestyles, alpha=alpha)
if line_labels:
nd = ndecs(clevs)
fmt = '%d'
if nd != 0:
fmt = '%1.' + str(nd) + 'f'
plotvars.plot.clabel(cs, fmt=fmt, colors=colors,
fontsize=text_fontsize,
fontweight=text_fontweight,
zorder=zorder)
# Thick zero contour line
if zero_thick:
cs = plotvars.plot.contour(x, y, field * fmult,
[-1e-32, 0], colors=colors,
linewidths=zero_thick,
linestyles=linestyles, alpha=alpha,
zorder=zorder)
# Color bar
if colorbar:
cfplot_colorbar(cfill=cfill, colorbar_labels=colorbar_labels,
colorbar_orientation=colorbar_orientation,
colorbar_position=colorbar_position,
colorbar_shrink=colorbar_shrink,
cb_orient=cb_orient,colorbar_title=colorbar_title,
text_fontsize=text_fontsize,
text_fontweight=text_fontweight,
colorbar_text_up_down=colorbar_text_up_down,
colorbar_text_down_up=colorbar_text_down_up,
colorbar_drawedges=colorbar_drawedges,
verbose=verbose)
# Title
plotvars.plot.set_title(
title,
y=1.03,
fontsize=title_fontsize,
fontweight=title_fontweight)
# reset plot limits if not a user plot
if user_gset == 0:
gset()
#############
# Rotated pole
#############
if ptype == 6:
# Extract x and y grid points
if plotvars.proj == 'cyl':
xpts = x
ypts = y
#if f.item('longitude'):
# xpts=f.item('longitude').array
#if f.item('latitude'):
# ypts=f.item('latitude').array
else:
xpts = np.arange(np.size(x))
ypts = np.arange(np.size(y))
if verbose:
print 'con - making a rotated pole plot'
user_gset = plotvars.user_gset
if plotvars.user_plot == 0:
gopen(user_plot=0)
# Set plot limits
if plotvars.proj == 'rotated':
plotargs={}
gset(xmin=0, xmax=np.size(xpts) - 1,
ymin=0, ymax=np.size(ypts) - 1,
user_gset=user_gset)
plot=plotvars.plot
if plotvars.proj == 'cyl':
rotated_pole = f.ref('rotated_latitude_longitude')
xpole = rotated_pole['grid_north_pole_longitude']
ypole = rotated_pole['grid_north_pole_latitude']
transform = ccrs.RotatedPole(pole_latitude=ypole,
pole_longitude=xpole)
plotargs={'transform' : transform}
if plotvars.user_mapset == 1:
set_map()
else:
# Extract lonmin, lonmax, latmin, latmax from rotated grid coordinates
#xvals, yvals = np.meshgrid(xpts, ypts)
points=ccrs.PlateCarree().transform_points(transform, xpts.flatten(),
ypts.flatten())
lons=np.array(points)[:, 0]
lats=np.array(points)[:, 1]
mapset(lonmin=np.min(lons), lonmax=np.max(lons),
latmin=np.min(lats), latmax=np.max(lats),
user_mapset=0, resolution=resolution_orig)
set_map()
plot=plotvars.mymap
# Get colour scale for use in contouring
# If colour bar extensions are enabled then the colour map goes
# from 1 to ncols-2. The colours for the colour bar extensions are
# then changed on the colourbar and plot after the plot is made
cscale_ncols = np.size(plotvars.cs)
colmap = cscale_get_map()
# Filled contours
if fill or blockfill:
colmap = cscale_get_map()
cmap = matplotlib.colors.ListedColormap(colmap)
if (plotvars.levels_extend ==
'min' or plotvars.levels_extend == 'both'):
cmap.set_under(plotvars.cs[0])
if (plotvars.levels_extend ==
'max' or plotvars.levels_extend == 'both'):
cmap.set_over(plotvars.cs[-1])
cfill = plot.contourf(xpts, ypts, field * fmult, clevs,
extend=plotvars.levels_extend,
cmap=cmap,
norm=plotvars.norm, alpha=alpha,
zorder=zorder, **plotargs)
# Block fill
if blockfill:
bfill(
f=field_orig *
fmult,
x=xpts,
y=ypts,
clevs=clevs,
lonlat=0,
bound=0,
alpha=alpha,
zorder=zorder)
# Contour lines and labels
if lines:
cs = plot.contour(
xpts, ypts, field * fmult, clevs, colors=colors,
linewidths=linewidths, linestyles=linestyles,
zorder=zorder, **plotargs)
if line_labels:
nd = ndecs(clevs)
fmt = '%d'
if nd != 0:
fmt = '%1.' + str(nd) + 'f'
plot.clabel(cs, fmt=fmt, colors=colors,
fontsize=text_fontsize,
fontweight=text_fontweight,
zorder=zorder)
# Thick zero contour line
if zero_thick:
cs = plot.contour(xpts, ypts, field * fmult,
[-1e-32, 0], colors=colors,
linewidths=zero_thick,
linestyles=linestyles, alpha=alpha,
zorder=zorder, **plotargs)
# Color bar
if colorbar:
cfplot_colorbar(cfill=cfill, colorbar_labels=colorbar_labels,
colorbar_orientation=colorbar_orientation,
colorbar_position=colorbar_position,
colorbar_shrink=colorbar_shrink,
cb_orient=cb_orient,colorbar_title=colorbar_title,
text_fontsize=text_fontsize,
text_fontweight=text_fontweight,
colorbar_text_up_down=colorbar_text_up_down,
colorbar_text_down_up=colorbar_text_down_up,
colorbar_drawedges=colorbar_drawedges,
verbose=verbose)
# Rotated grid axes
if axes:
if plotvars.proj == 'cyl':
plot_map_axes(axes=axes, xaxis=xaxis, yaxis=yaxis,
xticks=xticks, xticklabels=xticklabels,
yticks=yticks, yticklabels=yticklabels,
user_xlabel=user_xlabel, user_ylabel=user_ylabel,
verbose=verbose)
else:
rgaxes(xpole=xpole, ypole=ypole, xvec=x, yvec=y,
xticks=xticks, xticklabels=xticklabels,
yticks=yticks, yticklabels=yticklabels,
axes=axes, xaxis=xaxis, yaxis=yaxis,
xlabel=xlabel, ylabel=ylabel)
if plotvars.proj == 'rotated':
# Remove Matplotlib default axis labels
axes_plot(xticks=[100000000], xticklabels=[''],
yticks=[100000000], yticklabels=[''],
xlabel='', ylabel='')
# Add title and coastlines for cylindrical projection
if plotvars.proj == 'cyl':
# Coastlines
feature = cfeature.NaturalEarthFeature(
name='land', category='physical',
scale=plotvars.resolution,
facecolor='none')
plotvars.mymap.add_feature(feature, edgecolor=continent_color,
linewidth=continent_thickness,
linestyle=continent_linestyle, zorder=zorder)
# Title
if title != '':
map_title(title)
# Add title for native grid
if plotvars.proj == 'rotated':
# Title
plotvars.plot.set_title(title, y=1.03,
fontsize=title_fontsize,
fontweight=title_fontweight)
# reset plot limits if not a user plot
if plotvars.user_gset == 0:
gset()
############
# Other plots
############
if ptype == 0:
if verbose:
print 'con - making an other plot'
if plotvars.user_plot == 0:
gopen(user_plot=0)
user_gset = plotvars.user_gset
# Work out axes if none are supplied
if any(val is None for val in [
plotvars.xmin, plotvars.xmax, plotvars.ymin, plotvars.ymax]):
xmin = np.nanmin(x)
xmax = np.nanmax(x)
ymin = np.nanmin(y)
ymax = np.nanmax(y)
else:
xmin = plotvars.xmin
xmax = plotvars.xmax
ymin = plotvars.ymin
ymax = plotvars.ymax
xstep = (xmax - xmin) / 10.0
ystep = (ymax - ymin) / 10.0
# Set plot limits and set default plot labels
gset(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, user_gset=user_gset)
xaxisticks = gvals(
dmin=xmin,
dmax=xmax,
mystep=(
xmax -
xmin) /
10.0,
tight=1,
mod=0)[0]
xaxislabels = xaxisticks
if ymin < ymax:
yaxisticks = gvals(
dmin=ymin, dmax=ymax, mystep=(
ymax - ymin) / 10.0, tight=1, mod=0)[0]
if ymax < ymin:
yaxisticks = gvals(
dmin=ymax, dmax=ymin, mystep=(
ymin - ymax) / 10.0, tight=1, mod=0)[0]
yaxislabels = yaxisticks
if user_xlabel is not None:
xplotlabel = user_xlabel
else:
xplotlabel = ''
if user_ylabel is not None:
yplotlabel = user_ylabel
else:
yplotlabel = ''
# Draw axes
if axes:
if xaxis:
if xticks is not None:
xaxisticks = xticks
xaxislabels = xticks
if xticklabels is not None:
xaxislabels = xticklabels
else:
xaxisticks = [100000000]
xlabel = ''
if yaxis:
if yticks is not None:
yaxisticks = yticks
yaxislabels = yticks
if yticklabels is not None:
yaxislabels = yticklabels
else:
yaxisticks = [100000000]
ylabel = ''
else:
xplotticks = [100000000]
xplotticks = [100000000]
xlabel = ''
ylabel = ''
axes_plot(xticks=xaxisticks, xticklabels=xaxislabels,
yticks=yaxisticks, yticklabels=yaxislabels,
xlabel=xplotlabel, ylabel=yplotlabel)
# Get colour scale for use in contouring
# If colour bar extensions are enabled then the colour map goes
# then from 1 to ncols-2. The colours for the colour bar extensions
# are changed on the colourbar and plot after the plot is made
cscale_ncols = np.size(plotvars.cs)
colmap = cscale_get_map()
# Filled contours
if fill or blockfill:
colmap = cscale_get_map()
cmap = matplotlib.colors.ListedColormap(colmap)
if (plotvars.levels_extend ==
'min' or plotvars.levels_extend == 'both'):
cmap.set_under(plotvars.cs[0])
if (plotvars.levels_extend ==
'max' or plotvars.levels_extend == 'both'):
cmap.set_over(plotvars.cs[-1])
cfill = plotvars.plot.contourf(x, y, field * fmult, clevs,
extend=plotvars.levels_extend,
cmap=cmap,
norm=plotvars.norm, alpha=alpha,
zorder=zorder)
# Block fill
if blockfill:
bfill(
f=field_orig *
fmult,
x=x_orig,
y=y_orig,
clevs=clevs,
lonlat=0,
bound=0,
alpha=alpha,
zorder=zorder)
# Contour lines and labels
if lines:
cs = plotvars.plot.contour(
x, y, field * fmult, clevs, colors=colors,
linewidths=linewidths, linestyles=linestyles,
zorder=zorder)
if line_labels:
nd = ndecs(clevs)
fmt = '%d'
if nd != 0:
fmt = '%1.' + str(nd) + 'f'
plotvars.plot.clabel(cs, fmt=fmt, colors=colors,
fontsize=text_fontsize,
fontweight=text_fontweight,
zorder=zorder)
# Thick zero contour line
if zero_thick:
cs = plotvars.plot.contour(x, y, field * fmult, [-1e-32, 0],
colors=colors,
linewidths=zero_thick,
linestyles=linestyles, alpha=alpha,
zorder=zorder)
# Color bar
if colorbar:
cfplot_colorbar(cfill=cfill, colorbar_labels=colorbar_labels,
colorbar_orientation=colorbar_orientation,
colorbar_position=colorbar_position,
colorbar_shrink=colorbar_shrink,
cb_orient=cb_orient,colorbar_title=colorbar_title,
text_fontsize=text_fontsize,
text_fontweight=text_fontweight,
colorbar_text_up_down=colorbar_text_up_down,
colorbar_text_down_up=colorbar_text_down_up,
colorbar_drawedges=colorbar_drawedges,
verbose=verbose)
# Title
plotvars.plot.set_title(
title,
y=1.03,
fontsize=title_fontsize,
fontweight=title_fontweight)
# reset plot limits if not a user plot
if plotvars.user_gset == 0:
gset()
############################
# Set axis width if required
############################
if plotvars.axis_width is not None:
for axis in ['top','bottom','left','right']:
plotvars.plot.spines[axis].set_linewidth(plotvars.axis_width)
################################
# Add a master title if reqested
################################
if plotvars.master_title is not None:
location=plotvars.master_title_location
plotvars.master_plot.text(location[0], location[1],
plotvars.master_title,
horizontalalignment='center',
fontweight=plotvars.master_title_fontweight,
fontsize=plotvars.master_title_fontsize)
# Reset map resolution
if plotvars.user_mapset == 0:
mapset()
mapset(resolution=resolution_orig)
##################
# Save or view plot
##################
if plotvars.user_plot == 0:
if verbose:
print 'con - saving or viewing plot'
np.seterr(**old_settings) # reset to default numpy error settings
gclose()
[docs]def mapset(lonmin=None, lonmax=None, latmin=None, latmax=None, proj='cyl',
boundinglat=0, lon_0=0, lat_0=40, resolution='110m', user_mapset=1):
"""
| mapset sets the mapping parameters.
|
| lonmin=lonmin - minimum longitude
| lonmax=lonmax - maximum longitude
| latmin=latmin - minimum latitude
| latmax=latmax - maximum latitude
| proj=proj - 'cyl' for cylindrical projection. 'npstere' or 'spstere'
| for northern hemisphere or southern hemisphere polar stereographic.
| ortho, merc, moll, robin and lcc are abreviations for orthographic,
| mercator, mollweide, robinson and lambert conformal projections
| 'rotated' for contour plots on the native rotated grid.
|
| boundinglat=boundinglat - edge of the viewable latitudes in a
| stereographic plot
| lon_0=0 - longitude centre of desired map domain in polar
| stereographic and orthogrphic plots
| lat_0=40 - latitude centre of desired map domain in orthogrphic plots
| resolution='110m' - the map resolution - can be one of '110m',
| '50m' or '10m'. '50m' means 1:50,000,000 and not 50 metre.
| user_mapset=user_mapset - variable to indicate whether a user call
| to mapset has been made.
|
| The default map plotting projection is the cyclindrical equidistant
| projection from -180 to 180 in longitude and -90 to 90 in latitude.
| To change the map view in this projection to over the United Kingdom,
| for example, you would use
| mapset(lonmin=-6, lonmax=3, latmin=50, latmax=60)
| or
| mapset(-6, 3, 50, 60)
|
| The limits are -360 to 720 in longitude so to look at the equatorial
| Pacific you could use
| mapset(lonmin=90, lonmax=300, latmin=-30, latmax=30)
| or
| mapset(lonmin=-270, lonmax=-60, latmin=-30, latmax=30)
|
| The proj parameter accepts 'npstere' and 'spstere' for northern
| hemisphere or southern hemisphere polar stereographic projections.
| In addition to these the boundinglat parameter sets the edge of the
| viewable latitudes and lon_0 sets the centre of desired map domain.
|
|
|
| Map settings are persistent until a new call to mapset is made. To
| reset to the default map settings use mapset().
:Returns:
None
"""
# Set the continent resolution
plotvars.resolution = resolution
if all(val is None for val in [
lonmin, lonmax, latmin, latmax]) and proj == 'cyl':
plotvars.lonmin = -180
plotvars.lonmax = 180
plotvars.latmin = -90
plotvars.latmax = 90
plotvars.proj = 'cyl'
plotvars.user_mapset = 0
plotvars.plot_xmin=None
plotvars.plot_xmax=None
plotvars.plot_ymin=None
plotvars.plot_ymax=None
return
if lonmin is None:
lonmin = -180
if lonmax is None:
lonmax = 180
if latmin is None:
latmin = -90
if proj == 'merc':
latmin = -80
if latmax is None:
latmax = 90
if proj == 'merc':
latmax = 80
if proj == 'moll':
lonmin = lon_0 - 180
lonmax = lon_0 + 180
plotvars.lonmin = lonmin
plotvars.lonmax = lonmax
plotvars.latmin = latmin
plotvars.latmax = latmax
plotvars.proj = proj
plotvars.boundinglat = boundinglat
plotvars.lon_0 = lon_0
plotvars.lat_0 = lat_0
plotvars.user_mapset = user_mapset
#set_map()
[docs]def levs(min=None, max=None, step=None, manual=None, extend='both'):
"""
| The levs command manually sets the contour levels.
| min=min - minimum level
| max=max - maximum level
| step=step - step between levels
| manual= manual - set levels manually
| extend='neither', 'both', 'min', or 'max' - colour bar limit extensions
| Use the levs command when a predefined set of levels is required. The
| min, max and step parameters can be used to define a set of levels.
| These can take integer or floating point numbers. If just the step is
| defined then cf-plot will internally try to define a reasonable set
| of levels.
| If colour filled contours are plotted then the default is to extend
| the minimum and maximum contours coloured for out of range values
| - extend='both'.
| Once a user call is made to levs the levels are persistent.
| i.e. the next plot will use the same set of levels.
| Use levs() to reset to undefined levels.
:Returns:
None
"""
if all(val is None for val in [min, max, step, manual]):
plotvars.levels = None
plotvars.levels_min = None
plotvars.levels_max = None
plotvars.levels_step = None
plotvars.levels_extend = 'both'
plotvars.norm = None
plotvars.user_levs = 0
return
if manual is not None:
plotvars.levels = manual
plotvars.levels_min = None
plotvars.levels_max = None
plotvars.levels_step = None
# Set the normalization object as we are using potentially unevenly
# spaced levels
ncolors = np.size(plotvars.levels)
if extend == 'both' or extend == 'max':
ncolors = ncolors - 1
plotvars.norm = matplotlib.colors.BoundaryNorm(
boundaries=plotvars.levels, ncolors=ncolors)
plotvars.user_levs = 1
else:
if all(val is not None for val in [min, max, step]):
plotvars.levels_min = min
plotvars.levels_max = max
plotvars.levels_step = step
plotvars.norm = None
if all(isinstance(item, (int, long)) for item in [min, max, step]):
lstep = step * 1e-10
levs = (np.arange(min, max + lstep, step, dtype=np.float64))
levs = ((levs * 1e10).astype(np.int64)).astype(np.float64)
levs = (levs / 1e10).astype(np.int)
plotvars.levels = levs
else:
lstep = step * 1e-10
levs = np.arange(min, max + lstep, step, dtype=np.float64)
levs = (levs * 1e10).astype(np.int64).astype(np.float64)
levs = levs / 1e10
plotvars.levels = levs
plotvars.user_levs = 1
# Check for spurious decimal places due to numeric representation
# and fix if found
for pt in np.arange(np.size(plotvars.levels)):
ndecs = str(plotvars.levels[pt])[::-1].find('.')
if ndecs > 7:
plotvars.levels[pt] = round(plotvars.levels[pt], 7)
# If step only is set then reset user_levs to zero
if step is not None and all(val is None for val in [min, max]):
plotvars.user_levs = 0
plotvars.levels = None
plotvars.levels_step = step
# Check extend has a proper value
if extend not in ['neither', 'min', 'max', 'both']:
errstr = "\n\n extend must be one of 'neither', 'min', 'max', 'both'\n"
raise TypeError(errstr)
plotvars.levels_extend = extend
[docs]def mapaxis(min=None, max=None, type=None):
"""
| mapaxis is used to work out a sensible set of longitude and latitude
| tick marks and labels. This is an internal routine and is not used
| by the user.
| min=None - minimum axis value
| max=None - maximum axis value
| type=None - 1 = longitude, 2 = latitude
:Returns:
longtitude/latitude ticks and longitude/latitude tick labels
|
|
|
|
|
|
|
"""
import numpy as np
degsym=''
if plotvars.degsym:
degsym='$\degree$'
if type == 1:
lonmin = min
lonmax = max
lonrange = lonmax - lonmin
lonstep = 60
if lonrange <= 180:
lonstep = 30
if lonrange <= 90:
lonstep = 10
if lonrange <= 30:
lonstep = 5
if lonrange <= 10:
lonstep = 2
if lonrange <= 5:
lonstep = 1
lons = np.arange(-720, 720 + lonstep, lonstep)
lonticks = []
for lon in lons:
if lon >= lonmin and lon <= lonmax:
lonticks.append(lon)
lonlabels = []
for lon in lonticks:
lon2 = np.mod(lon + 180, 360) - 180
if lon2 < 0 and lon2 > -180:
if lon != 180:
lonlabels.append(str(abs(lon2)) + degsym + 'W')
if lon2 > 0 and lon2 <= 180:
lonlabels.append(str(lon2) + degsym + 'E')
if lon2 == 0:
lonlabels.append('0' + degsym )
if lon == 180 or lon == -180:
lonlabels.append('180' + degsym )
return(lonticks, lonlabels)
if type == 2:
latmin = min
latmax = max
latrange = latmax - latmin
latstep = 30
if latrange <= 90:
latstep = 10
if latrange <= 30:
latstep = 5
if latrange <= 10:
latstep = 2
if latrange <= 5:
latstep = 1
lats = np.arange(-90, 90 + latstep, latstep)
latticks = []
for lat in lats:
if lat >= latmin and lat <= latmax:
latticks.append(lat)
latlabels = []
for lat in latticks:
if lat < 0:
latlabels.append(str(abs(lat)) + degsym + 'S')
if lat > 0:
latlabels.append(str(lat) + degsym + 'N')
if lat == 0:
latlabels.append('0' + degsym)
return(latticks, latlabels)
def timeaxis(dtimes=None):
"""
| timeaxis is used to work out a sensible set of time labels and tick
| marks given a time span This is an internal routine and is not used
| by the user.
| dtimes=None - data times as a CF variable
:Returns:
time ticks and labels
|
|
|
|
|
|
|
"""
time_units = dtimes.Units
time_ticks = []
time_labels = []
axis_label = 'Time'
yearmin = min(dtimes.year.array)
yearmax = max(dtimes.year.array)
tmin = min(dtimes.dtarray)
tmax = max(dtimes.dtarray)
if plotvars.user_gset != 0:
if isinstance(plotvars.xmin, str):
t = cf.Data(cf.dt(plotvars.xmin), units=time_units)
yearmin = int(t.year)
t = cf.Data(cf.dt(plotvars.xmax), units=time_units)
yearmax = int(t.year)
tmin = cf.dt(plotvars.xmin)
tmax = cf.dt(plotvars.xmax)
if isinstance(plotvars.ymin, str):
t = cf.Data(cf.dt(plotvars.ymin), units=time_units)
yearmin = int(t.year)
t = cf.Data(cf.dt(plotvars.ymax), units=time_units)
yearmax = int(t.year)
tmin = cf.dt(plotvars.ymin)
tmax = cf.dt(plotvars.ymax)
# Years
span = yearmax - yearmin
if span > 4 and span < 3000:
axis_label = 'Time (year)'
tvals = []
if span <= 15:
step = 1
if span > 15:
step = 2
if span > 30:
step = 5
if span > 60:
step = 10
if span > 160:
step = 20
if span > 300:
step = 50
if span > 600:
step = 100
if span > 1300:
step = 200
if plotvars.tspace_year is not None:
step = plotvars.tspace_year
years = np.arange(yearmax / step + 2) * step
tvals = years[np.where((years >= yearmin) & (years <= yearmax))]
# Catch tvals if not properly defined and use gvals to generate some
# year tick marks
if np.size(tvals) < 2:
tvals = gvals(dmin=yearmin, dmax=yearmax, tight=1)[0]
for year in tvals:
time_ticks.append(np.min(
(cf.Data(cf.dt(str(int(year)) + '-01-01 00:00:00'),
units=time_units).array)))
time_labels.append(str(int(year)))
# Months
if yearmax - yearmin <= 4:
time_label = 'Time (month and year)'
months = [
'Jan',
'Feb',
'Mar',
'Apr',
'May',
'Jun',
'Jul',
'Aug',
'Sep',
'Oct',
'Nov',
'Dec']
# Check number of labels with 1 month steps
tsteps = 0
for year in np.arange(yearmax - yearmin + 1) + yearmin:
for month in np.arange(12):
mytime = cf.dt(str(year) + '-' +
str(month + 1) + '-01 00:00:00')
if mytime >= tmin and mytime <= tmax:
tsteps = tsteps + 1
if tsteps < 17:
mvals = np.arange(12)
if tsteps >= 17:
mvals = np.arange(4) * 3
for year in np.arange(yearmax - yearmin + 1) + yearmin:
for month in mvals:
mytime = cf.dt(str(year) + '-' +
str(month + 1) + '-01 00:00:00')
if mytime >= tmin and mytime <= tmax:
time_ticks.append(
np.min((cf.Data(mytime, units=time_units).array)))
time_labels.append(
str(months[month]) + ' ' + str(int(year)))
# Days and hours
if np.size(time_ticks) <= 2:
myday = cf.dt(int(tmin.year), int(tmin.month), int(tmin.day))
not_found = 0
hour_counter = 0
span = 0
while not_found <= 48:
mydate = cf.Data(myday, dtimes.Units) + \
cf.Data(hour_counter, 'hour')
if mydate >= tmin and mydate <= tmax:
span = span + 1
else:
not_found = not_found + 1
hour_counter = hour_counter + 1
step = 1
if span > 13:
step = 1
if span > 13:
step = 4
if span > 25:
step = 6
if span > 100:
step = 12
if span > 200:
step = 24
if span > 400:
step = 48
if span > 800:
step = 96
if plotvars.tspace_hour is not None:
step = plotvars.tspace_hour
if plotvars.tspace_day is not None:
step = plotvars.tspace_day * 24
not_found = 0
hour_counter = 0
span = 0
axis_label = 'Time (day)'
if span >= 24:
axis_label = 'Time (day and hour)'
time_ticks = []
time_labels = []
while not_found <= 48:
mytime = cf.Data(myday, dtimes.Units) + \
cf.Data(hour_counter, 'hour')
if mytime >= tmin and mytime <= tmax:
time_ticks.append(np.min(mytime.array))
label = str(mytime.year) + '-' + \
str(mytime.month) + '-' + str(mytime.day)
if step / 24 != step / 24.0:
label = label + ' ' + str(mytime.hour) + ':00:00'
time_labels.append(label)
else:
not_found = not_found + 1
hour_counter = hour_counter + step
return(time_ticks, time_labels, axis_label)
[docs]def ndecs(data=None):
"""
| ndecs finds the number of decimal places in an array. Needed to make the
| colour bar match the contour line labelling.
| data=data - input array of values
:Returns:
| maximum number of necimal places
|
|
|
|
|
|
|
|
"""
maxdecs = 0
for i in range(len(data)):
number = data[i]
a = str(number).split('.')
if np.size(a) == 2:
number_decs = len(a[1])
if number_decs > maxdecs:
maxdecs = number_decs
return maxdecs
[docs]def axes(xticks=None, xticklabels=None, yticks=None, yticklabels=None,
xstep=None, ystep=None, xlabel=None, ylabel=None, title=None):
"""
| axes is a function to set axes plotting parameters. The xstep and ystep
| parameters are used to label the axes starting at the left hand side and
| bottom of the plot respectively. For tighter control over labelling use
| xticks, yticks to specify the tick positions and xticklabels,
| yticklabels to specify the associated labels.
| xstep=xstep - x axis step
| ystep=ystep - y axis step
| xlabel=xlabel - label for the x-axis
| ylabel=ylabel - label for the y-axis
| xticks=xticks - values for x ticks
| xticklabels=xticklabels - labels for x tick marks
| yticks=yticks - values for y ticks
| yticklabels=yticklabels - labels for y tick marks
| title=None - set title
|
| Use axes() to reset all the axes plotting attributes to the default.
:Returns:
None
"""
if all(val is None for val in [
xticks, yticks, xticklabels, yticklabels, xstep, ystep, xlabel,
ylabel, title]):
plotvars.xticks = None
plotvars.yticks = None
plotvars.xticklabels = None
plotvars.yticklabels = None
plotvars.xstep = None
plotvars.ystep = None
plotvars.xlabel = None
plotvars.ylabel = None
plotvars.title = None
return
plotvars.xticks = xticks
plotvars.yticks = yticks
plotvars.xticklabels = xticklabels
plotvars.yticklabels = yticklabels
plotvars.xstep = xstep
plotvars.ystep = ystep
plotvars.xlabel = xlabel
plotvars.ylabel = ylabel
plotvars.title = title
def axes_plot(xticks=None, xticklabels=None, yticks=None, yticklabels=None,
xstep=None, ystep=None, xlabel=None, ylabel=None, title=None):
"""
| axes_plot is a system function to specify axes plotting parameters.
| The xstep and ystep parameters are used to label the axes starting at
| the left hand side and bottom of the plot respectively. For tighter
| control over labelling use xticks, yticks to specify the tick
| positions and xticklabels, yticklabels to specify the associated
| labels.
| xstep=xstep - x axis step
| ystep=ystep - y axis step
| xlabel=xlabel - label for the x-axis
| ylabel=ylabel - label for the y-axis
| xticks=xticks - values for x ticks
| xticklabels=xticklabels - labels for x tick marks
| yticks=yticks - values for y ticks
| yticklabels=yticklabels - labels for y tick marks
| title=None - set title
|
:Returns:
None
"""
if plotvars.plot_type == 1 or (plotvars.plot_type == 6 and plotvars.proj == 'cyl'):
xmin = plotvars.lonmin
xmax = plotvars.lonmax
ymin = plotvars.latmin
ymax = plotvars.latmax
else:
xmin = plotvars.xmin
xmax = plotvars.xmax
ymin = plotvars.ymin
ymax = plotvars.ymax
if plotvars.xstep is not None:
xstep = plotvars.xstep
xticks = None
xticklabels = None
if plotvars.ystep is not None:
ystep = plotvars.ystep
yticks = None
yticklabels = None
if plotvars.title is not None:
title = plotvars.title
title_fontsize = plotvars.title_fontsize
text_fontsize = plotvars.text_fontsize
axis_label_fontsize = plotvars.axis_label_fontsize
if title_fontsize is None:
title_fontsize = 15
if text_fontsize is None:
text_fontsize = 11
if axis_label_fontsize is None:
axis_label_fontsize = 11
axis_label_fontweight = plotvars.axis_label_fontweight
title_fontweight = plotvars.title_fontweight
text_fontweight = plotvars.text_fontweight
if (plotvars.plot_type == 1 or plotvars.plot_type == 6) and plotvars.proj == 'cyl':
plot = plotvars.mymap
lon_mid = plotvars.lonmin + (plotvars.lonmax - plotvars.lonmin) / 2.0
plotargs={'crs' : ccrs.PlateCarree()}
cyl = True
else:
plot = plotvars.plot
plotargs = {}
cyl = False
if xlabel is not None:
plotvars.plot.set_xlabel(xlabel, fontsize=axis_label_fontsize,
fontweight=axis_label_fontweight)
if ylabel is not None:
plotvars.plot.set_ylabel(ylabel, fontsize=axis_label_fontsize,
fontweight=axis_label_fontweight)
if xstep is not None:
ticks, mult = gvals(plotvars.xmin, plotvars.xmax,
tight=1, mystep=xstep)
plot_xticks = ticks * 10**mult
plot_xticks_labels = plot_xticks
if ystep is not None:
ticks, mult = gvals(plotvars.ymin, plotvars.ymax,
tight=1, mystep=ystep)
plot_yticks = ticks * 10**mult
plot_yticks_labels = plot_yticks
if xticks is not None:
plot_xticks = xticks
if xticklabels is not None:
plot_xticks_labels = xticklabels
else:
plot_xticks_labels = xticks
if yticks is not None:
plot_yticks = yticks
if yticklabels is not None:
plot_yticks_labels = yticklabels
else:
plot_yticks_labels = yticks
ticklen=np.min([(plotvars.lonmax - plotvars.lonmin)*0.014, (plotvars.latmax-plotvars.latmin)*0.014])
lonrange = plotvars.lonmax - plotvars.lonmin
latrange = plotvars.latmax - plotvars.latmin
lon_mid = plotvars.lonmin + (plotvars.lonmax - plotvars.lonmin) / 2.0
# Set the ticks and tick labels
if xticks is not None:
# fudge min and max longitude tick positions or the labels wrap
xticks_new = xticks
if lonrange >= 360:
xticks_new[0] = xticks_new[0] + 0.01
xticks_new[-1] = xticks_new[-1] - 0.01
plot.set_xticks(xticks_new, **plotargs)
plot.set_xticklabels(
xticklabels,
rotation=plotvars.xtick_label_rotation,
horizontalalignment=plotvars.xtick_label_align)
# Plot a corresponding tick on the top of the plot - **cartopy feature?
proj=ccrs.PlateCarree(central_longitude=lon_mid)
if plotvars.plot_type ==1:
for xval in xticks_new:
xpt, ypt = proj.transform_point(xval, plotvars.latmax, ccrs.PlateCarree())
ypt2 = ypt + ticklen
plot.plot([xpt, xpt], [ypt, ypt2], color='k',
linewidth=0.8, clip_on=False)
if yticks is not None:
plot.set_yticks(yticks, **plotargs)
plot.set_yticklabels(
yticklabels,
rotation=plotvars.ytick_label_rotation,
horizontalalignment=plotvars.ytick_label_align)
# Plot a corresponding tick on the right of the plot - **cartopy feature?
proj=ccrs.PlateCarree(central_longitude=lon_mid)
for ytick in yticks:
xpt, ypt = proj.transform_point(plotvars.lonmax-0.001, ytick, ccrs.PlateCarree())
xpt2 = xpt + ticklen * latrange / lonrange
plot.plot([xpt, xpt2], [ypt, ypt],
color='k', linewidth=0.8, clip_on=False)
# Set font size and weight
for label in plot.xaxis.get_ticklabels():
label.set_fontsize(axis_label_fontsize)
label.set_fontweight(axis_label_fontweight)
for label in plot.yaxis.get_ticklabels():
label.set_fontsize(axis_label_fontsize)
label.set_fontweight(axis_label_fontweight)
# Title
if title is not None:
plot.set_title(
title,
y=1.03,
fontsize=title_fontsize,
fontweight=title_fontweight)
[docs]def gset(xmin=None, xmax=None, ymin=None, ymax=None,
xlog=False, ylog=False, user_gset=1, twinx=None, twiny=None):
"""
| Set plot limits for all non longitude-latitide plots.
| xmin, xmax, ymin, ymax are all needed to set the plot limits.
| Set xlog/ylog to True or 1 to get a log axis.
|
| xmin=None - x minimum
| xmax=None - x maximum
| ymin=None - y minimum
| ymax=None - y maximum
| xlog=False - log x
| ylog=False - log y
| twinx=None - set to True to make a twin y axis plot
| twiny=None - set to True to make a twin x axis plot
|
| Once a user call is made to gset the plot limits are persistent.
| i.e. the next plot will use the same set of plot limits.
| Use gset() to reset to undefined plot limits i.e. the full range
| of the data.
|
| To set date axes use date strings i.e.
| cfp.gset(xmin = '1970-1-1', xmax = '1999-12-31', ymin = 285,
| ymax = 295)
|
| Note the correct date format is 'YYYY-MM-DD' or 'YYYY-MM-DD HH:MM:SS'
| anything else will give unexpected results.
:Returns:
None
|
|
|
|
"""
plotvars.user_gset = user_gset
if all(val is None for val in [xmin, xmax, ymin, ymax]):
plotvars.xmin = None
plotvars.xmax = None
plotvars.ymin = None
plotvars.ymax = None
plotvars.xlog = False
plotvars.ylog = False
plotvars.twinx = False
plotvars.twiny = False
plotvars.user_gset = 0
return
bcount=0
for val in [xmin, xmax, ymin, ymax]:
if val is None:
bcount=bcount + 1
if bcount !=0 and bcount != 4:
errstr = 'gset error\n'
errstr += 'xmin, xmax, ymin, ymax all need to be passed to gset\n'
errstr += 'to set the plot limits\n'
raise Warning(errstr)
plotvars.xmin = xmin
plotvars.xmax = xmax
plotvars.ymin = ymin
plotvars.ymax = ymax
plotvars.xlog = xlog
plotvars.ylog = ylog
# Check if any axes are time strings
time_xstr = False
time_ystr = False
try:
float(xmin)
except:
time_xstr = True
try:
float(ymin)
except:
time_ystr = True
# Set plot limits
if plotvars.plot is not None and twinx is None and twiny is None:
if not time_xstr and not time_ystr:
plotvars.plot.axis(
[plotvars.xmin, plotvars.xmax, plotvars.ymin, plotvars.ymax])
if plotvars.xlog:
plotvars.plot.set_xscale('log')
if plotvars.ylog:
plotvars.plot.set_yscale('log')
# Set twinx or twiny if requested
if twinx is not None:
plotvars.twinx=twinx
if twiny is not None:
plotvars.twiny=twiny
# Turn off default axis tick marks
#if plotvars.plot is not None:
# plotvars.plot.set_xticks([])
# plotvars.plot.set_yticks([])
[docs]def gopen(rows=1, columns=1, user_plot=1, file='python',
orientation='landscape', figsize=[11.7, 8.3],
left=0.12, right=0.92, top=0.92, bottom=0.08, wspace=0.2,
hspace=0.2, dpi=None, user_position=False):
"""
| gopen is used to open a graphic file.
|
| rows=1 - number of plot rows on the page
| columns=1 - number of plot columns on the page
| user_plot=1 - internal plot variable - do not use.
| file='python' - default file name
| orientation='landscape' - orientation - also takes 'portrait'
| figsize=[11.7, 8.3] - figure size in inches
| left=0.12 - left margin in normalised coordinates
| right=0.92 - right margin in normalised coordinates
| top=0.92 - top margin in normalised coordinates
| bottom=0.08 - bottom margin in normalised coordinates
| wspace=0.2 - width reserved for blank space between subplots
| hspace=0.2 - height reserved for white space between subplots
| dpi=None - resolution in dots per inch
| user_position=False - set to True to supply plot position via gpos
| xmin, xmax, ymin, ymax values
:Returns:
None
|
|
|
|
|
"""
# Set values in globals
plotvars.rows = rows
plotvars.columns = columns
if file != 'python':
plotvars.file = file
plotvars.orientation = orientation
plotvars.user_plot = user_plot
# Set user defined plot area to None
plotvars.plot_xmin = None
plotvars.plot_xmax = None
plotvars.plot_ymin = None
plotvars.plot_ymax = None
if orientation != 'landscape':
if orientation != 'portrait':
errstr = 'gopen error\n'
errstr += 'orientation incorrectly set\n'
errstr += 'input value was ' + orientation + '\n'
errstr += 'Valid options are portrait or landscape\n'
raise Warning(errstr)
# Set master plot size
if orientation == 'landscape':
plotvars.master_plot = plot.figure(figsize=(figsize[0], figsize[1]))
else:
plotvars.master_plot = plot.figure(figsize=(figsize[1], figsize[0]))
# Set margins
plotvars.master_plot.subplots_adjust(
left=left,
right=right,
top=top,
bottom=bottom,
wspace=wspace,
hspace=hspace)
# Set initial subplot
if user_position is False:
gpos(pos=1)
# Change tick length for plots > 2x2
if (columns > 2 or rows > 2):
matplotlib.rcParams['xtick.major.size'] = 2
matplotlib.rcParams['ytick.major.size'] = 2
# Set image resolution
if dpi is not None:
plotvars.dpi=dpi
[docs]def gclose(view=True):
"""
| gclose saves a graphics file. The default is to view the file as well
| - use view=0 to turn this off.
| view=True - view graphics file
:Returns:
None
|
|
|
|
|
|
|
|
|
"""
# Reset the user_plot variable to off
plotvars.user_plot = 0
file = plotvars.file
if file is not None:
type = 1
if file[-3:] == '.ps':
type = 1
if file[-4:] == '.eps':
type = 1
if file[-4:] == '.png':
type = 1
if file[-4:] == '.pdf':
type = 1
if type is None:
file = file + '.png'
plotvars.master_plot.savefig(
file, papertype='a4', orientation=plotvars.orientation, dpi=plotvars.dpi)
plot.close()
else:
if plotvars.viewer == 'display':
# Use Imagemagick display command if this exists
disp = which('display')
if disp is not None:
tfile = 'cfplot.png'
plotvars.master_plot.savefig(
tfile, papertype='a4', orientation=plotvars.orientation, dpi=plotvars.dpi)
subprocess.Popen([disp, tfile])
else:
plotvars.viewer = 'matplotlib'
if plotvars.viewer == 'matplotlib':
# Use Matplotlib viewer
plot.show()
plot.close()
# Reset plotting
plotvars.plot = None
plotvars.twinx = None
plotvars.twiny = None
plotvars.plot_xmin=None
plotvars.plot_xmax=None
plotvars.plot_ymin=None
plotvars.plot_ymax=None
def showplot(*args):
for data in args:
plot = data
plot.show()
[docs]def gpos(pos=1, xmin=None, xmax=None, ymin=None, ymax=None):
"""
| Set plot position. Plots start at top left and increase by one each plot
| to the right. When the end of the row has been reached then the next
| plot will be the leftmost plot on the next row down.
| pos=pos - plot position
|
| The following four parameters are used to get full user control
| over the plot position. In addition to these cfp.gopen
| must have the user_position=True parameter set.
| xmin=None xmin in normalised coordinates
| xmax=None xmax in normalised coordinates
| ymin=None ymin in normalised coordinates
| ymax=None ymax in normalised coordinates
|
|
:Returns:
None
|
|
|
|
|
|
|
|
"""
# Check inputs are okay
if pos < 1 or pos > plotvars.rows * plotvars.columns:
errstr = 'pos error - pos out of range:\n range = 1 - '
errstr = errstr + str(plotvars.rows * plotvars.columns)
errstr = errstr + '\n input pos was ' + str(pos)
errstr = errstr + '\n'
raise Warning(errstr)
user_pos = False
if all(val is not None for val in [xmin, xmax, ymin, ymax]):
user_pos = True
plotvars.plot_xmin = xmin
plotvars.plot_xmax = xmax
plotvars.plot_ymin = ymin
plotvars.plot_ymax = ymax
# Write user position values into plotvars as these need to
# be retrieved in set_map to add the mymap as a plot
plot
if user_pos is False:
plotvars.plot = plotvars.master_plot.add_subplot(
plotvars.rows, plotvars.columns, pos)
else:
delta_x = plotvars.plot_xmax - plotvars.plot_xmin
delta_y = plotvars.plot_ymax - plotvars.plot_ymin
plotvars.plot = plotvars.master_plot.add_axes([plotvars.plot_xmin,
plotvars.plot_ymin, delta_x, delta_y])
plotvars.plot.tick_params(which='both', direction='out', right=True, top=True)
# Set position in global variables
plotvars.pos = pos
#######################################
# pcon - convert mb to km and vice-versa
#######################################
[docs]def pcon(mb=None, km=None, h=7.0, p0=1000):
"""
| pcon is a function for converting pressure to height in kilometers and
| vice-versa. This function uses the equation P=P0exp(-z/H) to translate
| between pressure and height. In pcon the surface pressure P0 is set to
| 1000.0mb and the scale height H is set to 7.0. The value of H can vary
| from 6.0 in the polar regions to 8.5 in the tropics as well as
| seasonally. The value of P0 could also be said to be 1013.25mb rather
| than 1000.0mb.
| As this relationship is approximate:
| (i) Only use this for making the axis labels on y axis pressure plots
| (ii) Put the converted axis on the right hand side to indicate that
| this isn't the primary unit of measure
| print cfp.pcon(mb=[1000, 300, 100, 30, 10, 3, 1, 0.3])
| [0. 8.42780963 16.11809565 24.54590528 32.2361913
| 40.66400093 48.35428695, 56.78209658]
| mb=None - input pressure
| km=None - input height
| h=7.0 - default value for h
| p0=1000 - default value for p0
:Returns:
| pressure(mb) if height(km) input,
| height(km) if pressure(mb) input
"""
if all(val is None for val in [mb, km]) == 2:
errstr = 'pcon error - pcon must have mb or km input\n'
raise Warning(errstr)
if mb is not None:
return h * (np.log(p0) - np.log(mb))
if km is not None:
return np.exp(-1.0 * (np.array(km) / h - np.log(p0)))
[docs]def supscr(text=None):
"""
| supscr - add superscript text formatting for ** and ^
| This is an internal routine used in titles and colour bars
| and not used by the user.
| text=None - input text
:Returns:
Formatted text
|
|
|
|
|
|
|
"""
if text is None:
errstr = '\n supscr error - supscr must have text input\n'
raise Warning(errstr)
tform = ''
sup = 0
for i in text:
if (i == '^'):
sup = 2
if (i == '*'):
sup = sup + 1
if (sup == 0):
tform = tform + i
if (sup == 1):
if (i not in '*'):
tform = tform + '*' + i
sup = 0
if (sup == 3):
if i in '-0123456789':
tform = tform + i
else:
tform = tform + '}$' + i
sup = 0
if (sup == 2):
tform = tform + '$^{'
sup = 3
if (sup == 3):
tform = tform + '}$'
tform = tform.replace('m2', 'm$^{2}$')
tform = tform.replace('m3', 'm$^{3}$')
tform = tform.replace('m-2', 'm$^{-2}$')
tform = tform.replace('m-3', 'm$^{-3}$')
tform = tform.replace('s-1', 's$^{-1}$')
tform = tform.replace('s-2', 's$^{-2}$')
return tform
[docs]def gvals(dmin=None, dmax=None, tight=0, mystep=None, mod=1):
"""
| gvals - work out a sensible set of values between two limits
| This is an internal routine used for contour levels and axis
| labelling and is not generally used by the user.
| dmin=None - minimum
| dmax=None - maximum
| tight=0 - return values tight to input min and max
| mystep=None - use this step
| mod=1 - modify data to make use of a multipler
|
|
|
|
|
|
"""
if all(val is None for val in [dmin, dmax]) > 0:
errstr = '\n gvals error - gvals must have dmin and dmax input\n'
raise Warning(errstr)
# Return some values if dmin = dmax
if dmin == dmax:
vals = [dmin - 0.001, dmin, dmin + 0.001]
mult = 0
return vals, mult
# Copies of inputs
dmin1 = deepcopy(dmin)
dmax1 = deepcopy(dmax)
dmin2 = deepcopy(dmin)
dmax2 = deepcopy(dmax)
mult = 0 # field multiplyer
# Generate an inital step
step = (dmax - dmin) / 16.0
# Don't modify if dmin1 and dmax1 are both negative as this will create a
# race condition
if dmin1 < 0 and dmax1 <= 0:
mod = 0
# Small range code
# Return a linspace set of values if dmax1 - dmin1 is <= 3
if dmax1 - dmin1 <= 3:
mult = 0
mindecs = 40 # minimum number of decimal places - a high starting value
minval = 0 # spacing value of the minimum number of decimal spaces
for val in np.arange(19)+2:
vals = np.linspace(dmin, dmax, val)
nd = ndecs(vals)
if nd <= mindecs:
minval = val
mindecs = nd
vals = np.linspace(dmin, dmax, minval)
# Try to recalculate if too many decimal places
if mindecs >=4:
step=0.1
for val in np.arange(40)+1:
if dmax1-dmin1 <= 1.0/10**val:
step = 1.0/(10**(val+1))
newvals = [float(int(dmin1))]
while np.max(newvals) < dmax1:
newvals = np.append(newvals, newvals[-1] + step)
# Remove any points outside of limits
pts = np.where((newvals >= dmin1) & (newvals <= dmax1))
newvals = newvals[pts ]
# Remove upper and lower limits if tight=0 - i.e. a contour plot
if tight == 0 and np.size(vals) > 1:
if np.nanmax(newvals) >= dmax1:
newvals = newvals[0:-1]
if np.nanmin(newvals) <= dmin1:
newvals = newvals[1:]
if np.size(newvals) >2:
vals=newvals
return vals, mult
if mod == 1:
if (mystep is not None):
step = mystep
if step < 1:
while dmax1 <= 3:
step = step * 10.0
dmin1 = dmin1 * 10.0
dmax1 = dmax1 * 10.0
mult = mult - 1
if int(dmax2 - dmin2) > 3:
step = 1
mult = 0
for val in [2, 5, 10, 20, 25, 50, 100, 200, 250, 500, 1000, 2000,
2500, 5000, 10000, 20000, 25000, 50000, 100000,
200000, 250000, 500000, 1000000, 10000000 ]:
if int(dmax2 - dmin2) / val > 10:
step = val
mult = 0
if int(dmax2 - dmin2) / 1000000 > 10:
step = (dmax2 - dmin2) / 16.0
# Reset to user or non-zero values
if (mystep is not None):
step = mystep
if (step == 0):
step = 1
# Make integer step
if tight == 0:
vals = (int(dmin1) / step) * step
else:
vals = dmin1
while (np.nanmax(vals) + step) <= dmax1:
vals = np.append(vals, np.nanmax(vals) + step)
# Remove upper and lower limits if tight=0 - i.e. a contour plot
if tight == 0 and np.size(vals) > 1:
if np.nanmax(vals) >= dmax1:
vals = vals[0:-1]
if np.nanmin(vals) <= dmin1:
vals = vals[1:]
# Return values if there are five or more
if np.size(vals) >= 5:
return vals, mult
if mystep is not None:
if int(mystep) == mystep:
return vals, mult
if step == .9:
step = 1.0
if step == .8:
step = 1.0
if step == .7:
step = .5
if step == .6:
step = .5
if step == .3:
step = .25
if step == .09:
step = 0.1
if step == .08:
step = 0.1
if step == .07:
step = .1
if step == .06:
step = .1
if step == .04:
step = .1
if step == .03:
step = .02
if (dmax1 - dmin1 == step):
step = step / 10.
if step == 0.0:
step = (dmax1-dmin1)/16.
vals = float("%.2f" % (int(dmin1 / step) * step))
while (np.nanmax(vals) + step) <= dmax1:
vals = np.append(vals, float("%.2f" % (np.nanmax(vals) + step)))
if tight == 0:
if np.nanmax(vals) >= dmax1:
vals = vals[0:-1]
if np.nanmin(vals) <= dmin1:
vals = vals[1:]
return vals, mult
[docs]def cf_data_assign(f=None, colorbar_title=None, verbose=None, rotated_vect=False):
"""
| Check cf input data is okay and return data for contour plot.
| This is an internal routine not used by the user.
| f=None - input cf field
| colorbar_title=None - input colour bar title
| rotated vect=False - return 1D x and y for rotated plot vectors
| verbose=None - set to 1 to get a verbose idea of what the
| cf_data_assign is doing
:Returns:
| f - data for contouring
| x - x coordinates of data (optional)
| y - y coordinates of data (optional)
| ptype - plot type
| colorbar_title - colour bar title
| xlabel - x label for plot
| ylabel - y label for plot
|
|
|
|
|
"""
# Check input data has the correct number of dimensions
# Take into account rotated pole fields having extra dimensions
ndim = len(f.axes(size=cf.gt(1)))
if f.ref('rotated_latitude_longitude') is None:
if (ndim > 2 or ndim < 1):
print ''
if (ndim > 2):
errstr = 'cf_data_assign error - data has too many dimensions'
if (ndim < 1):
errstr = 'cf_data_assign error - data has too few dimensions'
errstr += '\n cf-plot requires one or two dimensional data\n'
for mydim in f.items():
sn = getattr(f.item(mydim), 'standard_name', False)
ln = getattr(f.item(mydim), 'long_name', False)
if sn:
errstr = errstr + \
str(mydim) + ',' + str(sn) + ',' + \
str(f.item(mydim).size) + '\n'
else:
if ln:
errstr = errstr + \
str(mydim) + ',' + str(ln) + ',' + \
str(f.item(mydim).size) + '\n'
raise Warning(errstr)
# Set up data arrays and variables
lons = None
lats = None
height = None
time = None
xlabel = ''
ylabel = ''
has_lons = None
has_lats = None
has_height = None
has_time = None
xpole = None
ypole = None
ptype = None
# Extract coordinate data if a matching CF standard_name or axis is found
for mydim in f.items():
sn = getattr(f.item(mydim), 'standard_name', 'NoName')
an = getattr(f.item(mydim), 'axis', 'NoName')
standard_name_x = ['longitude']
vs = 'cf_data_assign standard_name, axis - assigned '
if (sn in standard_name_x or an == 'X'):
if verbose:
print vs + 'lons -', sn, an
lons = np.squeeze(f.item(mydim).array)
standard_name_y = ['latitude']
if (sn in standard_name_y or an == 'Y'):
if verbose:
print vs + 'lats -', sn, an
lats = np.squeeze(f.item(mydim).array)
standard_name_z = ['pressure', 'air_pressure', 'height', 'depth']
if (sn in standard_name_z or an == 'Z'):
if verbose:
print vs + 'height -', sn, an
height = np.squeeze(f.item(mydim).array)
standard_name_t = ['time']
if (sn in standard_name_t or an == 'T'):
if verbose:
print vs + 'time -', sn, an
time = np.squeeze(f.item(mydim).array)
# CF defined units
lon_units = [
'degrees_east',
'degree_east',
'degree_E',
'degrees_E',
'degreeE',
'degreesE']
lat_units = [
'degrees_north',
'degree_north',
'degree_N',
'degrees_N',
'degreeN',
'degreesN']
height_units = ['mb', 'mbar', 'millibar', 'decibar', 'atmosphere',
'atm', 'pascal', 'Pa', 'hPa', 'meter', 'metre', 'm',
'kilometer', 'kilometre', 'km']
time_units = ['day', 'days', 'd', 'hour', 'hours', 'hr', 'h', 'minute',
'minutes', 'min', 'mins', 'second', 'seconds', 'sec',
'secs', 's']
# Extract coordinate data if a matching CF set of units is found
for mydim in f.items():
units = getattr(f.item(mydim), 'units', False)
if units in lon_units:
if lons is None:
if verbose:
print 'cf_data_assign units - assigned lons -', units
lons = np.squeeze(f.item(mydim).array)
if units in lat_units:
if lats is None:
if verbose:
print 'cf_data_assign units - assigned lats -', units
lats = np.squeeze(f.item(mydim).array)
if units in height_units:
if height is None:
if verbose:
print 'cf_data_assign units - assigned height -', units
height = np.squeeze(f.item(mydim).array)
if units in time_units:
if time is None:
if verbose:
print 'cf_data_assign units - assigned time -', units
time = np.squeeze(f.item(mydim).array)
# Extract coordinate data from variable name if not already assigned
for mydim in f.items():
if mydim[:3] == 'dim':
name = cf_var_name(field=f, dim=mydim)
vs = 'cf_data_assign dimension name - assigned '
if name[0:3] == 'lon':
if lons is None:
if verbose:
print vs + lons + '-', name
lons = np.squeeze(f.item(mydim).array)
if name[0:3] == 'lat':
if lats is None:
if verbose:
print vs+lats + '-', name
lats = np.squeeze(f.item(mydim).array)
if (name[0:5] == 'theta' or
name[0:1] == 'p' or
name == 'air_pressure'):
if height is None:
if verbose:
print vs+height + '-', name
height = np.squeeze(f.item(mydim).array)
if name[0:1] == 't':
if time is None:
if verbose:
print vs+time + '-', name
time = np.squeeze(f.item(mydim).array)
if np.size(lons) > 1:
has_lons = 1
if np.size(lats) > 1:
has_lats = 1
if np.size(height) > 1:
has_height = 1
if np.size(time) > 1:
has_time = 1
# assign field data
field = np.squeeze(f.array)
# Check what plot type is required.
# 0=simple contour plot, 1=map plot, 2=latitude-height plot,
# 3=longitude-time plot, 4=latitude-time plot.
if (np.size(lons) > 1 and np.size(lats) > 1):
ptype = 1
x = lons
y = lats
if (np.size(lats) > 1 and np.size(height) > 1):
ptype = 2
x = lats
y = height
for mydim in f.items():
name = cf_var_name(field=f, dim=mydim)
if name[0:3] == 'lat':
xunits = str(getattr(f.item(mydim), 'Units', ''))
if (xunits in lat_units):
xunits = 'degrees'
xlabel = name + ' (' + xunits + ')'
if name[0:1] == 'p' or \
name[0:5] == 'theta' or \
name[0:6] == 'height' or \
name[0:6] == 'hybrid' or \
name[0:5] == 'level' or \
name[0:5] == 'model':
yunits = str(getattr(f.item(mydim), 'Units', ''))
ylabel = name + ' (' + yunits + ')'
if (np.size(lons) > 1 and np.size(height) > 1):
ptype = 3
x = lons
y = height
for mydim in f.items():
name = cf_var_name(field=f, dim=mydim)
if name[0:3] == 'lon':
xunits = str(getattr(f.item(mydim), 'Units', ''))
if (xunits in lon_units):
xunits = 'degrees'
xlabel = name + ' (' + xunits + ')'
if name[0:1] == 'p' or \
name[0:5] == 'theta' or \
name[0:6] == 'height' or \
name[0:6] == 'hybrid' or \
name[0:5] == 'level' or \
name[0:5] == 'model':
yunits = str(getattr(f.item(mydim), 'Units', ''))
ylabel = name + ' (' + yunits + ')'
if (np.size(lons) > 1 and np.size(time) > 1):
ptype = 4
x = lons
y = time
if np.size(lats) > 1 and np.size(time) > 1:
ptype = 5
x = lats
y = time
# Rotated pole
if f.ref('rotated_latitude_longitude') is not None:
ptype = 6
rotated_pole = f.ref('rotated_latitude_longitude')
xpole = rotated_pole['grid_north_pole_longitude']
ypole = rotated_pole['grid_north_pole_latitude']
# Extract grid x and y coordinates
for mydim in f.items():
name = cf_var_name(field=f, dim=mydim)
if mydim[:3] == 'dim':
if name in ['grid_longitude', 'longitude', 'x']:
x = np.squeeze(f.item(mydim).array)
xunits = str(getattr(f.item(mydim), 'units', ''))
xlabel = cf_var_name(field=f, dim=mydim)
if name in ['grid_latitude', 'latitude', 'y']:
y = np.squeeze(f.item(mydim).array)
# Flip y and data if reversed
if y[0] > y[-1]:
y = y[::-1]
field = np.flipud(field)
yunits = str(getattr(f.item(mydim), 'Units', ''))
ylabel = cf_var_name(field=f, dim=mydim) + yunits
# Extract auxiliary lons and lats if they exist
if plotvars.proj != 'rotated' and not rotated_vect:
has_lons = False
has_lats = False
for mydim in f.items():
if mydim[:3] == 'aux':
name = cf_var_name(field=f, dim=mydim)
if name in ['longitude']:
x = np.squeeze(f.item(mydim).array)
has_lons = True
if name in ['latitude']:
y = np.squeeze(f.item(mydim).array)
has_lats = True
if not has_lons or not has_lats:
# Generate the grid using the x and y from above
xvals, yvals = np.meshgrid(x, y)
x = xvals
y = yvals
#rotated_pole = f.ref('rotated_latitude_longitude')
#xpole = rotated_pole['grid_north_pole_longitude']
#ypole = rotated_pole['grid_north_pole_latitude']
#transform = ccrs.RotatedPole(pole_latitude=ypole,
# pole_longitude=xpole)
#points=ccrs.PlateCarree().transform_points(transform,
# xvals.flatten(),
# yvals.flatten())
#x=np.array(points)[:, 0]
#y=np.array(points)[:, 1]
# time height plot
if has_height == 1 and has_time == 1:
ptype = 7
for mydim in f.items():
if mydim[:3] == 'dim':
if np.size(np.squeeze(f.item(mydim).array)
) == np.shape(np.squeeze(f.array))[0]:
x = np.squeeze(f.item(mydim).array)
xunits = str(getattr(f.item(mydim), 'units', ''))
xlabel = cf_var_name(field=f, dim=mydim) + xunits
if np.size(np.squeeze(f.item(mydim).array)
) == np.shape(np.squeeze(f.array))[1]:
y = np.squeeze(f.item(mydim).array)
yunits = '(' + str(getattr(f.item(mydim), 'Units', ''))
yunits += ')'
ylabel = cf_var_name(field=f, dim=mydim) + yunits
# Rotate array to get it as time vs height
field = np.rot90(field)
field = np.flipud(field)
# UKCP grid
if f.ref('transverse_mercator'):
ptype = 1
field = np.squeeze(f.array)
# Find the auxiliary lons and lats if provided
has_lons = False
has_lats = False
for mydim in f.items():
if mydim[:3] == 'aux':
name = cf_var_name(field=f, dim=mydim)
if name in ['longitude']:
x = np.squeeze(f.item(mydim).array)
has_lons = True
if name in ['latitude']:
y = np.squeeze(f.item(mydim).array)
has_lats = True
# Calculate lons and lats if no auxiliary data for these
if not has_lons or not has_lats:
xpts = f.item('X').array
ypts = f.item('Y').array
field = np.squeeze(f.array)
false_easting = f.ref('transverse_mercator')['false_easting']
false_northing = f.ref('transverse_mercator')['false_northing']
central_longitude = f.ref('transverse_mercator')['longitude_of_central_meridian']
central_latitude = f.ref('transverse_mercator')['latitude_of_projection_origin']
scale_factor = f.ref('transverse_mercator')['scale_factor_at_central_meridian']
# Set the transform
transform = ccrs.TransverseMercator(false_easting = false_easting,
false_northing = false_northing,
central_longitude = central_longitude,
central_latitude = central_latitude,
scale_factor = scale_factor)
# Calculate the longitude and latitude points
xvals, yvals = np.meshgrid(xpts, ypts)
points = ccrs.PlateCarree().transform_points(transform, xvals, yvals)
x = np.array(points)[:, :, 0]
y = np.array(points)[:, :, 1]
# None of the above
if ptype is None:
ptype = 0
for mydim in f.items():
if mydim[:3] == 'dim':
if np.size(np.squeeze(f.item(mydim).array)
) == np.shape(np.squeeze(f.array))[1]:
x = np.squeeze(f.item(mydim).array)
xunits = str(getattr(f.item(mydim), 'units', ''))
xlabel = cf_var_name(field=f, dim=mydim) + xunits
if np.size(np.squeeze(f.item(mydim).array)
) == np.shape(np.squeeze(f.array))[0]:
y = np.squeeze(f.item(mydim).array)
yunits = str(getattr(f.item(mydim), 'Units', ''))
ylabel = cf_var_name(field=f, dim=mydim) + yunits
# Assign colorbar_title
if (colorbar_title is None):
colorbar_title = 'No Name'
if hasattr(f, 'id'):
colorbar_title = f.id
if hasattr(f, 'ncvar'):
colorbar_title = f.ncvar
if hasattr(f, 'short_name'):
colorbar_title = f.short_name
if hasattr(f, 'long_name'):
colorbar_title = f.long_name
if hasattr(f, 'standard_name'):
colorbar_title = f.standard_name
if hasattr(f, 'Units'):
if str(f.Units) == '':
colorbar_title = colorbar_title + ''
else:
colorbar_title = colorbar_title + \
'(' + supscr(str(f.Units)) + ')'
# Return data
return(field, x, y, ptype, colorbar_title, xlabel, ylabel, xpole, ypole)
[docs]def check_data(field=None, x=None, y=None):
"""
| check_data - check user input contour data is correct.
| This is an internal routine and is not used by the user.
|
| field=None - field
| x=None - x points for field
| y=None - y points for field
|
|
|
|
|
|
"""
# Input error trapping
args = True
errstr = '\n'
if np.size(field) == 1:
if field is None:
errstr = errstr + 'con error - a field for contouring must be '
errstr += 'passed with the f= flag\n'
args = False
if np.size(x) == 1:
if x is None:
x = np.arange(np.shape(field)[1])
if np.size(y) == 1:
if y is None:
y = np.arange(np.shape(field)[0])
if not args:
raise Warning(errstr)
# Check input dimensions look okay.
# All inputs 2D
if np.ndim(field) == 2 and np.ndim(x) == 2 and np.ndim(y) == 2:
xpts = np.shape(field)[1]
ypts = np.shape(field)[0]
if xpts != np.shape(x)[1] or xpts != np.shape(y)[1]:
args = False
if ypts != np.shape(x)[0] or ypts != np.shape(y)[0]:
args = False
if args:
return
# Field x and y all 1D
if np.ndim(field) == 1 and np.ndim(x) == 1 and np.ndim(y) == 1:
if np.size(x) != np.size(field):
args = False
if np.size(y) != np.size(field):
args = False
if args:
return
# Field 2D, x and y 1D
if np.ndim(field) != 2:
args = False
if np.ndim(x) != 1:
args = False
if np.ndim(y) != 1:
args = False
if np.ndim(field) == 2:
if np.size(x) != np.shape(field)[1]:
args = False
if np.size(y) != np.shape(field)[0]:
args = False
if args is False:
errstr = errstr + 'Input arguments incorrectly shaped:\n'
errstr = errstr + 'x has shape:' + str(np.shape(x)) + '\n'
errstr = errstr + 'y has shape:' + str(np.shape(y)) + '\n'
errstr = errstr + 'field has shape' + str(np.shape(field)) + '\n\n'
errstr = errstr + 'Expected x=xpts, y=ypts, field=(ypts,xpts)\n'
errstr = errstr + 'x=npts, y=npts, field=npts\n'
errstr = errstr + \
'or x=[ypts, xpts], y=[ypts, xpts], field=[ypts, xpts]\n'
raise Warning(errstr)
[docs]def cscale(cmap=None, ncols=None, white=None, below=None,
above=None, reverse=False, uniform=False):
"""
| cscale - choose and manipulate colour maps. Around 200 colour scales are
| available - see the gallery section for more details.
|
| cmap=None - name of colour map
| ncols=None - number of colours for colour map
| white=None - change these colours to be white
| below=None - change the number of colours below the mid point of
| the colour scale to be this
| above=None - change the number of colours above the mid point of
| the colour scale to be this
| reverse=False - reverse the colour scale
| uniform=False - produce a uniform colour scale.
| For example: if below=3 and above=10 are specified
| then initially below=10 and above=10 are used. The
| colour scale is then cropped to use scale colours
| 6 to 19. This produces a more uniform intensity colour
| scale than one where all the blues are compressed into
| 3 colours.
|
|
| Personal colour maps are available by saving the map as red green blue
| to a file with a set of values on each line.
|
|
| Use cscale() To reset to the default settings.
|
:Returns:
None
|
|
|
|
"""
# If no map requested reset to default
if cmap is None:
cmap = 'scale1'
plotvars.cscale_flag = 0
return
else:
plotvars.cs_user = cmap
plotvars.cscale_flag = 1
vals = [ncols, white, below, above]
if any(val is not None for val in vals):
plotvars.cscale_flag = 2
if reverse is not False or uniform is not False:
plotvars.cscale_flag = 2
if cmap == 'scale1' or cmap == '':
if cmap == 'scale1':
myscale = cscale1
if cmap == 'viridis':
myscale = viridis
# convert cscale1 or viridis from hex to rgb
r = []
g = []
b = []
for myhex in myscale:
myhex = myhex.lstrip('#')
mylen = len(myhex)
rgb = tuple(int(myhex[i:i + mylen / 3], 16)
for i in range(0, mylen, mylen / 3))
r.append(rgb[0])
g.append(rgb[1])
b.append(rgb[2])
else:
import distutils.sysconfig as sysconfig
package_path = os.path.dirname(__file__)
file = os.path.join(package_path, 'colourmaps/' + cmap + '.rgb')
if os.path.isfile(file) is False:
if os.path.isfile(cmap) is False:
errstr = '\ncscale error - colour scale not found:\n'
errstr = errstr + 'File ' + file + ' not found\n'
errstr = errstr + 'File ' + cmap + ' not found\n'
raise Warning(errstr)
else:
file = cmap
# Read in rgb values and convert to hex
f = open(file, 'r')
lines = f.read()
lines = lines.splitlines()
r = []
g = []
b = []
hex = []
for line in lines:
vals = line.split()
r.append(int(vals[0]))
g.append(int(vals[1]))
b.append(int(vals[2]))
# Reverse the colour scale if requested
if reverse:
r = r[::-1]
g = g[::-1]
b = b[::-1]
# Interpolate to a new number of colours if requested
if ncols is not None:
x = np.arange(np.size(r))
xnew = np.linspace(0, np.size(r) - 1, num=ncols, endpoint=True)
f_red = interpolate.interp1d(x, r)
f_green = interpolate.interp1d(x, g)
f_blue = interpolate.interp1d(x, b)
r = f_red(xnew)
g = f_green(xnew)
b = f_blue(xnew)
# Change the number of colours below and above the mid-point if requested
if below is not None or above is not None:
# Mid-point of colour scale
npoints = np.size(r) / 2
# Below mid point x locations
x_below = []
lower = 0
if below == 1:
x_below = 0
if below is not None:
lower = below
if below is None:
lower = npoints
if below is not None and uniform:
lower = max(above, below)
if (lower > 1):
x_below = ((npoints - 1) / float(lower - 1)) * np.arange(lower)
# Above mid point x locations
x_above = []
upper = 0
if above == 1:
x_above = npoints * 2 - 1
if above is not None:
upper = above
if above is None:
upper = npoints
if above is not None and uniform:
upper = max(above, below)
if (upper > 1):
x_above = ((npoints - 1) / float(upper - 1)) * \
np.arange(upper) + npoints
# Append new colour positions
xnew = np.append(x_below, x_above)
# Interpolate to new colour scale
xpts = np.arange(np.size(r))
f_red = interpolate.interp1d(xpts, r)
f_green = interpolate.interp1d(xpts, g)
f_blue = interpolate.interp1d(xpts, b)
r = f_red(xnew)
g = f_green(xnew)
b = f_blue(xnew)
# Reset colours if uniform is set
if uniform:
mid_pt = max(below, above)
r = r[mid_pt - below:mid_pt + above]
g = g[mid_pt - below:mid_pt + above]
b = b[mid_pt - below:mid_pt + above]
# Convert to hex
hex = []
for col in np.arange(np.size(r)):
hex.append('#%02x%02x%02x' % (r[col], g[col], b[col]))
# White requested colour positions
if white is not None:
ccols = white
if np.size(white) == 1:
hex[white] = '#ffffff'
else:
for col in white:
hex[col] = '#ffffff'
# Set colour scale
plotvars.cs = hex
[docs]def cscale_get_map():
"""
| cscale_get_map - return colour map for use in contour plots.
| This depends on the colour bar extensions
| This is an internal routine and is not used by the user.
|
|
:Returns:
colour map
|
|
|
|
|
"""
cscale_ncols = np.size(plotvars.cs)
if (plotvars.levels_extend == 'both'):
colmap = plotvars.cs[1:cscale_ncols - 1]
if (plotvars.levels_extend == 'min'):
colmap = plotvars.cs[1:]
if (plotvars.levels_extend == 'max'):
colmap = plotvars.cs[:cscale_ncols - 1]
if (plotvars.levels_extend == 'neither'):
colmap = plotvars.cs
return (colmap)
[docs]def bfill(f=None, x=None, y=None, clevs=False, lonlat=False, bound=False,
alpha=1.0, single_fill_color=None, white=True, zorder=None):
"""
| bfill - block fill a field with colour rectangles
| This is an internal routine and is not generally used by the user.
|
| f=None - field
| x=None - x points for field
| y=None - y points for field
| clevs=None - levels for filling
| lonlat=False - lonlat data
| bound=False - x and y are cf data boundaries
| alpha=alpha - transparency setting 0 to 1
| white=True - colour unplotted areas white
| single_fill_color=None - colour for a blockfill between two levels
| - makes maplotlib named colours or
| - hexadecimal notation - '#d3d3d3' for grey
| zorder=None - plotting order
|
:Returns:
None
|
|
|
|
"""
# If single_fill_color is defined then turn off whiting out the background.
if single_fill_color is not None:
white = False
# Set the default map coordinates for the data to be PlateCarree
plotargs = {}
if lonlat:
plotargs = {'transform' : ccrs.PlateCarree()}
transverse_mercator = False
if isinstance(f, cf.Field):
if f.ref('transverse_mercator'):
transverse_mercator = True
lonlat = True
# Case of transverse mercator of which UKCP is an example
false_easting = f.ref('transverse_mercator')['false_easting']
false_northing = f.ref('transverse_mercator')['false_northing']
central_longitude = f.ref('transverse_mercator')['longitude_of_central_meridian']
central_latitude = f.ref('transverse_mercator')['latitude_of_projection_origin']
scale_factor = f.ref('transverse_mercator')['scale_factor_at_central_meridian']
transform = ccrs.TransverseMercator(false_easting = false_easting,
false_northing = false_northing,
central_longitude = central_longitude,
central_latitude = central_latitude,
scale_factor = scale_factor)
# Extract the axes and data
xpts = np.append(f.item('X').bounds.array[:,0], f.item('X').bounds.array[-1,1])
ypts = np.append(f.item('Y').bounds.array[:,0], f.item('Y').bounds.array[-1,1])
field = np.squeeze(f.array)
else:
# Assign f to field as this may be modified in lat-lon plots
field = f
if bound:
xpts = x
ypts = y
else:
# Find x box boundaries
xpts = x[0] - (x[1] - x[0]) / 2.0
for ix in np.arange(np.size(x) - 1):
xpts = np.append(xpts, x[ix] + (x[ix + 1] - x[ix]) / 2.0)
xpts = np.append(xpts, x[ix + 1] + (x[ix + 1] - x[ix]) / 2.0)
# Find y box boundaries
ypts = y[0] - (y[1] - y[0]) / 2.0
for iy in np.arange(np.size(y) - 1):
ypts = np.append(ypts, y[iy] + (y[iy + 1] - y[iy]) / 2.0)
ypts = np.append(ypts, y[iy + 1] + (y[iy + 1] - y[iy]) / 2.0)
# Shift lon grid if needed
if lonlat:
# Extract upper bound and original rhs of box longitude bounding points
upper_bound = ypts[-1]
xpts_orig = xpts
ypts_orig = ypts
# Reduce xpts and ypts by 1 or shifting of grid fails
# The last points are the right / upper bounds for the last data box
xpts = xpts[0:-1]
ypts = ypts[0:-1]
if plotvars.lonmin < np.nanmin(xpts):
xpts = xpts - 360
if plotvars.lonmin > np.nanmax(xpts):
xpts = xpts + 360
# Add cyclic information if missing.
lonrange = np.nanmax(xpts) - np.nanmin(xpts)
if lonrange < 360:
field, xpts = cartopy_util.add_cyclic_point(field, xpts)
right_bound = xpts[-1] + (xpts[-1] - xpts[-2])
# Add end x and y end points
xpts = np.append(xpts, right_bound)
ypts = np.append(ypts, upper_bound)
levels=np.array(deepcopy(clevs)).astype('float')
# Polar stereograpic
# Set points past plotting limb to be plotvars.boundinglat
# Also set any lats past the pole to be the pole
if plotvars.proj == 'npstere':
pts = np.where(ypts < plotvars.boundinglat)
if np.size(pts) > 0:
ypts[pts] = plotvars.boundinglat
pts=np.where(ypts > 90.0)
if np.size(pts) > 0:
ypts[pts] = 90.0
if plotvars.proj == 'spstere':
pts = np.where(ypts > plotvars.boundinglat)
if np.size(pts) > 0:
ypts[pts]=plotvars.boundinglat
pts=np.where(ypts < -90.0)
if np.size(pts) > 0:
ypts[pts] = -90.0
# Generate a Matplotlib colour map
if single_fill_color is None:
cols=plotvars.cs
else:
cols=single_fill_color
cmap = matplotlib.colors.ListedColormap(cols)
if single_fill_color is None:
if plotvars.levels_extend == 'both' or plotvars.levels_extend == 'min':
levels=np.insert(levels, 0, -1e30)
if plotvars.levels_extend == 'both' or plotvars.levels_extend == 'max':
levels=np.append(levels, 1e30)
if plotvars.levels_extend == 'both' or plotvars.levels_extend == 'min':
cmap.set_under(plotvars.cs[0])
cols=cols[1:]
if plotvars.levels_extend == 'both' or plotvars.levels_extend == 'max':
cmap.set_over(plotvars.cs[-1])
cols=cols[:-1]
# Colour array for storing the cell colour. Start with -1 as the default
# as the colours run from 0 to np.size(levels)-1
colarr=np.zeros([np.shape(field)[0], np.shape(field)[1]])-1
for i in np.arange(np.size(levels)-1):
lev=levels[i]
pts=np.where(np.logical_and(field > lev, field <= levels[i+1]))
colarr[pts]=int(i)
# Change points that are masked back to -1
if isinstance(field, np.ma.MaskedArray):
pts=np.ma.where(field.mask)
if np.size(pts) > 0:
colarr[pts]=-1
#print 'plotargs are ', **plotargs
if plotvars.plot_type == 1 and plotvars.proj != 'cyl':
for i in np.arange(np.size(levels)-1):
allverts = []
xy_stack=np.column_stack(np.where(colarr == i))
for pt in np.arange(np.shape(xy_stack)[0]):
ix = xy_stack[pt][1]
iy = xy_stack[pt][0]
lons = [xpts[ix], xpts[ix+1], xpts[ix+1], xpts[ix], xpts[ix]]
lats = [ypts[iy], ypts[iy], ypts[iy+1], ypts[iy+1], ypts[iy]]
#txpts, typts = plotvars.mymap(lons, lats)
txpts, typts = lons, lats
verts=[
(txpts[0], typts[0]),
(txpts[1], typts[1]),
(txpts[2], typts[2]),
(txpts[3], typts[3]),
(txpts[4], typts[4]),
]
allverts.append(verts)
# Make the collection and add it to the plot.
if single_fill_color is None:
color=plotvars.cs[i]
else:
color=single_fill_color
coll = PolyCollection(allverts, facecolor=color, edgecolors=color, alpha=alpha,
zorder=zorder, **plotargs)
if lonlat:
plotvars.mymap.add_collection(coll)
else:
plotvars.plot.add_collection(coll)
else:
for i in np.arange(np.size(levels)-1):
allverts = []
xy_stack=np.column_stack(np.where(colarr == i))
for pt in np.arange(np.shape(xy_stack)[0]):
ix = xy_stack[pt][1]
iy = xy_stack[pt][0]
verts=[
(xpts[ix], ypts[iy]),
(xpts[ix+1], ypts[iy]),
(xpts[ix+1], ypts[iy+1]),
(xpts[ix], ypts[iy+1]),
(xpts[ix], ypts[iy]),
]
allverts.append(verts)
# Make the collection and add it to the plot.
if single_fill_color is None:
color=plotvars.cs[i]
else:
color=single_fill_color
coll = PolyCollection(allverts, facecolor=color, edgecolors=color,
alpha=alpha, zorder=zorder, **plotargs)
if lonlat:
plotvars.mymap.add_collection(coll)
else:
plotvars.plot.add_collection(coll)
# Add white for undefined areas
if white:
allverts = []
xy_stack=np.column_stack(np.where(colarr == -1))
for pt in np.arange(np.shape(xy_stack)[0]):
ix = xy_stack[pt][1]
iy = xy_stack[pt][0]
verts=[
(xpts[ix], ypts[iy]),
(xpts[ix+1], ypts[iy]),
(xpts[ix+1], ypts[iy+1]),
(xpts[ix], ypts[iy+1]),
(xpts[ix], ypts[iy]),
]
allverts.append(verts)
# Make the collection and add it to the plot.
color=plotvars.cs[i]
coll = PolyCollection(allverts, facecolor='#ffffff', edgecolors='#ffffff',
alpha=alpha, zorder=zorder, **plotargs)
if lonlat:
plotvars.mymap.add_collection(coll)
else:
plotvars.plot.add_collection(coll)
[docs]def regrid(f=None, x=None, y=None, xnew=None, ynew=None):
"""
| regrid - bilinear interpolation of a grid to new grid locations
|
|
| f=None - original field
| x=None - original field x values
| y=None - original field y values
| xnew=None - new x points
| ynew=None - new y points
|
:Returns:
field values at requested locations
|
|
"""
# reassign input arrays
regrid_f = f
regrid_x = x
regrid_y = y
fieldout = []
# Reverse xpts and field if necessary
if regrid_x[0] > regrid_x[-1]:
regrid_x = regrid_x[::-1]
field = np.fliplr(regrid_f)
# Reverse ypts and field if necessary
if regrid_y[0] > regrid_y[-1]:
regrid_y = regrid_y[::-1]
regrid_f = np.flipud(regrid_f)
# Iterate over the new grid to get the new grid values.
for i in np.arange(np.size(xnew)):
xval = xnew[i]
yval = ynew[i]
# Find position of new grid point in the x and y arrays
myxpos = find_pos_in_array(vals=regrid_x, val=xval)
#print 'xval is ', xval
#print 'vals are ', regrid_x
#print 'myxpos is ', myxpos
myypos = find_pos_in_array(vals=regrid_y, val=yval)
myxpos2 = myxpos + 1
myypos2 = myypos + 1
if (myxpos2 != myxpos):
alpha = (xnew[i] - regrid_x[myxpos]) / \
(regrid_x[myxpos2] - regrid_x[myxpos])
else:
alpha = (xnew[i] - regrid_x[myxpos]) / 1E-30
newval1 = (regrid_f[myypos, myxpos] - regrid_f[myypos, myxpos2])
newval1 = newval1 * alpha
newval1 = regrid_f[myypos, myxpos] - newval1
newval2 = (regrid_f[myypos2, myxpos] - regrid_f[myypos2, myxpos2])
newval2 = newval2 * alpha
newval2 = regrid_f[myypos2, myxpos] - newval2
if (myypos2 != myypos):
alpha2 = (ynew[i] - regrid_y[myypos])
alpha2 = alpha2 / (regrid_y[myypos2] - regrid_y[myypos])
else:
alpha2 = (ynew[i] - regrid_y[myypos]) / 1E-30
newval3 = newval1 - (newval1 - newval2) * alpha2
fieldout = np.append(fieldout, newval3)
return fieldout
[docs]def stipple(f=None, x=None, y=None, min=None, max=None,
size=80, color='k', pts=50, marker='.', edgecolors='k',
alpha=1.0, ylog=False, zorder=None):
"""
| stipple - put markers on a plot to indicate value of interest
|
| f=None - cf field or field
| x=None - x points for field
| y=None - y points for field
| min=None - minimum threshold for stipple
| max=None - maximum threshold for stipple
| size=80 - default size for stipples
| color='k' - default colour for stipples
| pts=50 - number of points in the x direction
| marker='.' - default marker for stipples
| edegecolors='k' - outline colour
| alpha=1.0 - transparency setting - default is off
| ylog=False - set to True if a log pressure stipple plot
| is required
| zorder=None - plotting order
|
|
:Returns:
None
|
|
"""
if plotvars.plot_type not in [1,2,3]:
errstr = '\n stipple error - only X-Y, X-Z and Y-Z \n'
errstr = errstr + 'stipple supported at the present time\n'
errstr = errstr + 'Please raise a feature request if you see this error.\n'
raise Warning(errstr)
# Extract required data for contouring
# If a cf-python field
if isinstance(f, cf.Field):
colorbar_title = ''
field, xpts, ypts, ptype, colorbar_title, xlabel, ylabel, xpole, \
ypole = cf_data_assign(f, colorbar_title)
elif isinstance(f, cf.FieldList):
raise TypeError("Can't plot a field list")
else:
field = f # field data passed in as f
check_data(field, x, y)
xpts = x
ypts = y
if plotvars.plot_type == 1:
# Cylindrical projection
# Add cyclic information if missing.
lonrange = np.nanmax(xpts) - np.nanmin(xpts)
if lonrange < 360:
field, xpts = cartopy_util.add_cyclic_point(field, xpts)
if plotvars.proj == 'cyl':
# Calculate interpolation points
xnew, ynew = stipple_points(xmin=np.nanmin(xpts),
xmax=np.nanmax(xpts),
ymin=np.nanmin(ypts),
ymax=np.nanmax(ypts),
pts=pts, stype=2)
# Calculate points in map space
xnew_map=xnew
ynew_map=ynew
if plotvars.proj == 'npstere' or plotvars.proj == 'spstere':
# Calculate interpolation points
xnew, ynew, xnew_map, ynew_map = polar_regular_grid()
# Convert longitudes to be 0 to 360
# negative longitudes are incorrectly regridded in polar stereographic projection
xnew = np.mod(xnew + 360.0, 360.0)
if plotvars.plot_type >= 2 and plotvars.plot_type <= 3:
# Flip data if a lat-height plot and lats start at the north pole
if plotvars.plot_type == 2:
if xpts[0] > xpts[-1]:
xpts = xpts[::-1]
field = np.fliplr(field)
# Calculate interpolation points
ymin = np.nanmin(ypts)
ymax = np.nanmax(ypts)
if ylog:
ymin = np.log10(ymin)
ymax = np.log10(ymax)
xnew, ynew = stipple_points(xmin=np.nanmin(xpts),
xmax=np.nanmax(xpts),
ymin=ymin,
ymax=ymax,
pts=pts, stype=2)
if ylog:
ynew = 10**ynew
# Get values at the new points
vals = regrid(f=field, x=xpts, y=ypts, xnew=xnew, ynew=ynew)
# Work out which of the points are valid
valid_points = np.array([], dtype='int64')
for i in np.arange(np.size(vals)):
if vals[i] >= min and vals[i] <= max:
valid_points = np.append(valid_points, i)
if plotvars.plot_type == 1:
proj=ccrs.PlateCarree()
if np.size(valid_points) > 0:
plotvars.mymap.scatter(
xnew[valid_points],
ynew[valid_points],
s=size,
c=color,
marker=marker,
edgecolors=edgecolors,
alpha=alpha, transform=proj, zorder=zorder)
if plotvars.plot_type >= 2 and plotvars.plot_type <= 3:
plotvars.plot.scatter(
xnew[valid_points],
ynew[valid_points],
s=size,
c=color,
marker=marker,
edgecolors=edgecolors,
alpha=alpha, zorder=zorder)
[docs]def stipple_points(xmin=None, xmax=None, ymin=None,
ymax=None, pts=None, stype=None):
"""
| stipple_points - calculate interpolation points
|
| xmin=None - plot x minimum
| ymax=None - plot x maximum
| ymin=None - plot y minimum
| ymax=None - plot x maximum
| pts=None - number of points in the x and y directions
| one number gives the same in both directions
|
| stype=None - type of grid. 1=regular, 2=offset
|
|
|
:Returns:
stipple locations in x and y
|
|
"""
# Work out number of points in x and y directions
if np.size(pts) == 1:
pts_x = pts
pts_y = pts
if np.size(pts) == 2:
pts_x = pts[0]
pts_y = pts[1]
# Create regularly spaced points
xstep = (xmax - xmin) / float(pts_x)
x1 = [xmin + xstep / 4]
while (np.nanmax(x1) + xstep) < xmax - xstep / 10:
x1 = np.append(x1, np.nanmax(x1) + xstep)
nxpts = np.size(x1)
x2 = [xmin + xstep * 3 / 4]
while (np.nanmax(x2) + xstep) < xmax - xstep / 10:
x2 = np.append(x2, np.nanmax(x2) + xstep)
ystep = (ymax - ymin) / float(pts_y)
y1 = [ymin + ystep / 2]
while (np.nanmax(y1) + ystep) < ymax - ystep / 10:
y1 = np.append(y1, np.nanmax(y1) + ystep)
# Create interpolation points
xnew = []
ynew = []
iy = 0
for y in y1:
iy = iy + 1
if stype == 1:
xnew = np.append(xnew, x1)
y2 = np.zeros(np.size(x1))
y2.fill(y)
ynew = np.append(ynew, y2)
if stype == 2:
if iy % 2 == 0:
xnew = np.append(xnew, x1)
y2 = np.zeros(np.size(x1))
y2.fill(y)
ynew = np.append(ynew, y2)
if iy % 2 == 1:
xnew = np.append(xnew, x2)
y2 = np.zeros(np.size(x2))
y2.fill(y)
ynew = np.append(ynew, y2)
return xnew, ynew
[docs]def find_pos_in_array(vals=None, val=None, above=False):
"""
| find_pos_in_array - find the position of a point in an array
|
| vals - array values
| val - value to find position of
|
|
|
|
|
|
:Returns:
position in array
|
|
|
"""
pos = -1
if above is False:
for myval in vals:
if val > myval:
pos = pos + 1
if above is 1:
for myval in vals:
if val >= myval:
pos = pos + 1
if np.size(vals) - 1 > pos:
pos = pos + 1
return pos
[docs]def vect(u=None, v=None, x=None, y=None, scale=None, stride=None, pts=None,
key_length=None, key_label=None, ptype=None, title=None, magmin=None,
width=0.02, headwidth=3, headlength=5, headaxislength=4.5,
pivot='middle', key_location=[0.95, -0.06], key_show=True, axes=True,
xaxis=True, yaxis=True, xticks=None, xticklabels=None, yticks=None,
yticklabels=None, xlabel=None, ylabel=None, ylog=False, color='k',
zorder=None):
"""
| vect - plot vectors
|
| u=None - u wind
| v=None - v wind
| x=None - x locations of u and v
| y=None - y locations of u and v
| scale=None - data units per arrow length unit. A smaller values gives
| a larger vector. Generally takes one value but in the case
| of two supplied values the second vector scaling applies to
| the v field.
| stride=None - plot vector every stride points. Can take two values one
| for x and one for y.
| pts=None - use bilinear interpolation to interpolate vectors onto a new
| grid - takes one or two values.
| If one value is passed then this is used for both the x and
| y axes.
| magmin=None - don't plot any vects with less than this magnitude.
| Note this keyword doesn't work in the polar
stereographic projection.
| key_length=None - length of the key. Generally takes one value but in
| the case of two supplied values the second vector
| scaling applies to the v field.
| key_label=None - label for the key. Generally takes one value but in the
| case of two supplied values the second vector scaling
| applies to the v field.
| key_location=[0.9, -0.06] - location of the vector key relative to the
| plot in normalised coordinates.
| key_show=True - draw the key. Set to False if not required.
| ptype=0 - plot type - not needed for cf fields.
| 0 = no specific plot type,
| 1 = longitude-latitude,
| 2 = latitude - height,
| 3 = longitude - height,
| 4 = latitude - time,
| 5 = longitude - time
| 6 = rotated pole
|
| title=None - plot title
| width=0.005 - shaft width in arrow units; default is 0.005 times the
| width of the plot
| headwidth=3 - head width as multiple of shaft width, default is 3
| headlength=5 - head length as multiple of shaft width, default is 5
| headaxislength=4.5 - head length at shaft intersection, default is 4.5
| pivot='middle' - the part of the arrow that is at the grid point; the
| arrow rotates about this point
takes 'tail', 'middle', 'tip'
| axes=True - plot x and y axes
| xaxis=True - plot xaxis
| yaxis=True - plot y axis
| xticks=None - xtick positions
| xticklabels=None - xtick labels
| yticks=None - y tick positions
| yticklabels=None - ytick labels
| xlabel=None - label for x axis
| ylabel=None - label for y axis
| ylog=False - log y axis
| color='k' - colour for the vectors - default is black.
| zorder=None - plotting order
|
:Returns:
None
|
|
|
"""
# If the vector color is white set the quicker key colour to black
# so that it can be seen
qkey_color=color
if qkey_color == 'w' or qkey_color == 'white':
qkey_color = 'k'
colorbar_title = ''
text_fontsize = plotvars.text_fontsize
continent_thickness = plotvars.continent_thickness
continent_color = plotvars.continent_color
if text_fontsize is None:
text_fontsize = 11
if continent_thickness is None:
continent_thickness = 1.5
if continent_color is None:
continent_color = 'k'
# ylog=plotvars.ylog
title_fontsize = plotvars.title_fontsize
title_fontweight = plotvars.title_fontweight
if title_fontsize is None:
title_fontsize = 15
resolution_orig=plotvars.resolution
# Set potential user axis labels
user_xlabel = xlabel
user_ylabel = ylabel
rotated_vect=False
if u.ref('rotated_latitude_longitude') is not None:
rotated_vect = True
# Extract required data
# If a cf-python field
if isinstance(u, cf.Field):
u_data, u_x, u_y, ptype, colorbar_title, xlabel, ylabel, xpole, \
ypole = cf_data_assign(u, colorbar_title, rotated_vect=rotated_vect)
elif isinstance(u, cf.FieldList):
raise TypeError("Can't plot a field list")
else:
# field=f #field data passed in as f
check_data(u, x, y)
u_data = deepcopy(u)
u_x = deepcopy(x)
u_y = deepcopy(y)
xlabel = ''
ylabel = ''
if isinstance(v, cf.Field):
v_data, v_x, v_y, ptype, colorbar_title, xlabel, ylabel, xpole, \
ypole = cf_data_assign(v, colorbar_title, rotated_vect=rotated_vect)
elif isinstance(v, cf.FieldList):
raise TypeError("Can't plot a field list")
else:
# field=f #field data passed in as f
check_data(v, x, y)
v_data = deepcopy(v)
v_x = deepcopy(x)
v_y = deepcopy(y)
xlabel = ''
ylabel = ''
# If a minimum magnitude is specified mask these data points
if magmin is not None:
mag = np.sqrt(u_data**2 + v_data**2)
invalid = np.where(mag <= magmin)
if np.size(invalid) > 0:
u_data[invalid] = np.nan
v_data[invalid] = np.nan
# Reset xlabel and ylabel values with user defined labels in specified
if user_xlabel is not None:
xlabel = user_xlabel
if user_ylabel is not None:
ylabel = user_ylabel
# Retrieve any user defined axis labels
if xlabel == '' and plotvars.xlabel is not None:
xlabel = plotvars.xlabel
if ylabel == '' and plotvars.ylabel is not None:
ylabel = plotvars.ylabel
if xticks is None and plotvars.xticks is not None:
xticks = plotvars.xticks
if plotvars.xticklabels is not None:
xticklabels = plotvars.xticklabels
else:
xticklabels = map(str, xticks)
if yticks is None and plotvars.yticks is not None:
yticks = plotvars.yticks
if plotvars.yticklabels is not None:
yticklabels = plotvars.yticklabels
else:
yticklabels = map(str, yticks)
if scale is None:
scale = np.nanmax(u_data) / 4.0
if key_length is None:
key_length = scale
# Open a new plot if necessary
if plotvars.user_plot == 0:
gopen(user_plot=0)
# Set plot type if user specified
if (ptype is not None):
plotvars.plot_type = ptype
lonrange = np.nanmax(u_x) - np.nanmin(u_x)
latrange = np.nanmax(u_y) - np.nanmin(u_y)
if plotvars.plot_type == 1:
# Set up mapping
if (lonrange > 350 and latrange > 170) or plotvars.user_mapset == 1:
set_map()
else:
mapset(lonmin=np.nanmin(u_x), lonmax=np.nanmax(u_x),
latmin=np.nanmin(u_y), latmax=np.nanmax(u_y),
user_mapset=0, resolution=resolution_orig)
set_map()
mymap = plotvars.mymap
u_data, u_x = cartopy_util.add_cyclic_point(u_data, u_x)
v_data, v_x = cartopy_util.add_cyclic_point(v_data, v_x)
# stride data points to reduce vector density
if stride is not None:
if np.size(stride) == 1:
xstride = stride
ystride = stride
if np.size(stride) == 2:
xstride = stride[0]
ystride = stride[1]
u_x = u_x[0::xstride]
u_y = u_y[0::ystride]
u_data = u_data[0::ystride, 0::xstride]
v_data = v_data[0::ystride, 0::xstride]
# Map vectors
if plotvars.plot_type == 1:
lonmin=plotvars.lonmin
lonmax=plotvars.lonmax
proj=ccrs.PlateCarree()
if pts is None:
quiv = plotvars.mymap.quiver(u_x, u_y, u_data, v_data, scale=scale,
pivot=pivot, units='inches',
width=width, headwidth=headwidth,
headlength=headlength,
headaxislength=headaxislength,
color=color, transform=proj,
zorder=zorder)
else:
if plotvars.proj == 'cyl':
# **cartopy 0.16 fix for longitide points in cylindrical projection
# when regridding to a number of points
# Make points within the plotting region
for pt in np.arange(np.size(u_x)):
if u_x[pt] > lonmax:
u_x[pt]=u_x[pt]-360
quiv = plotvars.mymap.quiver(u_x, u_y, u_data, v_data, scale=scale,
pivot=pivot, units='inches',
width=width, headwidth=headwidth,
headlength=headlength,
headaxislength=headaxislength,
color=color,
regrid_shape=pts, transform=proj,
zorder=zorder)
# Make key_label if none exists
if key_label is None:
key_label = str(key_length)
if isinstance(u, cf.Field):
key_label = supscr(key_label + u.units)
if key_show:
quiv_key = plotvars.mymap.quiverkey(quiv, key_location[0],
key_location[1],
key_length,
key_label, labelpos='W',
color=qkey_color,
fontproperties={'size':str(plotvars.axis_label_fontsize)},
coordinates='axes')
# axes
plot_map_axes(axes=axes, xaxis=xaxis, yaxis=yaxis,
xticks=xticks, xticklabels=xticklabels,
yticks=yticks, yticklabels=yticklabels,
user_xlabel=user_xlabel, user_ylabel=user_ylabel,
verbose=False)
# Coastlines
continent_thickness = plotvars.continent_thickness
continent_color = plotvars.continent_color
continent_linestyle = plotvars.continent_linestyle
if continent_thickness is None:
continent_thickness = 1.5
if continent_color is None:
continent_color = 'k'
if continent_linestyle is None:
continent_linestyle = 'solid'
feature = cfeature.NaturalEarthFeature(
name='land', category='physical',
scale=plotvars.resolution,
facecolor='none')
mymap.add_feature(feature, edgecolor=continent_color,
linewidth=continent_thickness,
linestyle=continent_linestyle)
# Title
if title is not None:
map_title(title)
if plotvars.plot_type == 6:
if u.ref('rotated_latitude_longitude') is not None:
#rotated_pole = u.ref('rotated_latitude_longitude')
#xpole = float(rotated_pole['grid_north_pole_longitude'])
#ypole = float(rotated_pole['grid_north_pole_latitude'])
proj=ccrs.PlateCarree()
# Set up mapping
if (lonrange > 350 and latrange > 170) or plotvars.user_mapset == 1:
set_map()
else:
mapset(lonmin=np.nanmin(u_x), lonmax=np.nanmax(u_x),
latmin=np.nanmin(u_y), latmax=np.nanmax(u_y),
user_mapset=0, resolution=resolution_orig)
set_map()
quiv = plotvars.mymap.quiver(u_x, u_y, u_data, v_data, scale=scale*10, transform=proj,
pivot=pivot, units='inches',
width=width, headwidth=headwidth,
headlength=headlength,
headaxislength=headaxislength,
color=color, zorder=zorder)
# Make key_label if none exists
if key_label is None:
key_label = str(key_length)
if isinstance(u, cf.Field):
key_label = supscr(key_label + u.units)
if key_show:
quiv_key = plotvars.mymap.quiverkey(quiv, key_location[0],
key_location[1],
key_length,
key_label, labelpos='W',
color=qkey_color,
fontproperties={'size':str(plotvars.axis_label_fontsize)},
coordinates='axes')
# Axes on the native grid
if plotvars.plot == 'rotated':
rgaxes(xpole=xpole, ypole=ypole, xvec=x, yvec=y,
xticks=xticks, xticklabels=xticklabels,
yticks=yticks, yticklabels=yticklabels,
axes=axes, xaxis=xaxis, yaxis=yaxis,
xlabel=xlabel, ylabel=ylabel)
if plotvars.plot == 'cyl':
plot_map_axes(axes=axes, xaxis=xaxis, yaxis=yaxis,
xticks=xticks, xticklabels=xticklabels,
yticks=yticks, yticklabels=yticklabels,
user_xlabel=user_xlabel, user_ylabel=user_ylabel,
verbose=False)
# Title
if title is not None:
map_title(title)
######################################
# Latitude or longitude vs height plot
######################################
if plotvars.plot_type == 2 or plotvars.plot_type == 3:
user_gset = plotvars.user_gset
if user_gset == 0:
# Program selected data plot limits
xmin = np.nanmin(u_x)
xmax = np.nanmax(u_x)
if plotvars.plot_type == 2:
if xmin < -80 and xmin >= -90:
xmin = -90
if xmax > 80 and xmax <= 90:
xmax = 90
ymin = np.nanmin(u_y)
if ymin <= 10:
ymin = 0
ymax = np.nanmax(u_y)
else:
# User specified plot limits
xmin = plotvars.xmin
xmax = plotvars.xmax
if plotvars.ymin < plotvars.ymax:
ymin = plotvars.ymin
ymax = plotvars.ymax
else:
ymin = plotvars.ymax
ymax = plotvars.ymin
xstep = None
if ptype == 1:
if (xmin == -90 and xmax == 90):
xstep = 30
if ptype == 2:
if xmax - xmax >= 160:
xstep = 60
else:
xstep = 30
ystep = None
if (ymax == 1000):
ystep = 100
if (ymax == 100000):
ystep = 10000
ytype = 0 # pressure or similar y axis
if 'theta' in ylabel.split(' '):
ytype = 1
if 'height' in ylabel.split(' '):
ytype = 1
ystep = 100
if (ymax - ymin) > 5000:
ystep = 500.0
if (ymax - ymin) > 10000:
ystep = 1000.0
if (ymax - ymin) > 50000:
ystep = 10000.0
# Set plot limits and draw axes
if ylog != 1:
if ytype == 1:
gset(
xmin=xmin,
xmax=xmax,
ymin=ymin,
ymax=ymax,
user_gset=user_gset)
else:
gset(
xmin=xmin,
xmax=xmax,
ymin=ymax,
ymax=ymin,
user_gset=user_gset)
# Set default x-axis labels
lltype=1
if plotvars.plot_type == 2:
lltype = 2
llticks, lllabels = mapaxis(min=xmin, max=xmax, type=lltype)
heightticks = gvals(
dmin=ymin,
dmax=ymax,
tight=1,
mystep=ystep,
mod=0)[0]
heightlabels = heightticks
if axes:
if xaxis:
if xticks is not None:
llticks = xticks
lllabels = xticks
if xticklabels is not None:
lllabels = xticklabels
else:
llticks = [100000000]
xlabel = ''
if yaxis:
if yticks is not None:
heightticks = yticks
heightlabels = yticks
if yticklabels is not None:
heightlabels = yticklabels
else:
heightticks = [100000000]
ylabel = ''
else:
llticks = [100000000]
heightticks = [100000000]
xlabel = ''
ylabel = ''
axes_plot(xticks=llticks, xticklabels=lllabels,
yticks=heightticks, yticklabels=heightlabels,
xlabel=xlabel, ylabel=ylabel)
# Log y axis
if ylog:
if ymin == 0:
ymin = 1 # reset zero mb/height input to a small value
gset(
xmin=xmin,
xmax=xmax,
ymin=ymax,
ymax=ymin,
ylog=1,
user_gset=user_gset)
llticks, lllabels = mapaxis(min=xmin, max=xmax,
type=plotvars.plot_type)
if axes:
if xaxis:
if xticks is not None:
llticks = xticks
lllabels = xticks
if xticklabels is not None:
lllabels = xticklabels
else:
llticks = [100000000]
xlabel = ''
if yaxis:
if yticks is not None:
heightticks = yticks
heightlabels = yticks
if yticklabels is not None:
heightlabels = yticklabels
else:
heightticks = [100000000]
ylabel = ''
if yticks is None:
axes_plot(
xticks=llticks,
xticklabels=lllabels,
xlabel=xlabel,
ylabel=ylabel)
else:
axes_plot(xticks=llticks, xticklabels=lllabels,
yticks=heightticks, yticklabels=heightlabels,
xlabel=xlabel, ylabel=ylabel)
# Regrid the data if requested
if pts is not None:
xnew, ynew = stipple_points(xmin=np.min(u_x), xmax=np.max(u_x),
ymin=np.min(u_y), ymax=np.max(u_y),
pts=pts, stype=1)
if ytype == 0:
# Make y interpolation in log space as we have a pressure coordinate
u_vals = regrid(f=u_data, x=u_x, y=np.log10(u_y), xnew=xnew, ynew=np.log10(ynew))
v_vals = regrid(f=v_data, x=u_x, y=np.log10(u_y), xnew=xnew, ynew=np.log10(ynew))
else:
u_vals = regrid(f=u_data, x=u_x, y=u_y, xnew=xnew, ynew=ynew)
v_vals = regrid(f=v_data, x=u_x, y=u_y, xnew=xnew, ynew=ynew)
u_x = xnew
u_y = ynew
u_data = u_vals
v_data = v_vals
# set scale and key lengths
if np.size(scale) == 1:
scale_u = scale
scale_v = scale
else:
scale_u = scale[0]
scale_v = scale[1]
if np.size(key_length) == 2:
key_length_u = key_length[0]
key_length_v = key_length[1]
# scale v data
v_data = v_data * scale_u / scale_v
else:
key_length_u = key_length
# Plot the vectors
quiv = plotvars.plot.quiver(u_x, u_y, u_data, v_data, pivot=pivot,
units='inches', scale=scale_u,
width=width, headwidth=headwidth,
headlength=headlength,
headaxislength=headaxislength,
color=color, zorder=zorder)
# Plot single key
if np.size(scale) == 1:
# Single scale vector
if key_label is None:
key_label_u = str(key_length_u)
if isinstance(u, cf.Field):
key_label_u = supscr(key_label_u + ' (' + u.units + ')')
else:
key_label_u = key_label[0]
if key_show:
quiv_key = plotvars.plot.quiverkey(quiv, key_location[0],
key_location[1],
key_length_u, key_label_u,
labelpos='W',
color=qkey_color,
fontproperties={'size':str(plotvars.axis_label_fontsize)})
# Plot two keys
if np.size(scale) == 2:
# translate from normalised units to plot units
xpos = key_location[0] * \
(plotvars.xmax - plotvars.xmin) + plotvars.xmin
ypos = key_location[1] * \
(plotvars.ymax - plotvars.ymin) + plotvars.ymin
# horizontal and vertical spacings for offsetting vector reference
# text
xoffset = 0.01 * abs(plotvars.xmax - plotvars.xmin)
yoffset = 0.01 * abs(plotvars.ymax - plotvars.ymin)
# Assign key labels if necessary
if key_label is None:
key_label_u = str(key_length_u)
key_label_v = str(key_length_v)
if isinstance(u, cf.Field):
key_label_u = supscr(key_label_u + ' (' + u.units + ')')
if isinstance(v, cf.Field):
key_label_v = supscr(key_label_v + ' (' + v.units + ')')
else:
key_label_u = supscr(key_label[0])
key_label_v = supscr(key_label[1])
# Plot reference vectors and keys
if key_show:
quiv1 = plotvars.plot.quiver(xpos, ypos, key_length[0], 0,
pivot='tail', units='inches',
scale=scale[0],
headaxislength=headaxislength,
width=width, headwidth=headwidth,
headlength=headlength,
clip_on=False,
color=qkey_color)
quiv2 = plotvars.plot.quiver(xpos, ypos, 0, key_length[1],
pivot='tail', units='inches',
scale=scale[1],
headaxislength=headaxislength,
width=width, headwidth=headwidth,
headlength=headlength,
clip_on=False,
color=qkey_color)
plotvars.plot.text(
xpos,
ypos + yoffset,
key_label_u,
horizontalalignment='left',
verticalalignment='top')
plotvars.plot.text(
xpos - xoffset,
ypos,
key_label_v,
horizontalalignment='right',
verticalalignment='bottom')
if title is not None:
plotvars.plot.set_title(
title,
y=1.03,
fontsize=plotvars.title_fontsize,
fontweight=title_fontweight)
##########
# Save plot
##########
if plotvars.user_plot == 0:
gset()
cscale()
gclose()
if plotvars.user_mapset == 0:
mapset()
mapset(resolution=resolution_orig)
[docs]def set_map():
"""
| set_map - set map and write into plotvars.mymap
|
| No inputs
| This is an internal routine and not used by the user
|
|
|
|
|
:Returns:
None
|
|
|
"""
# Set up mapping
extent = True
lon_mid = plotvars.lonmin + (plotvars.lonmax - plotvars.lonmin) / 2.0
lat_mid = plotvars.latmin + (plotvars.latmax - plotvars.latmin) / 2.0
lonmin = plotvars.lonmin
lonmax = plotvars.lonmax
latmin = plotvars.latmin
latmax = plotvars.latmax
if plotvars.proj == 'cyl':
proj=ccrs.PlateCarree(central_longitude=lon_mid)
# Cartopy line plotting and identical left == right fix
if lonmax - lonmin == 360.0:
lonmax = lonmax + 0.01
if plotvars.proj == 'merc':
min_latitude=-80.0
if plotvars.lonmin > min_latitude:
min_latitude=plotvars.lonmin
max_latitude=84.0
if plotvars.lonmax < max_latitude:
max_latitude=plotvars.lonmax
proj=ccrs.Mercator(central_longitude=plotvars.lon_0,
min_latitude=min_latitude,
max_latitude=max_latitude)
if plotvars.proj == 'npstere':
proj = ccrs.NorthPolarStereo(central_longitude=plotvars.lon_0)
#mymap=plot.subplot(plotvars.rows, plotvars.columns, plotvars.pos, projection=proj)
# **cartopy 0.16 fix
# Here we add in 0.01 to the longitude extent as this helps with plotting
# lines and line labels
lonmin = plotvars.lon_0-180
lonmax = plotvars.lon_0+180.01
latmin = plotvars.boundinglat
latmax = 90
if plotvars.proj == 'spstere':
proj = ccrs.SouthPolarStereo(central_longitude=plotvars.lon_0)
#mymap=plot.subplot(plotvars.rows, plotvars.columns, plotvars.pos, projection=proj)
# **cartopy 0.16 fix
# Here we add in 0.01 to the longitude extent as this helps with plotting
# lines and line labels
lonmin = plotvars.lon_0-180
lonmax = plotvars.lonmax+180.01
latmin = -90
latmax = plotvars.boundinglat
if plotvars.proj == 'ortho':
proj=ccrs.Orthographic(central_longitude=plotvars.lon_0,
central_latitude=plotvars.lat_0)
lonmin = plotvars.lon_0-180.0
lonmax = plotvars.lon_0+180.01
extent = False
if plotvars.proj == 'moll':
proj=ccrs.Mollweide(central_longitude=plotvars.lon_0)
lonmin = plotvars.lon_0-180.0
lonmax = plotvars.lon_0+180.01
extent = False
if plotvars.proj == 'robin':
proj=ccrs.Robinson(central_longitude=plotvars.lon_0)
if plotvars.proj == 'lcc':
lon_0=plotvars.lonmin+(plotvars.lonmax-plotvars.lonmin)/2.0
lat_0=plotvars.latmin+(plotvars.latmax-plotvars.latmin)/2.0
proj=ccrs.LambertConformal(central_longitude=lon_0,
central_latitude=lat_0,
cutoff=-40)
if plotvars.proj == 'rotated':
proj=ccrs.PlateCarree(central_longitude=lon_mid)
if plotvars.proj == 'OSGB':
proj=ccrs.OSGB()
if plotvars.proj == 'EuroPP':
proj=ccrs.EuroPP()
if plotvars.proj == 'UKCP':
# Special case of TransverseMercator for UKCP
proj=ccrs.TransverseMercator()
if plotvars.proj == 'TransverseMercator':
proj=ccrs.TransverseMercator()
lonmin = plotvars.lon_0-180.0
lonmax = plotvars.lon_0+180.01
extent = False
# Add a plot containing the projection
if plotvars.plot_xmin:
delta_x = plotvars.plot_xmax - plotvars.plot_xmin
delta_y = plotvars.plot_ymax - plotvars.plot_ymin
mymap = plotvars.master_plot.add_axes([plotvars.plot_xmin,
plotvars.plot_ymin,
delta_x, delta_y],
projection=proj)
else:
mymap=plotvars.master_plot.add_subplot(plotvars.rows,
plotvars.columns,
plotvars.pos,
projection=proj)
# Set map extent
set_extent = True
if plotvars.proj == 'OSGB' or plotvars.proj == 'EuroPP' or plotvars.proj == 'UKCP':
set_extent = False
if extent and set_extent:
mymap.set_extent((lonmin, lonmax, latmin, latmax), crs=ccrs.PlateCarree())
if plotvars.proj == 'UKCP':
# Special case of TransverseMercator for UKCP
mymap.set_extent((-11, 3, 49, 61), crs=ccrs.PlateCarree())
if plotvars.proj == 'EuroPP':
# EuroPP somehow needs some limits setting.
mymap.set_extent((-12, 25, 30, 75), crs=ccrs.PlateCarree())
# Remove any plotvars.plot axes leaving just the plotvars.mymap axes
#if plotvars.plot is not None:
plotvars.plot.set_frame_on(False)
plotvars.plot.set_xticks([])
plotvars.plot.set_yticks([])
# Store map
plotvars.mymap = mymap
[docs]def polar_regular_grid(pts=50):
"""
| polar_regular_grid - return a regular grid over a polar
| stereographic area
|
| pts=50 - number of grid points in the x and y directions
|
|
|
|
|
|
:Returns:
lons, lats of grid in degrees
x, y locations of lons and lats
|
|
|
"""
mymap = plotvars.mymap
boundinglat = plotvars.boundinglat
lon_0 = plotvars.lon_0
if plotvars.proj == 'npstere':
thisproj=ccrs.NorthPolarStereo(central_longitude=lon_0)
else:
thisproj=ccrs.SouthPolarStereo(central_longitude=lon_0)
# Find min and max of plotting region in device coordinates
lons=np.array([lon_0-90, lon_0, lon_0+90, lon_0+180])
lats=np.array([boundinglat, boundinglat, boundinglat, boundinglat])
extent=thisproj.transform_points(ccrs.PlateCarree(), lons, lats)
xmin=np.min(extent[:,0])
xmax=np.max(extent[:,0])
ymin=np.min(extent[:,1])
ymax=np.max(extent[:,1])
# Make up a stipple of points for cover the pole
points_device = stipple_points(
xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, pts=pts, stype=2)
xnew = np.array(points_device)[0, :]
ynew = np.array(points_device)[1, :]
points_polar=ccrs.PlateCarree().transform_points(thisproj, xnew, ynew)
lons=np.array(points_polar)[:, 0]
lats=np.array(points_polar)[:, 1]
if plotvars.proj == 'npstere':
valid = np.where(lats >= boundinglat)
else:
valid = np.where(lats <= boundinglat)
return lons[valid], lats[valid], xnew[valid], ynew[valid]
[docs]def cf_var_name(field=None, dim=None):
"""
| cf_var_name - return the name from a supplied dimension
| in the following order
| ncvar
| short_name
| long_name
| standard_name
|
| field=None - field
| dim=None - dimension required - 'dim0', 'dim1' etc.
|
|
|
|
|
:Returns:
name
|
|
|
"""
id = getattr(field.item(dim), 'id', False)
ncvar = getattr(field.item(dim), 'ncvar', False)
short_name = getattr(field.item(dim), 'short_name', False)
long_name = getattr(field.item(dim), 'long_name', False)
standard_name = getattr(field.item(dim), 'standard_name', False)
name = 'No Name'
if id:
name = id
if ncvar:
name = ncvar
if short_name:
name = short_name
if long_name:
name = long_name
if standard_name:
name = standard_name
return name
[docs]def process_color_scales():
"""
| Process colour scales to generate images of them for the web
| documentation and the rst code for inclusion in the
| colour_scale.rst file.
|
|
| No inputs
| This is an internal routine and not used by the user
|
|
|
|
|
:Returns:
None
|
|
|
"""
# Define scale categories
uniform = ['viridis', 'magma', 'inferno', 'plasma', 'parula', 'gray']
ncl_large = ['amwg256', 'BkBlAqGrYeOrReViWh200', 'BlAqGrYeOrRe',
'BlAqGrYeOrReVi200', 'BlGrYeOrReVi200', 'BlRe', 'BlueRed',
'BlueRedGray', 'BlueWhiteOrangeRed', 'BlueYellowRed',
'BlWhRe', 'cmp_b2r', 'cmp_haxby', 'detail', 'extrema',
'GrayWhiteGray', 'GreenYellow', 'helix', 'helix1',
'hotres', 'matlab_hot', 'matlab_hsv', 'matlab_jet',
'matlab_lines', 'ncl_default', 'ncview_default',
'OceanLakeLandSnow', 'rainbow', 'rainbow_white_gray',
'rainbow_white', 'rainbow_gray', 'tbr_240_300',
'tbr_stdev_0_30', 'tbr_var_0_500', 'tbrAvg1', 'tbrStd1',
'tbrVar1', 'thelix', 'ViBlGrWhYeOrRe', 'wh_bl_gr_ye_re',
'WhBlGrYeRe', 'WhBlReWh', 'WhiteBlue',
'WhiteBlueGreenYellowRed', 'WhiteGreen',
'WhiteYellowOrangeRed', 'WhViBlGrYeOrRe', 'WhViBlGrYeOrReWh',
'wxpEnIR', '3gauss', '3saw', 'BrBG']
ncl_meteoswiss = ['hotcold_18lev', 'hotcolr_19lev', 'mch_default',
'perc2_9lev', 'percent_11lev', 'precip2_15lev',
'precip2_17lev', 'precip3_16lev', 'precip4_11lev',
'precip4_diff_19lev', 'precip_11lev',
'precip_diff_12lev', 'precip_diff_1lev',
'rh_19lev', 'spread_15lev']
ncl_color_blindness = ['StepSeq25', 'posneg_2', 'posneg_1',
'BlueDarkOrange18', 'BlueDarkRed18',
'GreenMagenta16', 'BlueGreen14', 'BrownBlue12',
'Cat12']
ncl_small = ['amwg', 'amwg_blueyellowred', 'BlueDarkRed18',
'BlueDarkOrange18', 'BlueGreen14', 'BrownBlue12', 'Cat12',
'cmp_flux', 'cosam12', 'cosam', 'GHRSST_anomaly',
'GreenMagenta16', 'hotcold_18lev', 'hotcolr_19lev',
'mch_default', 'nrl_sirkes', 'nrl_sirkes_nowhite',
'perc2_9lev', 'percent_11lev', 'posneg_2', 'prcp_1', 'prcp_2',
'prcp_3', 'precip_11lev', 'precip_diff_12lev',
'precip_diff_1lev', 'precip2_15lev', 'precip2_17lev',
'precip3_16lev', 'precip4_11lev', 'precip4_diff_19lev',
'radar', 'radar_1', 'rh_19lev', 'seaice_1', 'seaice_2',
'so4_21', 'spread_15lev', 'StepSeq25', 'sunshine_9lev',
'sunshine_diff_12lev', 'temp_19lev', 'temp_diff_18lev',
'temp_diff_1lev', 'topo_15lev', 'wgne15', 'wind_17lev']
orography = ['os250kmetres', 'wiki_1_0_2', 'wiki_1_0_3',
'wiki_2_0', 'wiki_2_0_reduced', 'arctic']
idl_guide = []
for i in np.arange(1, 45):
idl_guide.append('scale' + str(i))
for category in ['uniform', 'ncl_meteoswiss', 'ncl_small', 'ncl_large',
'ncl_color_blindness', 'orography', 'idl_guide']:
if category == 'uniform':
scales = uniform
div = '================== ====='
chars = 10
title = 'Perceptually uniform colour maps for use with continuous '
title += 'data'
print title
print '----------------------------------------------'
print ''
print div
print 'Name Scale'
print div
if category == 'ncl_meteoswiss':
scales = ncl_meteoswiss
div = '================== ====='
chars = 19
print 'NCAR Command Language - MeteoSwiss colour maps'
print '----------------------------------------------'
print ''
print div
print 'Name Scale'
print div
if category == 'ncl_small':
scales = ncl_small
div = '=================== ====='
chars = 20
print 'NCAR Command Language - small color maps (<50 colours)'
print '------------------------------------------------------'
print ''
print div
print 'Name Scale'
print div
if category == 'ncl_large':
scales = ncl_large
div = '======================= ====='
chars = 24
print 'NCAR Command Language - large colour maps (>50 colours)'
print '-------------------------------------------------------'
print ''
print div
print 'Name Scale'
print div
if category == 'ncl_color_blindness':
scales = ncl_color_blindness
div = '================ ====='
chars = 17
title = 'NCAR Command Language - Enhanced to help with colour'
title += 'blindness'
print title
title = '-----------------------------------------------------'
title += '---------'
print title
print ''
print div
print 'Name Scale'
print div
chars = 17
if category == 'orography':
scales = orography
div = '================ ====='
chars = 17
print 'Orography/bathymetry colour scales'
print '----------------------------------'
print ''
print div
print 'Name Scale'
print div
chars = 17
if category == 'idl_guide':
scales = idl_guide
div = '======= ====='
chars = 8
print 'IDL guide scales'
print '----------------'
print ''
print div
print 'Name Scale'
print div
chars = 8
for scale in scales:
# Make image of scale
fig = plot.figure(figsize=(8, 0.5))
ax1 = fig.add_axes([0.05, 0.1, 0.9, 0.2])
cscale(scale)
ncols = np.size(plotvars.cs)
cmap = matplotlib.colors.ListedColormap(plotvars.cs)
cb1 = matplotlib.colorbar.ColorbarBase(
ax1, cmap=cmap, orientation='horizontal', ticks=None)
cb1.set_ticks([0.0, 1.0])
cb1.set_ticklabels(['', ''])
file = '/home/andy/cf-docs/cfplot_sphinx/images/'
file += 'colour_scales/' + scale + '.png'
plot.savefig(file)
plot.close()
# Use covert to trim the png file to remove white space
subprocess.call(["convert", "-trim", file, file])
name_pad = scale
while len(name_pad) < chars:
name_pad = name_pad + ' '
fn = name_pad + '.. image:: images/colour_scales/' + scale + '.png'
print fn
print div
print ''
print ''
[docs]def reset():
"""
| reset all plotting variables
|
|
|
|
|
|
|
:Returns:
name
|
|
|
"""
axes()
cscale()
levs()
gset()
mapset()
setvars()
[docs]def setvars(file=None, title_fontsize=None, text_fontsize=None,
colorbar_fontsize=None, colorbar_fontweight=None,
axis_label_fontsize=None, title_fontweight=None,
text_fontweight=None, axis_label_fontweight=None, fontweight=None,
continent_thickness=None, continent_color=None,
continent_linestyle=None, viewer=None,
tspace_year=None, tspace_month=None, tspace_day=None,
tspace_hour=None, xtick_label_rotation=None,
xtick_label_align=None, ytick_label_rotation=None,
ytick_label_align=None, legend_text_weight=None,
legend_text_size=None, cs_uniform=None,
master_title=None, master_title_location=None,
master_title_fontsize=None, master_title_fontweight=None,
dpi=None, land_color=None, ocean_color=None,
lake_color=None,
rotated_grid_spacing=None, rotated_deg_spacing=None,
rotated_continents=None, rotated_grid=None,
rotated_labels=None, rotated_grid_thickness=None,
legend_frame=None,
legend_frame_edge_color=None, legend_frame_face_color=None,
degsym=None, axis_width=None, grid=None,
grid_spacing=None, grid_lons=None, grid_lats=None,
grid_colour=None, grid_linestyle=None, grid_thickness=None):
"""
| setvars - set plotting variables and their defaults
|
| file=None - output file name
| title_fontsize=None - title fontsize, default=15
| title_fontweight='normal' - title fontweight
| text_fontsize='normal' - text font size, default=11
| text_fontweight='normal' - text font weight
| axis_label_fontsize=None - axis label fontsize, default=11
| axis_label_fontweight='normal' - axis font weight
| legend_text_size='11' - legend text size
| legend_text_weight='normal' - legend text weight
| colorbar_textsize='11' - colorbar text size
| colorbar_fontweight='normal' - colorbar font weight
| legend_text_weight='normal' - legend text weight
| master_title_fontsize=30 - master title font size
| master_title_fontweight='normal' - master title font weight
| continent_thickness=1.5 - default=1.5
| continent_color='k' - default='k' (black)
| continent_linestyle='solid' - default='k' (black)
| viewer='display' - use ImageMagick display program to display
| the pictures. Set to 'matplotlib' to use the
| built in matplotlib viewer. display is non-blocking
| of the command prompt while the built in matplotlib
| viewer is blocking. Set to anything else to not
| view the picture.
| tspace_year=None - time axis spacing in years
| tspace_month=None - time axis spacing in months
| tspace_day=None - time axis spacing in days
| tspace_hour=None - time axis spacing in hours
| xtick_label_rotation=0 - rotation of xtick labels
| xtick_label_align='center' - alignment of xtick labels
| ytick_label_rotation=0 - rotation of ytick labels
| ytick_label_align='right' - alignment of ytick labels
| cs_uniform=True - make a uniform differential colour scale
| master_title=None - master title text
| master_title_location=[0.5,0.95] - master title location
| dpi=None - dots per inch setting
| land_color=None - land colour
| ocean_color=None - ocean colour
| lake_color=None - lake colour
| rotated_grid_spacing=10 - rotated grid spacing in degrees
| rotated_deg_spacing=0.75 - rotated grid spacing between graticule dots
| rotated_deg_tkickness=1.0 - rotated grid thickness for longitude and latitude lines
| rotated_continents=True - draw rotated continents
| rotated_grid=True - draw rotated grid
| rotated_labels=True - draw rotated grid labels
| legend_frame=True - draw a frame around a lineplot legend
| legend_frame_edge_color='k' - color for the legend frame
| legend_frame_face_color=None - color for the legend background
| degsym=True - add degree symbol to longitude and latitude axis labels
| axis_width=None - width of line for the axes
| grid=True - draw grid
| grid_spacing=1 - grid spacing in degrees
| grid_lons=None - grid longitudes
| grid_lats=None - grid latitudes
| grid_colour='grey' - grid colour
| grid_linestyle='--' - grid line style
| grid_thickness=1.0 - grid thickness
|
| Use setvars() to reset to the defaults
|
|
|
:Returns:
name
|
|
|
"""
vals = [file, title_fontsize, text_fontsize, axis_label_fontsize,
continent_thickness, title_fontweight, text_fontweight,
axis_label_fontweight, fontweight, continent_color,
continent_linestyle, tspace_year,
tspace_month, tspace_day, tspace_hour, xtick_label_rotation,
xtick_label_align, ytick_label_rotation, ytick_label_align,
legend_text_size, legend_text_weight, cs_uniform,
master_title, master_title_location,
master_title_fontsize, master_title_fontweight, dpi,
land_color, ocean_color, lake_color, rotated_grid_spacing,
rotated_deg_spacing, rotated_continents, rotated_grid,
rotated_grid_thickness,
rotated_labels, colorbar_fontsize, colorbar_fontweight,
legend_frame, legend_frame_edge_color, legend_frame_face_color,
degsym, axis_width, grid, grid_spacing, grid_lons, grid_lats,
grid_colour, grid_linestyle, grid_thickness]
if all(val is None for val in vals):
plotvars.file = None
plotvars.title_fontsize = 15
plotvars.text_fontsize = 11
plotvars.colorbar_fontsize = 11
plotvars.axis_label_fontsize = 11
plotvars.title_fontweight = 'normal'
plotvars.text_fontweight = 'normal'
plotvars.colorbar_fontweight = 'normal'
plotvars.axis_label_fontweight = 'normal'
plotvars.fontweight = 'normal'
plotvars.continent_thickness = None
plotvars.continent_color = None
plotvars.continent_linestyle = None
plotvars.tspace_year = None
plotvars.tspace_month = None
plotvars.tspace_day = None
plotvars.tspace_hour = None
plotvars.xtick_label_rotation = 0
plotvars.xtick_label_align = 'center'
plotvars.ytick_label_rotation = 0
plotvars.ytick_label_align = 'right'
plotvars.legend_text_size = 11
plotvars.legend_text_weight = 'normal'
plotvars.cs_uniform = True
plotvars.viewer = 'display'
plotvars.master_title = None
plotvars.master_title_location = [0.5, 0.95]
plotvars.master_title_fontsize = 30
plotvars.master_title_fontweight = 'normal'
plotvars.dpi=None
plotvars.land_color=None
plotvars.ocean_color=None
plotvars.lake_color=None
plotvars.rotated_grid_spacing=10
plotvars.rotated_deg_spacing=0.75
plotvars.rotated_grid_thickness=1.0
plotvars.rotated_continents=True
plotvars.rotated_grid=True
plotvars.rotated_labels=True
plotvars.legend_frame=True
plotvars.legend_frame_edge_color='k'
plotvars.legend_frame_face_color=None
plotvars.degsym=False
plotvars.axis_width=None
plotvars.grid=True
plotvars.grid_spacing=1
plotvars.grid_lons=None
plotvars.grid_lats=None
plotvars.grid_colour='grey'
plotvars.grid_linestyle='--'
plotvars.grid_thickness=1.0
if file is not None:
plotvars.file = file
if title_fontsize is not None:
plotvars.title_fontsize = title_fontsize
if axis_label_fontsize is not None:
plotvars.axis_label_fontsize = axis_label_fontsize
if continent_thickness is not None:
plotvars.continent_thickness = continent_thickness
if continent_color is not None:
plotvars.continent_color = continent_color
if continent_linestyle is not None:
plotvars.continent_linestyle = continent_linestyle
if text_fontsize is not None:
plotvars.text_fontsize = colorbar_fontsize
if colorbar_fontsize is not None:
plotvars.colorbar_fontsize = colorbar_fontsize
if text_fontweight is not None:
plotvars.text_fontweight = text_fontweight
if axis_label_fontweight is not None:
plotvars.axis_label_fontweight = axis_label_fontweight
if colorbar_fontweight is not None:
plotvars.colorbar_fontweight = colorbar_fontweight
if title_fontweight is not None:
plotvars.title_fontweight = title_fontweight
if viewer is not None:
plotvars.viewer = viewer
if tspace_year is not None:
plotvars.tspace_year = tspace_year
if tspace_month is not None:
plotvars.tspace_month = tspace_month
if tspace_day is not None:
plotvars.tspace_day = tspace_day
if tspace_hour is not None:
plotvars.tspace_hour = tspace_hour
if xtick_label_rotation is not None:
plotvars.xtick_label_rotation = xtick_label_rotation
if xtick_label_align is not None:
plotvars.xtick_label_align = xtick_label_align
if ytick_label_rotation is not None:
plotvars.ytick_label_rotation = ytick_label_rotation
if ytick_label_align is not None:
plotvars.ytick_label_align = ytick_label_align
if legend_text_size is not None:
plotvars.legend_text_size = legend_text_size
if legend_text_weight is not None:
plotvars.legend_text_weight = legend_text_weight
if cs_uniform is not None:
plotvars.cs_uniform = cs_uniform
if master_title is not None:
plotvars.master_title = master_title
if master_title_location is not None:
plotvars.master_title_location = master_title_location
if master_title_fontsize is not None:
plotvars.master_title_fontsize = master_title_fontsize
if master_title_fontweight is not None:
plotvars.master_title_fontweight = master_title_fontweight
if dpi is not None:
plotvars.dpi = dpi
if land_color is not None:
plotvars.land_color=land_color
if ocean_color is not None:
plotvars.ocean_color=ocean_color
if lake_color is not None:
plotvars.lake_color=lake_color
if rotated_grid_spacing is not None:
plotvars.rotated_grid_spacing=rotated_grid_spacing
if rotated_deg_spacing is not None:
plotvars.rotated_deg_spacing=rotated_deg_spacing
if rotated_grid_thickness is not None:
plotvars.rotated_grid_thickness=rotated_grid_thickness
if rotated_continents is not None:
plotvars.rotated_continents=rotated_continents
if rotated_grid is not None:
plotvars.rotated_grid=rotated_grid
if rotated_labels is not None:
plotvars.rotated_labels=rotated_labels
if legend_frame is not None:
plotvars.legend_frame=legend_frame
if legend_frame_edge_color is not None:
plotvars.legend_frame_edge_color=legend_frame_edge_color
if legend_frame_face_color is not None:
plotvars.legend_frame_face_color=legend_frame_face_color
if degsym is not None:
plotvars.degsym=degsym
if axis_width is not None:
plotvars.axis_width=axis_width
if grid is not None:
plotvars.grid=grid
if grid_spacing is not None:
plotvars.grid_spacing=grid_spacing
if grid_lons is not None:
plotvars.grid_lons=grid_lons
if grid_lats is not None:
plotvars.grid_lats=grid_lats
if grid_colour is not None:
plotvars.grid_colour=grid_colour
if grid_linestyle is not None:
plotvars.grid_linestyle=grid_linestyle
if grid_thickness is not None:
plotvars.grid_thickness=grid_thickness
[docs]def vloc(xvec=None, yvec=None, lons=None, lats=None):
"""
| vloc is used to locate the positions of a set of points in a vector
|
|
|
| xvec=None - data longitudes
| yvec=None - data latitudes
| lons=None - required longitude positions
| lats=None - required latitude positions
:Returns:
locations of user points in the longitude and latitude points
|
|
|
|
|
|
|
"""
import numpy as np
# Check input parameters
if any(val is None for val in [xvec, yvec, lons, lats]):
errstr = '\nvloc error\n'
errstr += 'xvec, yvec, lons, lats all need to be passed to vloc to\n'
errstr += 'generate a set of location points\n'
raise Warning(errstr)
xarr = np.zeros(np.size(lons))
yarr = np.zeros(np.size(lats))
# Convert longitudes to -180 to 180.
for i in np.arange(np.size(xvec)):
xvec[i] = ((xvec[i] + 180) % 360) - 180
for i in np.arange(np.size(lons)):
lons[i] = ((lons[i] + 180) % 360) - 180
# Centre around 180 degrees longitude if needed.
if (max(xvec) > 150):
for i in np.arange(np.size(xvec)):
xvec[i] = (xvec[i] + 360.0) % 360.0
pts = np.where(xvec < 0.0)
xvec[pts] = xvec[pts] + 360.0
for i in np.arange(np.size(lons)):
lons[i] = (lons[i] + 360.0) % 360.0
pts = np.where(lons < 0.0)
lons[pts] = lons[pts] + 360.0
# Find position in array
for i in np.arange(np.size(lons)):
if ((lons[i] < min(xvec)) or (lons[i] > max(xvec))):
xpt = -1
else:
xpts = np.where(lons[i] >= xvec)
xpt = np.nanmax(xpts)
if ((lats[i] < min(yvec)) or (lats[i] > max(yvec))):
ypt = -1
else:
ypts = np.where(lats[i] >= yvec)
ypt = np.nanmax(ypts)
if (xpt >= 0):
xarr[i] = xpt + (lons[i] - xvec[xpt]) / (xvec[xpt + 1] - xvec[xpt])
else:
xarr[i] = None
if (ypt >= 0) and ypt <= np.size(yvec) - 2:
yarr[i] = ypt + (lats[i] - yvec[ypt]) / (yvec[ypt + 1] - yvec[ypt])
else:
yarr[i] = None
return (xarr, yarr)
[docs]def rgaxes(xpole=None, ypole=None, xvec=None, yvec=None,
xticks=None, xticklabels=None, yticks=None, yticklabels=None,
axes=None, xaxis=None, yaxis=None, xlabel=None, ylabel=None):
"""
| rgaxes - label rotated grid plots
|
| xpole=None - location of xpole in degrees
| ypole=None - location of ypole in degrees
| xvec=None - location of x grid points
| yvec=None - location of y grid points
|
| axes=True - plot x and y axes
| xaxis=True - plot xaxis
| yaxis=True - plot y axis
| xticks=None - xtick positions
| xticklabels=None - xtick labels
| yticks=None - y tick positions
| yticklabels=None - ytick labels
| xlabel=None - label for x axis
| ylabel=None - label for y axis
|
:Returns:
name
|
|
|
|
|
|
"""
spacing = plotvars.rotated_grid_spacing
degspacing = plotvars.rotated_deg_spacing
continents = plotvars.rotated_continents
grid = plotvars.rotated_grid
labels = plotvars.rotated_labels
grid_thickness=plotvars.rotated_grid_thickness
# Invert y array if going from north to south
# Otherwise this gives nans for all output
yvec_orig = yvec
if (yvec[0] > yvec[np.size(yvec) - 1]):
yvec = yvec[::-1]
gset(xmin=0, xmax=np.size(xvec) - 1,
ymin=0, ymax=np.size(yvec) - 1, user_gset=0)
# Set continent thickness and color if not already set
if plotvars.continent_thickness is None:
continent_thickness = 1.5
if plotvars.continent_color is None:
continent_color = 'k'
# Draw continents
if continents:
import cartopy.io.shapereader as shpreader
import shapefile
shpfilename = shpreader.natural_earth(resolution=plotvars.resolution, category='physical', name='coastline')
reader = shapefile.Reader(shpfilename)
shapes = [s.points for s in reader.shapes()]
for shape in shapes:
lons, lats = zip(*shape)
lons=np.array(lons)
lats=np.array(lats)
rotated_transform = ccrs.RotatedPole(pole_latitude=ypole, pole_longitude=xpole)
points=rotated_transform.transform_points(ccrs.PlateCarree(), lons, lats)
xout = np.array(points)[:, 0]
yout = np.array(points)[:, 1]
xpts, ypts = vloc(lons=xout, lats=yout, xvec=xvec, yvec=yvec)
plotvars.plot.plot(xpts, ypts, linewidth=continent_thickness,
color=continent_color)
if xticks is None:
lons = -180 + np.arange(360 / spacing + 1) * spacing
else:
lons=xticks
if yticks is None:
lats = -90 + np.arange(180 / spacing + 1) * spacing
else:
lats=yticks
# Latitudes
if axes:
if xaxis:
for val in np.arange(np.size(lons)):
ipts = 179.0 / degspacing
lona = np.zeros(int(ipts)) + lons[val]
lata = -90 + np.arange(ipts - 1) * degspacing
rotated_transform = ccrs.RotatedPole(pole_latitude=ypole, pole_longitude=xpole)
points=rotated_transform.transform_points(ccrs.PlateCarree(), lona, lata)
xout = np.array(points)[:, 0]
yout = np.array(points)[:, 1]
xpts, ypts = vloc(lons=xout, lats=yout, xvec=xvec, yvec=yvec)
if grid:
plotvars.plot.plot(xpts, ypts, ':', linewidth=grid_thickness,
color='k')
if labels:
# Make a label unless the axis is all Nans
if (np.size(ypts[5:]) > np.sum(np.isnan(ypts[5:]))):
ymin = np.nanmin(ypts[5:])
loc = np.where(ypts == ymin)[0]
if np.size(loc) > 1:
loc = loc[1]
if loc > 0:
if np.isfinite(xpts[loc]):
line = matplotlib.lines.Line2D(
[xpts[loc], xpts[loc]], [0, -2], color='k')
plotvars.plot.add_line(line)
line.set_clip_on(False)
fw = plotvars.text_fontweight
if xticklabels is None:
xticklabel=mapaxis(lons[val], lons[val], type=1)[1][0]
else:
xticklabel=xticks[val]
plotvars.plot.text(xpts[loc], -5,
xticklabel,
horizontalalignment='center',
verticalalignment='top',
fontsize=plotvars.text_fontsize,
fontweight=fw)
# Longitudes
if axes:
if yaxis:
for val in np.arange(np.size(lats)):
ipts = 359.0 / degspacing
lata = np.zeros(int(ipts)) + lats[val]
lona = -180.0 + np.arange(ipts - 1) * degspacing
rotated_transform = ccrs.RotatedPole(pole_latitude=ypole, pole_longitude=xpole)
points=rotated_transform.transform_points(ccrs.PlateCarree(), lona, lata)
xout = np.array(points)[:, 0]
yout = np.array(points)[:, 1]
xpts, ypts = vloc(lons=xout, lats=yout, xvec=xvec, yvec=yvec)
if grid:
plotvars.plot.plot(xpts, ypts, ':', linewidth=grid_thickness,
color='k')
if labels:
# Make a label unless the axis is all Nans
if (np.size(xpts[5:]) > np.sum(np.isnan(xpts[5:]))):
xmin = np.nanmin(xpts[5:])
loc = np.where(xpts == xmin)[0]
if np.size(loc) == 1:
if loc > 0:
if np.isfinite(ypts[loc]):
line = matplotlib.lines.Line2D(
[0, -2], [ypts[loc], ypts[loc]], color='k')
plotvars.plot.add_line(line)
line.set_clip_on(False)
fw = plotvars.text_fontweight
if yticklabels is None:
yticklabel=mapaxis(lats[val], lats[val], type=2)[1][0]
else:
yticklabel=yticks[val]
plotvars.plot.text(-5, ypts[loc],
yticklabel,
horizontalalignment='right',
verticalalignment='center',
fontsize=plotvars.text_fontsize,
fontweight=fw)
# Reset yvec
yvec = yvec_orig
[docs]def lineplot(f=None, x=None, y=None, fill=True, lines=True, line_labels=True,
title=None, ptype=0, linestyle='-', linewidth=1.0, color='k',
xlog=False, ylog=False, verbose=None, swap_xy=False,
marker=None, markersize=5.0, markeredgecolor = 'k',
markeredgewidth=0.5, label=None,
legend_location='upper right', xunits=None, yunits=None,
xlabel=None, ylabel=None, xticks=None, yticks=None,
xticklabels=None, yticklabels=None, xname=None, yname=None,
axes=True, xaxis=True, yaxis=True,zorder=None):
"""
| lineplot is the interface to line plotting in cf-plot.
| The minimum use is lineplot(f) where f is a CF field.
| If x and y are passed then an appropriate plot is made allowing
| x vs data and y vs data plots.
| When making a labelled line plot:
| always have a label for each line
| always put the legend location as an option to the last call to lineplot
|
| f - CF data used to make a line plot
| x - x locations of data in y
| y - y locations of data in x
| linestyle='-' - line style
| color='k - line color
| linewidth=1.0 - line width
| marker=None - marker for points along the line
| markersize=5.0 - size of the marker
| markeredgecolor = 'k' - colour of edge around the marker
| markeredgewidth = 0.5 - width of edge around the marker
| xlog=False - log x-axis
| ylog=False - log y-axis
| label=None - line label - label for line
| legend_location='upper right' - default location of legend
| Other options are {'best': 0, 'center': 10, 'center left': 6,
| 'center right': 7, 'lower center': 8,
| 'lower left': 3, 'lower right': 4, 'right': 5,
| 'upper center': 9, 'upper left': 2, 'upper right': 1}
|
| verbose=None - change to 1 to get a verbose listing of what lineplot
| is doing
| zorder=None - plotting order
|
| The following parameters override any CF data defaults:
| title=None - plot title
| xunits=None - x units
| yunits=None - y units
| xlabel=None - x name
| ylabel=None - y name
| xname=None - depreciated keyword
| yname=None - depreciated keyword
| xticks=None - x ticks
| xticklabels=None - x tick labels
| yticks=None - y ticks
| yticklabels - y tick labels
| axes=True - plot x and y axes
| xaxis=True - plot xaxis
| yaxis=True - plot y axis
|
|
| When making a multiple lineplot the last call to lineplot is the one that
| any of the above axis overrides should be placed in.
|
|
"""
if verbose:
print 'lineplot - making a line plot'
# Catch depreciated keywords
if xname is not None or yname is not None:
print '\nlineplot error'
print 'xname and yname are now depreciated keywords'
print 'Please use xlabel and ylabel\n'
return
##################
# Open a new plot is necessary
##################
if plotvars.user_plot == 0:
gopen(user_plot=0)
##################
# Extract required data
# If a cf-python field
##################
cf_field = False
if f is not None:
if isinstance(f, cf.Field):
cf_field = True
elif isinstance(f, cf.FieldList):
raise TypeError("Can't plot a field list")
plot_xlabel = ''
plot_ylabel = ''
xlabel_units=''
ylabel_units=''
if cf_field:
# Extract data
if verbose:
print 'lineplot - CF field, extracting data'
has_count = 0
for mydim in f.items():
if np.size(np.squeeze(f.item(mydim).array)) > 1:
has_count = has_count + 1
x = np.squeeze(f.item(mydim).array)
#x label
xlabel_units = str(getattr(f.item(mydim), 'Units', ''))
plot_xlabel = cf_var_name(field=f, dim=mydim) + ' ('
plot_xlabel += xlabel_units + ')'
y = np.squeeze(f.array)
#y label
if hasattr(f, 'id'):
plot_ylabel = f.id
if hasattr(f, 'ncvar'):
plot_ylabel = f.ncvar
if hasattr(f, 'short_name'):
plot_ylabel = f.short_name
if hasattr(f, 'long_name'):
plot_ylabel = f.long_name
if hasattr(f, 'standard_name'):
plot_ylabel = f.standard_name
if hasattr(f, 'Units'):
ylabel_units = str(f.Units)
else:
ylabel_units = ''
plot_ylabel += ' (' + ylabel_units + ')'
if has_count != 1:
errstr = '\n lineplot error - passed field is not suitable '
errstr += 'for plotting as a line\n'
for mydim in f.items():
sn = getattr(f.item(mydim), 'standard_name', False)
ln = getattr(f.item(mydim), 'long_name', False)
if sn:
errstr = errstr + \
str(mydim) + ',' + str(sn) + ',' + \
str(f.item(mydim).size) + '\n'
else:
if ln:
errstr = errstr + \
str(mydim) + ',' + str(ln) + ',' + \
str(f.item(mydim).size) + '\n'
raise Warning(errstr)
else:
if verbose:
print 'lineplot - not a CF field, using passed data'
errstr = ''
if x is None or y is None:
errstr = 'lineplot error- must define both x and y'
if f is not None:
errstr += 'lineplot error- must define just x and y to make '
errstr += 'a lineplot'
if errstr != '':
raise Warning('\n' + errstr + '\n')
# Z on y-axis
ztype = None
if xlabel_units in ['mb', 'mbar', 'millibar', 'decibar',
'atmosphere', 'atm', 'pascal', 'Pa', 'hPa']:
ztype = 1
if xlabel_units in ['meter', 'metre', 'm', 'kilometer', 'kilometre', 'km']:
ztype = 2
if ztype is not None:
temporary_xpts = x
temporary_ypts = y
temporary_xlabel = plot_xlabel
temporary_ylabel = plot_ylabel
x = temporary_ypts
y = temporary_xpts
plot_xlabel = temporary_ylabel
plot_ylabel = temporary_xlabel
# Set data values
if verbose:
print 'lineplot - setting data values'
xpts = np.squeeze(x)
ypts = np.squeeze(y)
minx = np.min(x)
miny = np.min(y)
maxx = np.max(x)
maxy = np.max(y)
if cf_field:
taxis = f.item('T')
if ztype == 1:
miny = np.max(y)
maxy = np.min(y)
if ztype == 2:
if f.positive == 'down':
miny = np.max(y)
maxy = np.min(y)
# Use user set values if present
time_xstr = False
time_ystr = False
if plotvars.xmin is not None:
minx = plotvars.xmin
miny = plotvars.ymin
maxx = plotvars.xmax
maxy = plotvars.ymax
# Change from date string to a number if strings are passed
try:
float(minx)
except:
time_xstr = True
try:
float(miny)
except:
time_ystr = True
if cf_field:
taxis = f.item('T')
if time_xstr or time_ystr:
ref_time = f.item('T').units
ref_calendar = f.item('T').calendar
ref_time_origin = str(f.item('T').Units.reftime)
time_units = cf.Units(ref_time, ref_calendar)
if time_xstr:
t = cf.Data(cf.dt(minx), units=time_units)
minx = t.array
t = cf.Data(cf.dt(maxx), units=time_units)
maxx = t.array
taxis = cf.Data([cf.dt(plotvars.xmin),
cf.dt(plotvars.xmax)], units=time_units)
if time_ystr:
t = cf.Data(cf.dt(miny), units=time_units)
miny = t.array
t = cf.Data(cf.dt(maxy), units=time_units)
maxy = t.array
taxis = cf.Data([cf.dt(plotvars.ymin),
cf.dt(plotvars.ymax)], units=time_units)
# Set x and y labelling
# Retrieve any user defined axis labels
if plot_xlabel == '' and plotvars.xlabel is not None:
plot_xlabel = plotvars.xlabel
if plot_ylabel == '' and plotvars.ylabel is not None:
plot_ylabel = plotvars.ylabel
if xticks is None and plotvars.xticks is not None:
xticks = plotvars.xticks
if plotvars.xticklabels is not None:
xticklabels = plotvars.xticklabels
else:
xticklabels = map(str, xticks)
if yticks is None and plotvars.yticks is not None:
yticks = plotvars.yticks
if plotvars.yticklabels is not None:
yticklabels = plotvars.yticklabels
else:
yticklabels = map(str, yticks)
mod = 0
tight = 0
if plotvars.user_gset == 1:
tight = 1
if xticks is None:
if plot_xlabel[0:3].lower() == 'lon':
xticks, xticklabels = mapaxis(minx, maxx, type=1)
if plot_xlabel[0:3].lower() == 'lat':
xticks, xticklabels = mapaxis(minx, maxx, type=2)
if cf_field:
if xticks is None:
if f.item('T') is not None:
if np.size(f.item('T').array) > 1:
xticks, xticklabels, plot_xlabel = timeaxis(taxis)
if xticks is None:
xticks = gvals(dmin=minx, dmax=maxx, tight=tight, mod=mod)[0]
xticklabels = xticks
else:
if xticklabels is None:
xticklabels=[]
for val in xticks:
xticklabels.append('{}'.format(val))
if yticks is None:
if abs(maxy - miny) > 1:
if miny < maxy:
yticks = gvals(dmin=miny, dmax=maxy, tight=tight, mod=mod)[0]
if maxy < miny:
yticks = gvals(dmin=maxy, dmax=miny, tight=tight, mod=mod)[0]
else:
yticks = gvals(dmin=miny, dmax=maxy, tight=tight, mod=0)[0]
if yticklabels is None:
yticklabels=[]
for val in yticks:
yticklabels.append('{}'.format(val))
if xlabel is not None:
plot_xlabel = xlabel
if xunits is not None:
plot_xlabel += '('+xunits+')'
if ylabel is not None:
plot_ylabel = ylabel
if yunits is not None:
plot_ylabel += '('+yunits+')'
if swap_xy:
if verbose:
print 'lineplot - swapping x and y'
xpts, ypts = ypts, xpts
minx, miny = miny, minx
maxx, maxy = maxy, maxx
plot_xlabel, plot_ylabel = plot_ylabel, plot_xlabel
xticks, yticks = yticks, xticks
xticklabels, yticklabels = yticklabels, xticklabels
if plotvars.user_gset == 1:
if time_xstr is False and time_ystr is False:
minx = plotvars.xmin
maxx = plotvars.xmax
miny = plotvars.ymin
maxy = plotvars.ymax
if axes:
if xaxis is not True:
xticks = [100000000]
xticklabels = xticks
plot_xlabel = ''
if yaxis is not True:
yticks = [100000000]
yticklabels = yticks
plot_ylabel = ''
else:
xticks = [100000000]
xticklabels = xticks
yticks = [100000000]
yticklabels = yticks
plot_xlabel = ''
plot_ylabel = ''
# Make graph
if verbose:
print 'lineplot - making graph'
xlabelalignment=plotvars.xtick_label_align
ylabelalignment=plotvars.ytick_label_align
graph=plotvars.plot
if plotvars.twinx:
graph=graph.twinx()
ylabelalignment='left'
if plotvars.twiny:
graph=graph.twiny()
graph.axis([minx, maxx, miny, maxy])
graph.tick_params(direction='out', which='both', right=True, top=True)
graph.set_xlabel(plot_xlabel, fontsize=plotvars.axis_label_fontsize,
fontweight=plotvars.axis_label_fontweight)
graph.set_ylabel(plot_ylabel, fontsize=plotvars.axis_label_fontsize,
fontweight=plotvars.axis_label_fontweight)
if plotvars.xlog or xlog:
graph.set_xscale('log')
if plotvars.ylog or ylog:
graph.set_yscale('log')
if xticks is not None:
graph.set_xticks(xticks)
graph.set_xticklabels(
xticklabels,
rotation=plotvars.xtick_label_rotation,
horizontalalignment=xlabelalignment, fontsize=plotvars.axis_label_fontsize,
fontweight=plotvars.axis_label_fontweight)
if yticks is not None:
graph.set_yticks(yticks)
graph.set_yticklabels(
yticklabels,
rotation=plotvars.ytick_label_rotation,
horizontalalignment=ylabelalignment, fontsize=plotvars.axis_label_fontsize,
fontweight=plotvars.axis_label_fontweight)
graph.plot(xpts, ypts, color=color, linestyle=linestyle,
linewidth=linewidth, marker=marker,
markersize=markersize,
markeredgecolor=markeredgecolor,
markeredgewidth=markeredgewidth,
label=label, zorder=zorder)
# Time axes sometimes spill over the edge of the plot limits so
# use tight_layout() to adjust these plots
#if cf_field:
# if f.item('T') is not None:
# if np.size(f.item('T').array) > 1:
# plotvars.master_plot.tight_layout()
# Set axis width if required
if plotvars.axis_width is not None:
for axis in ['top','bottom','left','right']:
plotvars.plot.spines[axis].set_linewidth(plotvars.axis_width)
# Add a legend if needed
if label is not None:
legend_properties = {
'size': plotvars.legend_text_size,
'weight': plotvars.legend_text_weight}
graph.legend(loc=legend_location, prop=legend_properties,
frameon=plotvars.legend_frame,
edgecolor=plotvars.legend_frame_edge_color,
facecolor=plotvars.legend_frame_face_color)
# Set title
if title is not None:
graph.set_title(title, fontsize=plotvars.title_fontsize,
fontweight=plotvars.title_fontweight)
##################
# Save or view plot
##################
if plotvars.user_plot == 0:
if verbose:
print 'Saving or viewing plot'
gclose()
def regression_tests():
"""
| Test for cf-plot regressions
| Run through some standard levs, gvals, lon and lat labelling
| Make all the gallery plots and use Imagemaick to display them
| alongside a reference plot
|
|
|
|
|
"""
print '=================='
print 'Regression testing'
print '=================='
print ''
print '------------------'
print 'Testing for levels'
print '------------------'
ref_answer = [-35, -30, -25, -20, -15, -10, -5, 0, 5,
10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65]
compare_arrays(ref=ref_answer, levs_test=True, min=-35, max=65, step=5)
ref_answer = [-6., -4.8, -3.6, -2.4, -1.2, 0., 1.2, 2.4, 3.6, 4.8, 6.]
compare_arrays(ref=ref_answer, levs_test=True, min=-6, max=6, step=1.2)
ref_answer = [50000, 51000, 52000, 53000, 54000, 55000, 56000, 57000,
58000, 59000, 60000]
compare_arrays(ref=ref_answer, levs_test=True, min=50000, max=60000,
step=1000)
ref_answer = [-7000, -6500, -6000, -5500, -5000, -4500, -4000, -3500,
-3000, -2500, -2000, -1500, -1000, -500]
compare_arrays(
ref=ref_answer,
levs_test=True,
min=-7000,
max=-300,
step=500)
print ''
print '-----------------'
print 'Testing for gvals'
print '-----------------'
ref_answer = [281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292,
293]
compare_arrays(ref=ref_answer, min=280.50619506835938,
max=293.48431396484375, mult=0, gvals_test=True)
ref_answer = [0.356, 0.385, 0.414, 0.443, 0.472, 0.501, 0.53 , 0.559,
0.588, 0.617, 0.646, 0.675]
compare_arrays(ref=ref_answer, min=0.356, max=0.675, mult=0,
gvals_test=True)
ref_answer = [-45, -40, -35, -30, -25, -20, -15, -10, -5, 0, 5, 10, 15,
20, 25, 30, 35, 40, 45, 50]
compare_arrays(ref=ref_answer, min=-49.510975, max=53.206604, mult=0,
gvals_test=True)
ref_answer = [47000, 48000, 49000, 50000, 51000, 52000, 53000, 54000,
55000, 56000, 57000, 58000, 59000, 60000, 61000, 62000,
63000, 64000]
compare_arrays(ref=ref_answer, min=46956, max=64538, mult=0,
gvals_test=True)
ref_answer = [-1. , -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0. ,
0.1]
compare_arrays(ref=ref_answer, min=-1.0, max=0.1, mult=0,
gvals_test=True)
print ''
print '----------------------------------------'
print 'Testing for longitude/latitude labelling'
print '----------------------------------------'
ref_answer = ([-180, -120, -60, 0, 60, 120, 180],
['180', '120W', '60W', '0', '60E', '120E', '180'])
compare_arrays(ref=ref_answer, min=-180, max=180, type=1,
mapaxis_test=True)
ref_answer = ([150, 180, 210, 240, 270],
['150E', '180', '150W', '120W', '90W'])
compare_arrays(ref=ref_answer, min=135, max=280, type=1,
mapaxis_test=True)
ref_answer = ([0, 10, 20, 30, 40, 50, 60, 70, 80, 90], ['0', '10E', '20E',
'30E', '40E', '50E', '60E', '70E', '80E', '90E'])
compare_arrays(ref=ref_answer, min=0, max=90, type=1, mapaxis_test=True)
ref_answer = ([-90, -60, -30, 0, 30, 60, 90],
['90S', '60S', '30S', '0', '30N', '60N', '90N'])
compare_arrays(ref=ref_answer, min=-90, max=90, type=2, mapaxis_test=True)
ref_answer = ([0, 5, 10, 15, 20, 25, 30],
['0', '5N', '10N', '15N', '20N', '25N', '30N'])
compare_arrays(ref=ref_answer, min=0, max=30, type=2, mapaxis_test=True)
print ''
print '-----------------'
print 'Testing for plots'
print '-----------------'
# Run through gallery examples and compare to reference plots
# example1
reset()
setvars(file='fig1.png')
f = cf.read('/opt/graphics/cfplot_data/tas_A1.nc')[0]
con(f.subspace(time=15))
compare_images(1)
# example2
reset()
setvars(file='fig2.png')
f = cf.read('/opt/graphics/cfplot_data/tas_A1.nc')[0]
con(f.subspace(time=15), blockfill=True, lines=False)
compare_images(2)
# example3
reset()
setvars(file='fig3.png')
f = cf.read('/opt/graphics/cfplot_data/tas_A1.nc')[0]
mapset(lonmin=-15, lonmax=3, latmin=48, latmax=60)
levs(min=265, max=285, step=1)
con(f.subspace(time=15))
compare_images(3)
# example4
reset()
setvars(file='fig4.png')
f = cf.read('/opt/graphics/cfplot_data/ggap.nc')[1]
mapset(proj='npstere')
con(f.subspace(pressure=500))
compare_images(4)
# example5
reset()
setvars(file='fig5.png')
f = cf.read('/opt/graphics/cfplot_data/ggap.nc')[1]
mapset(proj='spstere', boundinglat=-30, lon_0=180)
con(f.subspace(pressure=500))
compare_images(5)
# example6
reset()
setvars(file='fig6.png')
f = cf.read('/opt/graphics/cfplot_data/ggap.nc')[3]
con(f.subspace(longitude=0))
compare_images(6)
# example7
reset()
setvars(file='fig7.png')
f = cf.read('/opt/graphics/cfplot_data/ggap.nc')[1]
con(f.collapse('mean', 'longitude'))
compare_images(7)
# example8
reset()
setvars(file='fig8.png')
f = cf.read('/opt/graphics/cfplot_data/ggap.nc')[1]
con(f.collapse('mean', 'longitude'), ylog=1)
compare_images(8)
# example9
reset()
setvars(file='fig9.png')
f = cf.read('/opt/graphics/cfplot_data/ggap.nc')[0]
con(f.collapse('mean', 'latitude'))
compare_images(9)
# example10
reset()
setvars(file='fig10.png')
f = cf.read('/opt/graphics/cfplot_data/tas_A1.nc')[0]
cscale('plasma')
con(f.subspace(longitude=0), lines=0)
compare_images(10)
# example11
reset()
setvars(file='fig11.png')
f = cf.read('/opt/graphics/cfplot_data/tas_A1.nc')[0]
gset(-30, 30, '1960-1-1', '1980-1-1')
levs(min=280, max=305, step=1)
cscale('plasma')
con(f.subspace(longitude=0), lines=0)
compare_images(11)
# example12
reset()
setvars(file='fig12.png')
f = cf.read('/opt/graphics/cfplot_data/tas_A1.nc')[0]
cscale('plasma')
con(f.subspace(latitude=0), lines=0)
compare_images(12)
# example13
reset()
setvars(file='fig13.png')
f = cf.read('/opt/graphics/cfplot_data/ggap.nc')
u = f[1].subspace(pressure=500)
v = f[3].subspace(pressure=500)
vect(u=u, v=v, key_length=10, scale=100, stride=5)
compare_images(13)
# example14
reset()
setvars(file='fig14.png')
f = cf.read('/opt/graphics/cfplot_data/ggap.nc')
u = f[1].subspace(pressure=500)
v = f[3].subspace(pressure=500)
t = f[0].subspace(pressure=500)
gopen()
mapset(lonmin=10, lonmax=120, latmin=-30, latmax=30)
levs(min=254, max=270, step=1)
con(t)
vect(u=u, v=v, key_length=10, scale=50, stride=2)
gclose()
compare_images(14)
# example15
reset()
setvars(file='fig15.png')
u = cf.read('/opt/graphics/cfplot_data/ggap.nc')[1]
v = cf.read('/opt/graphics/cfplot_data/ggap.nc')[3]
u = u.subspace(Z=500)
v = v.subspace(Z=500)
mapset(proj='npstere')
vect(u=u, v=v, key_length=10, scale=100, pts=40,
title='Polar plot with regular point distribution')
compare_images(15)
# example16
reset()
setvars(file='fig16.png')
c = cf.read('/opt/graphics/cfplot_data/vaAMIPlcd_DJF.nc')[0]
c = c.subspace(Y=cf.wi(-60, 60))
c = c.subspace(X=cf.wi(80, 160))
c = c.collapse('T: mean X: mean')
g = cf.read('/opt/graphics/cfplot_data/wapAMIPlcd_DJF.nc')[0]
g = g.subspace(Y=cf.wi(-60, 60))
g = g.subspace(X=cf.wi(80, 160))
g = g.collapse('T: mean X: mean')
vect(u=c, v=-g, key_length=[5, 0.05],
scale=[20, 0.2], title='DJF', key_location=[0.95, -0.05])
compare_images(16)
# example17
reset()
setvars(file='fig17.png')
f = cf.read('/opt/graphics/cfplot_data/tas_A1.nc')[0]
g = f.subspace(time=15)
gopen()
cscale('magma')
con(g)
stipple(f=g, min=220, max=260, size=100, color='#00ff00')
stipple(f=g, min=300, max=330, size=50, color='#0000ff', marker='s')
gclose()
compare_images(17)
# example18
reset()
setvars(file='fig18.png')
f = cf.read('/opt/graphics/cfplot_data/tas_A1.nc')[0]
g = f.subspace(time=15)
gopen()
cscale('magma')
mapset(proj='npstere')
con(g)
stipple(f=g, min=265, max=295, size=100, color='#00ff00')
gclose()
compare_images(18)
# example19
reset()
setvars(file='fig19.png')
f = cf.read('/opt/graphics/cfplot_data/ggap.nc')[1]
gopen(rows=2, columns=2, bottom=0.2)
gpos(1)
con(f.subspace(pressure=500), colorbar=None)
gpos(2)
mapset(proj='moll')
con(f.subspace(pressure=500), colorbar=None)
gpos(3)
mapset(proj='npstere', boundinglat=30, lon_0=180)
con(f.subspace(pressure=500), colorbar=None)
gpos(4)
mapset(proj='spstere', boundinglat=-30, lon_0=180)
con(f.subspace(pressure=500), colorbar_position=[
0.1, 0.1, 0.8, 0.02], colorbar_orientation='horizontal')
gclose()
compare_images(19)
# example20
reset()
setvars(file='fig20.png')
f = cf.read('/opt/graphics/cfplot_data/Geostropic_Adjustment.nc')[0]
con(f.subspace[9])
compare_images(20)
# example21
reset()
setvars(file='fig21.png')
f = cf.read('/opt/graphics/cfplot_data/Geostropic_Adjustment.nc')[0]
con(f.subspace[9], title='test data',
xticks=np.arange(5) * 100000 + 100000,
yticks=np.arange(7) * 2000 + 2000,
xlabel='x-axis', ylabel='z-axis')
compare_images(21)
# example22
reset()
setvars(file='fig22.png')
f = cf.read_field('/opt/graphics/cfplot_data/rgp.nc')
cscale('gray')
con(f)
compare_images(22)
# example23
reset()
setvars(file='fig23.png')
f = cf.read_field('/opt/graphics/cfplot_data/rgp.nc')
data = f.array
xvec = f.item('dim1').array
yvec = f.item('dim0').array
xpole = 160
ypole = 30
gopen()
cscale('plasma')
xpts = np.arange(np.size(xvec))
ypts = np.arange(np.size(yvec))
gset(xmin=0, xmax=np.size(xvec) - 1, ymin=0, ymax=np.size(yvec) - 1)
levs(min=980, max=1035, step=2.5)
con(data, xpts, ypts[::-1])
rgaxes(xpole=xpole, ypole=ypole, xvec=xvec, yvec=yvec)
gclose()
compare_images(23)
# example24
reset()
setvars(file='fig24.png')
from matplotlib.mlab import griddata
# Arrays for data
lons = []
lats = []
pressure = []
temp = []
# Read data
f = open('/opt/graphics/cfplot_data/synop_data.txt')
lines = f.readlines()
for line in lines:
mysplit = line.split()
lons = np.append(lons, float(mysplit[1]))
lats = np.append(lats, float(mysplit[2]))
pressure = np.append(pressure, float(mysplit[3]))
temp = np.append(temp, float(mysplit[4]))
# Linearly interpolate data to a regular grid
lons_new = np.arange(140) * 0.1 - 11.0
lats_new = np.arange(140) * 0.1 + 49.0
temp_new = griddata(lons, lats, temp, lons_new, lats_new, interp='linear')
cscale('parula')
con(x=lons_new, y=lats_new, f=temp_new, ptype=1)
compare_images(24)
# example25
reset()
setvars(file='fig25.png')
gopen()
con(x=lons_new, y=lats_new, f=temp_new, ptype=1)
for i in np.arange(len(lines)):
plotvars.plot.text(float(lons[i]), float(lats[i]), str(temp[i]),
horizontalalignment='center',
verticalalignment='center')
gclose()
compare_images(25)
# example26
reset()
setvars(file='fig26.png')
from netCDF4 import Dataset as ncfile
from matplotlib.mlab import griddata
# Get an Orca grid and flatten the arrays
nc = ncfile('/opt/graphics/cfplot_data/orca2.nc')
lons = np.array(nc.variables['longitude'])
lats = np.array(nc.variables['latitude'])
temp = np.array(nc.variables['sst'])
lons = lons.flatten()
lats = lats.flatten()
temp = temp.flatten()
# Add wrap around at both longitude limits
pts = np.squeeze(np.where(lons < -150))
lons = np.append(lons, lons[pts] + 360)
lats = np.append(lats, lats[pts])
temp = np.append(temp, temp[pts])
pts = np.squeeze(np.where(lons > 150))
lons = np.append(lons, lons[pts] - 360)
lats = np.append(lats, lats[pts])
temp = np.append(temp, temp[pts])
lons_new = np.arange(181 * 8) * 0.25 - 180.0
lats_new = np.arange(91 * 8) * 0.25 - 90.0
temp_new = griddata(lons, lats, temp, lons_new, lats_new, interp='linear')
con(x=lons_new, y=lats_new, f=temp_new, ptype=1)
compare_images(26)
# example27
reset()
setvars(file='fig27.png')
f = cf.read('/opt/graphics/cfplot_data/ggap.nc')[1]
g = f.collapse('X: mean')
lineplot(g.subspace(pressure=100), marker='o', color='blue',
title='Zonal mean zonal wind at 100mb')
compare_images(27)
# example28
reset()
setvars(file='fig28.png')
f = cf.read('/opt/graphics/cfplot_data/ggap.nc')[1]
g = f.collapse('X: mean')
xticks = [-90, -75, -60, -45, -30, -15, 0, 15, 30, 45, 60, 75, 90]
xticklabels = ['90S', '75S', '60S', '45S', '30S', '15S', '0', '15N',
'30N', '45N', '60N', '75N', '90N']
xpts = [-30, 30, 30, -30, -30]
ypts = [-8, -8, 5, 5, -8]
gset(xmin=-90, xmax=90, ymin=-10, ymax=50)
gopen()
lineplot(g.subspace(pressure=100), marker='o', color='blue',
title='Zonal mean zonal wind', label='100mb')
lineplot(g.subspace(pressure=200), marker='D', color='red',
label='200mb', xticks=xticks, xticklabels=xticklabels,
legend_location='upper right')
plotvars.plot.plot(xpts, ypts, linewidth=3.0, color='green')
plotvars.plot.text(35, -2, 'Region of interest',
horizontalalignment='left')
gclose()
compare_images(28)
# example29
reset()
setvars(file='fig29.png')
f = cf.read('/opt/graphics/cfplot_data/tas_A1.nc')[0]
temp = f.subspace(time=cf.wi(cf.dt('1900-01-01'), cf.dt('1980-01-01')))
temp_annual = temp.collapse('T: mean', group=cf.Y())
temp_annual_global = temp_annual.collapse('area: mean', weights='area')
temp_annual_global.Units -= 273.15
lineplot(
temp_annual_global,
title='Global average annual temperature',
color='blue')
compare_images(29)
def compare_images(example=None):
"""
| Compare images and return an error string if they don't match
|
|
|
|
|
|
|
"""
import hashlib
disp = which('display')
conv = which('convert')
comp = which('compare')
file = 'fig' + str(example) + '.png'
file_new = '/home/andy/cfplot.src/cfplot/' + file
file_ref = '/home/andy/regression/' + file
# Check md5 checksums are the same and display files if not
if hashlib.md5(open(file_new, 'rb').read()).hexdigest() != hashlib.md5(
open(file_ref, 'rb').read()).hexdigest():
print '***Failed example ' + str(example) + '**'
error_image = '/home/andy/cfplot.src/cfplot/' + 'error_' + file
diff_image = '/home/andy/cfplot.src/cfplot/' + 'difference_' + file
p = subprocess.Popen([comp, file_new, file_ref, diff_image])
(output, err) = p.communicate()
p_status = p.wait()
p = subprocess.Popen([conv, "+append", file_new,
file_ref, error_image])
(output, err) = p.communicate()
p_status = p.wait()
subprocess.Popen([disp, diff_image])
else:
print 'Passed example ' + str(example)
def compare_arrays(ref=None, levs_test=None, gvals_test=None,
mapaxis_test=None, min=None, max=None, step=None,
mult=None, type=None):
"""
| Compare arrays and return an error string if they don't match
|
|
|
|
|
|
|
"""
anom = 0
if levs_test:
levs(min, max, step)
if np.size(ref) != np.size(plotvars.levels):
anom = 1
else:
for val in np.arange(np.size(ref)):
if abs(ref[val] - plotvars.levels[val]) >= 1e-6:
anom = 1
if anom == 1:
print '***levs failure***'
print 'min, max, step are', min, max, step
print 'generated levels are:'
print plotvars.levels
print 'expected levels:'
print ref
else:
pass_str = 'Passed cfp.levs(min=' + str(min) + ', max='
pass_str += str(max) + ', step=' + str(step) + ')'
print pass_str
anom = 0
if gvals_test:
vals, testmult = gvals(min, max)
if np.size(ref) != np.size(vals):
anom = 1
else:
for val in np.arange(np.size(ref)):
if abs(ref[val] - vals[val]) >= 1e-6:
anom = 1
if mult != testmult:
anom = 1
if anom == 1:
print '***gvals failure***'
print 'cfp.gvals(' + str(min) + ', ' + str(max) + ')'
print ''
print 'generated values are:', vals
print 'with a multiplier of ', testmult
print ''
print 'expected values:', ref
print 'with a multiplier of ', mult
else:
pass_str = 'Passed cfp.gvals(' + str(min) + ', '+ str(max) + ')'
print pass_str
anom = 0
if mapaxis_test:
ref_ticks = ref[0]
ref_labels = ref[1]
test_ticks, test_labels = mapaxis(min=min, max=max, type=type)
if np.size(test_ticks) != np.size(ref_ticks):
anom = 1
else:
for val in np.arange(np.size(ref_ticks)):
if abs(ref_ticks[val] - test_ticks[val]) >= 1e-6:
anom = 1
if ref_labels[val] != test_labels[val]:
anom = 1
if anom == 1:
print '***mapaxis failure***'
print ''
print 'cfp.mapaxis(min=' + str(min) + ', max=' + str(max)
print ', type=' + str(type) + ')'
print 'generated values are:', test_ticks
print 'with labels:', test_labels
print ''
print 'expected ticks:', ref_ticks
print 'with labels:', ref_labels
else:
pass_str = 'Passed cfp.mapaxis(min=' + str(min) + ', max='
pass_str += str(max) + ', type=' + str(type) + ')'
print pass_str
[docs]def traj(f=None, title=None, ptype=0, linestyle='-', linewidth=1.0, linecolor='b',
marker='o', markersize=5.0, markerfacecolor='r',
markeredgecolor='g', latmax = None, latmin=None,
axes=True, xaxis=True, yaxis=True,
verbose=None, legend=False,
xlabel=None, ylabel=None, xticks=None, yticks=None,
xticklabels=None, yticklabels=None, colorbar=None,
colorbar_position=None, colorbar_orientation='horizontal',
colorbar_title=None,
vector=False, head_width=0.4, head_length=1.0,
fc='k', ec='k', zorder=None):
"""
| traj is the interface to trajectory plotting in cf-plot.
| The minimum use is traj(f) where f is a CF field.
|
| f - CF data used to make a line plot
| linestyle='-' - line style
| linecolor='b' - line colour
| linewidth=1.0 - line width
| marker='o' - marker for points along the line
| markersize=30 - size of the marker
| markerfacecolor='b' - colour of the marker face
| markeredgecolor='g' - colour of the marker edge
| legend=False - plot different colour markers based on a set of user levels
| zorder=None - order for plotting
| verbose=None - Set to True to get a verbose listing of what traj is doing
|
| The following parameters override any CF data defaults:
| title=None - plot title
| axes=True - plot x and y axes
| xaxis=True - plot xaxis
| yaxis=True - plot y axis
| xlabel=None - x name
| ylabel=None - y name
| xticks=None - x ticks
| xticklabels=None - x tick labels
| yticks=None - y ticks
| yticklabels=None - y tick labels
| colorbar=None - plot a colorbar
| colorbar_position=None - position of colorbar
| [xmin, ymin, x_extent,y_extent] in normalised
| coordinates. Use when a common colorbar
| is required for a set of plots. A typical set
| of values would be [0.1, 0.05, 0.8, 0.02]
| colorbar_orientation=None - orientation of the colorbar
| colorbar_title=None - title for the colorbar
|
| Vector options
| vector=False - Draw vectors
| head_width=2.0 - vector head width
| head_length=2.0 - vector head length
| fc='k' - vector face colour
| ec='k' - vector edge colour
"""
if verbose:
print 'traj - making a trajectory plot'
# Read in data
# Find the auxiliary lons and lats if provided
has_lons = False
has_lats = False
for mydim in f.items():
if mydim[:3] == 'aux':
name = cf_var_name(field=f, dim=mydim)
if name in ['longitude']:
lons = np.squeeze(f.item(mydim).array)
has_lons = True
if name in ['latitude']:
lats = np.squeeze(f.item(mydim).array)
has_lats = True
data=f.array
# Raise an error if lons and lats not found in the input data
if not has_lons or not has_lats:
message = '\n\n\ntraj error\n'
if not has_lons:
message += 'missing longitudes in the field auxiliary data\n'
if not has_lats:
message += 'missing latitudes in the field auxiliary data\n'
message += '\n\n\n'
raise TypeError(message)
if latmax is not None:
pts=np.where(lats >= latmax)
if np.size(pts) > 0:
lons[pts] = np.nan
lats[pts] = np.nan
if latmin is not None:
pts=np.where(lats <= latmin)
if np.size(pts) > 0:
lons[pts] = np.nan
lats[pts] = np.nan
# Set potential user axis labels
user_xlabel = xlabel
user_ylabel = ylabel
user_xlabel = ''
user_ylabel = ''
# Set plotting parameters
continent_thickness = 1.5
continent_color = 'k'
continent_linestyle = '-'
if plotvars.continent_thickness is not None:
continent_thickness = plotvars.continent_thickness
if plotvars.continent_color is not None:
continent_color = plotvars.continent_color
if plotvars.continent_linestyle is not None:
continent_linestyle = plotvars.continent_linestyle
land_color = plotvars.land_color
ocean_color = plotvars.ocean_color
lake_color = plotvars.lake_color
##################
# Open a new plot is necessary
##################
if plotvars.user_plot == 0:
gopen(user_plot=0)
# Set up mapping
if plotvars.user_mapset == 0:
plotvars.lonmin = -180
plotvars.lonmax = 180
plotvars.latmin = -90
plotvars.latmax = 90
set_map()
mymap = plotvars.mymap
user_mapset = plotvars.user_mapset
# Set the plot limits
gset(xmin=plotvars.lonmin, xmax=plotvars.lonmax,
ymin=plotvars.latmin, ymax=plotvars.latmax, user_gset=0)
# Set the map projection
if plotvars.proj == 'cyl':
lon_mid=(plotvars.lonmax - plotvars.lonmin)/2.0 + plotvars.lonmin
proj=ccrs.PlateCarree()
if plotvars.proj == 'npstere':
proj = ccrs.NorthPolarStereo(central_longitude=plotvars.lon_0)
if plotvars.proj == 'spstere':
proj = ccrs.SouthPolarStereo(central_longitude=plotvars.lon_0)
ntracks = np.shape(lons)[0]
##################################
# Line, symbol and vector plotting
##################################
for track in np.arange(ntracks):
xpts = lons[track, :]
ypts = lats[track, :]
xpts_orig = deepcopy(xpts)
xpts = np.mod(xpts + 180, 360) - 180
# Check if xpts are only within the remapped longitudes above
if np.min(xpts) < -170 or np.max(xpts) > 170:
xpts = xpts_orig
for ix in np.arange(np.size(xpts)-1):
diff = xpts[ix+1] - xpts[ix]
diff_limit = 60
if diff >= 60:
xpts[ix+1] = xpts[ix+1] - 360.0
if diff <= -60:
xpts[ix+1] = xpts[ix+1] + 360.0
# Plot lines and markers
plot_linewidth = linewidth
plot_markersize = markersize
if legend:
plot_markersize = 0.0
if plot_linewidth > 0.0 or plot_markersize > 0.0:
if verbose and track == 0 and linewidth > 0.0:
print 'plotting lines'
if verbose and track == 0 and markersize > 0.0:
print 'plotting markers'
mymap.plot(xpts, ypts, color=linecolor,
linewidth=plot_linewidth, linestyle=linestyle,
marker=marker, markersize=plot_markersize,
markerfacecolor=markerfacecolor, markeredgecolor=markeredgecolor,
zorder=zorder, clip_on=True, transform=ccrs.PlateCarree())
# Plot vectors
if vector:
if verbose and track == 0:
print 'plotting vectors'
if zorder is None:
plot_zorder = 101
else:
plot_zorder = zorder
if plotvars.proj == 'cyl':
pts = np.ma.MaskedArray.count(xpts)
for pt in np.arange(pts-1):
mymap.arrow(xpts[pt], ypts[pt],
xpts[pt+1] - xpts[pt],
ypts[pt+1] - ypts[pt],
head_width=head_width,
head_length=head_length,
fc=fc, ec=ec,
length_includes_head=True,
zorder=plot_zorder, clip_on=True,
transform=ccrs.PlateCarree())
# Plot different colour markers based on a user set of levels
if legend:
# Check levels are not None
levs = plotvars.levels
if plotvars.levels is not None:
if verbose:
print 'traj - plotting different colour markers based on a user set of levels'
levs = plotvars.levels
else:
# Automatic levels
if verbose:
print 'traj - generating automatic legend levels'
dmin = np.nanmin(data)
dmax = np.nanmax(data)
levs, mult = gvals(dmin=dmin, dmax=dmax, tight=0, mod=0)
# Add extend options to the levels if set
if plotvars.levels_extend == 'min' or plotvars.levels_extend == 'both':
levs = np.append(-1e-30, levs)
if plotvars.levels_extend == 'max' or plotvars.levels_extend == 'both':
levs = np.append(levs, 1e30)
# Set the default colour scale
if plotvars.cscale_flag == 0:
cscale('viridis', ncols=np.size(levs) + 1)
plotvars.cscale_flag = 0
# User selected colour map but no mods so fit to levels
if plotvars.cscale_flag == 1:
cscale(plotvars.cs_user, ncols=np.size(levs) + 1)
plotvars.cscale_flag = 1
# For polar stereographic plots mask any points outside the plotting limb
if plotvars.proj == 'npstere':
pts = np.where(lats < plotvars.boundinglat)
if np.size(pts) > 0:
lats[pts] = np.nan
if plotvars.proj == 'spstere':
pts = np.where(lats > plotvars.boundinglat)
if np.size(pts) > 0:
lats[pts] = np.nan
for track in np.arange(ntracks):
xpts = lons[track, :]
ypts = lats[track, :]
data2 = data[track, :]
for i in np.arange(np.size(levs)-1):
color = plotvars.cs[i]
pts = np.where(np.logical_and(data2 >= levs[i], data2 <= levs[i+1]))
if zorder is None:
plot_zorder = 101
else:
plot_zorder = zorder
if np.size(pts) > 0:
mymap.scatter(
xpts[pts], ypts[pts],
s=markersize*15,
c=color,
marker=marker,
edgecolors=markeredgecolor,
transform=ccrs.PlateCarree(), zorder=plot_zorder)
# Axes
plot_map_axes(axes=axes, xaxis=xaxis, yaxis=yaxis,
xticks=xticks, xticklabels=xticklabels,
yticks=yticks, yticklabels=yticklabels,
user_xlabel=user_xlabel, user_ylabel=user_ylabel,
verbose=verbose)
# Coastlines
feature = cfeature.NaturalEarthFeature(name='land', category='physical',
scale=plotvars.resolution, facecolor='none')
mymap.add_feature(feature, edgecolor=continent_color,
linewidth=continent_thickness,
linestyle=continent_linestyle)
if ocean_color is not None:
mymap.add_feature(cfeature.OCEAN, edgecolor='face', facecolor=ocean_color,
zorder=100)
if land_color is not None:
mymap.add_feature(cfeature.LAND, edgecolor='face', facecolor=land_color,
zorder=100)
if lake_color is not None:
mymap.add_feature(cfeature.LAKES, edgecolor='face', facecolor=lake_color,
zorder=100)
# Title
if title is not None:
map_title(title)
# Color bar
plot_colorbar = False
if colorbar is None and legend:
plot_colorbar = True
if colorbar:
plot_colorbar = True
if plot_colorbar:
if verbose:
print 'traj - adding colour bar'
if colorbar is None:
colorbar_orientation='horizontal'
if plotvars.proj == 'npstere' or plotvars.proj == 'spstere':
colorbar_orientation='vertical'
# Get plotting window coordinates and set colorbar coordinates
if colorbar_position is None:
bbox=plotvars.plot.get_position()
if colorbar_orientation == 'horizontal':
ax1 = plotvars.master_plot.add_axes([bbox.x0, bbox.y0, bbox.x1-bbox.x0, 0.02])
else:
ax1 = plotvars.master_plot.add_axes([bbox.x1, bbox.y0, 0.02, bbox.y1-bbox.y0])
else:
ax1 = plotvars.master_plot.add_axes(colorbar_position)
if (colorbar_title is None):
colorbar_title = 'No Name'
if hasattr(f, 'id'):
colorbar_title = f.id
if hasattr(f, 'ncvar'):
colorbar_title = f.ncvar
if hasattr(f, 'short_name'):
colorbar_title = f.short_name
if hasattr(f, 'long_name'):
colorbar_title = f.long_name
if hasattr(f, 'standard_name'):
colorbar_title = f.standard_name
if hasattr(f, 'Units'):
if str(f.Units) == '':
colorbar_title += ''
else:
colorbar_title += '(' + supscr(str(f.Units)) + ')'
# Get the colour map
colmap = cscale_get_map()
cmap = matplotlib.colors.ListedColormap(colmap)
extend=plotvars.levels_extend
if plotvars.levels is not None:
boundaries=plotvars.levels.astype(float)
else:
boundaries=levs[1:-1].astype(float)
# Add colorbar extensions if definded by levs. Using boundaries[0]-1
# for the lower and boundaries[-1]+1 is just for the colorbar and
# has no meaning for the plot.
if (extend == 'min' or extend == 'both'):
cmap.set_under(plotvars.cs[0])
boundaries = np.insert(boundaries, 0, boundaries[0]-1)
if (extend == 'max' or extend == 'both'):
cmap.set_over(plotvars.cs[-1])
boundaries = np.insert(boundaries, len(boundaries), boundaries[-1]+1)
cbar = matplotlib.colorbar.ColorbarBase(ax1, cmap=cmap,
norm=plotvars.norm,
extend=extend,
extendfrac='auto',
boundaries=boundaries,
ticks=plotvars.levels,
spacing='uniform',
orientation=colorbar_orientation)
# Reset plot position
gpos(pos=plotvars.pos)
cbar.set_label(colorbar_title, fontsize=plotvars.colorbar_fontsize,
fontweight=plotvars.colorbar_fontweight)
# Bug in Matplotlib colorbar labelling
# With clevs=[-1, 1, 10000, 20000, 30000, 40000, 50000, 60000]
# Labels are [0, 2, 10001, 20001, 30001, 40001, 50001, 60001]
# With a +1 near to the colorbar label
if plotvars.levels is not None:
cbar.set_ticks([i for i in plotvars.levels])
else:
cbar.set_ticks([i for i in levs[1:-1]])
for t in cbar.ax.get_xticklabels():
t.set_fontsize(plotvars.colorbar_fontsize)
t.set_fontweight(plotvars.colorbar_fontweight)
# reset plot limits if not a user plot
if plotvars.user_gset == 0:
gset()
##########
# Save plot
##########
if plotvars.user_plot == 0:
gclose()
def cfplot_colorbar(cfill=None, colorbar_labels=None,
colorbar_orientation=None,
colorbar_position=None,
colorbar_shrink=None,
cb_orient=None, colorbar_title=None,
text_fontsize=None,
text_fontweight=None,
colorbar_text_up_down=None,
colorbar_text_down_up=None,
colorbar_drawedges=None,
verbose=None):
"""
| cfplot_colorbar is the internal interface to the Matplotlib
| colorbar routine
|
| cfill=None
| colorbar_labels=None - colorbar labels
| colorbar_orientation=None - orientation 'horizontal' or 'vertical'
| colorbar_position=None
| colorbar_shrink=None
| cb_orient=None
| colorbar_title=None
| text_fontsize=None
| text_fontweight=None
| colorbar_text_up_down=None
| colorbar_text_down_up=None
| colorbar_drawedges=None - Draw internal delimeter lines in colorbar
| verbose=None
|
|
|
"""
if verbose:
print 'con - adding colour bar'
colorbar_fontsize = plotvars.colorbar_fontsize
colorbar_fontweight = plotvars.colorbar_fontweight
# Work out colour bar labeling
lbot=colorbar_labels
if colorbar_text_up_down:
lbot=colorbar_labels[1:][::2]
ltop=colorbar_labels[::2]
if colorbar_text_down_up:
lbot=colorbar_labels[::2]
ltop=colorbar_labels[1:][::2]
pad = 0.10
if plotvars.rows >= 3:
pad = 0.15
if plotvars.rows >= 5:
pad = 0.20
if colorbar_position is None:
cbar = plotvars.master_plot.colorbar(cfill,
ticks=lbot,
orientation=cb_orient,
aspect=75, pad=pad,
shrink=colorbar_shrink,
drawedges=colorbar_drawedges)
else:
position = plotvars.master_plot.add_axes(colorbar_position)
cbar = plotvars.master_plot.colorbar(
cfill, cax=position, ticks=lbot,
orientation=colorbar_orientation,
drawedges=colorbar_drawedges)
gpos(pos=plotvars.pos)
cbar.set_label(colorbar_title, fontsize=colorbar_fontsize,
fontweight=colorbar_fontweight)
#print dir(cbar.ax.xaxis)
#for tick in cbar.ax.xaxis.get_ticklines():
# tick.set_visible(False)
# Bug in Matplotlib colorbar labelling
# With clevs=[-1, 1, 10000, 20000, 30000, 40000, 50000, 60000]
# Labels are [0, 2, 10001, 20001, 30001, 40001, 50001, 60001]
# With a +1 near to the colorbar label
cbar.set_ticklabels([str(i) for i in colorbar_labels])
if colorbar_orientation == 'horizontal':
for tick in cbar.ax.xaxis.get_ticklines():
tick.set_visible(False)
for t in cbar.ax.get_xticklabels():
t.set_fontsize(colorbar_fontsize)
t.set_fontweight(colorbar_fontweight)
else:
for tick in cbar.ax.yaxis.get_ticklines():
tick.set_visible(False)
for t in cbar.ax.get_yticklabels():
t.set_fontsize(colorbar_fontsize)
t.set_fontweight(colorbar_fontweight)
# Alternate text top and bottom on a horizontal colorbar if requested
# Use method described at:
# https://stackoverflow.com/questions/37161022/matplotlib-colorbar-
# alternating-top-bottom-labels
if colorbar_text_up_down or colorbar_text_down_up:
vmin=cbar.norm.vmin
vmax=cbar.norm.vmax
if cbar.extend=='min':
shift_l=0.05
scaling=0.95
elif cbar.extend=='max':
shift_l=0.
scaling=0.95
elif cbar.extend=='both':
shift_l=0.05
scaling=0.9
else:
shift_l=0.
scaling=1.0
#-------------Print bottom tick labels-------------
cbar.ax.set_xticklabels(lbot)
#--------------Print top tick labels--------------
for ii in ltop:
cbar.ax.text(shift_l + scaling*(ii-vmin)/(vmax-vmin),
1.5, str(ii), transform=cbar.ax.transAxes,
va='bottom', ha='center', fontsize=colorbar_fontsize,
fontweight=colorbar_fontweight)
for t in cbar.ax.get_xticklabels():
t.set_fontsize(colorbar_fontsize)
t.set_fontweight(colorbar_fontweight)
def map_title(title=None):
"""
| map_title is an internal routine to draw a title on a map plot
|
| title=None - title to put on map plot
|
|
|
|
|
|
"""
boundinglat=plotvars.boundinglat
lon_0=plotvars.lon_0
lonmin=plotvars.lonmin
lonmax=plotvars.lonmax
latmin=plotvars.latmin
latmax=plotvars.latmax
polar_range=90-abs(boundinglat)
if plotvars.proj == 'npstere':
mylon=lon_0+180
mylat=boundinglat-polar_range/15.0
proj=ccrs.NorthPolarStereo(central_longitude=lon_0)
lon_mid=lon_0
xpt, ypt = proj.transform_point(mylon, mylat, ccrs.PlateCarree())
if plotvars.proj == 'spstere':
mylon=lon_0
mylat=boundinglat+polar_range/15.0
proj=ccrs.SouthPolarStereo(central_longitude=lon_0)
lon_mid=lon_0
xpt, ypt = proj.transform_point(mylon, mylat, ccrs.PlateCarree())
if plotvars.proj == 'cyl':
lon_mid = lonmin + (lonmax - lonmin) / 2.0
proj=ccrs.PlateCarree(central_longitude=lon_mid)
mylon=lon_mid
mylat=latmax
xpt, ypt = proj.transform_point(mylon, mylat, ccrs.PlateCarree())
ypt=ypt+(latmax-latmin)/40.0
if plotvars.proj=='lcc':
lon_mid = lonmin + (lonmax - lonmin) / 2.0
proj=ccrs.LambertConformal(central_longitude=plotvars.lon_0,
central_latitude=plotvars.lat_0,
cutoff=plotvars.latmin)
mylon=lon_mid
mylat=latmax
xpt, ypt = proj.transform_point(0.0, mylat, ccrs.PlateCarree())
plotvars.mymap.text(xpt, ypt, title, va='bottom',
ha='center',
rotation='horizontal', rotation_mode='anchor',
fontsize=plotvars.title_fontsize,
fontweight=plotvars.title_fontweight)
def plot_map_axes(axes=None, xaxis=None, yaxis=None,
xticks=None, xticklabels=None,
yticks=None, yticklabels=None,
user_xlabel=None, user_ylabel=None,
verbose=None):
"""
| plot_map_axes is an internal routine to draw the axes on a map plot
|
| axes=None - Boolean for drawing axes
| xaxis=None - Boolean for drawing x-axis
| yaxis=None - Boolean for drawing x-axis
| xticks=None - user defined xticks
| xticklabels=None - user defined xtick labels
| yticks=None - user defined yticks
| yticklabels=None - user defined ytick labels
| user_xlabel=None - user defined xlabel
| user_ylabel=None - user defined ylabel
| verbose=None
|
|
|
|
|
"""
# Font definitions
title_fontsize = plotvars.title_fontsize
text_fontsize = plotvars.text_fontsize
axis_label_fontsize = plotvars.axis_label_fontsize
text_fontweight = plotvars.text_fontweight
title_fontweight = plotvars.title_fontweight
axis_label_fontweight = plotvars.axis_label_fontweight
# Map parameters
boundinglat=plotvars.boundinglat
lon_0=plotvars.lon_0
lonmin=plotvars.lonmin
lonmax=plotvars.lonmax
latmin=plotvars.latmin
latmax=plotvars.latmax
polar_range=90-abs(boundinglat)
# Cylindrical
if plotvars.proj == 'cyl':
if verbose:
print 'con - adding cylindrical axes'
lonticks, lonlabels = mapaxis(
min=plotvars.lonmin, max=plotvars.lonmax, type=1)
latticks, latlabels = mapaxis(
min=plotvars.latmin, max=plotvars.latmax, type=2)
if axes:
if xaxis:
if xticks is None:
axes_plot(xticks=lonticks, xticklabels=lonlabels)
else:
if xticklabels is None:
axes_plot(xticks=xticks, xticklabels=xticks)
else:
axes_plot(xticks=xticks, xticklabels=xticklabels)
if yaxis:
if yticks is None:
axes_plot(yticks=latticks, yticklabels=latlabels)
else:
if yticklabels is None:
axes_plot(yticks=yticks, yticklabels=yticks)
else:
axes_plot(yticks=yticks, yticklabels=yticklabels)
if user_xlabel is not None:
plot.text(0.5, -0.10, user_xlabel, va='bottom',
ha='center',
rotation='horizontal', rotation_mode='anchor',
transform=plotvars.mymap.transAxes,
fontsize=axis_label_fontsize,
fontweight=axis_label_fontweight)
if user_ylabel is not None:
plot.text(-0.05, 0.50, user_ylabel, va='bottom',
ha='center',
rotation='vertical', rotation_mode='anchor',
transform=plotvars.mymap.transAxes,
fontsize=axis_label_fontsize,
fontweight=axis_label_fontweight)
# Polar stereographic
if plotvars.proj == 'npstere' or plotvars.proj == 'spstere':
if verbose:
print 'con - adding stereographic axes'
mymap=plotvars.mymap
latrange = 90-abs(boundinglat)
latstep = 30
if 90 - abs(boundinglat) <= 50:
latstep = 10
proj=ccrs.Geodetic()
# Add
if axes:
if xaxis:
if yticks is None:
latvals=np.arange(5)*30-60
else:
latvals=np.array(yticks)
if plotvars.proj == 'npstere':
latvals=latvals[np.where(latvals>=boundinglat)]
else:
latvals=latvals[np.where(latvals<=boundinglat)]
for lat in latvals:
if abs(lat - boundinglat) > 1:
lons=np.arange(361)
lats=np.zeros(361)+lat
mymap.plot(lons, lats, color='k',
linewidth=1, linestyle='--',
transform=proj)
if yaxis:
if xticks is None:
lonvals=np.arange(7)*60
else:
lonvals=xticks
for lon in lonvals:
label = mapaxis(lon, lon, 1)[1][0]
if plotvars.proj == 'npstere':
lats=np.arange(90-boundinglat)+boundinglat
else:
lats=np.arange(boundinglat+91)-90
lons=np.zeros(np.size(lats))+lon
mymap.plot(lons, lats, color='k',
linewidth=1, linestyle='--',
transform=proj)
# Add longitude labels
if plotvars.proj == 'npstere':
proj=ccrs.NorthPolarStereo(central_longitude=lon_0)
pole=90
latpt = boundinglat - latrange/40.0
else:
proj=ccrs.SouthPolarStereo(central_longitude=lon_0)
pole=-90
latpt = boundinglat + latrange/40.0
lon_mid, lat_mid = proj.transform_point(0,pole, ccrs.PlateCarree())
if axis_label_fontsize > 0.0:
for xtick in lonvals:
label = mapaxis(xtick, xtick, 1)[1][0]
lonr, latr = proj.transform_point(xtick,latpt, ccrs.PlateCarree())
v_align='center'
if lonr < 1:
h_align = 'right'
if lonr > 1:
h_align = 'left'
if abs(lonr) <= 1:
h_align = 'center'
if latr < 1:
v_align = 'top'
if latr > 1:
v_align = 'bottom'
mymap.text(lonr, latr, label, horizontalalignment=h_align,
verticalalignment=v_align,
fontsize=axis_label_fontsize,
fontweight=axis_label_fontweight, zorder=101)
# Make the plot circular by blanking off around the plot
# Find min and max of plotting region in map coordinates
lons=np.arange(360)
lats=np.zeros(np.size(lons))+boundinglat
device_coords=proj.transform_points(ccrs.PlateCarree(), lons, lats)
xmin=np.min(device_coords[:,0])
xmax=np.max(device_coords[:,0])
ymin=np.min(device_coords[:,1])
ymax=np.max(device_coords[:,1])
# blank off data past the bounding latitude
pts=np.where(device_coords[:,0] >= 0.0)
xpts=np.append(device_coords[:,0][pts], np.zeros(np.size(pts))+xmax)
ypts=np.append(device_coords[:,1][pts], device_coords[:,1][pts][::-1])
mymap.fill(xpts, ypts, alpha=1.0, color='w', zorder=100)
xpts=np.append(np.zeros(np.size(pts))+xmin, -1.0* device_coords[:,0][pts])
ypts=np.append(device_coords[:,1][pts], device_coords[:,1][pts][::-1])
mymap.fill(xpts, ypts, alpha=1.0, color='w', zorder=100)
# Turn off map outside the cicular plot area
mymap.outline_patch.set_visible(False)
# Draw a line around the bounding latitude
lons=np.arange(361)
lats=np.zeros(np.size(lons))+boundinglat
device_coords=proj.transform_points(ccrs.PlateCarree(), lons, lats)
mymap.plot(device_coords[:,0], device_coords[:,1], color='k',
zorder=100, clip_on=False)
# Modify xlim and ylim values as the default values clip the plot slightly
xmax = np.max(np.abs(mymap.set_xlim(None)))
mymap.set_xlim((-xmax, xmax), emit=False)
ymax = np.max(np.abs(mymap.set_ylim(None)))
mymap.set_ylim((-ymax, ymax), emit=False)
# Lambert conformal
if plotvars.proj == 'lcc':
lon_0=plotvars.lonmin+(plotvars.lonmax-plotvars.lonmin)/2.0
lat_0=plotvars.latmin+(plotvars.latmax-plotvars.latmin)/2.0
mymap=plotvars.mymap
proj=ccrs.LambertConformal(central_longitude=lon_0,
central_latitude=lat_0,
cutoff=-40)
lonmin=plotvars.lonmin
lonmax=plotvars.lonmax
latmin=plotvars.latmin
latmax=plotvars.latmax
# Modify xlim and ylim values as the default values clip the plot slightly
xmin = mymap.set_xlim(None)[0]
xmax = mymap.set_xlim(None)[1]
ymin = mymap.set_ylim(None)[0]
ymax = mymap.set_ylim(None)[1]
mymap.set_ylim(ymin*1.05, ymax, emit=False)
mymap.set_ylim(None)
# Mask off contours that appear because of the plot extention
#mymap.add_patch(mpatches.Polygon([[xmin, ymin], [xmax,ymin],
# [xmax, ymin*1.05], [xmin, ymin*1.05]],
# facecolor='red'))
#transform=ccrs.PlateCarree()))
lons=np.arange(lonmax-lonmin+1)+lonmin
lats=np.arange(latmax-latmin+1)+latmin
nlons=np.size(lons)
nlats=np.size(lats)
verts=[]
for lat in lats:
verts.append([lonmin, lat])
for lon in lons:
verts.append([lon, latmax])
for lat in lats[::-1]:
verts.append([lonmax, lat])
for lon in lons[::-1]:
verts.append([lon, latmin])
# Mask left and right of plot
lats=np.arange(latmax-latmin+1)+latmin
lons=np.zeros(np.size(lats))+lonmin
npts=np.size(lats)
device_coords=proj.transform_points(ccrs.PlateCarree(), lons, lats)
xmin=np.min(device_coords[:,0])
xmax=np.max(device_coords[:,0])
ymin=np.min(device_coords[:,1])
ymax=np.max(device_coords[:,1])
# Left
mymap.fill([xmin, xmin, xmax, xmin],
[ymin, ymax, ymax, ymin],
alpha=1.0, color='w', zorder=100)
mymap.plot([xmin, xmax], [ymin, ymax], color='k', zorder=101, clip_on=False)
# Right
mymap.fill([-xmin, -xmin, -xmax, -xmin],
[ymin, ymax, ymax, ymin],
alpha=1.0, color='w', zorder=100)
mymap.plot([-xmin, -xmax], [ymin, ymax], color='k', zorder=101, clip_on=False)
# Upper
lons=np.arange(lonmax-lonmin+1)+lonmin
lats=np.zeros(np.size(lons))+latmax
device_coords=proj.transform_points(ccrs.PlateCarree(), lons, lats)
ymax=np.max(device_coords[:,1])
xpts=np.append(device_coords[:,0], device_coords[:,0][::-1])
ypts=np.append(device_coords[:,1], np.zeros(np.size(lons))+ymax)
mymap.fill(xpts, ypts, alpha=1.0, color='w', zorder=100)
mymap.plot(device_coords[:,0], device_coords[:,1], color='k', zorder=101, clip_on=False)
# Lower
lons=np.arange(lonmax-lonmin+1)+lonmin
lats=np.zeros(np.size(lons))+latmin
device_coords=proj.transform_points(ccrs.PlateCarree(), lons, lats)
ymin=np.min(device_coords[:,1])*1.05
xpts=np.append(device_coords[:,0], device_coords[:,0][::-1])
ypts=np.append(device_coords[:,1], np.zeros(np.size(lons))+ymin)
mymap.fill(xpts, ypts, alpha=1.0, color='w', zorder=100)
mymap.plot(device_coords[:,0], device_coords[:,1], color='k', zorder=101, clip_on=False)
# Turn off drawing of the rectangular box around the plot
mymap.outline_patch.set_visible(False)
# Draw longitudes and latitudes if requested
fs = plotvars.axis_label_fontsize
fw = plotvars.axis_label_fontweight
lw=1.0
if axes and xaxis:
if xticks is None:
map_xticks, map_xticklabels=mapaxis(min=plotvars.lonmin, max=plotvars.lonmax, type=1)
else:
map_xticks=xticks
if xticklabels is None:
map_xticklabels=xticks
else:
map_xticklabels=xticklabels
if axes and xaxis:
lats=np.arange(latmax-latmin+1)+latmin
for tick in np.arange(np.size(map_xticks)):
lons=np.zeros(np.size(lats))+map_xticks[tick]
device_coords=proj.transform_points(ccrs.PlateCarree(), lons, lats)
mymap.plot(device_coords[:,0], device_coords[:,1],
linewidth=lw, linestyle='--', color='k', zorder=101)
device_coords=proj.transform_point(map_xticks[tick], latmin-3, ccrs.PlateCarree())
mymap.text(device_coords[0], device_coords[1],
map_xticklabels[tick],
horizontalalignment='center',
fontsize=fs,
fontweight=fw,
zorder=101)
if yticks is None:
map_yticks, map_yticklabels=mapaxis(min=plotvars.latmin, max=plotvars.latmax, type=2)
else:
map_yticks=yticks
if yticklabels is None:
map_yticklabels=yticks
else:
map_yticklabels=yticklabels
if axes and yaxis:
lons=np.arange(lonmax-lonmin+1)+lonmin
for tick in np.arange(np.size(map_yticks)):
lats=np.zeros(np.size(lons))+map_yticks[tick]
device_coords=proj.transform_points(ccrs.PlateCarree(), lons, lats)
mymap.plot(device_coords[:,0], device_coords[:,1],
linewidth=lw, linestyle='--', color='k', zorder=101)
device_coords=proj.transform_point(lonmin-1, map_yticks[tick], ccrs.PlateCarree())
mymap.text(device_coords[0], device_coords[1],
map_yticklabels[tick],
horizontalalignment='right',
verticalalignment='center',
fontsize=fs,
fontweight=fw,
zorder=101)
device_coords=proj.transform_point(lonmax+1, map_yticks[tick], ccrs.PlateCarree())
mymap.text(device_coords[0], device_coords[1],
map_yticklabels[tick],
horizontalalignment='left',
verticalalignment='center',
fontsize=fs,
fontweight=fw,
zorder=101)
# UKCP grid
if plotvars.proj == 'UKCP' and plotvars.grid:
lonmin = -11
lonmax = 3
latmin = 49
latmax = 61
spacing = plotvars.grid_spacing
if plotvars.grid_lons is None:
lons = np.arange(30 / spacing + 1) * spacing
lons = np.append((lons*-1)[::-1], lons[1:])
else:
lons=plotvars.grid_lons
if plotvars.grid_lats is None:
lats = np.arange(90.0 / spacing + 1) * spacing
else:
lats=plotvars.grid_lats
plotvars.mymap.gridlines(color=plotvars.grid_colour,
linewidth=plotvars.grid_thickness,
linestyle=plotvars.grid_linestyle,
xlocs=lons, ylocs=lats)