Source code for lsst.sims.maf.plots.spatialPlotters

from builtins import zip
import numbers
import copy
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, applyZPNorm

from lsst.sims.utils import _equatorialFromGalactic, _healbin
from .perceptual_rainbow import makePRCmap
perceptual_rainbow = makePRCmap()
import numpy.ma as ma

__all__ = ['setColorLims', 'setColorMap', 'HealpixSkyMap', 'HealpixPowerSpectrum',
           'HealpixHistogram', 'OpsimHistogram', 'BaseHistogram',
           'BaseSkyMap', 'HealpixSDSSSkyMap', 'LambertSkyMap']

baseDefaultPlotDict = {'title': None, 'xlabel': None, 'label': None,
                       'logScale': False, 'percentileClip': None, 'normVal': None, 'zp': None,
                       'cbarFormat': None, 'cmap': perceptual_rainbow, 'cbar_edge': True, 'nTicks': 10,
                       'colorMin': None, 'colorMax': None,
                       'xMin': None, 'xMax': None, 'yMin': None, 'yMax': None,
                       'labelsize': None, 'fontsize': None, 'figsize': None, 'subplot': 111,
                       'maskBelow': None}


[docs]def setColorLims(metricValue, plotDict): """Set up color bar limits.""" # Use plot dict if these values are set. colorMin = plotDict['colorMin'] colorMax = plotDict['colorMax'] # If not, try to use percentile clipping. if (plotDict['percentileClip'] is not None) & (np.size(metricValue.compressed()) > 0): pcMin, pcMax = percentileClipping(metricValue.compressed(), percentile=plotDict['percentileClip']) if colorMin is None: colorMin = pcMin if colorMax is None: colorMax = pcMax # If not, just use the data limits. if colorMin is None: colorMin = metricValue.compressed().min() if colorMax is None: colorMax = metricValue.compressed().max() # But make sure there is some range on the colorbar if colorMin == colorMax: colorMin = colorMin - 0.5 colorMax = colorMax + 0.5 return np.sort([colorMin, colorMax])
[docs]def setColorMap(plotDict): cmap = plotDict['cmap'] if cmap is None: cmap = 'perceptual_rainbow' if type(cmap) == str: cmap = getattr(cm, cmap) # Set background and masked pixel colors default healpy white and gray. cmap = copy.copy(cmap) cmap.set_over(cmap(1.0)) cmap.set_under('w') cmap.set_bad('gray') return cmap
[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 = {} self.defaultPlotDict.update(baseDefaultPlotDict) self.defaultPlotDict.update({'rot': (0, 0, 0), 'flip': 'astro', 'coord': 'C', 'nside': 8, 'reduceFunc': np.mean}) # Note: for alt/az sky maps using the healpix plotter, you can use # {'rot': (90, 90, 90), 'flip': 'geo'} self.healpy_visufunc = hp.mollview self.healpy_visufunc_params = {} self.ax = None self.im = 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. """ # Override the default plotting parameters with user specified values. plotDict = {} plotDict.update(self.defaultPlotDict) plotDict.update(userPlotDict) # Check if we have a valid HEALpix slicer if 'Heal' in slicer.slicerName: # Update the metric data with zeropoint or normalization. metricValue = applyZPNorm(metricValueIn, plotDict) else: # Bin the values up on a healpix grid. metricValue = _healbin(slicer.slicePoints['ra'], slicer.slicePoints['dec'], metricValueIn.filled(slicer.badval), nside=plotDict['nside'], reduceFunc=plotDict['reduceFunc'], fillVal=slicer.badval) mask = np.zeros(metricValue.size) mask[np.where(metricValue == slicer.badval)] = 1 metricValue = ma.array(metricValue, mask=mask) metricValue = applyZPNorm(metricValue, plotDict) if plotDict['maskBelow'] is not None: toMask = np.where(metricValue <= plotDict['maskBelow'])[0] metricValue.mask[toMask] = True badval = hp.UNSEEN else: badval = slicer.badval # Generate a Mollweide full-sky plot. fig = plt.figure(fignum, figsize=plotDict['figsize']) # Set up color bar limits. clims = setColorLims(metricValue, plotDict) cmap = setColorMap(plotDict) # Set log scale? norm = None if plotDict['logScale']: norm = 'log' # 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") if plotDict['coord'] == 'C': notext = True else: notext = False visufunc_params = {'title': plotDict['title'], 'cbar': False, 'min': clims[0], 'max': clims[1], 'rot': plotDict['rot'], 'flip': plotDict['flip'], 'coord': plotDict['coord'], 'cmap': cmap, 'norm': norm, 'sub': plotDict['subplot'], 'fig':fig.number, 'notext': notext} visufunc_params.update(self.healpy_visufunc_params) self.healpy_visufunc(metricValue.filled(badval), **visufunc_params) # Add a graticule (grid) over the globe. hp.graticule(dpar=30, dmer=30, verbose=False) # Add colorbar (not using healpy default colorbar because we want more tickmarks). self.ax = plt.gca() im = self.ax.get_images()[0] # Add label. if plotDict['label'] is not None: plt.figtext(0.8, 0.8, '%s' % (plotDict['label'])) # Make a color bar. Supress silly colorbar warnings. with warnings.catch_warnings(): warnings.simplefilter("ignore") cb = plt.colorbar(im, shrink=0.75, aspect=25, pad=0.1, 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 = {} self.defaultPlotDict.update(baseDefaultPlotDict) self.defaultPlotDict.update({'maxl': None, 'removeDipole': True, 'linestyle': '-'})
[docs] def __call__(self, metricValue, slicer, userPlotDict, fignum=None): """ Generate and plot the power spectrum of metricValue (calculated on a healpix grid). """ if 'Healpix' not in slicer.slicerName: raise ValueError('HealpixPowerSpectrum for use with healpix metricBundles.') plotDict = {} plotDict.update(self.defaultPlotDict) plotDict.update(userPlotDict) fig = plt.figure(fignum, figsize=plotDict['figsize']) ax = fig.add_subplot(plotDict['subplot']) # 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 = {} self.defaultPlotDict.update(baseDefaultPlotDict) self.defaultPlotDict.update({'ylabel': 'Area (1000s of square degrees)', 'bins': None, 'binsize': None, 'cumulative': False, 'scale': None, 'linestyle': '-'}) self.baseHist = BaseHistogram()
[docs] def __call__(self, metricValue, slicer, userPlotDict, fignum=None): """ Histogram metricValue for all healpix points. """ if 'Healpix' not in slicer.slicerName: 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 = {} self.defaultPlotDict.update(baseDefaultPlotDict) self.defaultPlotDict.update({'ylabel': 'Number of Fields', 'yaxisformat': '%d', 'bins': None, 'binsize': None, 'cumulative': False, 'scale': 1.0, 'linestyle': '-'}) 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 = {} self.defaultPlotDict.update(baseDefaultPlotDict) self.defaultPlotDict.update({'ylabel': 'Count', 'bins': None, 'binsize': None, 'cumulative': False, 'scale': 1.0, 'yaxisformat': '%.3f', 'linestyle': '-'})
[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) metricValue = applyZPNorm(metricValueIn, plotDict) metricValue = metricValue.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 cases of trying to histogram single-valued data. # First we try to use the range specified by a user, if there is one. Then use the data if not. # all of this only works if plotDict is not cumulative. 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: # If user did not specify bins or binsize, then we try to figure out a good number of bins. bins = optimalBins(metricValue) # Generate plots. fig = plt.figure(fignum, figsize=plotDict['figsize']) ax = fig.add_subplot(plotDict['subplot']) # 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']) hist_ylims = plt.ylim() if n.max() > hist_ylims[1]: plt.ylim(top = n.max()) if n.min() < hist_ylims[0] and not plotDict['logScale']: plt.ylim(bottom = n.min()) # Fill in axes labels and limits. # Option to use 'scale' to turn y axis into area or other value. def mjrFormatter(y, pos): if not isinstance(plotDict['scale'], numbers.Number): raise ValueError('plotDict["scale"] must be a number to scale the y axis.') return plotDict['yaxisformat'] % (y * plotDict['scale']) ax.yaxis.set_major_formatter(FuncFormatter(mjrFormatter)) # Set optional x, y limits. if 'xMin' in plotDict: plt.xlim(left=plotDict['xMin']) if 'xMax' in plotDict: plt.xlim(right=plotDict['xMax']) if 'yMin' in plotDict: plt.ylim(bottom=plotDict['yMin']) if 'yMax' in plotDict: plt.ylim(top=plotDict['yMax']) # 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 = {} self.defaultPlotDict.update(baseDefaultPlotDict) self.defaultPlotDict.update({'projection': 'aitoff', 'radius': np.radians(1.75), 'alpha': 1.0, 'plotMask': False, 'metricIsColor': False, 'cbar': True, 'raCen': 0.0, 'mwZone': True, 'bgcolor': 'gray'}) 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) metricValue = applyZPNorm(metricValueIn, plotDict) fig = plt.figure(fignum, figsize=plotDict['figsize']) # other projections available include # ['aitoff', 'hammer', 'lambert', 'mollweide', 'polar', 'rectilinear'] ax = fig.add_subplot(plotDict['subplot'], 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'. good = ~metricValue.mask # Add ellipses at RA/Dec locations - but don't add colors yet. lon = -(slicer.slicePoints['ra'][good] - plotDict['raCen'] - np.pi) % (np.pi * 2) - np.pi ellipses = self._plot_tissot_ellipse(lon, slicer.slicePoints['dec'][good], plotDict['radius'], rasterized=True, ax=ax) if plotDict['metricIsColor']: current = None for ellipse, mVal in zip(ellipses, metricValue.data[good]): 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: # Determine color min/max values. metricValue.compressed = non-masked points. clims = setColorLims(metricValue, plotDict) # Determine whether or not to use auto-log scale. if plotDict['logScale'] == 'auto': if clims[0] > 0: if np.log10(clims[1]) - np.log10(clims[0]) > 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. #clims[0] = 10 ** (int(np.log10(clims[0]))) #clims[1] = 10 ** (int(np.log10(clims[1]))) 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[good]) 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']) tick_locator = ticker.MaxNLocator(nbins=plotDict['nTicks']) cb.locator = tick_locator cb.update_ticks() # 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_facecolor(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 = {} self.defaultPlotDict.update(baseDefaultPlotDict) self.defaultPlotDict.update({'cbarFormat': '%.2f', 'raMin': -90, 'raMax': 90, 'raLen': 45, 'decMin': -2., 'decMax': 2.})
[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 """ plotDict = {} plotDict.update(self.defaultPlotDict) plotDict.update(userPlotDict) metricValue = applyZPNorm(metricValueIn, plotDict) norm = None if plotDict['logScale']: norm = 'log' clims = setColorLims(metricValue, plotDict) cmap = setColorMap(plotDict) racenters = np.arange(plotDict['raMin'], plotDict['raMax'], plotDict['raLen']) nframes = racenters.size fig = plt.figure(fignum) # Do not specify or use plotDict['subplot'] because this is done in each call to hp.cartview. 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). ax = 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") cb = mpl.colorbar.ColorbarBase(ax, cmap=cmap, norm=cnorm, orientation='horizontal', format=plotDict['cbarFormat']) cb.set_label(plotDict['xlabel']) 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") fig = plt.gcf() return fig.number
def project_lambert(longitude, latitude): """Project from RA,dec to plane https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection """ # flipping the sign on latitude goes north pole or south pole centered r_polar = 2*np.cos((np.pi/2+latitude)/2.) # Add pi/2 so north is up theta_polar = longitude + np.pi/2 x = r_polar * np.cos(theta_polar) y = r_polar * np.sin(theta_polar) return x, y def draw_grat(ax): """Draw some graticule lines on an axis """ decs = np.radians(90.-np.array([20, 40, 60, 80])) ra = np.radians(np.arange(0, 361, 1)) for dec in decs: temp_dec = ra*0+dec x, y = project_lambert(ra, temp_dec) ax.plot(x, y, 'k--', alpha=0.5) ras = np.radians(np.arange(0, 360+45, 45)) dec = np.radians(90.-np.arange(0, 81, 1)) for ra in ras: temp_ra = dec*0 + ra x, y = project_lambert(temp_ra, dec) ax.plot(x, y, 'k--', alpha=0.5) for dec in decs: x, y = project_lambert(np.radians(45.), dec) ax.text(x, y, '%i' % np.round(np.degrees(dec))) return ax
[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 arbitrary kwargs to use with the call to Basemap. """ def __init__(self): self.plotType = 'SkyMap' self.objectPlotter = False self.defaultPlotDict = {} self.defaultPlotDict.update(baseDefaultPlotDict) self.defaultPlotDict.update({'basemap': {'projection': 'nplaea', 'boundinglat': 1, 'lon_0': 180, 'resolution': None, 'celestial': False, 'round': False}, 'levels': 200, 'cbarFormat': '%i', 'norm': 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) metricValue = applyZPNorm(metricValueIn, plotDict) clims = setColorLims(metricValue, plotDict) # 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(plotDict['subplot']) x, y = project_lambert(slicer.slicePoints['ra'], slicer.slicePoints['dec']) # 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']), # np.degrees(slicer.slicePoints['dec']), # metricValue.filled(np.min(clims)-1), levels, tri=True, # cmap=plotDict['cmap'], ax=ax, latlon=True, # lw=1) # Set masked values to be below the lowest contour level. if plotDict['norm'] == 'log': z_val = metricValue.filled(np.min(clims)-0.9) norm = colors.LogNorm(vmin=z_val.min(), vmax=z_val.max()) else: norm = plotDict['norm'] tcf = ax.tricontourf(x, y, metricValue.filled(np.min(clims)-0.9), levels, cmap=plotDict['cmap'], norm=norm) ax = draw_grat(ax) ax.set_xticks([]) ax.set_yticks([]) alt_limit = 10. x, y = project_lambert(0, np.radians(alt_limit)) max_val = np.max(np.abs([x, y])) ax.set_xlim([-max_val, max_val]) ax.set_ylim([-max_val, max_val]) # Try to fix the ugly pdf contour problem for c in tcf.collections: c.set_edgecolor("face") cb = plt.colorbar(tcf, 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