Source code for astrobase.plotbase

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# plotbase.py - Waqas Bhatti (wbhatti@astro.princeton.edu) - Feb 2016
# License: MIT.

'''
Contains various useful functions for plotting light curves and associated data.

'''

#############
## LOGGING ##
#############

import logging
from astrobase import log_sub, log_fmt, log_date_fmt

DEBUG = False
if DEBUG:
    level = logging.DEBUG
else:
    level = logging.INFO
LOGGER = logging.getLogger(__name__)
logging.basicConfig(
    level=level,
    style=log_sub,
    format=log_fmt,
    datefmt=log_date_fmt,
)

LOGDEBUG = LOGGER.debug
LOGINFO = LOGGER.info
LOGWARNING = LOGGER.warning
LOGERROR = LOGGER.error
LOGEXCEPTION = LOGGER.exception


#############
## IMPORTS ##
#############

import os
import os.path
import pickle
from io import BytesIO as Strio

import numpy as np
from numpy import min as npmin, max as npmax

# FIXME: enforce no display for now
import matplotlib
matplotlib.use('Agg')
dispok = False

import matplotlib.axes
import matplotlib.pyplot as plt

# for convolving DSS stamps to simulate seeing effects
import astropy.convolution as aconv

from astropy.io import fits as pyfits
from astropy.wcs import WCS
from astropy.visualization import (
    ZScaleInterval,
    ImageNormalize,
    LinearStretch
)


###################
## LOCAL IMPORTS ##
###################

from .lcmath import phase_magseries, phase_magseries_with_errs, \
    phase_bin_magseries, phase_bin_magseries_with_errs, \
    time_bin_magseries, time_bin_magseries_with_errs, sigclip_magseries, \
    normalize_magseries, find_lc_timegroups

from .lcfit.nonphysical import spline_fit_magseries

from .services.skyview import get_stamp


#########################
## SIMPLE LIGHT CURVES ##
#########################

[docs]def plot_magseries(times, mags, magsarefluxes=False, errs=None, out=None, sigclip=30.0, normto='globalmedian', normmingap=4.0, timebin=None, yrange=None, segmentmingap=100.0, plotdpi=100): '''This plots a magnitude/flux time-series. Parameters ---------- times,mags : np.array The mag/flux time-series to plot as a function of time. magsarefluxes : bool Indicates if the input `mags` array is actually an array of flux measurements instead of magnitude measurements. If this is set to True, then the plot y-axis will be set as appropriate for mag or fluxes. In addition: - if `normto` is 'zero', then the median flux is divided from each observation's flux value to yield normalized fluxes with 1.0 as the global median. - if `normto` is 'globalmedian', then the global median flux value across the entire time series is multiplied with each measurement. - if `norm` is set to a `float`, then this number is multiplied with the flux value for each measurement. errs : np.array or None If this is provided, contains the measurement errors associated with each measurement of flux/mag in time-series. Providing this kwarg will add errbars to the output plot. out : str or StringIO/BytesIO object or None Sets the output type and target: - If `out` is a string, will save the plot to the specified file name. - If `out` is a StringIO/BytesIO object, will save the plot to that file handle. This can be useful to carry out additional operations on the output binary stream, or convert it to base64 text for embedding in HTML pages. - If `out` is None, will save the plot to a file called 'magseries-plot.png' in the current working directory. sigclip : float or int or sequence of two floats/ints or None If a single float or int, a symmetric sigma-clip will be performed using the number provided as the sigma-multiplier to cut out from the input time-series. If a list of two ints/floats is provided, the function will perform an 'asymmetric' sigma-clip. The first element in this list is the sigma value to use for fainter flux/mag values; the second element in this list is the sigma value to use for brighter flux/mag values. For example, `sigclip=[10., 3.]`, will sigclip out greater than 10-sigma dimmings and greater than 3-sigma brightenings. Here the meaning of "dimming" and "brightening" is set by *physics* (not the magnitude system), which is why the `magsarefluxes` kwarg must be correctly set. If `sigclip` is None, no sigma-clipping will be performed, and the time-series (with non-finite elems removed) will be passed through to the output. normto : {'globalmedian', 'zero'} or a float Sets the normalization target:: 'globalmedian' -> norms each mag to the global median of the LC column 'zero' -> norms each mag to zero a float -> norms each mag to this specified float value. normmingap : float This defines how much the difference between consecutive measurements is allowed to be to consider them as parts of different timegroups. By default it is set to 4.0 days. timebin : float or None The bin size to use to group together measurements closer than this amount in time. This is in seconds. If this is None, no time-binning will be performed. yrange : list of two floats or None This is used to provide a custom y-axis range to the plot. If None, will automatically determine y-axis range. segmentmingap : float or None This controls the minimum length of time (in days) required to consider a timegroup in the light curve as a separate segment. This is useful when the light curve consists of measurements taken over several seasons, so there's lots of dead space in the plot that can be cut out to zoom in on the interesting stuff. If `segmentmingap` is not None, the magseries plot will be cut in this way and the x-axis will show these breaks. plotdpi : int Sets the resolution in DPI for PNG plots (default = 100). Returns ------- str or BytesIO/StringIO object Returns based on the input: - If `out` is a str or None, the path to the generated plot file is returned. - If `out` is a StringIO/BytesIO object, will return the StringIO/BytesIO object to which the plot was written. ''' # sigclip the magnitude timeseries stimes, smags, serrs = sigclip_magseries(times, mags, errs, magsarefluxes=magsarefluxes, sigclip=sigclip) # now we proceed to binning if timebin and errs is not None: binned = time_bin_magseries_with_errs(stimes, smags, serrs, binsize=timebin) btimes, bmags, berrs = (binned['binnedtimes'], binned['binnedmags'], binned['binnederrs']) elif timebin and errs is None: binned = time_bin_magseries(stimes, smags, binsize=timebin) btimes, bmags, berrs = binned['binnedtimes'], binned['binnedmags'], None else: btimes, bmags, berrs = stimes, smags, serrs # check if we need to normalize if normto is not False: btimes, bmags = normalize_magseries(btimes, bmags, normto=normto, magsarefluxes=magsarefluxes, mingap=normmingap) btimeorigin = btimes.min() btimes = btimes - btimeorigin ################################## ## FINALLY PLOT THE LIGHT CURVE ## ################################## # if we're going to plot with segment gaps highlighted, then find the gaps if segmentmingap is not None: ntimegroups, timegroups = find_lc_timegroups(btimes, mingap=segmentmingap) # get the yrange for all the plots if it's given if yrange and isinstance(yrange,(list,tuple)) and len(yrange) == 2: ymin, ymax = yrange # if it's not given, figure it out else: # the plot y limits are just 0.05 mags on each side if mags are used if not magsarefluxes: ymin, ymax = (bmags.min() - 0.05, bmags.max() + 0.05) # if we're dealing with fluxes, limits are 2% of the flux range per side else: ycov = bmags.max() - bmags.min() ymin = bmags.min() - 0.02*ycov ymax = bmags.max() + 0.02*ycov # if we're supposed to make the plot segment-aware (i.e. gaps longer than # segmentmingap will be cut out) if segmentmingap and ntimegroups > 1: LOGINFO('%s time groups found' % ntimegroups) # our figure is now a multiple axis plot # the aspect ratio is a bit wider fig, axes = plt.subplots(1,ntimegroups,sharey=True) fig.set_size_inches(10,4.8) axes = np.ravel(axes) # now go through each axis and make the plots for each timegroup for timegroup, ax, axind in zip(timegroups, axes, range(len(axes))): tgtimes = btimes[timegroup] tgmags = bmags[timegroup] if berrs: tgerrs = berrs[timegroup] else: tgerrs = None LOGINFO('axes: %s, timegroup %s: JD %.3f to %.3f' % ( axind, axind+1, btimeorigin + tgtimes.min(), btimeorigin + tgtimes.max()) ) ax.errorbar(tgtimes, tgmags, fmt='go', yerr=tgerrs, markersize=2.0, markeredgewidth=0.0, ecolor='grey', capsize=0) # don't use offsets on any xaxis ax.get_xaxis().get_major_formatter().set_useOffset(False) # fix the ticks to use no yoffsets and remove right spines for first # axes instance if axind == 0: ax.get_yaxis().get_major_formatter().set_useOffset(False) ax.spines['right'].set_visible(False) ax.yaxis.tick_left() # remove the right and left spines for the other axes instances elif 0 < axind < (len(axes)-1): ax.spines['right'].set_visible(False) ax.spines['left'].set_visible(False) ax.tick_params(right='off', labelright='off', left='off',labelleft='off') # make the left spines invisible for the last axes instance elif axind == (len(axes)-1): ax.spines['left'].set_visible(False) ax.spines['right'].set_visible(True) ax.yaxis.tick_right() # set the yaxis limits if not magsarefluxes: ax.set_ylim(ymax, ymin) else: ax.set_ylim(ymin, ymax) # now figure out the xaxis ticklabels and ranges tgrange = tgtimes.max() - tgtimes.min() if tgrange < 10.0: ticklocations = [tgrange/2.0] ax.set_xlim(npmin(tgtimes) - 0.5, npmax(tgtimes) + 0.5) elif 10.0 < tgrange < 30.0: ticklocations = np.linspace(tgtimes.min()+5.0, tgtimes.max()-5.0, num=2) ax.set_xlim(npmin(tgtimes) - 2.0, npmax(tgtimes) + 2.0) elif 30.0 < tgrange < 100.0: ticklocations = np.linspace(tgtimes.min()+10.0, tgtimes.max()-10.0, num=3) ax.set_xlim(npmin(tgtimes) - 2.5, npmax(tgtimes) + 2.5) else: ticklocations = np.linspace(tgtimes.min()+20.0, tgtimes.max()-20.0, num=3) ax.set_xlim(npmin(tgtimes) - 3.0, npmax(tgtimes) + 3.0) ax.xaxis.set_ticks([int(x) for x in ticklocations]) # done with plotting all the sub axes # make the distance between sub plots smaller plt.subplots_adjust(wspace=0.07) # make the overall x and y labels fig.text(0.5, 0.00, 'JD - %.3f (not showing gaps > %.2f d)' % (btimeorigin, segmentmingap), ha='center') if not magsarefluxes: fig.text(0.02, 0.5, 'magnitude', va='center', rotation='vertical') else: fig.text(0.02, 0.5, 'flux', va='center', rotation='vertical') # make normal figure otherwise else: fig = plt.figure() fig.set_size_inches(7.5,4.8) plt.errorbar(btimes, bmags, fmt='go', yerr=berrs, markersize=2.0, markeredgewidth=0.0, ecolor='grey', capsize=0) # make a grid plt.grid(color='#a9a9a9', alpha=0.9, zorder=0, linewidth=1.0, linestyle=':') # fix the ticks to use no offsets plt.gca().get_yaxis().get_major_formatter().set_useOffset(False) plt.gca().get_xaxis().get_major_formatter().set_useOffset(False) plt.xlabel('JD - %.3f' % btimeorigin) # set the yaxis limits and labels if not magsarefluxes: plt.ylim(ymax, ymin) plt.ylabel('magnitude') else: plt.ylim(ymin, ymax) plt.ylabel('flux') is_Strio = isinstance(out, Strio) # write the plot out to a file if requested if out and not is_Strio: if out.endswith('.png'): plt.savefig(out,bbox_inches='tight',dpi=plotdpi) else: plt.savefig(out,bbox_inches='tight') plt.close() return os.path.abspath(out) elif out and is_Strio: plt.savefig(out, bbox_inches='tight', dpi=plotdpi, format='png') return out elif not out and dispok: plt.show() plt.close() return else: LOGWARNING('no output file specified and no $DISPLAY set, ' 'saving to magseries-plot.png in current directory') outfile = 'magseries-plot.png' plt.savefig(outfile,bbox_inches='tight',dpi=plotdpi) plt.close() return os.path.abspath(outfile)
######################### ## PHASED LIGHT CURVES ## #########################
[docs]def plot_phased_magseries(times, mags, period, epoch='min', fitknotfrac=0.01, errs=None, magsarefluxes=False, normto='globalmedian', normmingap=4.0, sigclip=30.0, phasewrap=True, phasesort=True, phasebin=None, plotphaselim=(-0.8,0.8), yrange=None, xtimenotphase=False, xaxlabel='phase', yaxlabel=None, modelmags=None, modeltimes=None, modelerrs=None, outfile=None, plotdpi=100): '''Plots a phased magnitude/flux time-series using the period provided. Parameters ---------- times,mags : np.array The mag/flux time-series to plot as a function of phase given `period`. period : float The period to use to phase-fold the time-series. Should be the same unit as `times` (usually in days) epoch : 'min' or float or None This indicates how to get the epoch to use for phasing the light curve: - If None, uses the `min(times)` as the epoch for phasing. - If epoch is the string 'min', then fits a cubic spline to the phased light curve using `min(times)` as the initial epoch, finds the magnitude/flux minimum of this phased light curve fit, and finally uses the that time value as the epoch. This is useful for plotting planetary transits and eclipsing binary phased light curves so that phase 0.0 corresponds to the mid-center time of primary eclipse (or transit). - If epoch is a float, then uses that directly to phase the light curve and as the epoch of the phased mag series plot. fitknotfrac : float If `epoch='min'`, this function will attempt to fit a cubic spline to the phased light curve to find a time of light minimum as phase 0.0. This kwarg sets the number of knots to generate the spline as a fraction of the total number of measurements in the input time-series. By default, this is set so that 100 knots are used to generate a spline for fitting the phased light curve consisting of 10000 measurements. errs : np.array or None If this is provided, contains the measurement errors associated with each measurement of flux/mag in time-series. Providing this kwarg will add errbars to the output plot. magsarefluxes : bool Indicates if the input `mags` array is actually an array of flux measurements instead of magnitude measurements. If this is set to True, then the plot y-axis will be set as appropriate for mag or fluxes. normto : {'globalmedian', 'zero'} or a float Sets the normalization target:: 'globalmedian' -> norms each mag to the global median of the LC column 'zero' -> norms each mag to zero a float -> norms each mag to this specified float value. normmingap : float This defines how much the difference between consecutive measurements is allowed to be to consider them as parts of different timegroups. By default it is set to 4.0 days. sigclip : float or int or sequence of two floats/ints or None If a single float or int, a symmetric sigma-clip will be performed using the number provided as the sigma-multiplier to cut out from the input time-series. If a list of two ints/floats is provided, the function will perform an 'asymmetric' sigma-clip. The first element in this list is the sigma value to use for fainter flux/mag values; the second element in this list is the sigma value to use for brighter flux/mag values. For example, `sigclip=[10., 3.]`, will sigclip out greater than 10-sigma dimmings and greater than 3-sigma brightenings. Here the meaning of "dimming" and "brightening" is set by *physics* (not the magnitude system), which is why the `magsarefluxes` kwarg must be correctly set. If `sigclip` is None, no sigma-clipping will be performed, and the time-series (with non-finite elems removed) will be passed through to the output. phasewrap : bool If this is True, the phased time-series will be wrapped around phase 0.0. phasesort : bool If this is True, the phased time-series will be sorted in phase. phasebin : float or None If this is provided, indicates the bin size to use to group together measurements closer than this amount in phase. This is in units of phase. The binned phased light curve will be overplotted on top of the phased light curve. Useful for when one has many measurement points and needs to pick out a small trend in an otherwise noisy phased light curve. plotphaselim : sequence of two floats or None The x-axis limits to use when making the phased light curve plot. By default, this is (-0.8, 0.8), which places phase 0.0 at the center of the plot and covers approximately two cycles in phase to make any trends clear. yrange : list of two floats or None This is used to provide a custom y-axis range to the plot. If None, will automatically determine y-axis range. xtimenotphase : bool If True, the x-axis gets units of time (multiplies phase by period). xaxlabel : str Sets the label for the x-axis. yaxlabel : str or None Sets the label for the y-axis. If this is None, the appropriate label will be used based on the value of the `magsarefluxes` kwarg. modeltimes,modelmags,modelerrs : np.array or None If all of these are provided, then this function will overplot the values of modeltimes and modelmags on top of the actual phased light curve. This is useful for plotting variability models on top of the light curve (e.g. plotting a Mandel-Agol transit model over the actual phased light curve. These arrays will be phased using the already provided period and epoch. outfile : str or StringIO/BytesIO or matplotlib.axes.Axes or None - a string filename for the file where the plot will be written. - a StringIO/BytesIO object to where the plot will be written. - a matplotlib.axes.Axes object to where the plot will be written. - if None, plots to 'magseries-phased-plot.png' in current dir. plotdpi : int Sets the resolution in DPI for PNG plots (default = 100). Returns ------- str or StringIO/BytesIO or matplotlib.axes.Axes This returns based on the input: - If `outfile` is a str or None, the path to the generated plot file is returned. - If `outfile` is a StringIO/BytesIO object, will return the StringIO/BytesIO object to which the plot was written. - If `outfile` is a matplotlib.axes.Axes object, will return the Axes object with the plot elements added to it. One can then directly include this Axes object in some other Figure. ''' # sigclip the magnitude timeseries stimes, smags, serrs = sigclip_magseries(times, mags, errs, magsarefluxes=magsarefluxes, sigclip=sigclip) # check if we need to normalize if normto is not False: stimes, smags = normalize_magseries(stimes, smags, normto=normto, magsarefluxes=magsarefluxes, mingap=normmingap) if ( isinstance(modelmags, np.ndarray) and isinstance(modeltimes, np.ndarray) ): stimes, smags = normalize_magseries(modeltimes, modelmags, normto=normto, magsarefluxes=magsarefluxes, mingap=normmingap) # figure out the epoch, if it's None, use the min of the time if epoch is None: epoch = stimes.min() # if the epoch is 'min', then fit a spline to the light curve phased # using the min of the time, find the fit mag minimum and use the time for # that as the epoch elif isinstance(epoch, str) and epoch == 'min': try: spfit = spline_fit_magseries(stimes, smags, serrs, period, knotfraction=fitknotfrac) epoch = spfit['fitinfo']['fitepoch'] if len(epoch) != 1: epoch = epoch[0] except Exception: LOGEXCEPTION('spline fit failed, using min(times) as epoch') epoch = npmin(stimes) # now phase the data light curve (and optionally, phase bin the light curve) if errs is not None: phasedlc = phase_magseries_with_errs(stimes, smags, serrs, period, epoch, wrap=phasewrap, sort=phasesort) plotphase = phasedlc['phase'] plotmags = phasedlc['mags'] ploterrs = phasedlc['errs'] # if we're supposed to bin the phases, do so if phasebin: binphasedlc = phase_bin_magseries_with_errs(plotphase, plotmags, ploterrs, binsize=phasebin) binplotphase = binphasedlc['binnedphases'] binplotmags = binphasedlc['binnedmags'] binploterrs = binphasedlc['binnederrs'] else: phasedlc = phase_magseries(stimes, smags, period, epoch, wrap=phasewrap, sort=phasesort) plotphase = phasedlc['phase'] plotmags = phasedlc['mags'] ploterrs = None # if we're supposed to bin the phases, do so if phasebin: binphasedlc = phase_bin_magseries(plotphase, plotmags, binsize=phasebin) binplotphase = binphasedlc['binnedphases'] binplotmags = binphasedlc['binnedmags'] binploterrs = None # phase the model light curve modelplotphase, modelplotmags = None, None if ( isinstance(modelerrs,np.ndarray) and isinstance(modeltimes,np.ndarray) and isinstance(modelmags,np.ndarray) ): modelphasedlc = phase_magseries_with_errs(modeltimes, modelmags, modelerrs, period, epoch, wrap=phasewrap, sort=phasesort) modelplotphase = modelphasedlc['phase'] modelplotmags = modelphasedlc['mags'] # note that we never will phase-bin the model (no point). elif ( not isinstance(modelerrs,np.ndarray) and isinstance(modeltimes,np.ndarray) and isinstance(modelmags,np.ndarray) ): modelphasedlc = phase_magseries(modeltimes, modelmags, period, epoch, wrap=phasewrap, sort=phasesort) modelplotphase = modelphasedlc['phase'] modelplotmags = modelphasedlc['mags'] # finally, make the plots # check if the outfile is actually an Axes object if isinstance(outfile, matplotlib.axes.Axes): ax = outfile # otherwise, it's just a normal file or StringIO/BytesIO else: fig = plt.figure() fig.set_size_inches(7.5,4.8) ax = plt.gca() if xtimenotphase: plotphase *= period if phasebin: ax.errorbar(plotphase, plotmags, fmt='o', color='#B2BEB5', yerr=ploterrs, markersize=3.0, markeredgewidth=0.0, ecolor='#B2BEB5', capsize=0) if xtimenotphase: binplotphase *= period ax.errorbar(binplotphase, binplotmags, fmt='bo', yerr=binploterrs, markersize=5.0, markeredgewidth=0.0, ecolor='#B2BEB5', capsize=0) else: ax.errorbar(plotphase, plotmags, fmt='ko', yerr=ploterrs, markersize=3.0, markeredgewidth=0.0, ecolor='#B2BEB5', capsize=0) if (isinstance(modelplotphase, np.ndarray) and isinstance(modelplotmags, np.ndarray)): if xtimenotphase: modelplotphase *= period ax.plot(modelplotphase, modelplotmags, zorder=5, linewidth=0.5, alpha=0.9, color='#181c19') # make a grid ax.grid(color='#a9a9a9', alpha=0.9, zorder=0, linewidth=1.0, linestyle=':') # make lines for phase 0.0, 0.5, and -0.5 ax.axvline(0.0,alpha=0.9,linestyle='dashed',color='g') if not xtimenotphase: ax.axvline(-0.5,alpha=0.9,linestyle='dashed',color='g') ax.axvline(0.5,alpha=0.9,linestyle='dashed',color='g') else: ax.axvline(-period*0.5,alpha=0.9,linestyle='dashed',color='g') ax.axvline(period*0.5,alpha=0.9,linestyle='dashed',color='g') # fix the ticks to use no offsets ax.get_yaxis().get_major_formatter().set_useOffset(False) ax.get_xaxis().get_major_formatter().set_useOffset(False) # get the yrange if yrange and isinstance(yrange,(list,tuple)) and len(yrange) == 2: ymin, ymax = yrange else: ymin, ymax = ax.get_ylim() # set the y axis labels and range if not yaxlabel: if not magsarefluxes: ax.set_ylim(ymax, ymin) yaxlabel = 'magnitude' else: ax.set_ylim(ymin, ymax) yaxlabel = 'flux' # set the x axis limit if not plotphaselim: ax.set_xlim((npmin(plotphase)-0.1, npmax(plotphase)+0.1)) else: if xtimenotphase: ax.set_xlim((period*plotphaselim[0],period*plotphaselim[1])) else: ax.set_xlim((plotphaselim[0],plotphaselim[1])) # set up the axis labels and plot title ax.set_xlabel(xaxlabel) ax.set_ylabel(yaxlabel) ax.set_title('period: %.6f d - epoch: %.6f' % (period, epoch)) LOGINFO('using period: %.6f d and epoch: %.6f' % (period, epoch)) # check if the output filename is actually an instance of StringIO is_Strio = isinstance(outfile, Strio) # make the figure if (outfile and not is_Strio and not isinstance(outfile, matplotlib.axes.Axes)): if outfile.endswith('.png'): fig.savefig(outfile, bbox_inches='tight', dpi=plotdpi) else: fig.savefig(outfile, bbox_inches='tight') plt.close() return period, epoch, os.path.abspath(outfile) elif outfile and is_Strio: fig.savefig(outfile, bbox_inches='tight', dpi=plotdpi, format='png') return outfile elif outfile and isinstance(outfile, matplotlib.axes.Axes): return outfile elif not outfile and dispok: plt.show() plt.close() return period, epoch else: LOGWARNING('no output file specified and no $DISPLAY set, ' 'saving to magseries-phased-plot.png in current directory') outfile = 'magseries-phased-plot.png' plt.savefig(outfile, bbox_inches='tight', dpi=plotdpi) plt.close() return period, epoch, os.path.abspath(outfile)
########################## ## PLOTTING FITS IMAGES ## ##########################
[docs]def skyview_stamp(ra, decl, survey='DSS2 Red', scaling='Linear', sizepix=300, flip=True, convolvewith=None, forcefetch=False, cachedir='~/.astrobase/stamp-cache', timeout=45.0, retry_failed=False, savewcsheader=True, verbose=False): '''This downloads a DSS FITS stamp centered on the coordinates specified. This wraps the function :py:func:`astrobase.services.skyview.get_stamp`, which downloads Digitized Sky Survey stamps in FITS format from the NASA SkyView service: https://skyview.gsfc.nasa.gov/current/cgi/query.pl Also adds some useful operations on top of the FITS file returned. Parameters ---------- ra,decl : float The center coordinates for the stamp in decimal degrees. survey : str The survey name to get the stamp from. This is one of the values in the 'SkyView Surveys' option boxes on the SkyView webpage. Currently, we've only tested using 'DSS2 Red' as the value for this kwarg, but the other ones should work in principle. scaling : str This is the pixel value scaling function to use. Can be any of the strings ("Log", "Linear", "Sqrt", "HistEq"). sizepix : int Size of the requested stamp, in pixels. (DSS scale is ~1arcsec/px). flip : bool Will flip the downloaded image top to bottom. This should usually be True because matplotlib and FITS have different image coord origin conventions. Alternatively, set this to False and use the `origin='lower'` in any call to `matplotlib.pyplot.imshow` when plotting this image. convolvewith : astropy.convolution Kernel object or None If `convolvewith` is an astropy.convolution Kernel object from: http://docs.astropy.org/en/stable/convolution/kernels.html then, this function will return the stamp convolved with that kernel. This can be useful to see effects of wide-field telescopes (like the HATNet and HATSouth lenses) degrading the nominal 1 arcsec/px of DSS, causing blending of targets and any variability. forcefetch : bool If True, will disregard any existing cached copies of the stamp already downloaded corresponding to the requested center coordinates and redownload the FITS from the SkyView service. cachedir : str This is the path to the astrobase cache directory. All downloaded FITS stamps are stored here as .fits.gz files so we can immediately respond with the cached copy when a request is made for a coordinate center that's already been downloaded. timeout : float Sets the timeout in seconds to wait for a response from the NASA SkyView service. retry_failed : bool If the initial request to SkyView fails, and this is True, will retry until it succeeds. savewcsheader : bool If this is True, also returns the WCS header of the downloaded FITS stamp in addition to the FITS image itself. Useful for projecting object coordinates onto image xy coordinates for visualization. verbose : bool If True, indicates progress. Returns ------- tuple or array or None This returns based on the value of `savewcsheader`: - If `savewcsheader=True`, returns a tuple: (FITS stamp image as a numpy array, FITS header) - If `savewcsheader=False`, returns only the FITS stamp image as numpy array. - If the stamp retrieval fails, returns None. ''' stampdict = get_stamp(ra, decl, survey=survey, scaling=scaling, sizepix=sizepix, forcefetch=forcefetch, cachedir=cachedir, timeout=timeout, retry_failed=retry_failed, verbose=verbose) # # DONE WITH FETCHING STUFF # if stampdict: # open the frame stampfits = pyfits.open(stampdict['fitsfile']) header = stampfits[0].header frame = stampfits[0].data stampfits.close() # finally, we can process the frame if flip: frame = np.flipud(frame) if verbose: LOGINFO('fetched stamp successfully for (%.3f, %.3f)' % (ra, decl)) if convolvewith: convolved = aconv.convolve(frame, convolvewith) if savewcsheader: return convolved, header else: return convolved else: if savewcsheader: return frame, header else: return frame else: LOGERROR('could not fetch the requested stamp for ' 'coords: (%.3f, %.3f) from survey: %s and scaling: %s' % (ra, decl, survey, scaling)) return None
[docs]def fits_finder_chart( fitsfile, outfile, fitsext=0, wcsfrom=None, scale=ZScaleInterval(), stretch=LinearStretch(), colormap=plt.cm.gray_r, findersize=None, finder_coordlimits=None, overlay_ra=None, overlay_decl=None, overlay_pltopts={'marker':'o', 'markersize':10.0, 'markerfacecolor':'none', 'markeredgewidth':2.0, 'markeredgecolor':'red'}, overlay_zoomcontain=False, grid=False, gridcolor='k' ): '''This makes a finder chart for a given FITS with an optional object position overlay. Parameters ---------- fitsfile : str `fitsfile` is the FITS file to use to make the finder chart. outfile : str `outfile` is the name of the output file. This can be a png or pdf or whatever else matplotlib can write given a filename and extension. fitsext : int Sets the FITS extension in `fitsfile` to use to extract the image array from. wcsfrom : str or None If `wcsfrom` is None, the WCS to transform the RA/Dec to pixel x/y will be taken from the FITS header of `fitsfile`. If this is not None, it must be a FITS or similar file that contains a WCS header in its first extension. scale : astropy.visualization.Interval object `scale` sets the normalization for the FITS pixel values. This is an astropy.visualization Interval object. See http://docs.astropy.org/en/stable/visualization/normalization.html for details on `scale` and `stretch` objects. stretch : astropy.visualization.Stretch object `stretch` sets the stretch function for mapping FITS pixel values to output pixel values. This is an astropy.visualization Stretch object. See http://docs.astropy.org/en/stable/visualization/normalization.html for details on `scale` and `stretch` objects. colormap : matplotlib Colormap object `colormap` is a matplotlib color map object to use for the output image. findersize : None or tuple of two ints If `findersize` is None, the output image size will be set by the NAXIS1 and NAXIS2 keywords in the input `fitsfile` FITS header. Otherwise, `findersize` must be a tuple with the intended x and y size of the image in inches (all output images will use a DPI = 100). finder_coordlimits : list of four floats or None If not None, `finder_coordlimits` sets x and y limits for the plot, effectively zooming it in if these are smaller than the dimensions of the FITS image. This should be a list of the form: [minra, maxra, mindecl, maxdecl] all in decimal degrees. overlay_ra, overlay_decl : np.array or None `overlay_ra` and `overlay_decl` are ndarrays containing the RA and Dec values to overplot on the image as an overlay. If these are both None, then no overlay will be plotted. overlay_pltopts : dict `overlay_pltopts` controls how the overlay points will be plotted. This a dict with standard matplotlib marker, etc. kwargs as key-val pairs, e.g. 'markersize', 'markerfacecolor', etc. The default options make red outline circles at the location of each object in the overlay. overlay_zoomcontain : bool `overlay_zoomcontain` controls if the finder chart will be zoomed to just contain the overlayed points. Everything outside the footprint of these points will be discarded. grid : bool `grid` sets if a grid will be made on the output image. gridcolor : str `gridcolor` sets the color of the grid lines. This is a usual matplotib color spec string. Returns ------- str or None The filename of the generated output image if successful. None otherwise. ''' # read in the FITS file if wcsfrom is None: hdulist = pyfits.open(fitsfile) img, hdr = hdulist[fitsext].data, hdulist[fitsext].header hdulist.close() frameshape = (hdr['NAXIS1'], hdr['NAXIS2']) w = WCS(hdr) elif os.path.exists(wcsfrom): hdulist = pyfits.open(fitsfile) img, hdr = hdulist[fitsext].data, hdulist[fitsext].header hdulist.close() frameshape = (hdr['NAXIS1'], hdr['NAXIS2']) w = WCS(wcsfrom) else: LOGERROR('could not determine WCS info for input FITS: %s' % fitsfile) return None # use the frame shape to set the output PNG's dimensions if findersize is None: fig = plt.figure(figsize=(frameshape[0]/100.0, frameshape[1]/100.0)) else: fig = plt.figure(figsize=findersize) # set the coord limits if zoomcontain is True # we'll leave 30 arcseconds of padding on each side if (overlay_zoomcontain and overlay_ra is not None and overlay_decl is not None): finder_coordlimits = [overlay_ra.min()-30.0/3600.0, overlay_ra.max()+30.0/3600.0, overlay_decl.min()-30.0/3600.0, overlay_decl.max()+30.0/3600.0] # set the coordinate limits if provided if finder_coordlimits and isinstance(finder_coordlimits, (list,tuple)): minra, maxra, mindecl, maxdecl = finder_coordlimits cntra, cntdecl = (minra + maxra)/2.0, (mindecl + maxdecl)/2.0 pixelcoords = w.all_world2pix([[minra, mindecl], [maxra, maxdecl], [cntra, cntdecl]],1) x1, y1, x2, y2 = (int(pixelcoords[0,0]), int(pixelcoords[0,1]), int(pixelcoords[1,0]), int(pixelcoords[1,1])) xmin = x1 if x1 < x2 else x2 xmax = x2 if x2 > x1 else x1 ymin = y1 if y1 < y2 else y2 ymax = y2 if y2 > y1 else y1 # create a new WCS with the same transform but new center coordinates whdr = w.to_header() whdr['CRPIX1'] = (xmax - xmin)/2 whdr['CRPIX2'] = (ymax - ymin)/2 whdr['CRVAL1'] = cntra whdr['CRVAL2'] = cntdecl whdr['NAXIS1'] = xmax - xmin whdr['NAXIS2'] = ymax - ymin w = WCS(whdr) else: xmin, xmax, ymin, ymax = 0, hdr['NAXIS2'], 0, hdr['NAXIS1'] # add the axes with the WCS projection # this should automatically handle subimages because we fix the WCS # appropriately above for these fig.add_subplot(111,projection=w) if scale is not None and stretch is not None: norm = ImageNormalize(img, interval=scale, stretch=stretch) plt.imshow(img[ymin:ymax,xmin:xmax], origin='lower', cmap=colormap, norm=norm) else: plt.imshow(img[ymin:ymax,xmin:xmax], origin='lower', cmap=colormap) # handle additional options if grid: plt.grid(color=gridcolor,ls='solid',lw=1.0) # handle the object overlay if overlay_ra is not None and overlay_decl is not None: our_pltopts = dict( transform=plt.gca().get_transform('fk5'), marker='o', markersize=10.0, markerfacecolor='none', markeredgewidth=2.0, markeredgecolor='red', rasterized=True, linestyle='none' ) if overlay_pltopts is not None and isinstance(overlay_pltopts, dict): our_pltopts.update(overlay_pltopts) plt.gca().set_autoscale_on(False) plt.gca().plot(overlay_ra, overlay_decl, **our_pltopts) plt.xlabel('Right Ascension [deg]') plt.ylabel('Declination [deg]') # get the x and y axes objects to fix the ticks xax = plt.gca().coords[0] yax = plt.gca().coords[1] yax.set_major_formatter('d.ddd') xax.set_major_formatter('d.ddd') # save the figure plt.savefig(outfile, dpi=100.0) plt.close('all') return outfile
################## ## PERIODOGRAMS ## ################## PLOTYLABELS = {'gls':'Generalized Lomb-Scargle normalized power', 'pdm':r'Stellingwerf PDM $\Theta$', 'aov':r'Schwarzenberg-Czerny AoV $\Theta$', 'mav':r'Schwarzenberg-Czerny AoVMH $\Theta$', 'bls':'Box Least-squared Search SR', 'acf':'Autocorrelation Function', 'win':'Lomb-Scargle normalized power', 'ext':'External period-finder power', 'tls':'Transit Least-Squares SDE'} METHODLABELS = {'gls':'Generalized Lomb-Scargle periodogram', 'pdm':'Stellingwerf phase-dispersion minimization', 'aov':'Schwarzenberg-Czerny AoV', 'mav':'Schwarzenberg-Czerny AoV multi-harmonic', 'bls':'Box Least-squared Search', 'acf':'McQuillan+ ACF Period Search', 'win':'Timeseries Sampling Lomb-Scargle periodogram', 'ext':'External period-finder periodogram', 'tls':'Transit Least-Squares periodogram'} METHODSHORTLABELS = {'gls':'Generalized L-S', 'pdm':'Stellingwerf PDM', 'aov':'Schwarzenberg-Czerny AoV', 'mav':'Schwarzenberg-Czerny AoVMH', 'acf':'McQuillan+ ACF', 'bls':'BLS', 'win':'Sampling L-S', 'ext':'External period-finder', 'tls':'TLS'}
[docs]def plot_periodbase_lsp(lspinfo, outfile=None, plotdpi=100): '''Makes a plot of periodograms obtained from `periodbase` functions. This takes the output dict produced by any `astrobase.periodbase` period-finder function or a pickle filename containing such a dict and makes a periodogram plot. Parameters ---------- lspinfo : dict or str If lspinfo is a dict, it must be a dict produced by an `astrobase.periodbase` period-finder function or a dict from your own period-finder function or routine that is of the form below with at least these keys:: {'periods': np.array of all periods searched by the period-finder, 'lspvals': np.array of periodogram power value for each period, 'bestperiod': a float value that is the period with the highest peak in the periodogram, i.e. the most-likely actual period, 'method': a three-letter code naming the period-finder used; must be one of the keys in the `METHODLABELS` dict above, 'nbestperiods': a list of the periods corresponding to periodogram peaks (`nbestlspvals` below) to annotate on the periodogram plot so they can be called out visually, 'nbestlspvals': a list of the power values associated with periodogram peaks to annotate on the periodogram plot so they can be called out visually; should be the same length as `nbestperiods` above} If lspinfo is a str, then it must be a path to a pickle file that contains a dict of the form described above. outfile : str or None If this is a str, will write the periodogram plot to the file specified by this string. If this is None, will write to a file called 'lsp-plot.png' in the current working directory. plotdpi : int Sets the resolution in DPI of the output periodogram plot PNG file. Returns ------- str Absolute path to the periodogram plot file created. ''' # get the lspinfo from a pickle file transparently if isinstance(lspinfo,str) and os.path.exists(lspinfo): LOGINFO('loading LSP info from pickle %s' % lspinfo) with open(lspinfo,'rb') as infd: lspinfo = pickle.load(infd) try: # get the things to plot out of the data periods = lspinfo['periods'] lspvals = lspinfo['lspvals'] bestperiod = lspinfo['bestperiod'] lspmethod = lspinfo['method'] # make the LSP plot on the first subplot plt.plot(periods, lspvals) plt.xscale('log',base=10) plt.xlabel('Period [days]') plt.ylabel(PLOTYLABELS[lspmethod]) plottitle = '%s best period: %.6f d' % (METHODSHORTLABELS[lspmethod], bestperiod) plt.title(plottitle) # show the best five peaks on the plot for bestperiod, bestpeak in zip(lspinfo['nbestperiods'], lspinfo['nbestlspvals']): plt.annotate('%.6f' % bestperiod, xy=(bestperiod, bestpeak), xycoords='data', xytext=(0.0,25.0), textcoords='offset points', arrowprops=dict(arrowstyle="->"), fontsize='x-small') # make a grid plt.grid(color='#a9a9a9', alpha=0.9, zorder=0, linewidth=1.0, linestyle=':') # make the figure if outfile and isinstance(outfile, str): if outfile.endswith('.png'): plt.savefig(outfile,bbox_inches='tight',dpi=plotdpi) else: plt.savefig(outfile,bbox_inches='tight') plt.close() return os.path.abspath(outfile) elif dispok: plt.show() plt.close() return else: LOGWARNING('no output file specified and no $DISPLAY set, ' 'saving to lsp-plot.png in current directory') outfile = 'lsp-plot.png' plt.savefig(outfile,bbox_inches='tight',dpi=plotdpi) plt.close() return os.path.abspath(outfile) except Exception: LOGEXCEPTION('could not plot this LSP, appears to be empty') return