What’s new in PyKE 3?

Developed since 2012, PyKE offers a user-friendly way to inspect and analyze the pixels and lightcurves obtained by NASA’s Kepler and K2.

The latest version of PyKE, v3.1, was released in January 2018 and adds a new object-oriented Python API which is intended to aid the development of custom pipelines and tools by the community.

In [47]:
%load_ext autoreload
%autoreload 2
%matplotlib inline
from matplotlib import rcParams
rcParams["figure.figsize"] = (14, 5)
The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Introducing a generic LightCurve class

The most notable change is the introduction of a generic LightCurve class which provides operations that are intended to suit time series data from any astronomical survey. A light curve is simply instantiated as follows:

In [45]:
from pyke import LightCurve
lc = LightCurve(time=[1, 2, 3], flux=[78.4, 79.6, 76.5])

A LightCurve object provides easy access to a range of common operations, such as fold(), flatten(), remove_outliers(), cdpp(), plot(), and more. To demonstrate these operations, let’s create a LightCurve object from a KeplerLightCurveFile we obtain from the data archive at MAST:

In [46]:
from pyke import KeplerLightCurveFile
lcfile = KeplerLightCurveFile("https://archive.stsci.edu/missions/kepler/lightcurves/0119/011904151/kplr011904151-2010009091648_llc.fits")
lc = lcfile.SAP_FLUX

Now lc is a LightCurve object on which you can run operations. For example, we can plot it:

In [43]:
lc.plot()
Out[43]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f4aadfd9b00>
../../_images/tutorials_ipython_notebooks_whatsnew31_9_1.png

We can access several of the metadata properties:

In [94]:
lc.keplerid
Out[94]:
11904151
In [95]:
lc.channel
Out[95]:
8
In [96]:
lc.quarter
Out[96]:
4

We can access the time and flux as arrays:

In [97]:
lc.time[:10]
Out[97]:
array([ 352.37710511,  352.39753829,  352.43840464,  352.45883781,
        352.47927099,  352.49970417,  352.52013735,  352.54057042,
        352.5610037 ,  352.58143688])
In [98]:
lc.flux[:10]
Out[98]:
array([ 562954.5625,  562910.3125,  562950.375 ,  562939.1875,
        562964.5625,  562967.    ,  562952.1875,  562950.375 ,
        562979.4375,  562973.4375], dtype=float32)

We don’t particularly care about the long-term trends, so let’s use a Savitzky-Golay filter to flatten the lightcurve:

In [86]:
detrended_lc, _ = lc.flatten(polyorder=1)
detrended_lc.plot()
Out[86]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f4aac1e1d30>
../../_images/tutorials_ipython_notebooks_whatsnew31_18_1.png
In [92]:
folded_lc = detrended_lc.fold(period=0.837495, phase=0.92)
folded_lc.plot();
../../_images/tutorials_ipython_notebooks_whatsnew31_19_0.png

We can also compute the CDPP noise metric:

In [88]:
lc.cdpp()
Out[88]:
34.399382449996828

Target Pixel File (TPF)

PyKE 3.1 includes class called KeplerTargetPixelFile which is used to handle target pixel files:

In [10]:
from pyke import KeplerTargetPixelFile

A KeplerTargetPixelFile can be instantiated either from a local file or a url:

In [11]:
tpf = KeplerTargetPixelFile('https://archive.stsci.edu/missions/k2/target_pixel_files/c14/'
                            '200100000/82000/ktwo200182949-c14_lpd-targ.fits.gz')

Additionally, we can mask out cadences that are flagged using the quality_bitmask argument in the constructor:

In [12]:
tpf = KeplerTargetPixelFile('https://archive.stsci.edu/missions/k2/target_pixel_files/c14/'
                            '200100000/82000/ktwo200182949-c14_lpd-targ.fits.gz',
                            quality_bitmask=KeplerQualityFlags.CONSERVATIVE_BITMASK)

Furthermore, we can mask out pixel values using the aperture_mask argument. The default behaviour is to use all pixels that have real values. This argument can also get a string value 'kepler-pipeline', in which case the default aperture used by Kepler’s pipeline is applied.

In [13]:
tpf = KeplerTargetPixelFile('https://archive.stsci.edu/missions/k2/target_pixel_files/c14/'
                            '200100000/82000/ktwo200182949-c14_lpd-targ.fits.gz',
                            aperture_mask='kepler-pipeline',
                            quality_bitmask=KeplerQualityFlags.CONSERVATIVE_BITMASK)
In [14]:
tpf.aperture_mask
Out[14]:
array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]], dtype=bool)

The TPF objects stores both data and a few metadata information, e.g., channel number, EPIC number, reference column and row, module, and shape. The whole header is also available:

In [15]:
tpf.header(ext=0)
Out[15]:
SIMPLE  =                    T / conforms to FITS standards
BITPIX  =                    8 / array data type
NAXIS   =                    0 / number of array dimensions
EXTEND  =                    T / file contains extensions
NEXTEND =                    2 / number of standard extensions
EXTNAME = 'PRIMARY '           / name of extension
EXTVER  =                    1 / extension version number (not format version)
ORIGIN  = 'NASA/Ames'          / institution responsible for creating this file
DATE    = '2017-10-30'         / file creation date.
CREATOR = '472097 TargetPixelExporterPipelineModule' / pipeline job and program
PROCVER = 'svn+ssh://murzim/repo/soc/tags/release/9.3.75 r62933' / SW version
FILEVER = '6.1     '           / file format version
TIMVERSN= 'OGIP/93-003'        / OGIP memo number for file format
TELESCOP= 'Kepler  '           / telescope
INSTRUME= 'Kepler Photometer'  / detector type
OBJECT  = 'EPIC 200182949'     / string version of target id
KEPLERID=            200182949 / unique Kepler target identifier
CHANNEL =                    1 / CCD channel
MODULE  =                    2 / CCD module
OUTPUT  =                    1 / CCD output
CAMPAIGN=                   14 / Observing campaign number
DATA_REL=                   20 / data release version number
OBSMODE = 'long cadence'       / observing mode
MISSION = 'K2      '           / Mission name
TTABLEID=                   92 / target table id
RADESYS = 'ICRS    '           / reference frame of celestial coordinates
RA_OBJ  =           165.706242 / [deg] right ascension
DEC_OBJ =             2.769444 / [deg] declination
EQUINOX =               2000.0 / equinox of celestial coordinate system
PMRA    =                      / [arcsec/yr] RA proper motion
PMDEC   =                      / [arcsec/yr] Dec proper motion
PMTOTAL =                      / [arcsec/yr] total proper motion
PARALLAX=                      / [arcsec] parallax
GLON    =                      / [deg] galactic longitude
GLAT    =                      / [deg] galactic latitude
GMAG    =                      / [mag] SDSS g band magnitude
RMAG    =                      / [mag] SDSS r band magnitude
IMAG    =                      / [mag] SDSS i band magnitude
ZMAG    =                      / [mag] SDSS z band magnitude
JMAG    =                      / [mag] J band magnitude from 2MASS
HMAG    =                      / [mag] H band magnitude from 2MASS
KMAG    =                      / [mag] K band magnitude from 2MASS
KEPMAG  =                      / [mag] Kepler magnitude (Kp)
GRCOLOR =                      / [mag] (g-r) color, SDSS bands
JKCOLOR =                      / [mag] (J-K) color, 2MASS bands
GKCOLOR =                      / [mag] (g-K) color, SDSS g - 2MASS K
TEFF    =                      / [K] Effective temperature
LOGG    =                      / [cm/s2] log10 surface gravity
FEH     =                      / [log10([Fe/H])]  metallicity
EBMINUSV=                      / [mag] E(B-V) reddening
AV      =                      / [mag] A_v extinction
RADIUS  =                      / [solar radii] stellar radius
TMINDEX =                      / unique 2MASS catalog ID
CHECKSUM= 'QkmmQijlQijlQijl'   / HDU checksum updated 2017-10-30T20:42:50Z

The pixel fluxes time series can be accessed using the flux property:

In [16]:
tpf.flux.shape
Out[16]:
(3209, 35, 35)

This shows that our TPF is a 35 x 35 image recorded over 3209 cadences.

One can visualize the pixel data at a given cadence using the plot method:

In [17]:
tpf.plot(frame=1)
../../_images/tutorials_ipython_notebooks_whatsnew31_38_0.png

We can perform aperture photometry using the method to_lightcurve:

In [18]:
lc = tpf.to_lightcurve()
In [19]:
plt.figure(figsize=[17, 4])
plt.plot(lc.time, lc.flux)
Out[19]:
[<matplotlib.lines.Line2D at 0x129c079e8>]
../../_images/tutorials_ipython_notebooks_whatsnew31_41_1.png

Let’s see how the previous light curve compares against the 'SAP_FLUX' produced by Kepler’s pipeline. For that, we are going to explore the KeplerLightCurveFile class:

In [20]:
from pyke.lightcurve import KeplerLightCurveFile
In [21]:
klc = KeplerLightCurveFile('https://archive.stsci.edu/missions/k2/lightcurves/'
                           'c14/200100000/82000/ktwo200182949-c14_llc.fits',
                           quality_bitmask=KeplerQualityFlags.CONSERVATIVE_BITMASK)
In [22]:
sap_lc = klc.SAP_FLUX
In [23]:
plt.figure(figsize=[17, 4])
plt.plot(lc.time, lc.flux)
plt.plot(sap_lc.time, sap_lc.flux)
plt.ylabel('Flux (e-/s)')
plt.xlabel('Time (BJD - 2454833)')
Out[23]:
<matplotlib.text.Text at 0x12a0af400>
../../_images/tutorials_ipython_notebooks_whatsnew31_46_1.png

Now, let’s correct this light curve using by fitting cotrending basis vectors. That can be achived either with the KeplerCBVCorrector class or the compute_cotrended_lightcurve in KeplerLightCurveFile. Let’s try the latter:

In [24]:
klc_corrected = klc.compute_cotrended_lightcurve(cbvs=range(1, 17))
In [25]:
plt.figure(figsize=[17, 4])
plt.plot(klc_corrected.time, klc_corrected.flux)
plt.ylabel('Flux (e-/s)')
plt.xlabel('Time (BJD - 2454833)')
Out[25]:
<matplotlib.text.Text at 0x12a1c1c50>
../../_images/tutorials_ipython_notebooks_whatsnew31_49_1.png
In [26]:
pdcsap_lc = klc.PDCSAP_FLUX
In [27]:
plt.figure(figsize=[17, 4])
plt.plot(klc_corrected.time, klc_corrected.flux)
plt.plot(pdcsap_lc.time, pdcsap_lc.flux)
plt.ylabel('Flux (e-/s)')
plt.xlabel('Time (BJD - 2454833)')
Out[27]:
<matplotlib.text.Text at 0x12a0bc0f0>
../../_images/tutorials_ipython_notebooks_whatsnew31_51_1.png

Utility functions

PyKE has included two convinience functions to convert between module.output to channel and vice-versa:

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
In [2]:
from pyke.utils import module_output_to_channel, channel_to_module_output
/Users/jvmirca/anaconda3/lib/python3.6/site-packages/matplotlib/__init__.py:1405: UserWarning:
This call to matplotlib.use() has no effect because the backend has already
been chosen; matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

  warnings.warn(_use_error_msg)
In [3]:
module_output_to_channel(module=19, output=3)
Out[3]:
67
In [4]:
channel_to_module_output(67)
Out[4]:
(19, 3)

PyKE 3.1 includes KeplerQualityFlags class which encodes the meaning of the Kepler QUALITY bitmask flags as documented in the Kepler Archive Manual (Table 2.3):

In [5]:
from pyke.utils import KeplerQualityFlags
In [6]:
KeplerQualityFlags.decode(1)
Out[6]:
['Attitude tweak']

It also can handle multiple flags:

In [7]:
KeplerQualityFlags.decode(1 + 1024 + 1048576)
Out[7]:
['Attitude tweak', 'Sudden sensitivity dropout', 'Thruster firing']

A few quality flags are already computed:

In [8]:
KeplerQualityFlags.decode(KeplerQualityFlags.DEFAULT_BITMASK)
Out[8]:
['Attitude tweak',
 'Safe mode',
 'Coarse point',
 'Earth point',
 'Desaturation event',
 'Cosmic ray in optimal aperture',
 'Manual exclude',
 'No data',
 'Thruster firing']
In [9]:
KeplerQualityFlags.decode(KeplerQualityFlags.CONSERVATIVE_BITMASK)
Out[9]:
['Attitude tweak',
 'Safe mode',
 'Coarse point',
 'Earth point',
 'Desaturation event',
 'Cosmic ray in optimal aperture',
 'Manual exclude',
 'Sudden sensitivity dropout',
 'Cosmic ray in collateral data',
 'No data',
 'Possible thruster firing',
 'Thruster firing']

Pixel Response Function (PRF) Photometry

PyKE 3.1 also includes tools to perform PRF Photometry:

In [28]:
from pyke.prf import PRFPhotometry, SceneModel, SimpleKeplerPRF

For that, let’s create a SceneModel which will be fitted to the object of the following TPF:

In [29]:
tpf = KeplerTargetPixelFile('https://archive.stsci.edu/missions/k2/target_pixel_files/c14/'
                            '201500000/43000/ktwo201543306-c14_lpd-targ.fits.gz',
                            quality_bitmask=KeplerQualityFlags.CONSERVATIVE_BITMASK)
In [30]:
tpf.plot(frame=100)
../../_images/tutorials_ipython_notebooks_whatsnew31_71_0.png
In [31]:
scene = SceneModel(prfs=[SimpleKeplerPRF(channel=tpf.channel, shape=tpf.shape[1:],
                                         column=tpf.column, row=tpf.row)])

We also need to define prior distributions on the parameters of our SceneModel model. Those parameters are the flux, center positions of the target, and a constant background level. We can do that with oktopus:

In [32]:
from oktopus.prior import UniformPrior
In [33]:
unif_prior = UniformPrior(lb=[0, 1090., 706., 0.],
                          ub=[1e5, 1096., 712., 1e5])
In [34]:
scene.plot(*unif_prior.mean)
../../_images/tutorials_ipython_notebooks_whatsnew31_76_0.png
In [35]:
prf_phot = PRFPhotometry(scene_model=scene, prior=unif_prior)
In [36]:
results = prf_phot.fit(tpf.flux + tpf.flux_bkg)
  0%|          | 0/3633 [00:00<?, ?it/s]/Users/jvmirca/anaconda3/lib/python3.6/site-packages/autograd/core.py:82: RuntimeWarning: invalid value encountered in log
  result_value = self.fun(*argvals, **kwargs)
/Users/jvmirca/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py:1927: RuntimeWarning: invalid value encountered in double_scalars
  tmp2 = (x - v) * (fx - fw)
/Users/jvmirca/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py:1928: RuntimeWarning: invalid value encountered in double_scalars
  p = (x - v) * tmp2 - (x - w) * tmp1
/Users/jvmirca/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py:1929: RuntimeWarning: invalid value encountered in double_scalars
  tmp2 = 2.0 * (tmp2 - tmp1)
100%|██████████| 3633/3633 [02:49<00:00, 21.46it/s]
In [37]:
plt.imshow(prf_phot.residuals[1], origin='lower')
plt.colorbar()
Out[37]:
<matplotlib.colorbar.Colorbar at 0x12aed6be0>
../../_images/tutorials_ipython_notebooks_whatsnew31_79_1.png
In [38]:
flux = results[:, 0]
xcenter = results[:, 1]
ycenter = results[:, 2]
bkg_density = results[:, 3]
In [39]:
plt.figure(figsize=[17, 4])
plt.plot(tpf.time, flux)
plt.ylabel('Flux (e-/s)')
plt.xlabel('Time (BJD - 2454833)')
Out[39]:
<matplotlib.text.Text at 0x12b2960f0>
../../_images/tutorials_ipython_notebooks_whatsnew31_81_1.png
In [40]:
plt.figure(figsize=[17, 4])
plt.plot(tpf.time, xcenter)
plt.ylabel('Column position')
plt.xlabel('Time (BJD - 2454833)')
Out[40]:
<matplotlib.text.Text at 0x128a809b0>
../../_images/tutorials_ipython_notebooks_whatsnew31_82_1.png
In [41]:
plt.figure(figsize=[17, 4])
plt.plot(tpf.time, ycenter)
plt.ylabel('Row position')
plt.xlabel('Time (BJD - 2454833)')
Out[41]:
<matplotlib.text.Text at 0x1299983c8>
../../_images/tutorials_ipython_notebooks_whatsnew31_83_1.png
In [42]:
plt.figure(figsize=[17, 4])
plt.plot(tpf.time, bkg_density)
plt.ylabel('Background density')
plt.xlabel('Time (BJD - 2454833)')
Out[42]:
<matplotlib.text.Text at 0x12b2bb278>
../../_images/tutorials_ipython_notebooks_whatsnew31_84_1.png

For more examples on PRF photometry, please see our tutorials page: http://pyke.keplerscience.org/tutorials/index.html