Intermediate: Finding Flares and Variable Stars in TASOC Light Curves#

This notebook demostrates MAST’s programmatic tools for accessing TESS time series data while exploring a flaring star from the literature (Günther et al 2019) and a nearby variable star.

The following topics will be covered:

  • Using the MAST API to get mission pipeline and TASOC light curves

  • Plotting TESS light curves in Python

  • Using the MAST API to make an full frame image (FFI) cutout

  • Creating a movie of TPF frames in Python

  • Using the MAST API to get a list of TESS Input Catalog (TIC) sources

  • Over plotting TIC sources on TESS images

  • Exploring the period of a variable star

Table of contents#

  1. Exploring a known stellar flare

    1. Querying MAST

    2. Downloading the light curves

    3. Plotting the light curves

    4. Making an animation

  2. Exploring a variable star

    1. Overlaying TESS Input Catalog sources

    2. Getting the variable star light curve

    3. Plotting the variable star light curve

    4. Finding the period of the variable star

  3. Additional Resources

  4. About this Notebook


  • TESS: The Transiting Exoplanet Survey Satellite

  • TASOC: The TESS Asteroseismic Science Operations Center

  • Sector: TESS observed the sky in regions of 24x96 degrees along the southern, then northern, ecliptic hemispheres. Each of these regions is referred to as a “sector”, starting with Sector 1.

  • TIC: The TESS input catalog.

  • FFI: TESS periodically reads out the entire frame of all four cameras, nominally every 30 minutes, and stores them as full frame images (FFIs).

  • TPF: Target Pixel File, a fits file containing stacks of small images centered on a target, one image for every timestamp the telescope took data.

  • HDU: Header Data Unit. A FITS file is made up of HDUs that contain data and metadata relating to the file. The first HDU is called the primary HDU, and anything that follows is considered an “extension”, e.g., “the first FITS extension”, “the second FITS extension”, etc.

  • HDUList: A list of HDUs that comprise a fits file.

  • BJD: Barycentric Julian Date, the Julian Date that has been corrected for differences in the Earth’s position with respect to the Solar System center of mass.

  • BTJD: Barycentric TESS Julian Date, the timestamp measured in BJD, but offset by 2457000.0. I.e., BTJD = BJD - 2457000.0

  • WCS: World Coordinate System, the coordinates that locate an astronomical object on the sky.


In this notbook we will use the MAST module of Astroquery (astroquery.mast) to query and download data, and various astropy and numpy functions to manipulate the data. We will use both the matplotlib and bokeh plotting packages to visualize our data as they have different strengths and weaknesses.

# For querying for data
from astroquery.mast import Tesscut, Observations, Catalogs

# For manipulating data
import numpy as np

from astropy.table import Table
from astropy.coordinates import SkyCoord
from import fits
from astropy.wcs import WCS
from astropy.timeseries import LombScargle
from astropy.time import Time
import astropy.units as u

# For matplotlib plotting
import matplotlib
%matplotlib inline

import matplotlib.pyplot as plt
import matplotlib.animation as animation

# For animation display
from matplotlib import rc
from IPython.display import HTML
rc('animation', html='jshtml')

# For bokeh plotting
from bokeh import plotting
from bokeh.models import Span
Loading BokehJS ...

Exploring a stellar flare#

Selecting the flare#

We will start with a known flare from the literature, in this case from Günther, M. N., Zhan, Z., Seager, S., et al. 2019, arXiv e-prints, arXiv:1901.00443. We picked a particularly long flare to give us the best chance of finding it in the 30 minute cadence data as well as the 2 minute cadence data.

We’ve made note of the TIC ID and sector for our flare of interest, as well as its peak time in BJD format:

tic_id = 141914082
sector = 1

tpeak = 2458341.89227 # Julian Day

Querying MAST#

Mission light curves#

Using the query_criteria function in the astroquery.mast.Observations class, we will specify that we are looking for TESS mission data using the obs_collection argument, that we want a specific TIC ID using the target_name argument, and that we want observations from a particular sector using the sequence_number argument.

mission_res = Observations.query_criteria(obs_collection="TESS", 
Table masked=True length=1
scienceTESSSPOCPhotometerTESSTESSOptical141914082--1tess2018206045859-s0001-0000000141914082-0120-s94.6175354047675-72.0448462219924G011175_G011180_G011176Ricker, Georgetimeseries358324.7932369675958352.67632608796120.0600.01000.0--58458.5833333--PUBLICFalse60829534110344410

TASOC light curves#

In addition to mission pipeline data, MAST also hosts a variety of community contributed High Level Science Products (HLSPs), all of which are given the mission (obs_collection) “HLSP”. In this case we will specifically search for HLSPs in the TESS project, which will return the light curves provided by the TASOC (note the provenance_name of “TASOC”). All other arguments remain the same.

tasoc_res = Observations.query_criteria(target_name=tic_id, 
Table length=6

This query returns two light curves. To understand the difference between the two light curves we look at the t_exptime column, and note the different values. These exposure times correspond to 2 minutes (short cadence) and 30 minutes (long cadence). We will explore both light curves.

Downloading the light curves#

For the rest of this notebook we will work with the TASOC light curves only, although we could do the same analysis on the mission light curves.

Querying for the list of associated data products#

Each observation may have one or more associated data products. In the case of the TASOC light curves, there is simply a single light curve file for each observation.

tasoc_prod = Observations.get_product_list(tasoc_res)
tasoc_prod["dataproduct_type", "description", "dataURI", "size"]
Table length=9

Downloading files#

We can choose to download some or all of the associated data files, in this case since we just have the two light curves, we will download all of the products.

tasoc_manifest = Observations.download_products(tasoc_prod)
Downloading URL to ./mastDownload/HLSP/hlsp_gsfc-eleanor-lite_tess_ffi_s0001-0000000141914082_tess_v1.0_lc/hlsp_gsfc-eleanor-lite_tess_ffi_s0001-0000000141914082_tess_v1.0_lc.fits ... [Done]
Downloading URL to ./mastDownload/HLSP/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc.fits ... [Done]
Downloading URL to ./mastDownload/HLSP/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc/hlsp_qlp_tess_ffi_s0001-0000000141914082_tess_v01_llc.txt ... [Done]
Downloading URL to ./mastDownload/HLSP/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_cbv-lc.fits ... [Done]
Downloading URL to ./mastDownload/HLSP/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05/hlsp_tasoc_tess_ffi_tic00141914082-s0001-cam4-ccd2-c1800_tess_v05_ens-lc.fits ... [Done]
Downloading URL to ./mastDownload/HLSP/hlsp_tasoc_tess_tpf_tic00141914082-s0001-cam4-ccd2-c0120_tess_v05/hlsp_tasoc_tess_tpf_tic00141914082-s0001-cam4-ccd2-c0120_tess_v05_cbv-lc.fits ... [Done]
Downloading URL to ./mastDownload/HLSP/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_tp/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_lc.fits ... [Done]
Downloading URL to ./mastDownload/HLSP/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_tp/hlsp_tess-spoc_tess_phot_0000000141914082-s0001_tess_v1_tp.fits ... [Done]
Downloading URL to ./mastDownload/HLSP/hlsp_tglc_tess_ffi_gaiaid-5266270443442455040-s0001-cam4-ccd2_tess_v1_llc/hlsp_tglc_tess_ffi_gaiaid-5266270443442455040-s0001-cam4-ccd2_tess_v1_llc.fits ... [Done]
Table length=9
Local PathStatusMessageURL

Plotting the light curves#

We will use the bokeh plotting library so that we can have interactivity, and will plot both the 2 minute and 30 minute cadence data together.

We can tell which file corresponds to which cadence length by examining the filenames and noting that one contains c0120 (2 minutes) and the other c1800 (30 minutes).

# Loading the short cadence light curve
hdu =["Local Path"][0])
short_cad_lc = Table(hdu[1].data)

# Loading the long cadence light curve
hdu =["Local Path"][1])
long_cad_lc = Table(hdu[1].data)
Table length=1281
bfig = plotting.figure(width=850, height=300, title=f"Detrended Lightcurve (TIC{tic_id})")

# Short cadence["TIME"], short_cad_lc["RAW_FLUX"]/np.median(short_cad_lc["RAW_FLUX"]), fill_color="black",size=2, line_color=None)
bfig.line(short_cad_lc["TIME"], short_cad_lc["RAW_FLUX"]/np.median(short_cad_lc["RAW_FLUX"]), line_color='black')

# Long cadence["TIME"],long_cad_lc["SAP_FLUX"], fill_color="#0384f7",size=6, line_color=None)
bfig.line(long_cad_lc["TIME"],long_cad_lc["SAP_FLUX"], line_color='#0384f7')

# Marking the flare (tpeak is in BJD, while the time column in the light curve is BTJD, so we must convert)
vline = Span(location=(tpeak - 2457000), dimension='height', line_color='#bf006e', line_width=3)

# Labeling the axes
bfig.xaxis.axis_label = "Time (BTJD)"
bfig.yaxis.axis_label = "Flux"

Experiment with the controls on the right to zoom in and and pan around on this light curve to look at the marked flare and other features (some sort of periodic feature maybe?)

Making an animation#

Looking at the above plot we can see the flare event in both the long and short cadence light curves. Since we can see it even in the 30 minute data, we should be able to make an animation of the area around the flaring star and see the flare happen.

We will use TESScut, the MAST cutout tool for full-frame images to cutout the area around the flaring star across the entire sector, and then make a movie that shows how it changes over time.

We will use the astroquery.mast Tesscut class to make this cutout.
We will use two functions:

  • Find the sky coordinate of our flare star: Observations._resolve_object*

  • Query for cutouts and get the result as a list of HDUList objects: Tesscut.get_cutouts **

* Observations._resolve_object is a private (not documented) function which is being removed in favor of the public function Observations.resolve_object in the next Astroquery release.

** We must start by finding the sky coordinate of our star, however starting with the next Astroquery release, Tesscut functions will be able to take an object name such as a TIC ID as well.

coord = Observations.resolve_object(f"TIC {tic_id}")

**Requesting a cutout target pixel file. **

This query will return a list of HDUList objects, each of which is the cutout target pixel file (TPF) for a single sector. In this case, because we specified a single sector we know that the resulting list will only have one element and can pull it out directly.

cutout_hdu = Tesscut.get_cutouts(coordinates=coord, size=40, sector=1)[0]
Filename: <class '_io.BytesIO'>
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU      57   ()      
  1  PIXELS        1 BinTableHDU    281   1282R x 12C   [D, E, J, 1600J, 1600E, 1600E, 1600E, 1600E, J, E, E, 38A]   
  2  APERTURE      1 ImageHDU        82   (40, 40)   int32   

A TESScut TPF has three extensions:

  • No. 0 (Primary): This HDU contains meta-data related to the entire file.

  • No. 1 (Pixels): This HDU contains a binary table that holds data like cutout image arrays and times. We will extract information from here to get our set of cutout images.

  • No. 2 (Aperture): This HDU contains the image extension with data collected from the aperture. We will use this extension to get the WCS associated with out cutout.

cutout_table = Table(cutout_hdu[1].data)

The Pixels extension contains the cutout data table, which has a number of columns. For our puposes we care about the “TIME” column which has the observation times in BTJD, and the “FLUX” column which contains the cutout images (they are calibrated, but not background subtracted).

Exploring the cutout time series#

We want to explore what is happening within our cutout area over the time that the flare occurs, so we will make an animated plot of the cutout frames.

We can’t make a movie of the whole sector (it would take too long), so we will choose only the time range around the flare.

Use the light curve plot to figure out what time range you want to animate, or use our selections.

start_btjd = 1341.5
end_btjd = 1342.5

start_index = (np.abs(cutout_table['TIME'] - start_btjd)).argmin()
end_index = (np.abs(cutout_table['TIME'] - end_btjd)).argmin()

print(f"Frames {start_index}-{end_index} ({end_index-start_index} frames)")
Frames 721-769 (48 frames)

Looking at the animated cutout#

def make_animation(data_array, start_frame=0, end_frame=None, vmin=None, vmax=None, delay=50):
    Function that takes an array where each frame is a 2D image array and make an animated plot
    that runs through the frames.
    Note: This can take a long time to run if you have a lot of frames.    
    data_array : array
        Array of 2D images.
    start_frame : int
        The index of the initial frame to show. Default is the first frame.
    end_frame : int
        The index of the final frame to show. Default is the last frame.
    vmin : float
        Data range min for the colormap. Defaults to data minimum value.
    vmax : float
        Data range max for the colormap. Defaults to data maximum value.
        Delay before the next frame is shown in milliseconds.

    response : `animation.FuncAnimation`
    if not vmin:
        vmin = np.min(data_array)
    if not vmax:
        vmax = np.max(data_array)
    if not end_frame:
        end_frame = len(data_array) - 1 # set to the end of the array
    num_frames = end_frame - start_frame + 1 # include the end frame
    def animate(i, fig, ax, binarytab, start=0):
        """Function used to update the animation"""
        ax.set_title("Epoch #" + str(i+start))
        im = ax.imshow(binarytab[i+start],, vmin=vmin, vmax=vmax)
        return im,
    # Create initial plot.
    fig, ax = plt.subplots(figsize=(10,10))
    ax.imshow(data_array[start_frame],, vmin=vmin, vmax=vmax)

    ani = animation.FuncAnimation(fig, animate, fargs=(fig, ax, data_array, start_frame), frames=num_frames, 
                                  interval=delay, repeat_delay=1000)
    return ani
make_animation(cutout_table['FLUX'], start_index, end_index, vmax=700, delay=150)