eleanor/eleanor_tess_workshop.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from lightkurve.lightcurve import LightCurve as LC\n",
    "from bls import BLS\n",
    "from matplotlib import gridspec\n",
    "import eleanor\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord, Angle\n",
    "from IPython.display import YouTubeVideo\n",
    "\n",
    "print(eleanor.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Goal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "id = 'evHF_mnIdj4'\n",
    "YouTubeVideo(id=id,width=900,height=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to create light curves for stars (or other objects!) found within the TESS Full-Frame Images (FFIs). Grabbing pixels from the FFI and doing a bit of photometry shouldn't be a problem, right? The big issue is that the FFIs are big (~35Mb) and there are a lot of them for each sector (~21,000) and it's safe to assume most people don't have access to that kind of computing space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "data = np.load('ffi_2d.npy')\n",
    "\n",
    "plt.figure(figsize=(8,8))\n",
    "ax = plt.gca()\n",
    "im = ax.imshow(data, origin='lower', vmin=0, vmax=600)\n",
    "plt.title('Full-Frame Image', fontsize=18)\n",
    "\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax).set_label(label='Pixel Intensity', fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's start making light curves!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![SegmentLocal](eleanor.gif \"ellie\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pass in a TIC ID $\\rightarrow$ Get a Light Curve!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a simple first example, we will use $\\texttt{eleanor}$ to create a target pixel file (TPF) and a light curve for a given TESS Input Catalog (TIC) target.\n",
    "\n",
    "We will load in the data for TIC 234503282. Calling $\\texttt{eleanor.Source()}$ will assign a $\\textit{Gaia}$ ID, coordinates, TESS magnitude, sector, camera, CCD, and associated postcard."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "star = eleanor.Source(tic=234503282, sector=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is just some of the information we receive when we call $\\texttt{eleanor.Source()}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false
   },
   "outputs": [],
   "source": [
    "print('TIC ID: ',star.tic)\n",
    "print('Coordinates: ', star.coords)\n",
    "print('Camera, Chip Location: ', star.camera, ',', star.chip)\n",
    "print('TESS mag: ', star.tess_mag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned previously, a postcard is assigned when calling $\\texttt{eleanor.Source()}$, but what exactly is a postcard? A postcard is a 148 $\\times$ 104 pixel cutout of the FFIs that are time-stacked and background subtracted. For the postcard background, we subtract a constant from the frame using $\\texttt{photutils.MMMBackground()}$. The background at each cadence is passed through with the postcard object. This is the postcard that is assigned to our source TIC234503282.\n",
    "\n",
    "Now that we have our Source information, we simply call the $\\texttt{eleanor.TargetData()}$ function, which will extracted a TPF, perform aperture photometry with an automatically optimized aperture choice, and complete some simple systematics corrections.\n",
    "\n",
    "For first time users, or users looking at this part of the TESS FOV for the first time, a file will download. Don't be alarmed! We promise it's not a virus. The downloading file is the postcard that we just discussed.\n",
    "\n",
    "If you call do_psf=True (which enables you to complete a simple point-spred function model), a loading bar will appear to let you know which cadence the model is on and when. If you do not call do_psf, no loading bar will appear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false
   },
   "outputs": [],
   "source": [
    "data = eleanor.TargetData(star, height=15, width=15, do_psf=True, do_pca=True)\n",
    "\n",
    "plt.figure(figsize=(14,6))\n",
    "plt.imshow(data.post_obj.flux[:,:,0], origin='lower', vmax=200)\n",
    "plt.title('Postcard', fontsize=18)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "star.postcard"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yay! We have our light curve object! First thing's first: let's take a peak at the pixels themselves and the \"best\" aperture automatically selected by $\\texttt{eleanor}$. If all goes as planned, your target should be in or very close to the center of the TPF frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(14,6))\n",
    "ax1.imshow(data.tpf[0], origin='lower', vmax=200)\n",
    "ax1.set_title('TPF', fontsize=20)\n",
    "ax2.imshow(data.aperture, origin='lower')\n",
    "ax2.set_title('Best Aperture', fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At first glance, this aperture may be a little big... But we'll deal with that in a little bit. First, let's see what the background looks like. There are two steps to the background subtraction: on the postcard level and on the TPF level. We can check how much background was subtracted by calling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.bkg_type"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"PC_LEVEL\" means Postcard Level, so no additional background was subtracted at the TPF level. Let's take a look at what the background actually looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,6))\n",
    "plt.plot(data.time, data.flux_bkg, linewidth=4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is pretty close to what we expect for the background. What's causing this structure? Earth-shine! In addition to background, we also create Quality Flags for the user. These will alert the user of spacecraft issues such as thruster fires, where the data at that time might not be good/useful. We can plot the Quality Flags:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,6))\n",
    "plt.plot(data.time, data.quality, '.', linewidth=4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any cadence where the Quality Flag != 0 is potentially not a good one. I think it's time to plot our first light curve! We're only going to plot the cadences with good quality flags. We provide FOUR different light curve options:\n",
    "    1. Raw Flux -- Just the aperture photometry.\n",
    "    2. Corrected Flux -- Minimal correction pointing variations, background, and removal of long term trends.\n",
    "    3. PCA Flux -- Principal Component Flux is used to help remove any common systematics seen in light curves on \n",
    "                   the same chip.\n",
    "    4. PSF Flux -- the Point-Spread Function modeled flux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We provide quality flags! You should use them\n",
    "q = data.quality == 0\n",
    "\n",
    "plt.figure(figsize=(14,7))\n",
    "plt.plot(data.time[q], data.raw_flux[q]/np.median(data.raw_flux[q])-0.015, 'k')\n",
    "plt.plot(data.time[q], data.corr_flux[q]/np.median(data.corr_flux[q]) + 0.01, 'r')\n",
    "plt.plot(data.time[q], data.pca_flux[q]/np.median(data.pca_flux[q]) + 0.035, 'y')\n",
    "plt.plot(data.time[q], data.psf_flux[q]/np.median(data.psf_flux[q]) + 0.06, 'b')\n",
    "\n",
    "plt.ylabel('Normalized Flux')\n",
    "plt.xlabel('Time')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yay light curves! Yay a signal! What could it be?? First, we should try creating a light curve with a smaller aperture, to make sure there isn't a background source causing this signal. TESS's pixels are large, so this is a good test to do. LUCKILY! Every $\\texttt{eleanor.TargetData()}$ object comes equipt with all of the apertures we tested and all of the raw and corrected fluxes for these additional apertures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(14,4))\n",
    "gs = gridspec.GridSpec(1, 2, width_ratios=[1, 3], height_ratios=[2]) \n",
    "ax0 = plt.subplot(gs[0])\n",
    "ax0.imshow(data.all_apertures[13], origin='lower')\n",
    "ax1 = plt.subplot(gs[1])\n",
    "ax1.plot(data.time[q], data.all_raw_lc[13][q])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the 21 apertures $\\texttt{eleanor.TargetData()}$ provides you with isn't enough, or you want to try a different shape or weighting scheme, you can create your own 2D array of the same shape as the TPF and pass that in as an aperture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = np.zeros(np.shape(data.tpf[0]))\n",
    "mask[5,7] = 1\n",
    "mask[6:8,6:9] = 1\n",
    "mask[8,7] = 1\n",
    "\n",
    "data.get_lightcurve(aperture=mask)\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(14,4))\n",
    "gs = gridspec.GridSpec(1, 2, width_ratios=[1, 3], height_ratios=[2]) \n",
    "ax0 = plt.subplot(gs[0])\n",
    "ax0.imshow(mask, origin='lower')\n",
    "ax1 = plt.subplot(gs[1])\n",
    "ax1.plot(data.time[q], data.raw_flux[q])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our potentially new planet candidate transit signal is still there! This is a promising sign. \n",
    "\n",
    "Before we move on to convincing ourselves that this transit signal may be real, what else can we do within $\\texttt{eleanor.TargetData()}$?\n",
    "\n",
    "We can run PCA with the flux from this new mask and choose the number of components we want to use. The default number of modes within $\\texttt{eleanor}$ is 8, but we have 16 components created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.pca(flux=data.corr_flux, modes=10)\n",
    "plt.figure(figsize=(14,5))\n",
    "plt.plot(data.time[q], data.pca_flux[q])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pass in Coordinates $\\rightarrow$ Get Light Curve for a Target Observed in Multiple Sectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! Now we know how to create a light curve using $\\texttt{eleanor.Source()}$ and $\\texttt{eleanor.TargetData()}$ with a TIC ID. Now, let's do it all over again, but this time we'll highlight two additional features:\n",
    "    1. Pass in coordinates to get a light curve.\n",
    "    2. Create a light curve for a source observed in multiple sectors (in this case, we only have Sectors 1 \n",
    "       and 2 ready for use).\n",
    "       \n",
    "First, we have to assign our coordinates with a unit. This allows users to pass in coordinates in degrees or hourangle. Then, instead of calling $\\texttt{eleanor.Source()}$ as we did last time, we'll call $\\texttt{eleanor.multi_sectors()}$, which additionally requires a list of sectors you want. If we call $\\texttt{eleanor.Source()}$, it will return the light curve for the most recent sector the target was observed in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ra  = Angle(11.595497, u.deg)\n",
    "dec = Angle(-63.473076,u.deg)\n",
    "\n",
    "stars = eleanor.multi_sectors(coords=SkyCoord(ra,dec), sectors=[1,2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\texttt{eleanor.multi_sectors()}$ returns a list of $\\texttt{eleanor.Source()}$ objects! We can then pass these into $\\texttt{eleanor.TargetData()}$ individually to create light curves for each sector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data0 = eleanor.TargetData(stars[0], height=15, width=15, bkg_size=31, do_psf=True, do_pca=True)\n",
    "data1 = eleanor.TargetData(stars[1], height=15, width=15, bkg_size=31, do_psf=True, do_pca=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's plot the light curve for this source from both sectors..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q0 = data0.quality == 0\n",
    "q1 = data1.quality == 0\n",
    "\n",
    "plt.figure(figsize=(14,6))\n",
    "plt.plot(data0.time[q0], data0.psf_flux[q0]/np.median(data0.psf_flux[q0]) + 0.05, 'k.')\n",
    "plt.plot(data1.time[q1], data1.psf_flux[q1]/np.median(data1.psf_flux[q1]) + 0.04, 'r.')\n",
    "\n",
    "plt.ylabel('Normalized Flux')\n",
    "plt.xlabel('Time')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oh my! Another new potential planet candidate signal! Just kidding -- This is the same source as we used in the first part :) We just used the coordinates rather than the TIC ID to create the light curve. I leave the rest of the steps we looked at before as an exercise to the user."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.save(output_fn='test.fits', directory='.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's fit it!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've all been super curious as to what this signal could be. Let's do a bit of initial vetting! The first trick to try is a Box-Least Squares analysis using $\\texttt{bls.py}$, an open-source packge created by Dan Foreman-Mackey. Planet transits, unlike eclipsing binaries, are quite boxy in shape. We can fit a box to a phase-folded light curve and see if the data matches well. The rest of this notebook will go over how to call $\\texttt{bls.py}$ and its outputs.\n",
    "\n",
    "We'll be calling $\\texttt{bls.py}$ a few times so, as is good practice, let's make this a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bls_fitting(t, lc):  \n",
    "    # Make a lightkurve LightCurve object. This will allow us to easily flatten our light curve.\n",
    "    lk = LC(time=t, flux=lc)\n",
    "    lc_filt = lk.flatten(window_length=51).flux\n",
    "\n",
    "\n",
    "    # Create BLS Object:\n",
    "    #     BLS will fit an array of transit duration times.\n",
    "    #     You also can pass in a minimum and maximum period, given in days.\n",
    "    durations = np.linspace(0.05, 0.2, 10)\n",
    "    model     = BLS(t, lc_filt)\n",
    "    results   = model.autopower(durations, frequency_factor=5.0,\n",
    "                                maximum_period=15)\n",
    "    \n",
    "    # Find the period and duration from the maximum power peak\n",
    "    period = results.period[np.argmax(results.power)]\n",
    "    dur = results.duration[np.argmax(results.power)]\n",
    "    \n",
    "    # Find the transit time\n",
    "    t0 = results.transit_time[np.argmax(results.power)]\n",
    "    \n",
    "    # Fold the light curve over the transit time\n",
    "    t_fold = (t - t0 + 0.5*period) % period - 0.5*period\n",
    "\n",
    "    \n",
    "    # EVERYTHING BEYOND HERE IN THE FUNCTION IS FOR PLOTTING PURPOSES\n",
    "    fig, (ax1, ax2, ax3) = plt.subplots(ncols=1,nrows=3, figsize=(16,11))\n",
    "\n",
    "    ax2.axvline(period, color='goldenrod', alpha=0.4, lw=5)\n",
    "    for n in range(2, 10):\n",
    "        ax2.axvline(n*period, alpha=0.4, lw=1, linestyle=\"dashed\")\n",
    "        ax2.axvline(period / n, alpha=0.4, lw=1, linestyle=\"dashed\")\n",
    "\n",
    "    ax2.plot(results.period, results.power, \"k\", lw=1)\n",
    "\n",
    "    ax2.set_xlim(np.min(results.period), np.max(results.period))\n",
    "    ax2.set_xlabel(\"period [days]\")\n",
    "    ax2.set_ylabel(\"log likelihood\")\n",
    "\n",
    "    ax1.plot(t, lc_filt, '.', color=\"midnightblue\")\n",
    "\n",
    "    ax3.plot(t_fold, lc_filt, 'k.')\n",
    "    ax3.set_xlim(-4*dur, 4*dur)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's pass in a time and a light curve into our $\\texttt{bls_fitting}$ function and see what plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bls_fitting(data0.time[q0], data0.raw_flux[q0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first subplot shows our light curve and the deep potential transits we know and love. The second subplot shows the periodogram, where the x-axis is the period in days and the y-axis is the log likelihood. Highlighted in yellow is the peak of the function and the dashed blue lines represent harmonics of the transit. We see that all of the lines line up nicely with peaks in the function. The third subplot is our phase-folded light curve. It looks really good!\n",
    "\n",
    "We can do the same analysis with the light curves from both sectors as well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now you try!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TIC 139771134\n",
    "\n",
    "Gaia 6582115549524788480 (Hint: You do not need to look up the associated TIC ID or coordinates for this source.)\n",
    "\n",
    "RA, Dec = 04:35:54.594, -66.133666\n",
    "\n",
    "TIC 260304296\n",
    "\n",
    "RA, Dec = 113.593357, -64.936218"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Guiding questions:\n",
    "\n",
    "    1. Which have been observed in both Sectors 1 & 2?\n",
    "    2. Which have potential new planet candidate transits?\n",
    "    3. Which option of flux gives the clearest signal of a transit?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Currently, only Sectors 1 and 2 are available with $\\texttt{eleanor}$ and we're working to get all available sectors working soon!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}