from builtins import zip
import numpy as np
import warnings
import healpy as hp
from matplotlib import colors
from matplotlib import ticker
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.ticker import FuncFormatter
import matplotlib as mpl
from matplotlib.patches import Ellipse
from matplotlib.collections import PatchCollection
from lsst.sims.maf.utils import optimalBins, percentileClipping
from .plotHandler import BasePlotter
from lsst.sims.utils import _equatorialFromGalactic
import inspect
from .perceptual_rainbow import makePRCmap
perceptual_rainbow = makePRCmap()
__all__ = ['HealpixSkyMap', 'HealpixPowerSpectrum', 'HealpixHistogram', 'OpsimHistogram',
'BaseHistogram', 'BaseSkyMap', 'HealpixSDSSSkyMap', 'LambertSkyMap']
[docs]class HealpixSkyMap(BasePlotter):
"""
Generate a sky map of healpix metric values using healpy's mollweide view.
"""
def __init__(self):
super(HealpixSkyMap, self).__init__()
# Set the plotType
self.plotType = 'SkyMap'
self.objectPlotter = False
# Set up the default plotting parameters.
self.defaultPlotDict = {'title': None, 'xlabel': None, 'label': None,
'logScale': False, 'cbarFormat': None, 'cmap': perceptual_rainbow,
'percentileClip': None, 'colorMin': None, 'colorMax': None,
'zp': None, 'normVal': None, 'labelsize': None, 'fontsize': None,
'cbar_edge': True, 'nTicks': 15, 'rot': (0, 0, 0), 'figsize': None}
[docs] def __call__(self, metricValueIn, slicer, userPlotDict, fignum=None):
"""
Parameters
----------
metricValue : numpy.ma.MaskedArray
slicer : lsst.sims.maf.slicers.HealpixSlicer
userPlotDict: dict
Dictionary of plot parameters set by user (overrides default values).
fignum : int
Matplotlib figure number to use (default = None, starts new figure).
Returns
-------
int
Matplotlib figure number used to create the plot.
"""
# Check that the slicer is a HealpixSlicer, or subclass thereof
# Using the names rather than just comparing the classes themselves
# to avoid circular dependency between slicers and plots
classes = inspect.getmro(slicer.__class__)
cnames = [cls.__name__ for cls in classes]
if 'HealpixSlicer' not in cnames:
raise ValueError('HealpixSkyMap is for use with healpix slicers')
# Override the default plotting parameters with user specified values.
plotDict = {}
plotDict.update(self.defaultPlotDict)
plotDict.update(userPlotDict)
# Generate a Mollweide full-sky plot.
fig = plt.figure(fignum, figsize=plotDict['figsize'])
norm = None
if plotDict['logScale']:
norm = 'log'
cmap = plotDict['cmap']
if type(cmap) == str:
cmap = getattr(cm, cmap)
# Set background and masked pixel colors default healpy white and gray.
cmap.set_over(cmap(1.0))
cmap.set_under('w')
cmap.set_bad('gray')
if plotDict['zp'] is not None:
metricValue = metricValueIn - plotDict['zp']
elif plotDict['normVal'] is not None:
metricValue = metricValueIn / plotDict['normVal']
else:
metricValue = metricValueIn
# Set up color bar limits.
colorMin = plotDict['colorMin']
colorMax = plotDict['colorMax']
if plotDict['percentileClip'] is not None:
pcMin, pcMax = percentileClipping(metricValue.compressed(), percentile=plotDict['percentileClip'])
if colorMin is None and plotDict['percentileClip']:
colorMin = pcMin
if colorMax is None and plotDict['percentileClip']:
colorMax = pcMax
if (colorMin is not None) or (colorMax is not None):
clims = [colorMin, colorMax]
else:
clims = None
# Make sure there is some range on the colorbar
if clims is None:
if metricValue.compressed().size > 0:
clims = [metricValue.compressed().min(), metricValue.compressed().max()]
else:
clims = [-1, 1]
if clims[0] == clims[1]:
clims[0] = clims[0] - 1
clims[1] = clims[1] + 1
# Avoid trying to log scale when zero is in the range.
if (norm == 'log') & ((clims[0] <= 0 <= clims[1]) or (clims[0] >= 0 >= clims[1])):
# Try something simple
above = metricValue[np.where(metricValue > 0)]
if len(above) > 0:
clims[0] = above.max()
# If still bad, give up and turn off norm
if ((clims[0] <= 0 <= clims[1]) or (clims[0] >= 0 >= clims[1])):
norm = None
warnings.warn("Using norm was set to log, but color limits pass through 0. "
"Adjusting so plotting doesn't fail")
hp.mollview(metricValue.filled(slicer.badval), title=plotDict['title'], cbar=False,
min=clims[0], max=clims[1], rot=plotDict['rot'], flip='astro',
cmap=cmap, norm=norm, fig=fig.number)
# This graticule call can fail with old versions of healpy and matplotlib 1.4.0.
# Make sure the latest version of healpy in the stack is setup
hp.graticule(dpar=20, dmer=20, verbose=False)
# Add colorbar (not using healpy default colorbar because we want more tickmarks).
ax = plt.gca()
im = ax.get_images()[0]
# Add label.
if plotDict['label'] is not None:
plt.figtext(0.8, 0.8, '%s' % (plotDict['label']))
# supress silly colorbar warnings
with warnings.catch_warnings():
warnings.simplefilter("ignore")
cb = plt.colorbar(im, shrink=0.75, aspect=25, orientation='horizontal',
format=plotDict['cbarFormat'], extendrect=True)
cb.set_label(plotDict['xlabel'], fontsize=plotDict['fontsize'])
if plotDict['labelsize'] is not None:
cb.ax.tick_params(labelsize=plotDict['labelsize'])
if norm == 'log':
tick_locator = ticker.LogLocator(numticks=plotDict['nTicks'])
cb.locator = tick_locator
cb.update_ticks()
if (plotDict['nTicks'] is not None) & (norm != 'log'):
tick_locator = ticker.MaxNLocator(nbins=plotDict['nTicks'])
cb.locator = tick_locator
cb.update_ticks()
# If outputing to PDF, this fixes the colorbar white stripes
if plotDict['cbar_edge']:
cb.solids.set_edgecolor("face")
return fig.number
[docs]class HealpixPowerSpectrum(BasePlotter):
def __init__(self):
self.plotType = 'PowerSpectrum'
self.objectPlotter = False
self.defaultPlotDict = {'title': None, 'label': None,
'maxl': None, 'removeDipole': True,
'logScale': True, 'linestyle': '-',
'fontsize': None, 'labelsize': None, 'figsize': None}
[docs] def __call__(self, metricValue, slicer, userPlotDict, fignum=None):
"""
Generate and plot the power spectrum of metricValue (calculated on a healpix grid).
"""
if slicer.slicerName != 'HealpixSlicer':
raise ValueError('HealpixPowerSpectrum for use with healpix metricBundles.')
plotDict = {}
plotDict.update(self.defaultPlotDict)
plotDict.update(userPlotDict)
fig = plt.figure(fignum, figsize=plotDict['figsize'])
# If the mask is True everywhere (no data), just plot zeros
if False not in metricValue.mask:
return None
if plotDict['removeDipole']:
cl = hp.anafast(hp.remove_dipole(metricValue.filled(slicer.badval)), lmax=plotDict['maxl'])
else:
cl = hp.anafast(metricValue.filled(slicer.badval), lmax=plotDict['maxl'])
ell = np.arange(np.size(cl))
if plotDict['removeDipole']:
condition = (ell > 1)
else:
condition = (ell > 0)
ell = ell[condition]
cl = cl[condition]
# Plot the results.
plt.plot(ell, (cl * ell * (ell + 1)) / 2.0 / np.pi,
color=plotDict['color'], linestyle=plotDict['linestyle'], label=plotDict['label'])
if cl.max() > 0 and plotDict['logScale']:
plt.yscale('log')
plt.xlabel(r'$l$', fontsize=plotDict['fontsize'])
plt.ylabel(r'$l(l+1)C_l/(2\pi)$', fontsize=plotDict['fontsize'])
if plotDict['labelsize'] is not None:
plt.tick_params(axis='x', labelsize=plotDict['labelsize'])
plt.tick_params(axis='y', labelsize=plotDict['labelsize'])
if plotDict['title'] is not None:
plt.title(plotDict['title'])
# Return figure number (so we can reuse/add onto/save this figure if desired).
return fig.number
[docs]class HealpixHistogram(BasePlotter):
def __init__(self):
self.plotType = 'Histogram'
self.objectPlotter = False
self.defaultPlotDict = {'title': None, 'xlabel': None,
'ylabel': 'Area (1000s of square degrees)', 'label': None,
'bins': None, 'binsize': None, 'cumulative': False,
'scale': None, 'xMin': None, 'xMax': None,
'logScale': False, 'linestyle': '-',
'fontsize': None, 'labelsize': None, 'figsize': None}
self.baseHist = BaseHistogram()
[docs] def __call__(self, metricValue, slicer, userPlotDict, fignum=None):
"""
Histogram metricValue for all healpix points.
"""
if slicer.slicerName != 'HealpixSlicer':
raise ValueError('HealpixHistogram is for use with healpix slicer.')
plotDict = {}
plotDict.update(self.defaultPlotDict)
plotDict.update(userPlotDict)
if plotDict['scale'] is None:
plotDict['scale'] = (hp.nside2pixarea(slicer.nside, degrees=True) / 1000.0)
fignum = self.baseHist(metricValue, slicer, plotDict, fignum=fignum)
return fignum
[docs]class OpsimHistogram(BasePlotter):
def __init__(self):
self.plotType = 'Histogram'
self.objectPlotter = False
self.defaultPlotDict = {'title': None, 'xlabel': None, 'label': None,
'ylabel': 'Number of Fields', 'yaxisformat': '%d',
'bins': None, 'binsize': None, 'cumulative': False,
'scale': 1.0, 'xMin': None, 'xMax': None,
'logScale': False, 'linestyle': '-',
'fontsize': None, 'labelsize': None, 'figsize': None}
self.baseHist = BaseHistogram()
[docs] def __call__(self, metricValue, slicer, userPlotDict, fignum=None):
"""
Histogram metricValue for all healpix points.
"""
if slicer.slicerName != 'OpsimFieldSlicer':
raise ValueError('OpsimHistogram is for use with OpsimFieldSlicer.')
plotDict = {}
plotDict.update(self.defaultPlotDict)
plotDict.update(userPlotDict)
fignum = self.baseHist(metricValue, slicer, plotDict, fignum=fignum)
return fignum
[docs]class BaseHistogram(BasePlotter):
def __init__(self):
self.plotType = 'Histogram'
self.objectPlotter = False
self.defaultPlotDict = {'title': None, 'xlabel': None, 'ylabel': 'Count', 'label': None,
'bins': None, 'binsize': None, 'cumulative': False,
'scale': 1.0,
'xMin': None, 'xMax': None,
'yMin': None, 'yMax': None,
'logScale': False,
'yaxisformat': '%.3f', 'linestyle': '-',
'zp': None, 'normVal': None, 'percentileClip': None,
'fontsize': None, 'labelsize': None, 'figsize': None}
[docs] def __call__(self, metricValueIn, slicer, userPlotDict, fignum=None):
"""
Plot a histogram of metricValues (such as would come from a spatial slicer).
"""
# Adjust metric values by zeropoint or normVal, and use 'compressed' version of masked array.
plotDict = {}
plotDict.update(self.defaultPlotDict)
plotDict.update(userPlotDict)
if plotDict['zp'] is not None:
metricValue = metricValueIn.compressed() - plotDict['zp']
elif plotDict['normVal'] is not None:
metricValue = metricValueIn.compressed() / plotDict['normVal']
else:
metricValue = metricValueIn.compressed()
# Toss any NaNs or infs
metricValue = metricValue[np.isfinite(metricValue)]
# Determine percentile clipped X range, if set. (and xmin/max not set).
if plotDict['xMin'] is None and plotDict['xMax'] is None:
if plotDict['percentileClip']:
plotDict['xMin'], plotDict['xMax'] = percentileClipping(metricValue,
percentile=plotDict['percentileClip'])
# Set the histogram range values, to avoid None/Number comparisons.
histRange = [plotDict['xMin'], plotDict['xMax']]
if histRange[0] is None:
histRange[0] = metricValue.min()
if histRange[1] is None:
histRange[1] = metricValue.max()
# Need to have some range of values on the histogram, or it will fail.
if histRange[0] == histRange[1]:
warnings.warn('Histogram range was single-valued; expanding default range.')
histRange[1] = histRange[0] + 1.0
# Set up the bins for the histogram. User specified 'bins' overrides 'binsize'.
# Note that 'bins' could be a single number or an array, simply passed to plt.histogram.
if plotDict['bins'] is not None:
bins = plotDict['bins']
elif plotDict['binsize'] is not None:
# If generating a cumulative histogram, want to use full range of data (but with given binsize).
# .. but if user set histRange to be wider than full range of data, then
# extend bins to cover this range, so we can make prettier plots.
if plotDict['cumulative']:
if plotDict['xMin'] is not None:
# Potentially, expand the range for the cumulative histogram.
bmin = np.min([metricValue.min(), plotDict['xMin']])
else:
bmin = metricValue.min()
if plotDict['xMax'] is not None:
bmax = np.max([metricValue.max(), plotDict['xMax']])
else:
bmax = metricValue.max()
bins = np.arange(bmin, bmax + plotDict['binsize'] / 2.0, plotDict['binsize'])
# Otherwise, not cumulative so just use metric values, without potential expansion.
else:
bins = np.arange(histRange[0], histRange[1] + plotDict['binsize'] / 2.0, plotDict['binsize'])
# Catch edge-case where there is only 1 bin value
if bins.size < 2:
bins = np.arange(bins.min() - plotDict['binsize'] * 2.0,
bins.max() + plotDict['binsize'] * 2.0, plotDict['binsize'])
else:
bins = optimalBins(metricValue)
# Generate plots.
fig = plt.figure(fignum, figsize=plotDict['figsize'])
# Check if any data falls within histRange, because otherwise histogram generation will fail.
if isinstance(bins, np.ndarray):
condition = ((metricValue >= bins.min()) & (metricValue <= bins.max()))
else:
condition = ((metricValue >= histRange[0]) & (metricValue <= histRange[1]))
plotValue = metricValue[condition]
if len(plotValue) == 0:
# No data is within histRange/bins. So let's just make a simple histogram anyway.
n, b, p = plt.hist(metricValue, bins=50, histtype='step', cumulative=plotDict['cumulative'],
log=plotDict['logScale'], label=plotDict['label'],
color=plotDict['color'])
else:
# There is data to plot, and we've already ensured histRange/bins are more than single value.
n, b, p = plt.hist(metricValue, bins=bins, range=histRange,
histtype='step', log=plotDict['logScale'],
cumulative=plotDict['cumulative'],
label=plotDict['label'], color=plotDict['color'])
# Fill in axes labels and limits.
# Option to use 'scale' to turn y axis into area or other value.
def mjrFormatter(y, pos):
return plotDict['yaxisformat'] % (y * plotDict['scale'])
ax = plt.gca()
ax.yaxis.set_major_formatter(FuncFormatter(mjrFormatter))
# Set optional x, y limits.
plt.ylim([plotDict['yMin'], plotDict['yMax']])
plt.xlim([plotDict['xMin'], plotDict['xMax']])
# Set/Add various labels.
plt.xlabel(plotDict['xlabel'], fontsize=plotDict['fontsize'])
plt.ylabel(plotDict['ylabel'], fontsize=plotDict['fontsize'])
plt.title(plotDict['title'])
if plotDict['labelsize'] is not None:
plt.tick_params(axis='x', labelsize=plotDict['labelsize'])
plt.tick_params(axis='y', labelsize=plotDict['labelsize'])
# Return figure number
return fig.number
[docs]class BaseSkyMap(BasePlotter):
def __init__(self):
self.plotType = 'SkyMap'
self.objectPlotter = False # unless 'metricIsColor' is true..
self.defaultPlotDict = {'title': None, 'xlabel': None, 'label': None,
'projection': 'aitoff', 'radius': np.radians(1.75),
'logScale': 'auto', 'cbar': True, 'cbarFormat': None,
'cmap': perceptual_rainbow, 'alpha': 1.0,
'zp': None, 'normVal': None,
'colorMin': None, 'colorMax': None, 'percentileClip': None,
'cbar_edge': True, 'plotMask': False, 'metricIsColor': False,
'raCen': 0.0, 'mwZone': True, 'bgcolor': 'gray',
'labelsize': None, 'fontsize': None, 'figsize': None}
def _plot_tissot_ellipse(self, lon, lat, radius, ax=None, **kwargs):
"""Plot Tissot Ellipse/Tissot Indicatrix
Parameters
----------
lon : float or array_like
longitude-like of ellipse centers (radians)
lat : float or array_like
latitude-like of ellipse centers (radians)
radius : float or array_like
radius of ellipses (radians)
ax : Axes object (optional)
matplotlib axes instance on which to draw ellipses.
Other Parameters
----------------
other keyword arguments will be passed to matplotlib.patches.Ellipse.
# The code in this method adapted from astroML, which is BSD-licensed.
# See http: //github.com/astroML/astroML for details.
"""
# Code adapted from astroML, which is BSD-licensed.
# See http: //github.com/astroML/astroML for details.
ellipses = []
if ax is None:
ax = plt.gca()
for l, b, diam in np.broadcast(lon, lat, radius * 2.0):
el = Ellipse((l, b), diam / np.cos(b), diam, **kwargs)
ellipses.append(el)
return ellipses
def _plot_ecliptic(self, raCen=0, ax=None):
"""
Plot a red line at location of ecliptic.
"""
if ax is None:
ax = plt.gca()
ecinc = 23.439291 * (np.pi / 180.0)
ra_ec = np.arange(0, np.pi * 2., (np.pi * 2. / 360.))
dec_ec = np.sin(ra_ec) * ecinc
lon = -(ra_ec - raCen - np.pi) % (np.pi * 2) - np.pi
ax.plot(lon, dec_ec, 'r.', markersize=1.8, alpha=0.4)
def _plot_mwZone(self, raCen=0, peakWidth=np.radians(10.), taperLength=np.radians(80.), ax=None):
"""
Plot blue lines to mark the milky way galactic exclusion zone.
"""
if ax is None:
ax = plt.gca()
# Calculate galactic coordinates for mw location.
step = 0.02
galL = np.arange(-np.pi, np.pi + step / 2., step)
val = peakWidth * np.cos(galL / taperLength * np.pi / 2.)
galB1 = np.where(np.abs(galL) <= taperLength, val, 0)
galB2 = np.where(np.abs(galL) <= taperLength, -val, 0)
# Convert to ra/dec.
# Convert to lon/lat and plot.
ra, dec = _equatorialFromGalactic(galL, galB1)
lon = -(ra - raCen - np.pi) % (np.pi * 2) - np.pi
ax.plot(lon, dec, 'b.', markersize=1.8, alpha=0.4)
ra, dec = _equatorialFromGalactic(galL, galB2)
lon = -(ra - raCen - np.pi) % (np.pi * 2) - np.pi
ax.plot(lon, dec, 'b.', markersize=1.8, alpha=0.4)
[docs] def __call__(self, metricValueIn, slicer, userPlotDict, fignum=None):
"""
Plot the sky map of metricValue for a generic spatial slicer.
"""
if 'ra' not in slicer.slicePoints or 'dec' not in slicer.slicePoints:
errMessage = 'SpatialSlicer must contain "ra" and "dec" in slicePoints metadata.'
errMessage += ' SlicePoints only contains keys %s.' % (slicer.slicePoints.keys())
raise ValueError(errMessage)
plotDict = {}
plotDict.update(self.defaultPlotDict)
plotDict.update(userPlotDict)
if plotDict['zp'] is not None:
metricValue = metricValueIn - plotDict['zp']
elif plotDict['normVal'] is not None:
metricValue = metricValueIn / plotDict['normVal']
else:
metricValue = metricValueIn
fig = plt.figure(fignum, figsize=plotDict['figsize'])
# other projections available include
# ['aitoff', 'hammer', 'lambert', 'mollweide', 'polar', 'rectilinear']
ax = fig.add_subplot(111, projection=plotDict['projection'])
# Set up valid datapoints and colormin/max values.
if plotDict['plotMask']:
# Plot all data points.
mask = np.ones(len(metricValue), dtype='bool')
else:
# Only plot points which are not masked. Flip numpy ma mask where 'False' == 'good'.
mask = ~metricValue.mask
# Determine color min/max values. metricValue.compressed = non-masked points.
if not plotDict['metricIsColor']:
if plotDict['percentileClip'] is not None:
pcMin, pcMax = percentileClipping(metricValue.compressed(),
percentile=plotDict['percentileClip'])
if plotDict['colorMin'] is None:
if plotDict['percentileClip']:
plotDict['colorMin'] = pcMin
else:
plotDict['colorMin'] = metricValue.compressed().min()
if plotDict['colorMax'] is None:
if plotDict['percentileClip']:
plotDict['colorMax'] = pcMax
else:
plotDict['colorMax'] = metricValue.compressed().max()
# Avoid colorbars with no range.
if plotDict['colorMax'] == plotDict['colorMin']:
plotDict['colorMax'] = plotDict['colorMax'] + 1
plotDict['colorMin'] = plotDict['colorMin'] - 1
# Combine to make clims:
clims = [plotDict['colorMin'], plotDict['colorMax']]
# Determine whether or not to use auto-log scale.
if plotDict['logScale'] == 'auto':
if plotDict['colorMin'] > 0:
if np.log10(plotDict['colorMax']) - np.log10(plotDict['colorMin']) > 3:
plotDict['logScale'] = True
else:
plotDict['logScale'] = False
else:
plotDict['logScale'] = False
if plotDict['logScale']:
# Move min/max values to things that can be marked on the colorbar.
plotDict['colorMin'] = 10**(int(np.log10(plotDict['colorMin'])))
plotDict['colorMax'] = 10**(int(np.log10(plotDict['colorMax'])))
# Add ellipses at RA/Dec locations
lon = -(slicer.slicePoints['ra'][mask] - plotDict['raCen'] - np.pi) % (np.pi * 2) - np.pi
ellipses = self._plot_tissot_ellipse(lon, slicer.slicePoints['dec'][mask],
plotDict['radius'], rasterized=True, ax=ax)
if plotDict['metricIsColor']:
current = None
for ellipse, mVal in zip(ellipses, metricValue.data[mask]):
if mVal[3] > 1:
ellipse.set_alpha(1.0)
ellipse.set_facecolor((mVal[0], mVal[1], mVal[2]))
ellipse.set_edgecolor('k')
current = ellipse
else:
ellipse.set_alpha(mVal[3])
ellipse.set_color((mVal[0], mVal[1], mVal[2]))
ax.add_patch(ellipse)
if current:
ax.add_patch(current)
else:
if plotDict['logScale']:
norml = colors.LogNorm()
p = PatchCollection(ellipses, cmap=plotDict['cmap'], alpha=plotDict['alpha'],
linewidth=0, edgecolor=None, norm=norml, rasterized=True)
else:
p = PatchCollection(ellipses, cmap=plotDict['cmap'], alpha=plotDict['alpha'],
linewidth=0, edgecolor=None, rasterized=True)
p.set_array(metricValue.data[mask])
p.set_clim(clims)
ax.add_collection(p)
# Add color bar (with optional setting of limits)
if plotDict['cbar']:
cb = plt.colorbar(p, aspect=25, extendrect=True, orientation='horizontal',
format=plotDict['cbarFormat'])
# If outputing to PDF, this fixes the colorbar white stripes
if plotDict['cbar_edge']:
cb.solids.set_edgecolor("face")
cb.set_label(plotDict['xlabel'], fontsize=plotDict['fontsize'])
cb.ax.tick_params(labelsize=plotDict['labelsize'])
# Add ecliptic
self._plot_ecliptic(plotDict['raCen'], ax=ax)
if plotDict['mwZone']:
self._plot_mwZone(plotDict['raCen'], ax=ax)
ax.grid(True, zorder=1)
ax.xaxis.set_ticklabels([])
if plotDict['bgcolor'] is not None:
ax.set_axis_bgcolor(plotDict['bgcolor'])
# Add label.
if plotDict['label'] is not None:
plt.figtext(0.75, 0.9, '%s' % plotDict['label'], fontsize=plotDict['fontsize'])
if plotDict['title'] is not None:
plt.text(0.5, 1.09, plotDict['title'], horizontalalignment='center',
transform=ax.transAxes, fontsize=plotDict['fontsize'])
return fig.number
[docs]class HealpixSDSSSkyMap(BasePlotter):
def __init__(self):
self.plotType = 'SkyMap'
self.objectPlotter = False
self.defaultPlotDict = {'title': None, 'xlabel': None, 'logScale': False,
'cbarFormat': '%.2f', 'cmap': perceptual_rainbow,
'percentileClip': None, 'colorMin': None,
'colorMax': None, 'zp': None, 'normVal': None,
'cbar_edge': True, 'label': None, 'raMin': -90,
'raMax': 90, 'raLen': 45, 'decMin': -2., 'decMax': 2.,
'labelsize': None, 'fontsize': None}
[docs] def __call__(self, metricValueIn, slicer, userPlotDict, fignum=None, ):
"""
Plot the sky map of metricValue using healpy cartview plots in thin strips.
raMin: Minimum RA to plot (deg)
raMax: Max RA to plot (deg). Note raMin/raMax define the centers that will be plotted.
raLen: Length of the plotted strips in degrees
decMin: minimum dec value to plot
decMax: max dec value to plot
metricValueIn: metric values
"""
fig = plt.figure(fignum)
plotDict = {}
plotDict.update(self.defaultPlotDict)
plotDict.update(userPlotDict)
norm = None
if plotDict['logScale']:
norm = 'log'
if plotDict['cmap'] is None:
cmap = cm.cubehelix
else:
cmap = plotDict['cmap']
if type(cmap) == str:
cmap = getattr(cm, cmap)
# Make colormap compatible with healpy
cmap = colors.LinearSegmentedColormap('cmap', cmap._segmentdata, cmap.N)
cmap.set_over(cmap(1.0))
cmap.set_under('w')
cmap.set_bad('gray')
if plotDict['zp'] is not None:
metricValue = metricValueIn - plotDict['zp']
elif plotDict['normVal'] is not None:
metricValue = metricValueIn / plotDict['normVal']
else:
metricValue = metricValueIn
if plotDict['percentileClip'] is not None:
pcMin, pcMax = percentileClipping(metricValue.compressed(),
percentile=plotDict['percentileClip'])
if plotDict['colorMin'] is None:
plotDict['colorMin'] = pcMin
if plotDict['colorMax'] is None:
plotDict['colorMax'] = pcMax
if (plotDict['colorMin'] is not None) or (plotDict['colorMax'] is not None):
clims = [plotDict['colorMin'], plotDict['colorMax']]
else:
clims = None
# Make sure there is some range on the colorbar
if clims is None:
if metricValue.compressed().size > 0:
clims = [metricValue.compressed().min(), metricValue.compressed().max()]
else:
clims = [-1, 1]
if clims[0] == clims[1]:
clims[0] = clims[0] - 1
clims[1] = clims[1] + 1
racenters = np.arange(plotDict['raMin'], plotDict['raMax'], plotDict['raLen'])
nframes = racenters.size
for i, racenter in enumerate(racenters):
if i == 0:
useTitle = plotDict['title'] + ' /n' + '%i < RA < %i' % (racenter - plotDict['raLen'],
racenter + plotDict['raLen'])
else:
useTitle = '%i < RA < %i' % (racenter - plotDict['raLen'], racenter + plotDict['raLen'])
hp.cartview(metricValue.filled(slicer.badval), title=useTitle, cbar=False,
min=clims[0], max=clims[1], flip='astro', rot=(racenter, 0, 0),
cmap=cmap, norm=norm, lonra=[-plotDict['raLen'], plotDict['raLen']],
latra=[plotDict['decMin'], plotDict['decMax']], sub=(nframes + 1, 1, i + 1), fig=fig)
hp.graticule(dpar=20, dmer=20, verbose=False)
# Add colorbar (not using healpy default colorbar because want more tickmarks).
fig = plt.gcf()
ax1 = fig.add_axes([0.1, .15, .8, .075]) # left, bottom, width, height
# Add label.
if plotDict['label'] is not None:
plt.figtext(0.8, 0.9, '%s' % plotDict['label'])
# Make the colorbar as a seperate figure,
# from http: //matplotlib.org/examples/api/colorbar_only.html
cnorm = colors.Normalize(vmin=clims[0], vmax=clims[1])
# supress silly colorbar warnings
with warnings.catch_warnings():
warnings.simplefilter("ignore")
cb1 = mpl.colorbar.ColorbarBase(ax1, cmap=cmap, norm=cnorm,
orientation='horizontal', format=plotDict['cbarFormat'])
cb1.set_label(plotDict['xlabel'])
cb1.ax.tick_params(labelsize=plotDict['labelsize'])
# If outputing to PDF, this fixes the colorbar white stripes
if plotDict['cbar_edge']:
cb1.solids.set_edgecolor("face")
fig = plt.gcf()
return fig.number
[docs]class LambertSkyMap(BasePlotter):
"""
Use basemap and contour to make a Lambertian projection.
Note that the plotDict can include a 'basemap' key with a dictionary of
kwargs to use with the call to Basemap.
"""
def __init__(self):
self.plotType = 'SkyMap'
self.objectPlotter = False
self.defaultPlotDict = {'basemap': {'projection': 'nplaea', 'boundinglat': 20,
'lon_0': 0., 'resolution': 'l', 'celestial': False},
'cbar': True, 'cmap': perceptual_rainbow, 'levels': 200,
'cbarFormat': '%i', 'cbar_edge': True, 'zp': None,
'normVal': None, 'percentileClip': None, 'colorMin': None,
'colorMax': None, 'linewidths': 0,
'fontsize': None, 'labelsize': None, 'figsize': None}
def __call__(self, metricValueIn, slicer, userPlotDict, fignum=None):
if 'ra' not in slicer.slicePoints or 'dec' not in slicer.slicePoints:
errMessage = 'SpatialSlicer must contain "ra" and "dec" in slicePoints metadata.'
errMessage += ' SlicePoints only contains keys %s.' % (slicer.slicePoints.keys())
raise ValueError(errMessage)
plotDict = {}
plotDict.update(self.defaultPlotDict)
plotDict.update(userPlotDict)
if plotDict['zp'] is not None:
metricValue = metricValueIn - plotDict['zp']
elif plotDict['normVal'] is not None:
metricValue = metricValueIn / plotDict['normVal']
else:
metricValue = metricValueIn
if plotDict['percentileClip'] is not None:
pcMin, pcMax = percentileClipping(metricValue.compressed(),
percentile=plotDict['percentileClip'])
if plotDict['colorMin'] is None:
plotDict['colorMin'] = pcMin
if plotDict['colorMax'] is None:
plotDict['colorMax'] = pcMax
if (plotDict['colorMin'] is not None) or (plotDict['colorMax'] is not None):
clims = [plotDict['colorMin'], plotDict['colorMax']]
else:
clims = None
# Make sure there is some range on the colorbar
if clims is None:
if metricValue.compressed().size > 0:
clims = [metricValue.compressed().min(), metricValue.compressed().max()]
else:
clims = [-1, 1]
if clims[0] == clims[1]:
clims[0] = clims[0] - 1
clims[1] = clims[1] + 1
# Calculate the levels to use for the contour
if np.size(plotDict['levels']) > 1:
levels = plotDict['levels']
else:
step = (clims[1] - clims[0]) / plotDict['levels']
levels = np.arange(clims[0], clims[1] + step, step)
fig = plt.figure(fignum, figsize=plotDict['figsize'])
ax = fig.add_subplot(111)
# Hide this extra dependency down here for now
# if using anaconda, to get basemap:
# conda install basemap
# Note, this should be possible without basemap, but there are
# matplotlib bugs:
# http: //stackoverflow.com/questions/31975303/matplotlib-tricontourf-with-an-axis-projection
try:
from mpl_toolkits.basemap import Basemap
except ModuleNotFoundError:
raise('To use this plotting function, please install Basemap into your python distribution')
m = Basemap(**plotDict['basemap'])
good = np.where(metricValue != slicer.badval)
# Contour the plot first to remove any anti-aliasing artifacts. Doesn't seem to work though. See:
# http: //stackoverflow.com/questions/15822159/aliasing-when-saving-matplotlib\
# -filled-contour-plot-to-pdf-or-eps
# tmpContour = m.contour(np.degrees(slicer.slicePoints['ra'][good]),
# np.degrees(slicer.slicePoints['dec'][good]),
# metricValue[good], levels, tri=True,
# cmap=plotDict['cmap'], ax=ax, latlon=True,
# lw=1)
CS = m.contourf(np.degrees(slicer.slicePoints['ra'][good]),
np.degrees(slicer.slicePoints['dec'][good]),
metricValue[good], levels, tri=True,
cmap=plotDict['cmap'], ax=ax, latlon=True)
# Try to fix the ugly pdf contour problem
for c in CS.collections:
c.set_edgecolor("face")
para = np.arange(0, 89, 20)
m.drawparallels(para, labels=para)
m.drawmeridians(np.arange(-180, 181, 60))
cb = plt.colorbar(CS, format=plotDict['cbarFormat'])
cb.set_label(plotDict['xlabel'])
if plotDict['labelsize'] is not None:
cb.ax.tick_params(labelsize=plotDict['labelsize'])
# Pop in an extra line to raise the title a bit
ax.set_title(plotDict['title']+'\n ')
# If outputing to PDF, this fixes the colorbar white stripes
if plotDict['cbar_edge']:
cb.solids.set_edgecolor("face")
return fig.number