How to recover the first TESS planet candidate with Lightkurve?#
Data from the TESS mission are available from the data archive at MAST. This tutorial demonstrates how the Lightkurve Python package can be used to read in these data and create your own TESS light curves with different aperture masks.
Below is a quick tutorial on how to get started using Lightkurve and TESS data. We’ll use the nearby, bright target Pi Mensae (ID 261136679), around which the mission team recently discovered a short period planet candidate on a 6.27 day orbit. See the pre-print paper by Huang et al (2018) for more details.
TESS data is stored in a binary file format which is documented in the TESS Science Data Products Description Document. Lightkurve provides a TessTargetPixelFile class which allows you to interact with the data easily.
[1]:
import lightkurve as lk
[2]:
search_result = lk.search_targetpixelfile('Pi Mensae', mission='TESS', sector=1)
[3]:
search_result
[3]:
[4]:
tpf = search_result.download(quality_bitmask='default')
/Users/tapritc2/tessgi/lightkurve_party/lightkurve/src/lightkurve/search.py:424: LightkurveWarning: Warning: 2 files available to download. Only the first file has been downloaded. Please use `download_all()` or specify additional criteria (e.g. quarter, campaign, or sector) to limit your search.
warnings.warn(
[5]:
tpf
[5]:
TessTargetPixelFile(TICID: 261136679)
TessTargetPixelFile’s have many helpful methods and attributes. For example, you can access basic meta data on the target easily:
[6]:
tpf.mission
[6]:
'TESS'
[7]:
tpf.targetid # TESS Input Catalog (TIC) Identifier
[7]:
261136679
[8]:
tpf.sector # TESS Observation Sector
[8]:
1
[9]:
tpf.camera # TESS Camera Number
[9]:
4
[10]:
tpf.ccd # TESS CCD Number
[10]:
2
We might want to plot the data, we can do this with the plot() method. You can add the keyword aperture_mask to plot an aperture on top of the image. In this case we’ve used the pipeline_mask which is stored in the original .fits file, but you can use any aperture you like.
[11]:
%matplotlib inline
tpf.plot(aperture_mask=tpf.pipeline_mask);
If you want to access the original fits file that generated the data you can use the hdu attribute of the tpf. This will return an astropy.io.fits object, for example
[12]:
tpf.hdu
[12]:
[<astropy.io.fits.hdu.image.PrimaryHDU object at 0x1774b2b90>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x1774f51d0>, <astropy.io.fits.hdu.image.ImageHDU object at 0x1774b4350>, <astropy.io.fits.hdu.table.BinTableHDU object at 0x1776c7510>]
You can access each extension and the data inside it in the same way you’d use astropy.io.fits. If you want to access data held in the TPF, such as the time of the observations, you can do that easily by using
[13]:
tpf.time
[13]:
<Time object: scale='tdb' format='btjd' value=[1325.2969605 1325.29834936 1325.29973823 ... 1353.17428819 1353.17567704
1353.1770659 ]>
This returns the time in units of days counted since Julian Day 2457000.
You can access the corresponding flux values using
[14]:
tpf.flux
[14]:
Flux is a numpy.ndarray with a shape of (TIME x PIXELS x PIXELS). If you want to access just the first frame you can use
[15]:
tpf.flux[0]
[15]:
These values are in units electrons per second.
Building Light Curves from TPFs#
We can use the to_lightcurve() method to turn this TPF into a light curve using Simple Aperture Photometry. This will put an aperture on the target, and sum up the flux in all the pixels inside the aperture.
The default for to_lightcurve() is to use the mask generated by the TESS pipeline.
[16]:
lc = tpf.to_lightcurve()
Now we can use the plot function to take a look at the data.
[17]:
lc.errorbar();
This looks pretty good, but maybe we can improve things by creating a new aperture.
[18]:
aperture_mask = tpf.create_threshold_mask(threshold=10)
# Plot that aperture
tpf.plot(aperture_mask=aperture_mask);
[19]:
lc = tpf.to_lightcurve(aperture_mask=aperture_mask)
[20]:
lc.errorbar();
There’s a long term trend in this dataset, which we can remove with a simple smoothing filter. You can use the lc.flatten() method to apply and divide the Savitzky-Golay smoothing filter. Here we’ll use a window_length of 1001 cadences, which is roughly a 5% of the full length of the light curve.
[21]:
# Number of cadences in the full light curve
print(lc.time.shape)
(18279,)
[22]:
flat_lc = lc.flatten(window_length=1001)
flat_lc.errorbar();
The light curve looks much flatter. Unfortunately there is a portion of the light curve that is very noisy, due to a jitter in the TESS spacecraft. We can remove this simply by masking the light curve. First we’ll select the times that had the jitter.
[23]:
# Flag the times that are good quality
mask = (flat_lc.time.value < 1346) | (flat_lc.time.value > 1350)
Then we can just clip those times out.
[24]:
masked_lc = flat_lc[mask]
masked_lc.errorbar();
We can use Lightkurve to plot these two light curves over each other to see the difference.
[25]:
# First define the `matplotlib.pyplot.axes`
ax = flat_lc.errorbar()
# Pass that axis to the next plot
masked_lc.errorbar(ax=ax, label='masked');
This looks much better. Now we might want to clip out some outliers from the light curve. We can do that with a simple lightkurve function remove_outliers().
[26]:
clipped_lc = masked_lc.remove_outliers(sigma=6)
clipped_lc.errorbar();
It’s a little hard to see these data because of the plotting style. Let’s use a scatter plot instead. We can do this with the lc.scatter() method. This method works in the same way that matplotlib.pyplot.scatter works, and takes in the same keyword arguments.
[27]:
clipped_lc.scatter(s=0.1);
We can also add errorbars using the lc.errorbar() method.
[28]:
ax = clipped_lc.scatter(s=0.1)
clipped_lc.errorbar(ax=ax, alpha=0.2); # alpha determines the transparency
Finally let’s use lightkurve to fold the data at the exoplanet orbital period and see if we can see the transit.
[29]:
folded_lc = clipped_lc.fold(period=6.27, epoch_time=1325.504)
folded_lc.errorbar();
It looks like there’s something there, but it’s hard to see. Let’s bin the light curve to reduce the number of points, but also reduce the uncertainty of those points.
[30]:
import astropy.units as u
binned_lc = folded_lc.bin(time_bin_size=5*u.minute)
binned_lc.errorbar();
And now we can see the transit of Pi Men c!
Note that you can actually do all these steps in just a few lines:
[31]:
lc = tpf.to_lightcurve(aperture_mask=aperture_mask).flatten(window_length=1001)
lc = lc[(lc.time.value < 1346) | (lc.time.value > 1350)]
lc.remove_outliers(sigma=6).fold(period=6.27, epoch_time=1325.504).bin(time_bin_size=5*u.minute).errorbar();
Comparing two apertures#
In the above tutorial we used our own aperture instead of the pipeline aperture. Let’s compare the results from using these two different apertures.
[32]:
# Use the default
lc = tpf.to_lightcurve(aperture_mask=tpf.pipeline_mask).flatten(window_length=1001)
lc = lc[(lc.time.value < 1346) | (lc.time.value > 1350)].remove_outliers(6).fold(period=6.27, epoch_time=1325.504).bin(5*u.minute)
# Use a custom aperture
custom_lc = tpf.to_lightcurve(aperture_mask=aperture_mask).flatten(window_length=1001)
custom_lc = custom_lc[(custom_lc.time.value < 1346) | (custom_lc.time.value > 1350)].remove_outliers(6).fold(period=6.27, epoch_time=1325.504).bin(5*u.minute)
[33]:
ax = lc.errorbar(label='Default aperture')
custom_lc.errorbar(ax=ax, label='Custom aperture');
The importance of using different aperture masks is clearly visible in the figure above. Note however that the data archive at MAST also contains lightcurve products which have more advanced systematics removal methods applied. We will explore those in a future tutorial!