Intro to PSFMachine Tutorial 1.0¶

In [1]:

import numpy as np
import psfmachine as pm
import lightkurve as lk
import matplotlib.pyplot as plt

In [2]:

hidecode

import warnings
from scipy import sparse
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=sparse.SparseEfficiencyWarning)

Target Pixel Files¶

In this tutorial we will see how PSFMachine can be used to do PSF photmetry and extract light curves from Kepler/K2 data.

TPFs are dta products created by Kepler, K2 and TESS missions. It contain custom cutouts of the full telescope CCD around selected targets. Each file contains about 5x5 pixels around the target and all the available cadences. But many of these TPFs also contain background sources that can be easily isolated or blended. Kepler and TESS pipelines only created light curves for main targets but there are many sources (~270,000 background sources in Kepler TPFs down to 18th g mag) that can be extracted to create new light curves. PSFMachine provides the tools to create these new light curves.

We can download these files using the LightKurve package by specifying the target name, mission, quarter of observation, etc. Let's download 200 TPFs around Kepler-16 observed during quarter 5

In [3]:

tpfs = lk.search_targetpixelfile('Kepler-16', mission='Kepler', quarter=5, 
                                 radius=1000, limit=200, cadence='long'
                                ).download_all(quality_bitmask=None)

TPFMachine¶

As explained in the LFD paper, PSFMachine performs PSF photometry of all sources (using Gaia as source catalog) in the TPFs. To create a PSFMachine object from TPFs, we use the submodule TPFMachine that includes useful methods to parse and work with TPF data. We'll do this using pm.TPFMachine.from_TPFs() which inputs a collection of TPFs.

Printing the new machine object will tell us how many sources are in the TPFs pixels, the number of cadences, and total number of pixel data.

In [4]:

machine = pm.TPFMachine.from_TPFs(tpfs)
machine

Out[4]:

TPFMachine (N sources, N times, N pixels): (359, 4486, 8057)

To see the source catalog, we can use the attribute machine.sources

In [5]:

machine.sources

Out[5]:

	designation	ra	dec	parallax	parallax_error	pmra	pmdec	phot_g_mean_flux	phot_g_mean_mag	phot_bp_mean_mag	phot_rp_mean_mag	phot_bp_mean_flux	phot_rp_mean_flux	clean_flag	tpf_id
0	Gaia DR3 2133262530247167616	289.324552	51.582686	1.616466	0.032562	-2.298416	-11.012768	7062.031645	16.065042	16.800795	15.232625	2600.753065	6398.939581	0	12554848
1	Gaia DR3 2133263080002987392	289.207900	51.586092	0.504336	0.050699	-3.558866	-7.815627	2738.845449	17.093449	17.545729	16.499073	1309.559992	1993.098815	0	None
2	Gaia DR3 2133263080003584128	289.205602	51.582282	0.984947	0.017882	-7.581524	-6.793829	20097.434517	14.929515	15.299093	14.393923	10370.022892	13854.449111	0	12554673
3	Gaia DR3 2133263084302086272	289.194951	51.583371	0.883457	0.019191	-5.343587	-7.699694	16423.362618	15.148711	15.525050	14.608281	8421.642868	11372.234414	0	12554659
4	Gaia DR3 2133263084302086784	289.196992	51.586797	0.309922	0.017201	-4.139091	-5.795632	20483.462439	14.908858	15.399604	14.256479	9453.125928	15724.123963	0	None
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
354	Gaia DR3 2133492779152324736	288.838971	51.893487	0.821284	0.040056	-6.652471	-2.473718	3941.997142	16.698076	17.246286	15.988624	1725.451307	3189.398445	0	None
355	Gaia DR3 2133492916591300864	288.889496	51.894211	0.144796	0.069053	-1.658150	-3.073686	1581.248367	17.689867	18.231911	17.057388	696.069188	1191.798859	0	None
356	Gaia DR3 2133492916591304320	288.894867	51.896601	0.648541	0.050593	-5.933056	-0.122113	2980.270501	17.001728	17.525383	16.370127	1334.331618	2244.435653	0	None
357	Gaia DR3 2133492912292382976	288.893595	51.890141	2.585587	0.015439	-18.046519	-29.085882	530062.957277	11.376548	11.654636	10.933944	297551.702144	335413.728047	0	12689094
358	Gaia DR3 2133493088390017792	288.909527	51.930597	0.841148	0.021900	-3.869092	5.411347	11906.104779	15.497943	15.880922	14.951333	6068.033664	8291.345107	0	12734396

359 rows × 15 columns

Create Light Curves¶

To create light curves, we just have to do machine.fit_lightcurves(). This method will fit a PSF model to all the sources available simultaneously, use the pos_corr to remove velocity aberration from the scene, and will use these two model to obtain the PSF flux. Also, the method will do aperture photometry (SAP) on all sources where each aperture is defined accordingly to the PSF profile. Setting plot=True will show diagnostic plots for the PSF model and the time model.

For more details on how the PSF shape mode and time model fitting work, see the LFD paper and the more detailed tutorial for PSMachine Tutorial_11_TPFs.

machine.fit_lightcurves() will create the new attribute machine.lcs with the new light curves as lk.LightCurveCollection object.

In [6]:

machine.fit_lightcurves(plot=True)

WARNING: Input data contains invalid values (NaNs or infs), which were automatically clipped. [astropy.stats.sigma_clipping]
2023-06-09 15:13:21,845 - astroquery - WARNING - Input data contains invalid values (NaNs or infs), which were automatically clipped.
Optimizing apertures per source: 100%|███████████████████████████████████████████████| 359/359 [00:15<00:00, 22.54it/s]
SAP: 100%|█████████████████████████████████████████████████████████████████████████| 359/359 [00:00<00:00, 1052.80it/s]
Fitting 359 Sources (w. VA): 100%|█████████████████████████████████████████████████| 4486/4486 [01:12<00:00, 61.87it/s]
Fitting 359 Sources (w. VA): 100%|█████████████████████████████████████████████████| 4486/4486 [01:05<00:00, 68.42it/s]

Let's see our new light curves:

In [7]:

machine.lcs[0]

Out[7]:

KeplerLightCurve length=4486 LABEL="KIC 12554848" QUARTER=5 CAMPAIGN=None

time	flux	flux_err	psf_flux_NVA	psf_flux_err_NVA	sap_flux	sap_flux_err
	electron / s	electron / s	electron / s	electron / s	electron / s	electron / s
Time	float64	float64	float64	float64	float64	float64
2455276.4906453583	3183.8373152452436	2.6301012934246994	3186.0180004965164	2.638509325265348	3809.272216796875	3.4859983921051025
...	...	...	...	...	...	...
2455371.1620332897	3068.1642802787674	2.660412461705618	3011.3345861487705	2.5855987035461587	3669.705322265625	3.437415361404419

The lightkurve object has 7 columns:

time: the observing time in MJD
flux and lux_err: have the PSF photometry corrected for velocity aberration
flux_NVA and lux_err_NVA: have the PSF photometry without correction for velocity aberration
sap_flux and sap_lux_err: have the aperture photometry

Because each element in machine.lcs is a lightkurve.LightKurve object, we can use all the functionalities the library offers, like plotting, saving to FITS, doing periodograms, etc.

Plot Light Curves¶

Kepler-16 (KIC 12644769) is a binary system with a K-type MS star and a M-type red dwarf. Let's see how its light curve looks like:

In [8]:

machine.sources[machine.sources.tpf_id == 12644769]

Out[8]:

	designation	ra	dec	parallax	parallax_error	pmra	pmdec	phot_g_mean_flux	phot_g_mean_mag	phot_bp_mean_mag	phot_rp_mean_mag	phot_bp_mean_flux	phot_rp_mean_flux	clean_flag	tpf_id
270	Gaia DR3 2133476355197071616	289.075798	51.757299	13.280588	0.028397	14.09245	-48.939436	379684.909773	11.738809	12.534245	10.867677	132348.418185	356522.977737	0	12644769

In [11]:

machine.lcs[270].plot()

Out[11]:

<AxesSubplot: xlabel='Time [JD]', ylabel='Flux [$\\mathrm{e^{-}\\,s^{-1}}$]'>

What about other light curve in the data:

In [12]:

for lc in machine.lcs[::25]:
    if np.isnan(lc.flux).all():
        continue
    ax = lc.plot()
    plt.show()

In [ ]: