<jupyter_start><jupyter_text>Sector1This notebook downloads and pre-processes the data for Sector 1 targets.<jupyter_code>%matplotlib inline
from astropy.io import fits
from astroquery.mast import Observations
import matplotlib.pyplot as plt
import numpy as np
from tqdm import tqdm
# Global settings
sector = 1
ntargets = 1000
earth_in_view = "((time > 1325.5) & (time < 1334)) | ((time > 1340) & (time < 1349))"
no_earth = "(time > 1336) & (time < 1337)"
ngroups = 4
nreg_per_group = 25
thresh = [1.0, 1.0, 1.0]
# Get the data
obs_id = obs_id="tess*-s00%02d*" % sector
obsTable = Observations.query_criteria(obs_id=obs_id,
                                       dataproduct_type="timeseries",
                                       obs_collection="TESS",
                                       t_exptime=[119, 121])
obsTable.sort("obs_id")
np.random.seed(41)
inds = np.random.choice(np.arange(len(obsTable), dtype=int), size=ntargets, replace=False)
times = []
fluxes = []
flux_errors = []
for i, ind in tqdm(enumerate(inds), total=ntargets):
    fits_file = "https://mast.stsci.edu/api/v0.1/Download/file/?uri=" + obsTable[ind]["dataURL"]
    with fits.open(fits_file, mode="readonly") as hdulist:
        t = hdulist[1].data['TIME']
        f = hdulist[1].data['SAP_BKG']
        e = hdulist[1].data['SAP_BKG_ERR']
        q = hdulist[1].data['QUALITY']
        times.append(t[q==0])
        fluxes.append(f[q==0])
        flux_errors.append(e[q==0])
# Interpolate to a common time array
ind = np.argmax([len(t) for t in times])
time = np.linspace(times[ind][0], times[ind][-1], len(times[ind]))
flux = np.empty((len(times), len(time)))
flux_err = np.empty((len(times), len(time)))
fnorm = np.empty((len(times), len(time)))
rngs = np.empty(len(times))
for i in range(len(fluxes)):
    flux[i] = np.interp(time, times[i], fluxes[i])
    flux_err[i] = np.interp(time, times[i], flux_errors[i])
    fsort = (flux[i])[np.argsort(flux[i])]
    maxf = np.nanmax(fsort[:int(0.95 * len(fsort))])
    minf = np.nanmin(fsort[int(0.05 * len(fsort)):])
    rngs[i] = maxf - minf
    fnorm[i] = (flux[i] - minf) / (maxf - minf)
# Group them by features
group = np.zeros(len(flux), dtype=int)
score = np.zeros(len(flux), dtype=float) * np.nan
med = np.zeros((ngroups, len(time)), dtype=float) * np.nan
sig = np.zeros((ngroups, len(time)), dtype=float) * np.nan

# Loop over groups
for g in tqdm(range(ngroups - 1), total=ngroups - 1):
    # Compute the median absolute deviation of the group
    med[g] = np.nanmedian(fnorm[group==g], axis=0)
    mad = np.nanmedian(np.abs(fnorm[group==g] - med[g]), axis=0)
    sig[g] = 1.4826 * mad

    # Compute reduced chi^2 "score"
    for i in range(len(fluxes)):
        if group[i] == g:
            score[i] = np.nansum((fnorm[i] - med[g]) ** 2 / sig[g] ** 2) / len(time)

    # Remove the outlier light curves from this group
    group[score > thresh[g]] = g + 1
# Sort the light curves by (robust) min-max range
inds = np.argsort(rngs)[::-1]
flux = flux[inds]
flux_err = flux_err[inds]
fnorm = fnorm[inds]
group = group[inds]
score = score[inds]
# Plot the first several
nrows = 10
fig, ax = plt.subplots(nrows, ngroups, figsize=(16, 12))
lo = np.nanmin(med[0] - 2 * sig[0])
hi = np.nanmax(med[0] + 2 * sig[0])
for g in range(ngroups):
    for n in range(nrows):
        if n < len(fnorm[group==g]):
            ax[n, g].plot(time, fnorm[group==g][n])
            ax[n, g].plot(time, med[g], "r-", lw=1, alpha=0.25)
            ax[n, g].fill_between(time, med[g] - sig[g], med[g] + sig[g], color="r", alpha=0.25)
            if np.isfinite(score[group==g][n]):
                ax[n, g].annotate("%.2f" % score[group==g][n], xy=(-0.075, 0.5), xycoords="axes fraction")
        ax[n, g].axis('off')
        ax[n, g].set_ylim(lo, hi)
# Mask indices where Earth is not in view (~eyeballed)
inds = eval(earth_in_view)
bkg_inds = eval(no_earth)

# Keep only the first few in each group
times = time[inds]
bkg_times = time[bkg_inds]
fluxes = np.empty((nreg_per_group * (ngroups - 1), len(times)))
bkg_fluxes = np.empty((nreg_per_group * (ngroups - 1), len(bkg_times)))
bkg = np.empty(nreg_per_group * (ngroups - 1))
flux_errs = np.empty((nreg_per_group * (ngroups - 1), len(times)))
i = 0
for g in range(ngroups - 1):
    for n in range(nreg_per_group):
        fluxes[i] = flux[group==g][n][inds]
        norm = np.nanmedian(fluxes[i])
        fluxes[i] /= norm
        flux_errs[i] = flux_err[group==g][n][inds] / norm
        bkg_fluxes[i] = flux[group==g][n][bkg_inds]
        bkg_fluxes[i] /= norm
        bkg[i] = np.nanmedian(bkg_fluxes[i])
        i += 1
        
# Plot our data
fig, ax = plt.subplots(1, figsize=(16, 7))
for i, f in enumerate(fluxes):
    plt.plot(times, f, lw=0.5)
    plt.plot(bkg_times, bkg_fluxes[i], 'k-')
    plt.plot(1336.5, np.nanmedian(bkg_fluxes[i]), 'ro')
# Save the data
np.savez("../data/sector%02d-data.npz" % sector, time=times, flux=fluxes, flux_err=flux_errs, bkg=bkg)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17746478873239438, 'fraction_numerical': 0.026156941649899398, 'mean_word_length': 3.3952254641909816, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '305819', 'license': 'MIT', 'n_tokens_mistral': 2110, 'n_tokens_neox': 1980, 'n_words': 426, 'path': 'notebooks/Sector1.ipynb', 'repo_name': 'rodluger/earthshine', 'size': '4970'}