"""
Helpers for plotting of maps etc
--------------------------------
.. note::
BETA stage, code might undergo revisions
"""
import numpy as np
import matplotlib.pyplot as plt
from cartopy import crs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from geonum.helpers import order_of_magnitude
def _get_tick_formatter(value):
"""
Get tick formatting string for map plots
Parameters
----------
value : float
value used to estimate formatter (e.g. one of the longitude
coordinates).
Returns
-------
str
formatting string
"""
digits = 2 - order_of_magnitude(value)
digits = 0 if digits < 0 else digits
return f'.{digits:d}f'
def _calculate_map_ticks(ax):
"""
Calculate longitude (x) and latitude (y) ticks for map plots.
Parameters
----------
ax : cartopy.mpl.geoaxes.GeoAxes
map axes instances.
Returns
-------
tuple
2-element tuple containing x and y ticks
"""
lonleft, lonright = ax.get_xlim()
num_lonticks = 7 if lonleft == -lonright else 6
xticks = np.linspace(lonleft, lonright, num_lonticks)
latleft, latright = ax.get_ylim()
num_latticks = 7 if latleft == - latright else 6
yticks = np.linspace(latleft, latright, num_latticks)
return (xticks, yticks)
[docs]
def set_map_ticks(ax, xticks=None, yticks=None, tick_format=None):
"""Set or update ticks in instance of GeoAxes object (cartopy)
Note
----
The input GeoAxes must be setup with :class:`cartopy.crs.PlateCarree`
projection.
Parameters
----------
ax : cartopy.GeoAxes
map axes instance (needs to be set up with PlateCarree projection).
xticks : iterable, optional
ticks of x-axis (longitudes). If None, it is retrieved automatically
based on xlim of input axes. Defaults to None.
yticks : iterable, optional
ticks of y-axis (latitudes). If None (or if xticks is None),
it is retrieved automatically based on ylim of input axes. Defaults
to None.
tick_formatter : str, optional
string formatter for axes labels in the plot.
Returns
-------
cartopy.GeoAxes
modified axes instance
"""
if not isinstance(ax.projection, crs.PlateCarree):
raise AttributeError('Map ticks can only be added to GeoAxes with '
'PlateCarree projection.')
if xticks is None:
xticks, yticks = _calculate_map_ticks(ax)
if tick_format is None:
val = np.mean(xticks)
if val == 0:
val = xticks[next((i for i, x in enumerate(xticks) if x), None)]
if val != 0:
tick_format = _get_tick_formatter(val)
else:
tick_format = '.1f'
ax.set_xticks(xticks, crs=crs.PlateCarree())
ax.set_yticks(yticks, crs=crs.PlateCarree())
lon_formatter = LongitudeFormatter(number_format=tick_format,
degree_symbol='°',
dateline_direction_label=True)
lat_formatter = LatitudeFormatter(number_format=tick_format,
degree_symbol='°')
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
return ax
def _infer_text_location(ax,p):
info = dict(ha='center')
left,right = ax.get_xlim()
dx = right - left
bottom,top = ax.get_ylim()
dy = top - bottom
lat, lon = p.latitude, p.longitude
distleft = abs(left - lon) / dx
if distleft < 0.3:
info['ha'] = 'left'
elif distleft > 0.7:
info['ha'] = 'right'
yoffs = dy*0.08 # put text 5% of range above or below the points latitude
distbottom = abs(bottom - lat) / dx
if distbottom > 0.6: # point is in upper 40% of map, put text below
yoffs *= -1
info['xytext'] = (lon, lat+yoffs)
return info
[docs]
def plot_topo_contourf(ax, topo, oceans_separate=True, levels=50, cmap=None):
if cmap is None:
cmap = 'Oranges'
# Create meshgrid for domain for contour plot
X,Y = topo.init_mesh()
# Separate land areas from sea areas for plotting below
if oceans_separate:
seamask = np.ma.masked_less_equal(topo.data, 0)
data = np.ma.MaskedArray(topo.data, mask=seamask.mask)
if np.any(seamask):
seadata = np.ma.MaskedArray(topo.data, mask=~seamask.mask)
# Plot sea in a light blue
ax.contourf(X, Y, seadata, 50,
colors='#e6f7ff')
else:
data = topo.data
# Plot topographic land data
pdata = ax.contourf(X, Y, data, levels=levels,
cmap=cmap)
return ax, pdata
[docs]
def plot_geopoint_into_map(ax, p, annotate=True, annot_kwargs=None, **kwargs):
if annot_kwargs is None:
annot_kwargs = {}
ax.scatter(x=[p.longitude],y=[p.latitude], **kwargs)
if annotate:
annot_loc = _infer_text_location(ax,p)
try:
c = kwargs['color']
except KeyError:
c = 'k'
annot = dict(
text = p.name,
xy = (p.longitude,
p.latitude), # location of observatory in plot
arrowprops = dict(color=c, lw=1, arrowstyle='->', shrinkB=4),
size=7, color=c, fontweight='bold', **annot_loc
)
annot.update(annot_kwargs)
ax.annotate(**annot)
return ax
[docs]
def plot_geovector3d_into_map(ax, vec, **kwargs):
if vec.anchor is None:
raise AttributeError(
f'input vector {vec} does not have an anchor point assigned. '
f'Please set anchor location using vec.add_anchor')
start = vec.anchor
end = start + vec
ax.plot([start.longitude, end.longitude],
[start.latitude, end.latitude], **kwargs)
return ax
[docs]
def rotate_xtick_labels(ax, deg=30, ha="right"):
"""Rotate xtick labels in matplotlib axes object"""
lbls = ax.get_xticklabels()
lbls = [lbl.get_text() for lbl in lbls]
ax.set_xticklabels(lbls, rotation=deg, ha=ha)
return ax
[docs]
def shifted_color_map(vmin, vmax, cmap=None):
"""Shift center of a diverging colormap to value 0
Note
----
This method was found `here <http://stackoverflow.com/questions/
7404116/defining-the-midpoint-of-a-colormap-in-matplotlib>`__
(last access: 17/01/2017). Thanks to `Paul H <http://stackoverflow.com/
users/1552748/paul-h>`_ who provided it.
Function to offset the "center" of a colormap. Useful for
data with a negative min and positive max and if you want the
middle of the colormap's dynamic range to be at zero level
Parameters
----------
vmin : float
lower end of data value range
vmax : float
upper end of data value range
cmap
matplotlib colormap to be shifted (if None, use "seismic")
Returns
-------
matplotlib.colors.LinearSegmentedColormap
shifted colormap
"""
import matplotlib.cm as colormaps
import matplotlib.colors as colors
if cmap is None:
cmap = colormaps.seismic
midpoint = 1 - np.abs(vmax) / (np.abs(vmax) + np.abs(vmin))
cdict = {
'red': [],
'green': [],
'blue': [],
'alpha': []
}
# regular index to compute the colors
reg_index = np.linspace(0, 1, 257)
# shifted index to match the data
shift_index = np.hstack([
np.linspace(0.0, midpoint, 128, endpoint=False),
np.linspace(midpoint, 1.0, 129, endpoint=True)
])
for ri, si in zip(reg_index, shift_index):
r, g, b, a = cmap(ri)
cdict['red'].append((si, r, r))
cdict['green'].append((si, g, g))
cdict['blue'].append((si, b, b))
cdict['alpha'].append((si, a, a))
return colors.LinearSegmentedColormap('shiftedcmap', cdict)