# -*- coding: utf-8 -*-
r"""
Extensions for matplotlib's pyplot with plot types typically
needed for particle data.
Author: Martin Wieser
Module: irfpy.ica.pyplot
Various extensions for matplotlib's pyplot with plotting
functions suitable for plotting of ICA data.
"""
import matplotlib.pylab as plt
import matplotlib as mpl
from matplotlib.pyplot import * # noqa: F401,F403
# this makes it possible to share the namespace of pyplot
import copy
import matplotlib.dates as mdates
from matplotlib.ticker import FuncFormatter, MaxNLocator
from mpl_toolkits import axes_grid1
# from matplotlib.projections.polar import PolarAxes
# import matplotlib.axes as maxes
try:
if __name__ == '__main__':
import tools as icatools
import colormaps as icacolors
else:
from . import tools as icatools
from . import colormaps as icacolors
except:
import irfpy.ica.tools as icatools
import irfpy.ica.colormaps as icacolors
import numpy as np
import gc
[docs]def setfont(font="Helvetica", unicode=True):
r"""
Set Matplotlibs rcParams to use LaTeX for font rendering.
Calling this makes sure there will be identical looking
fonts independent of the output device and wether latex
is used or not.
Revert all changes by calling rcdefault() from matplotlib.
Parameters:
font (string):
"Helvetica"
"Times"
"Computer Modern"
unicode (boolean):
Use unicode. Default: False.
"""
# Use TeX for all figure text!
plt.rc('text', usetex=True)
font = font.lower().replace(" ", "")
if font == 'times':
# Times
font = {'family': 'serif', 'serif': ['Times']}
preamble = r"""
\usepackage{color}
\usepackage{mathptmx}
"""
elif font == 'helvetica':
# Helvetica
# set serif, too. Otherwise setting to times and then
# Helvetica causes an error.
font = {'family': 'sans-serif', 'sans-serif': ['Helvetica'],
'serif': ['cm10']}
preamble = r"""
\usepackage{color}
\usepackage[tx]{sfmath}
\usepackage{helvet}
"""
else:
# Computer modern serif
font = {'family': 'serif', 'serif': ['cm10']}
preamble = r"""
\usepackage{color}
"""
if unicode:
# Unicode for Tex
# preamble = r"""\usepackage[utf8]{inputenc}""" + preamble
# inputenc should be set automatically
plt.rcParams['text.latex.unicode'] = True
# print font, preamble
plt.rc('font', **font)
plt.rcParams['text.latex.preamble'] = preamble
[docs]def removefig(fig=None):
r'''
Removes fig or the current figure reference from pyplot and garbage collects
it (effectively removes it from memory).
'''
if fig is None:
plt.close(plt.gcf())
else:
plt.close(fig)
gc.collect()
[docs]def removeallfigs():
r'''
Removes all figures from pyplot and garbage collects
it (effectively removes it from memory).
'''
for f in plt.get_fignums():
removefig(fig=f)
#def format_exponenta(ax, axis='y'):
# r'''
# Change the ticklabel format to scientific format
# https://stackoverflow.com/questions/31517156/adjust-exponent-text-after-setting-scientific-limits-on-matplotlib-axis::
#
# plt.format_exponent(ax,axis='y')
#
# '''
#
# ax.ticklabel_format(axis=axis, style='sci', scilimits=(-2, 2))
#
# # Get the appropriate axis
# if axis == 'y':
# ax_axis = ax.yaxis
# x_pos = 0.0
# y_pos = 1.0
# horizontalalignment='left'
# verticalalignment='bottom'
# else:
# ax_axis = ax.xaxis
# x_pos = 1.0
# y_pos = -0.05
# horizontalalignment='right'
# verticalalignment='top'
#
# # Run plt.tight_layout() because otherwise the offset text doesn't update
# plt.tight_layout()
# ##### THIS IS A BUG
# ##### Well, at least it's sub-optimal because you might not
# ##### want to use tight_layout(). If anyone has a better way of
# ##### ensuring the offset text is updated appropriately
# ##### please comment!
#
# # Get the offset value
# offset = ax_axis.get_offset_text().get_text()
#
# if len(offset) > 0:
# # Get that exponent value and change it into latex format
# minus_sign = u'\u2212'
# expo = np.float(offset.replace(minus_sign, '-').split('e')[-1])
# offset_text = r'x$\mathregular{10^{%d}}$' %expo
#
# # Turn off the offset text that's calculated automatically
# ax_axis.offsetText.set_visible(False)
#
# # Add in a text box at the top of the y axis
# ax.text(x_pos, y_pos, offset_text, transform=ax.transAxes,
# horizontalalignment=horizontalalignment,
# verticalalignment=verticalalignment)
# return ax
[docs]def add_colorbar(im, aspect=20, pad_fraction=0.5, **kwargs):
r"""Add a vertical color bar to an image plot with the height
of the colorbar equal the height of the plot. Otherwise works the same way
as plt.colorbar().::
im = plt.imshow(np.arange(200).reshape((20, 10)))
plt.add_colorbar(im)
"""
divider = axes_grid1.make_axes_locatable(im.axes)
width = axes_grid1.axes_size.AxesY(im.axes, aspect=1./aspect)
pad = axes_grid1.axes_size.Fraction(pad_fraction, width)
current_ax = plt.gca()
cax = divider.append_axes("right", size=width, pad=pad)
plt.sca(current_ax)
return im.axes.figure.colorbar(im, cax=cax, **kwargs)
[docs]def legend_ordered(order, **kwargs):
r"""Adds a legend with the element order given by the parameter order.
The length of the list order must correspond to the number of
legend entries n and the elements must be consecutievly numbered
from 0 to n-1. Otherwiese works the same as plt.legend()
plt.legend_ordered([2,0,1])
"""
#get handles and labels
handles, labels = plt.gca().get_legend_handles_labels()
#add legend to plot
return plt.legend([handles[idx] for idx in order],[labels[idx] for idx in order], **kwargs)
[docs]def plotbinnedmatrix(*args, **kwargs):
r'''
Plots a color matrix (also known as a quadrilateral mesh)
using matplotlib.pyplot.pcolormesh(). It takes the similar parameters
as pcolormesh(). In addition to pcolormesh(), this function makes sure
the lables are centered on each bin on both axes and that all bins
are plotted, including the last one.
Call signatures::
plotbinnedmatrix(C)
plotbinnedmatrix(X, Y, C)
plotbinnedmatrix(C, **kwargs)
plotbinnedmatrix(X, Y, C, **kwargs)
If X and Y are given, then the values in them are used as axis labels.
The shape of C must be (Y,X), otherwise an error is generated. C can
be a numpy.masked array.
Example: Plot a 16x16 bin matrix::
import numpy as np
import matplotlib.pyplot as plt
import irfpy.ica.plot as icaplt
C=np.arange(128).reshape((16,8)) # (y,x)
icaplt.plotbinnedmatrix(np.arange(8),np.arange(16), C)
plt.title('A binned 8x16 matrix')
Or also shorter::
plotbinnedmatrix(C)
'''
# make my default colormap
# cmap = plt.get_cmap('viridis')
# cmap = plt.get_cmap('jet')
# cmap = icacolors.comet
mycmap = kwargs.pop('cmap', None)
if mycmap is None:
mycmap = copy.copy(icacolors.viridis)
mycmap.set_bad(color=(0.2, 0.2, 0.2), alpha=1.0)
funcname = 'irfpy.ica.pyplot.plotbinnedmatrix'
if len(args) == 1:
X = np.arange(np.shape(args[0])[1])
Y = np.arange(np.shape(args[0])[0])
C = args[0]
elif len(args) == 3:
X = args[0]
Y = args[1]
C = args[2]
else:
raise TypeError(
'Incompatible X, Y or C inputs to %s; see help(%s)' % (
funcname, funcname))
# more error checking:
if len(np.shape(X)) != 1:
raise TypeError(
'X must be one-dimensional but is of shape ' +
str(np.shape(X))
)
if len(np.shape(Y)) != 1:
raise TypeError(
'Y must be one-dimensional but is of shape ' +
str(np.shape(Y))
)
if np.shape(C)[1] != np.shape(X)[0]:
raise TypeError(
'Dimension of C (y,x)=' + str(np.shape(C)) +
' is incompatible with the dimension of X = ' +
str(np.shape(X)[0]) + ', it should be ' +
str(np.shape(C)[1])
)
if np.shape(C)[0] != np.shape(Y)[0]:
raise TypeError(
'Dimension of C (y,x)=' + str(np.shape(C)) +
' is incompatible with the dimension of Y = ' +
str(np.shape(Y)[0]) + ', it should be ' +
str(np.shape(C)[0])
)
def format_fnx(val, pos):
k = int(round(val))
if k >= 0 and k < X.size:
return str(X[k])
else:
return ''
def format_fny(val, pos):
k = int(round(val))
if k >= 0 and k < Y.size:
return str(Y[k])
else:
return ''
r = plt.pcolormesh(np.arange(len(X) + 1) - 0.5,
np.arange(len(Y) + 1) - 0.5,
C,
cmap=mycmap, **kwargs)
ax = plt.gca()
ax.xaxis.set_major_formatter(FuncFormatter(format_fnx))
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
ax.yaxis.set_major_formatter(FuncFormatter(format_fny))
ax.yaxis.set_major_locator(MaxNLocator(integer=True))
plt.xlim(-0.5, len(X) - 0.5)
plt.ylim(-0.5, len(Y) - 0.5)
return r
[docs]def plotenergytime(Time, Energy, Data, **kwargs):
r'''
Plots a rainbow energy time spectrogram using matplotlib.pyplot.pcolormesh().
It takes the same parameters as pcolormesh except for:
Time:
The X coordinate, a numpy array of datetime objects as in time_instances.
Alternatively, a numpy array of Matplotlib dates can be used.
Energy:
The Y coordinate, centers of energy bins, sorted in ascending order
Energy may be filled with NaN above a certain index up to index 95.
Internally the energy centers are converted to energy intervals using
irfpy.ica.tools.energyintervals(Energy).
Energy may be of the shape (nE,) or (nE, nT) with nE the number of
energy bins and nT the number of timeinstances. If the shape (nE,nT)
is used, the each timeinstance can have a different energy vector.
This is e.g. useful if the energy scale is compensated for the
spacecraft potential.
Data:
The C values of the shape (Energy,Time).
linear=False:
If set to True, the interval boundaries are constructed using
arithmetic averages of neighboring energy centers, otherwise
geometric intervall boundaries are used.
rectangular=False:
If set to True, the data is plotted as rectangular elements.
dateformat='':
Date formatting string, as passed to matplotlib.dates.DateFormatter()
If set to None, then no date formatting takes place and the x-axis is
in units of days.
maxtimebinwidth=None:
Maximum width in seconds of the rectangle of data associated with
one timeinstance.
Call signatures::
plotenergytime(Time, Energy, Data)
plotenergytime(Time, Energy, Data, **kwargs)
A complete example:
Note the use of norm=LogNorm(...) to get a logarithmic color scale.
The data itself does not need to be processed with np.log10().::
import irfpy.ica.io as icaio
import matplotlib as mpl
import irfpy.ica.pyplot as plt
import irfpy.ica.tools as icatools
from irfpy.util.irfpyrc import Rc
rc = Rc()
path = rc.get('ica', 'dataPath')
starttime="20150605T121500"
stoptime="20150605T134000"
mat = icaio.readproc(path,starttime,stoptime,
variables=['orig_ionspectra',
'cur_pacc',
'mode',
'sw_version',
'E',
'time_instances',
'sum_orig_ionspectra']
)
orig_ionspectra=mat['orig_ionspectra']
sw_version=mat['sw_version']
time_instances=mat['time_instances']
D = mat['sum_orig_ionspectra']
plt.clf()
for start,stop,version in icatools.timeslices(sw_version):
# This loop allows different sw_versions in one plot
# get the expected length in time of one measuremnt in seconds
# +0.1s is added to avoid plotting artefacts
deltaT = icatools.deltaTfromversion(version)+0.1
# get an energy vector for each time instance
E = icatools.getE(mat['E'],sw_version[start:stop])
# epsilon is added to D to allow to plot the value 0.0 in log-scale:
epsilon = 1e-10
plt.plotenergytime(time_instances[start:stop],
E,
D[:,start:stop]+epsilon,
norm=mpl.colors.LogNorm(vmin=0.2, vmax=10000.0),
dateformat='%H:%M', # optional
maxtimebinwidth=deltaT, # optional
cmap='viridis') # optional
plt.xlim(time_instances[0],time_instances[-1])
plt.colorbar(label='Counts')
plt.yscale('log')
plt.ylim(3,4e4)
plt.ylabel('E (eV)')
plt.title('ICA Energy-Time Spectrogram')
plt.show()
'''
# make my default colormap
# cmap = plt.get_cmap('viridis')
# cmap = plt.get_cmap('jet')
# cmap = icacolors.comet
mycmap = kwargs.pop('cmap', None)
if mycmap is None:
mycmap = copy.copy(icacolors.viridis)
mycmap.set_bad(color=(0.2, 0.2, 0.2), alpha=1.0)
myedgecolors = kwargs.pop('edgecolors', 'None')
mydateformat = kwargs.pop('dateformat', '')
mymaxtimebinwidth = kwargs.pop('maxtimebinwidth', -1.0)
makerectangular = kwargs.pop('rectangular', False)
# #find largest element in E that is not nan
# Emax = np.where(np.isnan(Energy)==False)[0][-1] +1
# if Emax < 1:
# funcname = 'irfpy.ica.pyplot.plotenergytime'
# raise TypeError(
# '%s: Energy does not contain enough non-NaN elements '+
# '(has: %d, but >1 required); see help(%s)' % (
# funcname, Emax, 'irfpy.ica.tools.energyintervals'))
#
linear = kwargs.pop('linear', False)
#print('DEBUG: linear='+str(linear))
Einterval = icatools.energyintervals(Energy, linear=linear)
#print('DEBUG: Einterval='+str(Einterval))
Emax = Einterval.shape[0]-1
if Time.dtype != float:
T2 = mdates.date2num(Time)
else:
T2 = Time.copy()
if np.shape(T2)[0] == np.shape(Data[0:Emax, :])[1]:
if len(T2) > 1:
T2 = np.append(T2, [T2[-1] + (T2[-1] - T2[-2])])
elif len(T2) == 1:
if mymaxtimebinwidth > 0.0:
T2 = np.append(T2, [T2[0] + mymaxtimebinwidth / 86400.0])
# assume mymaxtimebinwidth
else:
# assume arbitrarily 12 seconds as we do not know
# how wide the time bin should be
T2 = np.append(T2, [T2[0] + 12.0 / 86400.0])
if len(Einterval.shape) == 2:
# add energies for the extra time instance
Einterval = np.concatenate((Einterval, Einterval[:, -1][:, None]), axis=1)
# Check if somewhere the time difference is larger than maxtimebinwith
# and if so insert a dummy time with a NaN data vector
if mymaxtimebinwidth > 0.0:
T2diff = np.diff(T2) * 86400.0
longpos = np.where(T2diff > mymaxtimebinwidth)[0]
T2 = np.insert(T2, longpos + 1, T2[longpos] + mymaxtimebinwidth / 86400.0)
NanData = np.zeros_like(Data[:, longpos]) * np.nan
Data = np.insert(Data, longpos + 1, NanData, axis=1)
if len(Einterval.shape) == 2:
# add energies for the extra time instances if Einterval is 2D
Einterval = np.insert(Einterval, longpos + 1,
Einterval[:, longpos], axis=1)
thedata = np.ma.masked_invalid(Data[0:Emax, :])
if len(np.shape(T2)) == 1 and len(Einterval.shape) == 2:
# now we have a 1D time axis but a 2D energy axis.
# the time axis needs to be expanded to match
T2 = np.tile(T2, (Einterval.shape[0], 1))
if makerectangular:
# duplicate elements such that the plotted mesh-elements
# get rectangular
thedata = np.repeat(thedata, 2, axis=1)
Einterval = np.repeat(Einterval, 2, axis=1)[:,:-1]
T2 = np.repeat(T2, 2, axis=1)[:,1:]
T2[:,1::2] -= 0.001*mymaxtimebinwidth/86400.0
#thedata[:,1::2] = np.ma.masked
r = plt.pcolormesh(T2,
Einterval,
thedata,
cmap=mycmap,
edgecolors=myedgecolors,
**kwargs)
if mydateformat != None:
if mydateformat != '':
hms = mdates.DateFormatter(mydateformat)
plt.gca().xaxis.set_major_formatter(hms)
plt.gca().xaxis_date()
return r
[docs]def plottime(*args, **kwargs):
r'''
Plots a lineplot with time on the xaxis..
It takes the same parameters as plot() except for:
Time:
The X coordinate, a numpy array of datetime objects as in time_instances.
Alternatively, a numpy array of Matplotlib dates can be used.
Data:
The Y values of the shape (Time).
Other parameters
dateformat='':
Date formatting string, as passed to matplotlib.dates.DateFormatter()
If set to None, then no date formatting takes place and the x-axis is
in units of days.
Call signatures::
plottime(Time, Data)
plottime(Time, Data, *args, **kwargs)
'''
args = list(args)
Time = args.pop(0)
mydateformat = kwargs.pop('dateformat', '')
if Time.dtype != float:
mdat = mdates.date2num(Time)
else:
mdat = Time
r = plt.plot(mdat, *args, **kwargs)
if mydateformat != None:
if mydateformat != '':
hms = mdates.DateFormatter(mydateformat)
plt.gca().xaxis.set_major_formatter(hms)
plt.gca().xaxis_date()
return r
[docs]def filltime_between(*args, **kwargs):
r'''
Plots a fill plot with time on the xaxis..
It takes the same parameters as filltime() except for:
Time:
the X coordinate, a numpy array of datetime objects as in time_instances
Data:
the Y values of the shape (Time).
Other parameters
dateformat='':
Date formatting string, as passed to matplotlib.dates.DateFormatter()
Call signatures
baseline = -10
filltime_between(Time, Data, baseline,where=Data>baseline, interpolate=True)
'''
args = list(args)
Time = args.pop(0)
mydateformat = kwargs.pop('dateformat', '')
if Time.dtype != float:
mdat = mdates.date2num(Time)
else:
mdat = Time
r = plt.fill_between(mdat, *args, **kwargs)
if mydateformat != None:
if mydateformat != '':
hms = mdates.DateFormatter(mydateformat)
plt.gca().xaxis.set_major_formatter(hms)
plt.gca().xaxis_date()
return r
[docs]def filltime(*args, **kwargs):
r'''
Plots a fill plot with time on the xaxis..
It takes the same parameters as fill() except for:
Time:
the X coordinate, a numpy array of datetime objects as in time_instances
Data:
the Y values of the shape (Time).
Other parameters
dateformat='':
Date formatting string, as passed to matplotlib.dates.DateFormatter()
Call signatures
filltime(Time, Data)
filltime(Time, Data, *args, **kwargs)
'''
args = list(args)
Time = args.pop(0)
mydateformat = kwargs.pop('dateformat', '')
if Time.dtype != float:
mdat = mdates.date2num(Time)
else:
mdat = Time
r = plt.fill(mdat, *args, **kwargs)
if mydateformat != None:
if mydateformat != '':
hms = mdates.DateFormatter(mydateformat)
plt.gca().xaxis.set_major_formatter(hms)
plt.gca().xaxis_date()
return r
#%%
[docs]def xlimtime(t1,t2):
r""" Same as plt.xlim but accepts datetime objects that are then converted
using matplotlib.datesdate2num() to ordinary plt.xlim(...) parameters."""
plt.xlim(mdates.date2num(t1),mdates.date2num(t2))
#%% Control plot boundaries
[docs]def copy_plot_width(ax_src, ax_dest):
r'''
copies the wifth of the plot using ax_src to the plot using ax_dest
'''
ll, bb, _, hh = ax_dest.get_position().bounds
_, _, ww, _ = ax_src.get_position().bounds
ax_dest.set_position([ll, bb, ww, hh])
[docs]def set_plot_left(ax_dest,left):
r'''
sets the left position of a plot while keeping the right side
'''
ll, bb, ww, hh = ax_dest.get_position().bounds
ax_dest.set_position([left, bb, ww-(left-ll), hh])
[docs]def copy_plot_height(ax_src, ax_dest):
r'''
copies the wifth of the plot using ax_src to the plot using ax_dest
'''
ll, bb, ww, _ = ax_dest.get_position().bounds
_, _, _, hh = ax_src.get_position().bounds
ax_dest.set_position([ll, bb, ww, hh])
[docs]def copy_plot_position(ax_src, ax_dest):
r'''
copies the position of the plot using ax_src to the plot using ax_dest
'''
ll, bb, ww, hh = ax_src.get_position().bounds
ax_dest.set_position([ll, bb, ww, hh])
#%% Circular plots
[docs]class radialscaleobject:
def __init__(self):
self.axis = None
# self.label=None
self.tick_labels = None
[docs] def get_ticklabels(self):
if self.tick_labels is None:
return list()
else:
return list(self.tick_labels)
[docs]def radialscale(ax=None, pos=0.1, labelpos=0.23, label="", location='left', linewidth=1, **kwargs):
r'''
Plots a radial scale besides a of a plot made with plotsectorpolar().
Call signatures::
radialscale(label="Energy [eV]")
pos:
Horizontal pos where the scale is drawn in axis units of
plotsectorpolar radii (0.0 = left or right axis border. Default is 0.1
which draws the scale outside of the circle of plotsectorpolar()
on the left.
labelpos:
Horizontal pos where the scale label is drawn in axis units
Default is 0.23 which draws the axis labels on the left outside
of the circle of plotsectorpolar() on the
outside of the radial scale bar.
label:
The label text to be shown. The value is assigned to the
yaxis-label of the plot, and the label is moved to an suitable position
defined by labelpos.
location='left':
The axis is drawn on the left side of the polar plot. 'right' the
axis is drawn on the right side.
returns: axis of the scalebar.
Example:
See plotsectorpolar()
'''
if ax is None:
ax = plt.gca()
plt.gcf().canvas.draw() # this renders the tickmark labels
(minradius, maxradius) = ax.get_ylim()
linear = plt.gca().get_yscale() == 'linear'
rrr = list()
aaa = list()
ttt = list()
locs = [minradius] + [x for x in ax.yaxis.get_ticklocs()]
labl = [""] + [t.get_text() for t in ax.yaxis.get_ticklabels()]
if maxradius < locs[-1]:
locs[-1] = maxradius
labl[-1] = ""
elif maxradius > locs[-1]:
locs.append(maxradius)
labl.append("")
if location == 'right':
ax_data = (1+pos, 0.5)
labelpos=1+labelpos
else:
ax_data = (0+pos, 0.5)
labelpos=0+labelpos
disp_coords = plt.gca().transAxes.transform(ax_data)
data_coords = plt.gca().transData.inverted().transform(disp_coords)
# print(data_coords)
if linear:
ppos = -data_coords[1] # /maxradius
else:
ppos = -np.log10(data_coords[1])
if location == 'right':
ppos = -ppos
# print(ppos)
for r, t in zip(locs, labl):
if linear:
y1 = r
x1 = ppos # *maxradius
rr1 = np.sqrt(x1 * x1 + y1 * y1)
theta1 = np.arctan2(y1, x1)
else:
y1 = np.log10(r)
x1 = ppos # *np.log10(maxradius)
rr1 = 10**np.sqrt(x1 * x1 + y1 * y1)
theta1 = np.arctan2(y1, x1)
rrr.append(rr1)
aaa.append(theta1)
ttt.append(" "+ t + " ")
ax.tick_params(axis='y', labelsize=0.001)
# hides the tick labels without removing the label text
# plot the 'spine'
spine = ax.plot(aaa, rrr, '-_k', clip_on=False, linewidth=linewidth)
spine[0].set_linewidth(linewidth)
spine[0].set_markeredgewidth(linewidth)
sticklabels = list()
for ii in range(len(aaa)):
if location=='right':
sticklabels.append(ax.text(aaa[ii], rrr[ii], str(ttt[ii]),
ha="left", va="center",
clip_on=False, **kwargs))
else:
sticklabels.append(ax.text(aaa[ii], rrr[ii], str(ttt[ii]),
ha="right", va="center",
clip_on=False, **kwargs))
# plot the 'axis label'
plt.gca().yaxis.set_label_coords(labelpos, 0.75)
plt.ylabel(str(label))
# slabel= ax.text(thetal, rrl, str(label), ha='right', va='center',
# clip_on=False, rotation=90, **kwargs)
rax = radialscaleobject()
rax.axis = spine[0]
# rax.label= slabel
rax.tick_labels = sticklabels
return rax
[docs]def plotsectorpolar(*args, **kwargs):
r"""
Plots a color matrix as polar plot with matplotlib.pyplot.pcolormesh()
using a quadrilateral mesh. It takes the similar parameters
as pcolormesh(). In addition to pcolormesh(), this function makes sure
the lables are centered on each bin on both axes and that all bins
are plotted, including the last one.
Call signatures::
plotsectorpolar(C)
plotsectorpolar(X, Y, C)
plotsectorpolar(C, **kwargs)
plotsectorpolar(X, Y, C, **kwargs)
If X and Y are given, then the values in them are used as axis labels.
The shape of C must be (Y,X), otherwise an error is generated. C can
be a numpy.masked array.
Allowed for kwargs is everything that pcolormesh accepts plus the ones
listed below:
linear:
True : Energy intervals are interpolated in a linear fashion
False : Energy intervals are interpolated in a geometric fashion
zeroangle:
The azimuth angle that should be plotted in the +X direction of the
plot in degrees.
cmap:
A color map.
Plot a polar sector plot::
import numpy as np
import irfpy.ica.io as icaio
import matplotlib as mpl
import irfpy.ica.pyplot as plt
import irfpy.ica.tools as icatools
from irfpy.util.irfpyrc import Rc
rc = Rc()
path = rc.get('icadds', 'dataroot')
starttime="20150605T121500"
stoptime="20150605T134000"
cmin = 0.1
cmax = 1e4
def sumorigionspacc(orig_ionspectra,mode,version):
epsilon=1e-10
D3 = np.zeros((16, 96))
for start,stop,(vm,vv) in icatools.timeslices(zip(mode,version)):
D3 = D3 + np.nansum(np.nansum(orig_ionspectra[:,:,:,start:stop],axis=2),axis=2)
# sect x energy
D3[D3 < epsilon] = epsilon
return D3.T
mat = icaio.readproc(path,starttime,stoptime,
variables=['orig_ionspectra',
'cur_pacc','mode','sw_version','E'])
orig_ionspectra=mat['orig_ionspectra']
mode=mat['mode']
sw_version=mat['sw_version']
E = icatools.getE(mat['E'],sw_version)
D = sumorigionspacc(orig_ionspectra,mode,sw_version)
A=np.array(['S'+str(n) for n in np.arange(16)])
plt.clf()
plt.subplot(1,2,1,projection='polar')
pax = plt.plotsectorpolar(A,E[:,0],D,
norm=mpl.colors.LogNorm(vmin=cmin, vmax=cmax))
plt.colorbar(label='Counts')
plt.title('Sector/Energy linear')
plt.yscale("linear")
# Example: remove the radial labels:
# This hides the tick labels without removing the label text
pax.tick_params(axis='y', labelsize=0)
plt.subplot(1,2,2,projection='polar',axisbg='#333333')
pax=plt.plotsectorpolar(A,E[:,0],D,
norm=mpl.colors.LogNorm(vmin=cmin, vmax=cmax))
plt.title('Sector/Energy log')
cbar = plt.colorbar(label='Counts')
rax = plt.radialscale(label='Energy [eV]',linewidth=2)
# Example: change font sizes everywhere to 12:
for item in ([pax.title, # the title
pax.xaxis.label, # azimuth axis label
pax.yaxis.label, # radial axis label
cbar.ax.yaxis.label] + # color bar label
pax.get_xticklabels() + # azimuth tick labels
rax.get_ticklabels() + # radial scale tick labels
cbar.ax.yaxis.get_ticklabels() # color bar tick labels
):
item.set_fontsize(12)
plt.tight_layout()
plt.show()
Or also shorter::
plotsectorpolar(C)
"""
# make my default colormap
# cmap = plt.get_cmap('viridis')
# cmap = plt.get_cmap('jet')
# cmap = icacolors.comet
mycmap = kwargs.pop('cmap', None)
if mycmap is None:
mycmap = copy.copy(icacolors.viridis)
mycmap.set_bad(color=(0.2, 0.2, 0.2), alpha=1.0)
funcname = 'irfpy.ica.pyplot.plotsectorpolar'
if len(args) == 1:
X = np.arange(np.shape(args[0])[1])
Y = np.arange(np.shape(args[0])[0]).astype(float)
C = np.array(args[0]).astype(float).astype(float)
elif len(args) == 3:
X = np.array(args[0])
Y = np.array(args[1]).astype(float).copy()
C = np.array(args[2]).astype(float)
else:
raise TypeError(
'Incompatible X, Y or C inputs to %s; see help(%s)' % (
funcname, funcname))
# more error checking:
if len(np.shape(X)) != 1:
raise TypeError(
'X must be one-dimensional but is of shape ' +
str(np.shape(X))
)
if len(np.shape(Y)) != 1:
raise TypeError(
'Y must be one-dimensional but is of shape ' +
str(np.shape(Y))
)
if np.shape(C)[1] != np.shape(X)[0]:
raise TypeError(
'Dimension of C (y,x)=' + str(np.shape(C)) +
' is incompatible with the dimension of X = ' +
str(np.shape(X)[0]) + ', it should be ' +
str(np.shape(C)[1])
)
if np.shape(C)[0] != np.shape(Y)[0]:
raise TypeError(
'Dimension of C (y,x)=' + str(np.shape(C)) +
' is incompatible with the dimension of Y = ' +
str(np.shape(Y)[0]) + ', it should be ' +
str(np.shape(C)[0])
)
# check if the Y-vector contains nan. If so the nan are removed together
# with the corresponding elements in C
notnanmask = np.invert(np.isnan(Y))
Y = Y[notnanmask]
C = C[notnanmask, :]
def format_fnx(val, pos):
k = int(np.floor(val / (2 * np.pi) * X.size))
if k >= 0 and k < X.size:
return str(X[k])
else:
return ''
# subdivide each sector in this number of subsectors
# but do not make more than 256 sub sectors in total
multiplier = 256 // np.shape(X)[0]
if multiplier < 1:
multiplier = 1
thetazerodirection = kwargs.pop('zeroangle', 0)
centersectoroffset = 360.0/len(X)/2.0
shift = len(X) / 360.0 * centersectoroffset
theta = (np.arange(multiplier * len(X) + 1) - shift * multiplier) / (multiplier * len(X)) * 2 * np.pi
linear = kwargs.pop('linear', len(Y) <= 33)
radius = icatools.energyintervals(Y, linear=linear)
# the polar axis:
# plt.axes(polar=True,axisbg='#333333')
# this must be called before anything is plotted
D = np.repeat(C, multiplier, axis=1) # smooth
ax = plt.gca()
if '.PolarAxes' not in str(type(ax)):
raise ValueError("Expected the plot axis to be: \n" +
"<class 'matplotlib.axes._subplots.PolarAxesSubplot'> or <class 'matplotlib.axes._subplots.PolarAxes'>,\n" +
"but it is instead: \n" +
str(type(ax)) + ".\n\n" +
"plotsectorpolar() requires that the polar " +
"projection is activated:\n" +
"If you do not use plt.subplot() then add:\n" +
" plt.subplot(1,1,1,projection='polar')\n" +
"to your code. If you already have a plt.subplot() " +
"statement, then add the\n" +
"projection='polar' parameter to the subplot() " +
"statement that defines the \n" +
"subplot for plotsectorpolar().")
# matplotlib.axes._subplots.PolarAxesSubplot': ok
# matplotlib.axes._subplots.AxesSubplot': fail
#def scatter_logpolar(ax, theta, r_, bullseye=0.3, **kwargs):
# min10 = np.log10(np.min(r_))
# max10 = np.log10(np.max(r_))
# r = np.log10(r_) - min10 + bullseye
# ax.scatter(theta, r, **kwargs)
# l = np.arange(np.floor(min10), max10)
# ax.set_rticks(l - min10 + bullseye)
# ax.set_yticklabels(["1e%d" % x for x in l])
# ax.set_rlim(0, max10 - min10 + bullseye)
# ax.set_title('log-polar manual')
# return ax
plt.pcolormesh(theta,
radius,
D,
cmap=mycmap,
**kwargs)
ax.xaxis.set_major_formatter(FuncFormatter(format_fnx))
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
startangle = 360 / len(X) * 0.5 - centersectoroffset
thetaticks = np.arange(startangle, 360 + startangle, 360 / len(X))
# set ticklabels location at 1.3 times the axes' radius
maxradius = np.nanmax(radius)
ax.set_thetagrids(thetaticks)
ax.set_theta_zero_location('E',offset=thetazerodirection)
if linear:
minradius = 0.001
ax.set_rscale('linear')
ax.set_rmin(minradius)
ax.set_rmax(maxradius)
else:
minradius = 1.0
ax.set_rscale('log')
ax.set_rmin(minradius)
ax.set_rmax(maxradius*0.1)
ax.xaxis.set_tick_params(pad=25.0)
ax.set_rlabel_position(28)
ax.grid(True, color='w', linestyle=':', linewidth=0.5)
return ax
[docs]def plotenergypolar(*args, **kwargs):
r"""
Plots a color matrix as polar plot with matplotlib.pyplot.pcolormesh()
using a quadrilateral mesh. It takes the similar parameters
as pcolormesh(). In addition to pcolormesh(), this function makes sure
the lables are centered on each bin on both axes and that all bins
are plotted, including the last one.
This function superseeds plotsectorpolar() and fixes a number of issues
when the radial axis is logarithmic. This funciton always plots internally the
data with a linear scale but takes the log10 of the radial coordinate in case
linear=False.
Call signatures::
plotenergypolar(C)
plotenergypolar(X, Y, C)
plotenergypolar(C, **kwargs)
plotenergypolar(X, Y, C, **kwargs)
If X and Y are given, then the values in them are used as axis labels.
The shape of C must be (Y,X), otherwise an error is generated. C can
be a numpy.masked array.
Allowed for kwargs is everything that pcolormesh accepts plus the ones
listed below:
linear:
True : Energy intervals are interpolated in a linear fashion
The radial axis is linear
False : Energy intervals are interpolated in a geometric fashion
The radial axis is the log10(Y).
zeroangle :
The azimuth angle that should be plotted in the +X direction of the
plot in degrees.
rmin : the center of the circle corresponds to this value.
Default: 0.0 if linear == True, 1.0 otherwise
cmap:
A color map.
Plot a polar sector plot::
import numpy as np
import irfpy.ica.io as icaio
import matplotlib as mpl
import irfpy.ica.pyplot as plt
import irfpy.ica.tools as icatools
from irfpy.util.irfpyrc import Rc
rc = Rc()
path = rc.get('icadds', 'dataroot')
starttime="20150605T121500"
stoptime="20150605T134000"
cmin = 0.1
cmax = 1e4
def sumorigionspacc(orig_ionspectra,mode,version):
epsilon=1e-10
D3 = np.zeros((16, 96))
for start,stop,(vm,vv) in icatools.timeslices(zip(mode,version)):
D3 = D3 + np.nansum(np.nansum(orig_ionspectra[:,:,:,start:stop],axis=2),axis=2)
# sect x energy
D3[D3 < epsilon] = epsilon
return D3.T
mat = icaio.readproc(path,starttime,stoptime,
variables=['orig_ionspectra',
'cur_pacc','mode','sw_version','E'])
orig_ionspectra=mat['orig_ionspectra']
mode=mat['mode']
sw_version=mat['sw_version']
E = icatools.getE(mat['E'],sw_version)
D = sumorigionspacc(orig_ionspectra,mode,sw_version)
A=np.array(['S'+str(n) for n in np.arange(16)])
plt.clf()
plt.subplot(1,2,1,projection='polar')
pax = plt.plotsectorpolar(A,E[:,0],D,
norm=mpl.colors.LogNorm(vmin=cmin, vmax=cmax))
plt.colorbar(label='Counts')
plt.title('Sector/Energy linear')
plt.yscale("linear")
# Example: remove the radial labels:
# This hides the tick labels without removing the label text
pax.tick_params(axis='y', labelsize=0)
plt.subplot(1,2,2,projection='polar',axisbg='#333333')
pax=plt.plotsectorpolar(A,E[:,0],D,
norm=mpl.colors.LogNorm(vmin=cmin, vmax=cmax))
plt.title('Sector/Energy log')
cbar = plt.colorbar(label='Counts')
rax = plt.radialscale(label='Energy [eV]',linewidth=2)
# Example: change font sizes everywhere to 12:
for item in ([pax.title, # the title
pax.xaxis.label, # azimuth axis label
pax.yaxis.label, # radial axis label
cbar.ax.yaxis.label] + # color bar label
pax.get_xticklabels() + # azimuth tick labels
rax.get_ticklabels() + # radial scale tick labels
cbar.ax.yaxis.get_ticklabels() # color bar tick labels
):
item.set_fontsize(12)
plt.tight_layout()
plt.show()
Or also shorter::
plotenergypolar(C)
"""
# make my default colormap
# cmap = plt.get_cmap('viridis')
# cmap = plt.get_cmap('jet')
# cmap = icacolors.comet
mycmap = kwargs.pop('cmap', None)
if mycmap is None:
mycmap = copy.copy(icacolors.viridis)
mycmap.set_bad(color=(0.2, 0.2, 0.2), alpha=1.0)
funcname = 'irfpy.ica.pyplot.plotsectorpolar'
if len(args) == 1:
X = np.arange(np.shape(args[0])[1])
Y = np.arange(np.shape(args[0])[0]).astype(float)
C = np.array(args[0]).astype(float).astype(float)
elif len(args) == 3:
X = np.array(args[0]).copy()
Y = np.array(args[1]).astype(float).copy()
C = np.array(args[2]).astype(float)
else:
raise TypeError(
'Incompatible X, Y or C inputs to %s; see help(%s)' % (
funcname, funcname))
# more error checking:
if len(np.shape(X)) != 1:
raise TypeError(
'X must be one-dimensional but is of shape ' +
str(np.shape(X))
)
if len(np.shape(Y)) != 1:
raise TypeError(
'Y must be one-dimensional but is of shape ' +
str(np.shape(Y))
)
if np.shape(C)[1] != np.shape(X)[0]:
raise TypeError(
'Dimension of C (y,x)=' + str(np.shape(C)) +
' is incompatible with the dimension of X = ' +
str(np.shape(X)[0]) + ', it should be ' +
str(np.shape(C)[1])
)
if np.shape(C)[0] != np.shape(Y)[0]:
raise TypeError(
'Dimension of C (y,x)=' + str(np.shape(C)) +
' is incompatible with the dimension of Y = ' +
str(np.shape(Y)[0]) + ', it should be ' +
str(np.shape(C)[0])
)
# check if the Y-vector contains nan. If so the nan are removed together
# with the corresponding elements in C
notnanmask = np.invert(np.isnan(Y))
Y = Y[notnanmask]
C = C[notnanmask, :]
def format_fnx(val, pos):
k = int(np.floor(val / (2 * np.pi) * X.size))
if k >= 0 and k < X.size:
return str(X[k])
else:
return ''
multiplier = 256 // np.shape(X)[0]
if multiplier < 1:
multiplier = 1
thetazerodirection = kwargs.pop('zeroangle', 0)
centersectoroffset = 360.0/len(X)/2.0
shift = len(X) / 360.0 * centersectoroffset
theta = (np.arange(multiplier * len(X) + 1) - shift * multiplier) / (multiplier * len(X)) * 2 * np.pi
linear = kwargs.pop('linear', len(Y) <= 33)
minradius = kwargs.pop('minradius', 0.0)
radius = icatools.energyintervals(Y, linear=linear)
if not linear:
radius = np.log10(radius)
Y = np.log10(Y)
# the polar axis:
# plt.axes(polar=True,axisbg='#333333')
# this must be called before anything is plotted
D = np.repeat(C, multiplier, axis=1) # smooth
ax = plt.gca()
if '.PolarAxesSubplot' not in str(type(ax)):
raise ValueError("Expected the plot axis to be: \n" +
"<class 'matplotlib.axes._subplots.PolarAxesSubplot'>,\n" +
"but it is instead: \n" +
str(type(ax)) + ".\n\n" +
"plotenergypolar() requires that the polar " +
"projection is activated:\n" +
"If you do not use plt.subplot() then add:\n" +
" plt.subplot(1,1,1,projection='polar')\n" +
"to your code. If you already have a plt.subplot() " +
"statement, then add the\n" +
"projection='polar' parameter to the subplot() " +
"statement that defines the \n" +
"subplot for plotsectorpolar().")
# matplotlib.axes._subplots.PolarAxesSubplot': ok
# matplotlib.axes._subplots.AxesSubplot': fail
plt.pcolormesh(theta,
radius,
D,
cmap=mycmap,
**kwargs)
ax.xaxis.set_major_formatter(FuncFormatter(format_fnx))
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
startangle = 360 / len(X) * 0.5 - centersectoroffset
thetaticks = np.arange(startangle, 360 + startangle, 360 / len(X))
# set ticklabels location at 1.3 times the axes' radius
maxradius = np.nanmax(radius)
ax.set_thetagrids(thetaticks)
ax.set_theta_zero_location('E',offset=thetazerodirection)
if linear:
minradius = 0.001
ax.set_rscale('linear')
ax.set_rmin(minradius)
ax.set_rmax(maxradius)
else:
minradius = 1.0
ax.set_rscale('linear')
ax.set_rmin(minradius)
ax.set_rmax(maxradius)
ax.set_rlabel_position(28)
ax.grid(True, color='w', linestyle=':', linewidth=0.5)
return ax
#%%
if __name__ == '__main__':
import irfpy.ica.io as icaio
from irfpy.util.irfpyrc import Rc
import matplotlib as mpl
import irfpy.ica.tools as icatools
path = '/Users/wieser/ica/processed'
starttime="20150605T121500"
stoptime="20150605T134000"
cmin = 0.1
cmax = 1e4
def sumorigionspacc(orig_ionspectra,mode,version):
epsilon=1e-10
D3 = np.zeros((16, 96))
for start,stop,(vm,vv) in icatools.timeslices(zip(mode,version)):
D3 = D3 + np.nansum(np.nansum(orig_ionspectra[:,:,:,start:stop],axis=2),axis=2)
# sect x energy
D3[D3 < epsilon] = epsilon
return D3.T
mat = icaio.readproc(path,starttime,stoptime,
variables=['orig_ionspectra',
'cur_pacc','mode','sw_version','E'])
orig_ionspectra=mat['orig_ionspectra']
mode=mat['mode']
sw_version=mat['sw_version']
E = icatools.getE(mat['E'],sw_version)
D = sumorigionspacc(orig_ionspectra,mode,sw_version)
A=np.array(['S'+str(n) for n in np.arange(16)])
# A=np.array(['' for n in np.arange(16)])
plt.figure(99)
plt.clf()
plt.subplot(1,2,1,projection='polar')
pax = plotenergypolar(A,E[:,0],D,
norm=mpl.colors.LogNorm(vmin=cmin, vmax=cmax),
linear=True)
plt.colorbar(label='Counts')
plt.title('Sector/Energy linear')
# Example: remove the radial labels:
# This hides the tick labels without removing the label text
pax.tick_params(axis='y', labelsize=0)
plt.subplot(1,2,2,projection='polar')
pax=plotenergypolar(A,E[:,0],D,
norm=mpl.colors.LogNorm(vmin=cmin, vmax=cmax),
linear = False)
plt.title('Sector/Energy log')
cbar = plt.colorbar(label='Counts')
# rax = radialscale(label='Energy [eV]',linewidth=2)
# Example: change font sizes everywhere to 12:
pax.tick_params(axis='y', labelsize=0)
plt.tight_layout()
plt.show()
plt.figure(100)
starttime="20150605T121500"
stoptime="20150605T134000"
mat = icaio.readproc(path,starttime,stoptime,
variables=['orig_ionspectra',
'cur_pacc',
'mode',
'sw_version',
'E',
'time_instances',
'sum_orig_ionspectra']
)
orig_ionspectra=mat['orig_ionspectra']
sw_version=mat['sw_version']
time_instances=mat['time_instances']
D = mat['sum_orig_ionspectra']
plt.clf()
for start,stop,version in icatools.timeslices(sw_version):
# This loop allows different sw_versions in one plot
# get the expected length in time of one measuremnt in seconds
# +0.1s is added to avoid plotting artefacts
deltaT = icatools.deltaTfromversion(version)+0.1
# get an energy vector for each time instance
E = icatools.getE(mat['E'],sw_version[start:stop])
Eofs = np.arange(start,stop)
E= E + Eofs
# epsilon is added to D to allow to plot the value 0.0 in log-scale:
epsilon = 1e-10
plotenergytime(time_instances[start:stop],
E,
D[:,start:stop]+epsilon,
norm=mpl.colors.LogNorm(vmin=0.2, vmax=10000.0),
dateformat='%H:%M', # optional
maxtimebinwidth=deltaT, # optional
cmap=plt.get_cmap('viridis'),
rectangular=True) # optional
plt.xlim(time_instances[0],time_instances[-1])
plt.colorbar(label='Counts')
plt.yscale('log')
plt.ylim(3,4e4)
plt.ylabel('E (eV)')
plt.title('ICA Energy-Time Spectrogram')
plt.show()