Viewing Kepler Light Curves via Matplotlib and PyFITS

With all the attention that the Curiosity rover is getting, I thought I'd take a look back at the Kepler mission.  The objective of Kepler is, essentially, to discover planets outside our own Solar System.  The Kepler instrument is a space-based, wide-view telescope/photometer (see here for details).

It observes over 100,000 stars in Cygnus and Lyra, recording their brightness over time.  If a star has a planet and its orbit is aligned in the line of sight, the planet will pass across the face of the star, partially blocking its light.  The Kepler instrument is sensitive enough to detect these changes for planets of sufficient size.

The public has access to the data generated by the project.  The "meat" of the data are the light curves.  These are available via a web interface here, but I wanted to see if I could get the same curves from the raw data myself.

The raw data is available in FITS-formatted files.  FITS is a standard format used for image data but can contain non-image data and metadata as well.  For Kepler, the light curves are stored as electron fluxes read from the sensor over each observation.  Higher electron fluxes are brighter observations, and vice versa.

It turns out that there's at least one module that makes dealing with FITS files a lot easier: PyFITS.  Using PyFITS, I was able to put together a simple program to show a light curve stored in a given file.  Once you have this much, it can be a starting point to analyze the data from that source and from other sources.  My code is below: 

import matplotlib.pyplot as plt
import pyfits

filename=raw_input("\nWhat is the FITS file name? ")
FITSfile=pyfits.open(filename)

dataheader=FITSfile[1].header
topheader=FITSfile[0].header
jdadj=str(dataheader["BJDREFI"])
obsobject=str(dataheader["OBJECT"])

lightdata=FITSfile[1].data

flux=lightdata.field("PDCSAP_FLUX")     #Start with the PDCSAP_FLUX if
                                        #you're not sure which one to use.
                                        #This is the corrected electron
                                        #flux (electrons/sec) in the
                                        #optimal aperture after correction.

time=lightdata.field("TIME")            #Barycenter corrected Julian date

fig1=plt.figure()
sub1=fig1.add_subplot(111)
sub1.plot(time,flux,color="black",marker=",",linestyle="None")
xlab="Time (days, Julian date - %s)"%jdadj
sub1.set_xlabel(xlab)
sub1.set_ylabel("Brightness (electron flux)")
plottitle="Light Curve for %s"%obsobject
sub1.set_title(plottitle)

plt.show()

FITSfile.close()


Here is some sample output:

 

This is a light curve for KIC 8191672, a.k.a. Kepler-5.  The planet Kepler-5b orbits this star with a period of about 3.55 days and is roughly twice as massive as Jupiter.  The spots where the reading drops significantly from the average are the times when the planet is transiting across the face of the star, blocking a bit of the star's light.  See here and here for more details on Kepler-5b. 

The raw data by sources can be found here, listed by Kepler ID.