1624 lines
60 KiB
Python
1624 lines
60 KiB
Python
from astropy import units as u
|
|
from astropy.io import fits
|
|
from astropy.io import registry
|
|
from astropy.io.fits import HDUList
|
|
from astropy.table import Table
|
|
from astropy.timeseries import TimeSeries
|
|
from astropy.time import TimeBase, Time, TimeDelta
|
|
from astropy.units import Quantity, UnitsWarning
|
|
from scipy.signal import savgol_filter
|
|
from scipy.interpolate import interp1d
|
|
from copy import deepcopy
|
|
import numpy as np
|
|
import logging
|
|
import warnings
|
|
|
|
from astropy.time.formats import TimeFromEpoch
|
|
|
|
log = logging.getLogger(__name__)
|
|
|
|
class LightkurveWarning(Warning):
|
|
"""Class for all Lightkurve warnings."""
|
|
pass
|
|
|
|
class TimeBKJD(TimeFromEpoch):
|
|
name = 'bkjd'
|
|
unit = 1.0
|
|
epoch_val = 2454833
|
|
epoch_val2 = None
|
|
epoch_scale = 'tdb'
|
|
epoch_format = 'jd'
|
|
|
|
|
|
class TimeBTJD(TimeFromEpoch):
|
|
name = 'btjd'
|
|
unit = 1.0
|
|
epoch_val = 2457000
|
|
epoch_val2 = None
|
|
epoch_scale = 'tdb'
|
|
epoch_format = 'jd'
|
|
|
|
def detect_filetype(hdulist: HDUList) -> str:
|
|
# Is it a MIT/QLP TESS FFI Quicklook Pipeline light curve?
|
|
# cf. http://archive.stsci.edu/hlsp/qlp
|
|
if "mit/qlp" in hdulist[0].header.get("origin", "").lower():
|
|
return "QLP"
|
|
|
|
# Is it a vanilla eleanor or GSFC-ELEANOR-LITE light curve?
|
|
if (
|
|
hdulist[0].header.get("LITE") is not None
|
|
and hdulist[0].header.get("PCORIGIN") is not None
|
|
):
|
|
return "ELEANOR"
|
|
|
|
# Is it a PATHOS TESS light curve?
|
|
# cf. http://archive.stsci.edu/hlsp/pathos
|
|
# The 'pathos' name doesn't stay in the filename if when we use fits.open
|
|
# to download a file, so we have to check for all the important columns
|
|
# This will cause problems if another HLSP has the exact same colnames...
|
|
if all(
|
|
x in hdulist[1].columns.names
|
|
for x in [
|
|
"PSF_FLUX_RAW",
|
|
"PSF_FLUX_COR",
|
|
"AP4_FLUX_RAW",
|
|
"AP4_FLUX_COR",
|
|
"SKY_LOCAL",
|
|
]
|
|
):
|
|
return "PATHOS"
|
|
|
|
# Is it a TASOC TESS light curve?
|
|
# cf. https://tasoc.dk and https://archive.stsci.edu/hlsp/tasoc
|
|
if hdulist[0].header.get("ORIGIN") == "TASOC/Aarhus":
|
|
return "TASOC"
|
|
|
|
# Is it a CDIPS TESS light curve?
|
|
# cf. http://archive.stsci.edu/hlsp/cdips
|
|
if "cdips" in hdulist[0].header.get("ORIGIN", "").lower():
|
|
return "CDIPS"
|
|
|
|
# Is it a K2VARCAT file?
|
|
# There are no self-identifying keywords in the header, so go by filename.
|
|
fname = hdulist.filename() if callable(hdulist.filename) else hdulist.filename
|
|
if fname is not None:
|
|
if "hlsp_k2varcat" in fname:
|
|
return "K2VARCAT"
|
|
|
|
# Is it a K2SC file?
|
|
if "k2sc" in hdulist[0].header.get("creator", "").lower():
|
|
return "K2SC"
|
|
|
|
# Is it a K2SFF file?
|
|
try:
|
|
# There are no metadata keywords identifying K2SFF FITS files,
|
|
# so we go by structure.
|
|
if (
|
|
hdulist[1].header.get("EXTNAME") == "BESTAPER"
|
|
and hdulist[1].header.get("TTYPE4") == "ARCLENGTH"
|
|
):
|
|
return "K2SFF"
|
|
except Exception:
|
|
pass
|
|
|
|
# Is it an EVEREST file?
|
|
try:
|
|
if "EVEREST" in str(hdulist[0].header.get("COMMENT")):
|
|
return "EVEREST"
|
|
except Exception:
|
|
pass
|
|
|
|
# Is it a KEPSEISMIC file?
|
|
if hdulist[0].header.get("ORIGIN") == "CEA & SSI":
|
|
return "KEPSEISMIC"
|
|
|
|
# Is it a TGLC file?
|
|
if hdulist[0].header.get("ORIGIN") == "UCSB/TGLC":
|
|
return "TGLC"
|
|
|
|
# Is it an official data product?
|
|
header = hdulist[0].header
|
|
try:
|
|
# use `telescop` keyword to determine mission
|
|
# and `creator` to determine tpf or lc
|
|
if "TELESCOP" in header.keys():
|
|
telescop = header["telescop"].lower()
|
|
else:
|
|
# Some old custom TESS data did not define the `TELESCOP` card
|
|
telescop = header["mission"].lower()
|
|
creator = header["creator"].lower()
|
|
origin = header["origin"].lower()
|
|
if telescop == "kepler":
|
|
# Kepler TPFs will contain "TargetPixelExporterPipelineModule"
|
|
if "targetpixel" in creator:
|
|
return "KeplerTargetPixelFile"
|
|
# Kepler LCFs will contain "FluxExporter2PipelineModule"
|
|
elif (
|
|
"fluxexporter" in creator
|
|
or "lightcurve" in creator
|
|
or "lightcurve" in creator
|
|
):
|
|
return "KeplerLightCurve"
|
|
elif telescop == "tess":
|
|
# TESS TPFs will contain "TargetPixelExporterPipelineModule"
|
|
if "targetpixel" in creator:
|
|
return "TessTargetPixelFile"
|
|
# TESS LCFs will contain "LightCurveExporterPipelineModule"
|
|
elif "lightcurve" in creator:
|
|
return "TessLightCurve"
|
|
# Early versions of TESScut did not set a good CREATOR keyword
|
|
elif "stsci" in origin:
|
|
return "TessTargetPixelFile"
|
|
# If the TELESCOP or CREATOR keywords don't exist we expect a KeyError;
|
|
# if one of them is Undefined we expect `.lower()` to yield an AttributeError.
|
|
except (KeyError, AttributeError):
|
|
return None
|
|
|
|
def read(path_or_url, **kwargs):
|
|
# pass header into `detect_filetype()`
|
|
try:
|
|
with fits.open(path_or_url) as temp:
|
|
filetype = detect_filetype(temp)
|
|
except OSError as e:
|
|
filetype = None
|
|
# Raise an explicit FileNotFoundError if file not found
|
|
if "No such file" in str(e):
|
|
raise e
|
|
|
|
try:
|
|
if filetype == "KeplerLightCurve":
|
|
return KeplerLightCurve.read(path_or_url, format="kepler2", **kwargs)
|
|
elif filetype == "TessLightCurve":
|
|
return TessLightCurve.read(path_or_url, format="tess2", **kwargs)
|
|
elif filetype == "QLP":
|
|
return TessLightCurve.read(path_or_url, format="qlp2", **kwargs)
|
|
elif filetype == "ELEANOR":
|
|
return TessLightCurve.read(path_or_url, format="eleanor2", **kwargs)
|
|
elif filetype == "PATHOS":
|
|
return TessLightCurve.read(path_or_url, format="pathos2", **kwargs)
|
|
elif filetype == "CDIPS":
|
|
return TessLightCurve.read(path_or_url, format="cdips2", **kwargs)
|
|
elif filetype == "TASOC":
|
|
return TessLightCurve.read(path_or_url, format="tasoc2", **kwargs)
|
|
elif filetype == "K2SFF":
|
|
return KeplerLightCurve.read(path_or_url, format="k2sff2", **kwargs)
|
|
elif filetype == "EVEREST":
|
|
return KeplerLightCurve.read(path_or_url, format="everest2", **kwargs)
|
|
elif filetype == "KEPSEISMIC":
|
|
return KeplerLightCurve.read(path_or_url, format="kepseismic2", **kwargs)
|
|
elif filetype == "TGLC":
|
|
return TessLightCurve.read(path_or_url, format="tglc2", **kwargs)
|
|
except BaseException as exc:
|
|
# ensure path_or_url is in the error
|
|
print(exc)
|
|
print("End of read")
|
|
|
|
class QualityFlags(object):
|
|
"""Abstract class"""
|
|
|
|
STRINGS = {}
|
|
OPTIONS = {}
|
|
|
|
@classmethod
|
|
def decode(cls, quality):
|
|
# If passed an astropy quantity object, get the value
|
|
if isinstance(quality, Quantity):
|
|
quality = quality.value
|
|
result = []
|
|
for flag in cls.STRINGS.keys():
|
|
if quality & flag > 0:
|
|
result.append(cls.STRINGS[flag])
|
|
return result
|
|
|
|
@classmethod
|
|
def create_quality_mask(cls, quality_array, bitmask=None):
|
|
# Return an array filled with `True` by default (i.e. ignore nothing)
|
|
if bitmask is None:
|
|
return np.ones(len(quality_array), dtype=bool)
|
|
if isinstance(quality_array, u.Quantity):
|
|
quality_array = quality_array.value
|
|
# A few pre-defined bitmasks can be specified as strings
|
|
if isinstance(bitmask, str):
|
|
try:
|
|
bitmask = cls.OPTIONS[bitmask]
|
|
except KeyError:
|
|
valid_options = tuple(cls.OPTIONS.keys())
|
|
raise ValueError(
|
|
"quality_bitmask='{}' is not supported, "
|
|
"expected one of {}"
|
|
"".format(bitmask, valid_options)
|
|
)
|
|
# The bitmask is applied using the bitwise AND operator
|
|
quality_mask = (quality_array & bitmask) == 0
|
|
# Log the quality masking as info or warning
|
|
n_cadences = len(quality_array)
|
|
n_cadences_masked = (~quality_mask).sum()
|
|
percent_masked = 100.0 * n_cadences_masked / n_cadences
|
|
logmsg = (
|
|
"{:.0f}% ({}/{}) of the cadences will be ignored due to the "
|
|
"quality mask (quality_bitmask={})."
|
|
"".format(percent_masked, n_cadences_masked, n_cadences, bitmask)
|
|
)
|
|
if percent_masked > 20:
|
|
log.warning("Warning: " + logmsg)
|
|
else:
|
|
log.info(logmsg)
|
|
return quality_mask
|
|
|
|
|
|
class KeplerQualityFlags(QualityFlags):
|
|
"""
|
|
This class encodes the meaning of the various Kepler QUALITY bitmask flags,
|
|
as documented in the Kepler Archive Manual (Ref. [1], Table 2.3).
|
|
|
|
References
|
|
----------
|
|
.. [1] Kepler: A Search for Terrestrial Planets. Kepler Archive Manual.
|
|
http://archive.stsci.edu/kepler/manuals/archive_manual.pdf
|
|
"""
|
|
|
|
AttitudeTweak = 1
|
|
SafeMode = 2
|
|
CoarsePoint = 4
|
|
EarthPoint = 8
|
|
ZeroCrossing = 16
|
|
Desat = 32
|
|
Argabrightening = 64
|
|
ApertureCosmic = 128
|
|
ManualExclude = 256
|
|
# Bit 2**10 = 512 is unused by Kepler
|
|
SensitivityDropout = 1024
|
|
ImpulsiveOutlier = 2048
|
|
ArgabrighteningOnCCD = 4096
|
|
CollateralCosmic = 8192
|
|
DetectorAnomaly = 16384
|
|
NoFinePoint = 32768
|
|
NoData = 65536
|
|
RollingBandInAperture = 131072
|
|
RollingBandInMask = 262144
|
|
PossibleThrusterFiring = 524288
|
|
ThrusterFiring = 1048576
|
|
|
|
#: DEFAULT bitmask identifies all cadences which are definitely useless.
|
|
DEFAULT_BITMASK = (
|
|
AttitudeTweak
|
|
| SafeMode
|
|
| CoarsePoint
|
|
| EarthPoint
|
|
| Desat
|
|
| ManualExclude
|
|
| DetectorAnomaly
|
|
| NoData
|
|
| ThrusterFiring
|
|
)
|
|
#: HARD bitmask is conservative and may identify cadences which are useful.
|
|
HARD_BITMASK = (
|
|
DEFAULT_BITMASK
|
|
| SensitivityDropout
|
|
| ApertureCosmic
|
|
| CollateralCosmic
|
|
| PossibleThrusterFiring
|
|
)
|
|
#: HARDEST bitmask identifies cadences with any flag set. Its use is not recommended.
|
|
HARDEST_BITMASK = 2096639
|
|
|
|
#: Dictionary which provides friendly names for the various bitmasks.
|
|
OPTIONS = {
|
|
"none": 0,
|
|
"default": DEFAULT_BITMASK,
|
|
"hard": HARD_BITMASK,
|
|
"hardest": HARDEST_BITMASK,
|
|
}
|
|
|
|
#: Pretty string descriptions for each flag
|
|
STRINGS = {
|
|
1: "Attitude tweak",
|
|
2: "Safe mode",
|
|
4: "Coarse point",
|
|
8: "Earth point",
|
|
16: "Zero crossing",
|
|
32: "Desaturation event",
|
|
64: "Argabrightening",
|
|
128: "Cosmic ray in optimal aperture",
|
|
256: "Manual exclude",
|
|
1024: "Sudden sensitivity dropout",
|
|
2048: "Impulsive outlier",
|
|
4096: "Argabrightening on CCD",
|
|
8192: "Cosmic ray in collateral data",
|
|
16384: "Detector anomaly",
|
|
32768: "No fine point",
|
|
65536: "No data",
|
|
131072: "Rolling band in optimal aperture",
|
|
262144: "Rolling band in full mask",
|
|
524288: "Possible thruster firing",
|
|
1048576: "Thruster firing",
|
|
}
|
|
|
|
|
|
class TessQualityFlags(QualityFlags):
|
|
"""
|
|
This class encodes the meaning of the various TESS QUALITY bitmask flags,
|
|
as documented in the TESS Data Products Description Document (Ref. [1], Table 28).
|
|
|
|
References
|
|
----------
|
|
.. [1] TESS Science Data Products Description Document (EXP-TESS-ARC-ICD-0014)
|
|
https://archive.stsci.edu/missions/tess/doc/EXP-TESS-ARC-ICD-TM-0014.pdf
|
|
"""
|
|
|
|
AttitudeTweak = 1
|
|
SafeMode = 2
|
|
CoarsePoint = 4
|
|
EarthPoint = 8
|
|
Argabrightening = 16
|
|
Desat = 32
|
|
ApertureCosmic = 64
|
|
ManualExclude = 128
|
|
Discontinuity = 256
|
|
ImpulsiveOutlier = 512
|
|
CollateralCosmic = 1024
|
|
#: The first stray light flag is set manually by MIT based on visual inspection.
|
|
Straylight = 2048
|
|
#: The second stray light flag is set automatically by Ames/SPOC based on background level thresholds.
|
|
Straylight2 = 4096
|
|
# See TESS Science Data Products Description Document
|
|
PlanetSearchExclude = 8192
|
|
BadCalibrationExclude = 16384
|
|
# Set in the sector 20 data release notes
|
|
InsufficientTargets = 32768
|
|
|
|
#: DEFAULT bitmask identifies all cadences which are definitely useless.
|
|
DEFAULT_BITMASK = (
|
|
AttitudeTweak | SafeMode | CoarsePoint | EarthPoint | Desat | ManualExclude
|
|
)
|
|
#: HARD bitmask is conservative and may identify cadences which are useful.
|
|
HARD_BITMASK = (
|
|
DEFAULT_BITMASK | ApertureCosmic | CollateralCosmic | Straylight | Straylight2
|
|
)
|
|
#: HARDEST bitmask identifies cadences with any flag set. Its use is not recommended.
|
|
HARDEST_BITMASK = 65535
|
|
|
|
#: Dictionary which provides friendly names for the various bitmasks.
|
|
OPTIONS = {
|
|
"none": 0,
|
|
"default": DEFAULT_BITMASK,
|
|
"hard": HARD_BITMASK,
|
|
"hardest": HARDEST_BITMASK,
|
|
}
|
|
|
|
#: Pretty string descriptions for each flag
|
|
STRINGS = {
|
|
1: "Attitude tweak",
|
|
2: "Safe mode",
|
|
4: "Coarse point",
|
|
8: "Earth point",
|
|
16: "Argabrightening",
|
|
32: "Desaturation event",
|
|
64: "Cosmic ray in optimal aperture",
|
|
128: "Manual exclude",
|
|
256: "Discontinuity corrected",
|
|
512: "Impulsive outlier",
|
|
1024: "Cosmic ray in collateral data",
|
|
2048: "Straylight",
|
|
4096: "Straylight2",
|
|
8192: "Planet Search Exclude",
|
|
16384: "Bad Calibration Exclude",
|
|
32768: "Insufficient Targets for Error Correction Exclude",
|
|
}
|
|
|
|
def validate_method(method, supported_methods):
|
|
method = method.lower()
|
|
if method in supported_methods:
|
|
return method
|
|
raise ValueError(
|
|
"method '{}' is not supported; "
|
|
"must be one of {}".format(method, supported_methods)
|
|
)
|
|
|
|
def _is_dict_like(data1):
|
|
return hasattr(data1, "keys") and callable(getattr(data1, "keys"))
|
|
|
|
class LightCurve(TimeSeries):
|
|
|
|
# The constructor of the `TimeSeries` base class will enforce the presence
|
|
# of these columns:
|
|
_required_columns = ["time", "flux", "flux_err"]
|
|
|
|
# The following keywords were removed in Lightkurve v2.0.
|
|
# Their use will trigger a warning.
|
|
_deprecated_keywords = (
|
|
"targetid",
|
|
"label",
|
|
"time_format",
|
|
"time_scale",
|
|
"flux_unit",
|
|
)
|
|
_deprecated_column_keywords = [
|
|
"centroid_col",
|
|
"centroid_row",
|
|
"cadenceno",
|
|
"quality",
|
|
]
|
|
|
|
# If an iterable is passed for ``time``, we will initialize an AstroPy
|
|
# ``Time`` object using the following format and scale:
|
|
_default_time_format = "jd"
|
|
_default_time_scale = "tdb"
|
|
|
|
# To emulate pandas, we do not support creating new columns or meta data
|
|
# fields via attribute assignment, and raise a warning in __setattr__ when
|
|
# a new attribute is created. We need to relax this warning during the
|
|
# initial construction of the object using `_new_attributes_relax`.
|
|
_new_attributes_relax = True
|
|
|
|
# cf. issue #925
|
|
__array_priority__ = 100_000
|
|
|
|
def __init__(self, data=None, *args, time=None, flux=None, flux_err=None, **kwargs):
|
|
|
|
# the ` {has,get,set}_time_in_data()`: helpers to handle `data` of different types
|
|
# in some cases, they also need to access kwargs["names"] as well
|
|
|
|
def get_time_idx_in(names):
|
|
time_indices = np.argwhere(np.asarray(names) == "time")
|
|
if len(time_indices) > 0:
|
|
return time_indices[0][0]
|
|
else:
|
|
return None
|
|
|
|
def get_time_in_data_list():
|
|
if len(data) < 1:
|
|
return None
|
|
names = kwargs.get("names")
|
|
if names is None:
|
|
# the first item MUST be time if no names specified
|
|
if isinstance(data[0], TimeBase): # Time or TimeDelta
|
|
return data[0]
|
|
else:
|
|
return None
|
|
else:
|
|
time_idx = get_time_idx_in(names)
|
|
if time_idx is not None:
|
|
return data[time_idx]
|
|
else:
|
|
return None
|
|
|
|
def set_time_in_data_list(value):
|
|
if len(data) < 1:
|
|
raise AssertionError("data should be non-empty")
|
|
names = kwargs.get("names")
|
|
if names is None:
|
|
# the first item MUST be time if no names specified
|
|
# this is to support base Table's select columns
|
|
# in __getitem__()
|
|
# https://github.com/astropy/astropy/blob/326435449ad8d859f1abf36800c3fb88d49c27ea/astropy/table/table.py#L1888
|
|
data[0] = value
|
|
else:
|
|
time_idx = get_time_idx_in(names)
|
|
if time_idx is not None:
|
|
data[time_idx] = value
|
|
else:
|
|
raise AssertionError("data should have time column")
|
|
|
|
def get_time_in_data_np_structured_array():
|
|
if data.dtype.names is None: # no labeled filed, not a structured array
|
|
return None
|
|
if "time" not in data.dtype.names:
|
|
return None
|
|
return data["time"]
|
|
|
|
def remove_time_from_data_np_structured_array():
|
|
if data.dtype.names is None:
|
|
raise AssertionError("data should be a numpy structured array")
|
|
if "time" not in data.dtype.names:
|
|
raise AssertionError("data should have a time field")
|
|
filtered_names = [n for n in data.dtype.names if n != "time"]
|
|
return data[filtered_names]
|
|
|
|
def has_time_in_data():
|
|
"""Check if the data has a column with the name"""
|
|
if data is None:
|
|
return False
|
|
elif _is_dict_like(data):
|
|
# data is a dict-like object with keys
|
|
return "time" in data.keys()
|
|
elif _is_list_like(data):
|
|
# case data is a list-like object (a list of columns, etc.)
|
|
return get_time_in_data_list() is not None
|
|
elif _is_np_structured_array(data):
|
|
# case numpy structured array (supported by base TimeSeries)
|
|
# https://numpy.org/doc/stable/user/basics.rec.html
|
|
return get_time_in_data_np_structured_array() is not None
|
|
else:
|
|
raise ValueError(f"Unsupported type for time in data: {type(data)}")
|
|
|
|
def get_time_in_data():
|
|
if _is_dict_like(data):
|
|
# data is a dict-like object with keys
|
|
return data["time"]
|
|
elif _is_list_like(data):
|
|
return get_time_in_data_list()
|
|
elif _is_np_structured_array(data):
|
|
return get_time_in_data_np_structured_array()
|
|
else:
|
|
# should never reach here. It'd have been caught by `has_time_in()``
|
|
raise AssertionError("Unsupported type for time in data")
|
|
|
|
def set_time_in_data(value):
|
|
if _is_dict_like(data):
|
|
# data is a dict-like object with keys
|
|
data["time"] = value
|
|
elif _is_list_like(data):
|
|
set_time_in_data_list(value)
|
|
elif _is_np_structured_array(data):
|
|
# astropy Time cannot be assigned to a column in np structured array
|
|
# we have special codepath handling it outside this function
|
|
raise AssertionError("Setting Time instances to np structured array is not supported")
|
|
else:
|
|
# should never reach here. It'd have been caught by `has_time_in()``
|
|
raise AssertionError("Unsupported type for time in data")
|
|
|
|
# Delay checking for required columns until the end
|
|
self._required_columns_relax = True
|
|
|
|
# Lightkurve v1.x supported passing time, flux, and flux_err as
|
|
# positional arguments. We support it here for backwards compatibility.
|
|
if len(args) in [1, 2]:
|
|
warnings.warn(
|
|
"passing flux as a positional argument is deprecated"
|
|
", please use ``flux=...`` instead.",
|
|
LightkurveDeprecationWarning,
|
|
)
|
|
time = data
|
|
flux = args[0]
|
|
data = None
|
|
if len(args) == 2:
|
|
flux_err = args[1]
|
|
|
|
# For backwards compatibility with Lightkurve v1.x,
|
|
# we support passing deprecated keywords via **kwargs.
|
|
deprecated_kws = {}
|
|
for kw in self._deprecated_keywords:
|
|
if kw in kwargs:
|
|
deprecated_kws[kw] = kwargs.pop(kw)
|
|
|
|
deprecated_column_kws = {}
|
|
for kw in self._deprecated_column_keywords:
|
|
if kw in kwargs:
|
|
deprecated_column_kws[kw] = kwargs.pop(kw)
|
|
|
|
# If `time` is passed as keyword argument, we populate it with integer numbers
|
|
if data is None or not has_time_in_data():
|
|
if time is None and flux is not None:
|
|
time = np.arange(len(flux))
|
|
# We are tolerant of missing time format
|
|
if time is not None and not isinstance(time, (Time, TimeDelta)):
|
|
# Lightkurve v1.x supported specifying the time_format
|
|
# as a constructor kwarg
|
|
time = Time(
|
|
time,
|
|
format=deprecated_kws.get("time_format", self._default_time_format),
|
|
scale=deprecated_kws.get("time_scale", self._default_time_scale),
|
|
)
|
|
|
|
# Also be tolerant of missing time format if time is passed via `data`
|
|
if data is not None and has_time_in_data():
|
|
if not isinstance(get_time_in_data(), (Time, TimeDelta)):
|
|
tmp_time = Time(
|
|
get_time_in_data(),
|
|
format=deprecated_kws.get("time_format", self._default_time_format),
|
|
scale=deprecated_kws.get("time_scale", self._default_time_scale),
|
|
)
|
|
if _is_np_structured_array(data):
|
|
# special case for np structured array
|
|
# one cannot set a `Time` instance to it
|
|
# so we set the time to the `time` param, and take it out of data
|
|
time = tmp_time
|
|
data = remove_time_from_data_np_structured_array()
|
|
else:
|
|
set_time_in_data(tmp_time)
|
|
|
|
# Allow overriding the required columns
|
|
self._required_columns = kwargs.pop("_required_columns", self._required_columns)
|
|
|
|
# Call the SampledTimeSeries constructor.
|
|
# Disable required columns for now; we'll check those later.
|
|
tmp = self._required_columns
|
|
self._required_columns = []
|
|
super().__init__(data=data, time=time, **kwargs)
|
|
self._required_columns = tmp
|
|
|
|
# For some operations, an empty time series needs to be created, then
|
|
# columns added one by one. We should check that when columns are added
|
|
# manually, time is added first and is of the right type.
|
|
if data is None and time is None and flux is None and flux_err is None:
|
|
self._required_columns_relax = True
|
|
return
|
|
|
|
# Load `time`, `flux`, and `flux_err` from the table as local variable names
|
|
time = self.columns["time"] # super().__init__() guarantees this is a column
|
|
if "flux" in self.colnames:
|
|
if flux is None:
|
|
flux = self.columns["flux"]
|
|
else:
|
|
raise TypeError(
|
|
f"'flux' has been given both in the `data` table and as a keyword argument"
|
|
)
|
|
if "flux_err" in self.colnames:
|
|
if flux_err is None:
|
|
flux_err = self.columns["flux_err"]
|
|
else:
|
|
raise TypeError(
|
|
f"'flux_err' has been given both in the `data` table and as a keyword argument"
|
|
)
|
|
|
|
# Ensure `flux` and `flux_err` are populated with NaNs if missing
|
|
if flux is None and time is not None:
|
|
flux = np.empty(len(time))
|
|
flux[:] = np.nan
|
|
if not isinstance(flux, Quantity):
|
|
flux = Quantity(flux, deprecated_kws.get("flux_unit"))
|
|
|
|
if flux_err is None:
|
|
flux_err = np.empty(len(flux))
|
|
flux_err[:] = np.nan
|
|
if not isinstance(flux_err, Quantity):
|
|
flux_err = Quantity(flux_err, flux.unit)
|
|
|
|
# Backwards compatibility with Lightkurve v1.x
|
|
# Ensure attributes are set if passed via deprecated kwargs
|
|
for kw in deprecated_kws:
|
|
if kw not in self.meta:
|
|
self.meta[kw.upper()] = deprecated_kws[kw]
|
|
|
|
# Ensure all required columns are in the right order
|
|
with self._delay_required_column_checks():
|
|
for idx, col in enumerate(self._required_columns):
|
|
if col in self.colnames:
|
|
self.remove_column(col)
|
|
self.add_column(locals()[col], index=idx, name=col)
|
|
|
|
# Ensure columns are set if passed via deprecated kwargs
|
|
for kw in deprecated_column_kws:
|
|
if kw not in self.meta and kw not in self.columns:
|
|
self.add_column(deprecated_column_kws[kw], name=kw)
|
|
|
|
# Ensure flux and flux_err have the same units
|
|
if self["flux"].unit != self["flux_err"].unit:
|
|
raise ValueError("flux and flux_err must have the same units")
|
|
|
|
self._new_attributes_relax = False
|
|
self._required_columns_relax = False
|
|
self._check_required_columns()
|
|
|
|
def __getattr__(self, name, **kwargs):
|
|
"""Expose all columns and meta keywords as attributes."""
|
|
if name in self.__dict__:
|
|
return self.__dict__[name]
|
|
elif name in self.__class__.__dict__:
|
|
return self.__class__.__dict__[name].__get__(self)
|
|
elif name in self.columns:
|
|
return self[name]
|
|
elif "_meta" in self.__dict__:
|
|
if name in self.__dict__["_meta"]:
|
|
return self.__dict__["_meta"][name]
|
|
elif name.upper() in self.__dict__["_meta"]:
|
|
return self.__dict__["_meta"][name.upper()]
|
|
raise AttributeError(f"object has no attribute {name}")
|
|
|
|
def __setattr__(self, name, value, **kwargs):
|
|
"""To get copied, attributes have to be stored in the meta dictionary!"""
|
|
to_set_as_attr = False
|
|
if name in self.__dict__:
|
|
to_set_as_attr = True
|
|
elif name == "time":
|
|
self["time"] = value # astropy will convert value to Time if needed
|
|
elif ("columns" in self.__dict__) and (name in self.__dict__["columns"]):
|
|
self.replace_column(name, value)
|
|
elif "_meta" in self.__dict__:
|
|
if name in self.__dict__["_meta"]:
|
|
self.__dict__["_meta"][name] = value
|
|
elif name.upper() in self.__dict__["_meta"]:
|
|
self.__dict__["_meta"][name.upper()] = value
|
|
else:
|
|
to_set_as_attr = True
|
|
else:
|
|
to_set_as_attr = True
|
|
if to_set_as_attr:
|
|
if (
|
|
name not in self.__dict__
|
|
and not name.startswith("_")
|
|
and not self._new_attributes_relax
|
|
and name != 'meta'
|
|
):
|
|
warnings.warn(
|
|
(
|
|
"Lightkurve doesn't allow columns or meta values to be created via a new attribute name."
|
|
"A new attribute is created. It will not be carried over when the object is copied."
|
|
" - see https://docs.lightkurve.org/reference/api/lightkurve.LightCurve.html"
|
|
),
|
|
UserWarning,
|
|
stacklevel=2,
|
|
)
|
|
super().__setattr__(name, value, **kwargs)
|
|
|
|
def _repr_simple_(self) -> str:
|
|
"""Returns a simple __repr__.
|
|
|
|
Used by `LightCurveCollection`.
|
|
"""
|
|
result = f"<{self.__class__.__name__}"
|
|
if "LABEL" in self.meta:
|
|
result += f" LABEL=\"{self.meta.get('LABEL')}\""
|
|
for kw in ["QUARTER", "CAMPAIGN", "SECTOR", "AUTHOR", "FLUX_ORIGIN"]:
|
|
if kw in self.meta:
|
|
result += f" {kw}={self.meta.get(kw)}"
|
|
result += ">"
|
|
return result
|
|
|
|
def _base_repr_(self, html=False, descr_vals=None, **kwargs):
|
|
"""Defines the description shown by `__repr__` and `_html_repr_`."""
|
|
if descr_vals is None:
|
|
descr_vals = [self.__class__.__name__]
|
|
if self.masked:
|
|
descr_vals.append("masked=True")
|
|
descr_vals.append("length={}".format(len(self)))
|
|
if "LABEL" in self.meta:
|
|
descr_vals.append(f"LABEL=\"{self.meta.get('LABEL')}\"")
|
|
for kw in ["QUARTER", "CAMPAIGN", "SECTOR", "AUTHOR", "FLUX_ORIGIN"]:
|
|
if kw in self.meta:
|
|
descr_vals.append(f"{kw}={self.meta.get(kw)}")
|
|
return super()._base_repr_(html=html, descr_vals=descr_vals, **kwargs)
|
|
|
|
# Define `time`, `flux`, `flux_err` as class attributes to enable IDE
|
|
# of these required columns auto-completion.
|
|
|
|
@property
|
|
def time(self) -> Time:
|
|
"""Time values stored as an AstroPy `~astropy.time.Time` object."""
|
|
return self["time"]
|
|
|
|
@time.setter
|
|
def time(self, time):
|
|
self["time"] = time
|
|
|
|
@property
|
|
def flux(self) -> Quantity:
|
|
"""Brightness values stored as an AstroPy `~astropy.units.Quantity` object."""
|
|
return self["flux"]
|
|
|
|
@flux.setter
|
|
def flux(self, flux):
|
|
self["flux"] = flux
|
|
|
|
@property
|
|
def flux_err(self) -> Quantity:
|
|
"""Brightness uncertainties stored as an AstroPy `~astropy.units.Quantity` object."""
|
|
return self["flux_err"]
|
|
|
|
@flux_err.setter
|
|
def flux_err(self, flux_err):
|
|
self["flux_err"] = flux_err
|
|
|
|
def flatten(
|
|
self,
|
|
window_length=101,
|
|
polyorder=2,
|
|
return_trend=False,
|
|
break_tolerance=5,
|
|
niters=3,
|
|
sigma=3,
|
|
mask=None,
|
|
**kwargs,
|
|
):
|
|
if mask is None:
|
|
mask = np.ones(len(self.time), dtype=bool)
|
|
else:
|
|
# Deep copy ensures we don't change the original.
|
|
mask = deepcopy(~mask)
|
|
# Add NaNs & outliers to the mask
|
|
extra_mask = np.isfinite(self.flux)
|
|
extra_mask &= np.nan_to_num(np.abs(self.flux - np.nanmedian(self.flux))) <= (
|
|
np.nanstd(self.flux) * sigma
|
|
)
|
|
# In astropy>=5.0, extra_mask is a masked array
|
|
if hasattr(extra_mask, 'mask'):
|
|
mask &= extra_mask.filled(False)
|
|
else: # support astropy<5.0
|
|
mask &= extra_mask
|
|
|
|
for iter in np.arange(0, niters):
|
|
if break_tolerance is None:
|
|
break_tolerance = np.nan
|
|
if polyorder >= window_length:
|
|
polyorder = window_length - 1
|
|
# Split the lightcurve into segments by finding large gaps in time
|
|
dt = self.time.value[mask][1:] - self.time.value[mask][0:-1]
|
|
with warnings.catch_warnings(): # Ignore warnings due to NaNs
|
|
warnings.simplefilter("ignore", RuntimeWarning)
|
|
cut = np.where(dt > break_tolerance * np.nanmedian(dt))[0] + 1
|
|
low = np.append([0], cut)
|
|
high = np.append(cut, len(self.time[mask]))
|
|
# Then, apply the savgol_filter to each segment separately
|
|
trend_signal = Quantity(np.zeros(len(self.time[mask])), unit=self.flux.unit)
|
|
for l, h in zip(low, high):
|
|
# Reduce `window_length` and `polyorder` for short segments;
|
|
# this prevents `savgol_filter` from raising an exception
|
|
# If the segment is too short, just take the median
|
|
if np.any([window_length > (h - l), (h - l) < break_tolerance]):
|
|
trend_signal[l:h] = np.nanmedian(self.flux[mask][l:h])
|
|
else:
|
|
# Scipy outputs a warning here that is not useful, will be fixed in version 1.2
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("ignore", FutureWarning)
|
|
trsig = savgol_filter(
|
|
x=self.flux.value[mask][l:h],
|
|
window_length=window_length,
|
|
polyorder=polyorder,
|
|
**kwargs,
|
|
)
|
|
trend_signal[l:h] = Quantity(trsig, trend_signal.unit)
|
|
# Ignore outliers; note we add `1e-14` below to avoid detecting
|
|
# outliers which are merely caused by numerical noise.
|
|
mask1 = np.nan_to_num(np.abs(self.flux[mask] - trend_signal)) < (
|
|
np.nanstd(self.flux[mask] - trend_signal) * sigma
|
|
+ Quantity(1e-14, self.flux.unit)
|
|
)
|
|
f = interp1d(
|
|
self.time.value[mask][mask1],
|
|
trend_signal[mask1],
|
|
fill_value="extrapolate",
|
|
)
|
|
trend_signal = Quantity(f(self.time.value), self.flux.unit)
|
|
# In astropy>=5.0, mask1 is a masked array
|
|
if hasattr(mask1, 'mask'):
|
|
mask[mask] &= mask1.filled(False)
|
|
else: # support astropy<5.0
|
|
mask[mask] &= mask1
|
|
|
|
flatten_lc = self.copy()
|
|
with warnings.catch_warnings():
|
|
# ignore invalid division warnings
|
|
warnings.simplefilter("ignore", RuntimeWarning)
|
|
flatten_lc.flux = flatten_lc.flux / trend_signal
|
|
flatten_lc.flux_err = flatten_lc.flux_err / trend_signal
|
|
|
|
flatten_lc.meta["NORMALIZED"] = True
|
|
|
|
if return_trend:
|
|
trend_lc = self.copy()
|
|
trend_lc.flux = trend_signal
|
|
return flatten_lc, trend_lc
|
|
return flatten_lc
|
|
|
|
def to_periodogram(self, method="lombscargle", **kwargs):
|
|
supported_methods = ["ls", "bls", "lombscargle", "boxleastsquares"]
|
|
method = validate_method(method.replace(" ", ""), supported_methods)
|
|
if method in ["bls", "boxleastsquares"]:
|
|
from .MinimalPeriodogram import BoxLeastSquaresPeriodogram
|
|
|
|
return BoxLeastSquaresPeriodogram.from_lightcurve(lc=self, **kwargs)
|
|
else:
|
|
from .MinimalPeriodogram import LombScarglePeriodogram
|
|
|
|
return LombScarglePeriodogram.from_lightcurve(lc=self, **kwargs)
|
|
|
|
def fold(
|
|
self,
|
|
period=None,
|
|
epoch_time=None,
|
|
epoch_phase=0,
|
|
wrap_phase=None,
|
|
normalize_phase=False,
|
|
):
|
|
# Lightkurve v1.x assumed that `period` was given in days if no unit
|
|
# was specified. We maintain this behavior for backwards-compatibility.
|
|
if period is not None and not isinstance(period, Quantity):
|
|
period *= u.day
|
|
if epoch_time is not None and not isinstance(epoch_time, Time):
|
|
epoch_time = Time(
|
|
epoch_time, format=self.time.format, scale=self.time.scale
|
|
)
|
|
if (
|
|
epoch_phase is not None
|
|
and not isinstance(epoch_phase, Quantity)
|
|
and not normalize_phase
|
|
):
|
|
epoch_phase *= u.day
|
|
if wrap_phase is not None and not isinstance(wrap_phase, Quantity):
|
|
wrap_phase *= u.day
|
|
|
|
# Warn if `epoch_time` appears to use the wrong format
|
|
if epoch_time is not None and epoch_time.value > 2450000:
|
|
if self.time.format == "bkjd":
|
|
warnings.warn(
|
|
"`epoch_time` appears to be given in JD, "
|
|
"however the light curve time uses BKJD "
|
|
"(i.e. JD - 2454833).",
|
|
LightkurveWarning,
|
|
)
|
|
elif self.time.format == "btjd":
|
|
warnings.warn(
|
|
"`epoch_time` appears to be given in JD, "
|
|
"however the light curve time uses BTJD "
|
|
"(i.e. JD - 2457000).",
|
|
LightkurveWarning,
|
|
)
|
|
|
|
ts = super().fold(
|
|
period=period,
|
|
epoch_time=epoch_time,
|
|
epoch_phase=epoch_phase,
|
|
wrap_phase=wrap_phase,
|
|
normalize_phase=normalize_phase,
|
|
)
|
|
|
|
# The folded time would pass the `TimeSeries` validation check if
|
|
# `normalize_phase=True`, so creating a `FoldedLightCurve` object
|
|
# requires the following three-step workaround:
|
|
# 1. Give the folded light curve a valid time column again
|
|
with ts._delay_required_column_checks():
|
|
folded_time = ts.time.copy()
|
|
ts.remove_column("time")
|
|
ts.add_column(self.time, name="time", index=0)
|
|
# 2. Create the folded object
|
|
lc = FoldedLightCurve(data=ts)
|
|
# 3. Restore the folded time
|
|
with lc._delay_required_column_checks():
|
|
lc.remove_column("time")
|
|
lc.add_column(folded_time, name="time", index=0)
|
|
|
|
# Add extra column and meta data specific to FoldedLightCurve
|
|
lc.add_column(
|
|
self.time.copy(), name="time_original", index=len(self._required_columns)
|
|
)
|
|
lc.meta["PERIOD"] = period
|
|
lc.meta["EPOCH_TIME"] = epoch_time
|
|
lc.meta["EPOCH_PHASE"] = epoch_phase
|
|
lc.meta["WRAP_PHASE"] = wrap_phase
|
|
lc.meta["NORMALIZE_PHASE"] = normalize_phase
|
|
lc.sort("time")
|
|
|
|
return lc
|
|
|
|
def normalize(self, unit="unscaled"):
|
|
validate_method(unit, ["unscaled", "percent", "ppt", "ppm"])
|
|
median_flux = np.nanmedian(self.flux)
|
|
std_flux = np.nanstd(self.flux)
|
|
|
|
# If the median flux is within half a standard deviation from zero, the
|
|
# light curve is likely zero-centered and normalization makes no sense.
|
|
if (median_flux == 0) or (
|
|
np.isfinite(std_flux) and (np.abs(median_flux) < 0.5 * std_flux)
|
|
):
|
|
warnings.warn(
|
|
"The light curve appears to be zero-centered "
|
|
"(median={:.2e} +/- {:.2e}); `normalize()` will divide "
|
|
"the light curve by a value close to zero, which is "
|
|
"probably not what you want."
|
|
"".format(median_flux, std_flux),
|
|
LightkurveWarning,
|
|
)
|
|
# If the median flux is negative, normalization will invert the light
|
|
# curve and makes no sense.
|
|
if median_flux < 0:
|
|
warnings.warn(
|
|
"The light curve has a negative median flux ({:.2e});"
|
|
" `normalize()` will therefore divide by a negative "
|
|
"number and invert the light curve, which is probably"
|
|
"not what you want".format(median_flux),
|
|
LightkurveWarning,
|
|
)
|
|
|
|
# Create a new light curve instance and normalize its values
|
|
lc = self.copy()
|
|
lc.flux = lc.flux / median_flux
|
|
lc.flux_err = lc.flux_err / median_flux
|
|
if not lc.flux.unit:
|
|
lc.flux *= u.dimensionless_unscaled
|
|
if not lc.flux_err.unit:
|
|
lc.flux_err *= u.dimensionless_unscaled
|
|
|
|
# Set the desired relative (dimensionless) units
|
|
if unit == "percent":
|
|
lc.flux = lc.flux.to(u.percent)
|
|
lc.flux_err = lc.flux_err.to(u.percent)
|
|
elif unit in ("ppt", "ppm"):
|
|
lc.flux = lc.flux.to(unit)
|
|
lc.flux_err = lc.flux_err.to(unit)
|
|
|
|
lc.meta["NORMALIZED"] = True
|
|
return lc
|
|
|
|
def remove_nans(self, column: str = "flux"):
|
|
return self[~np.isnan(self[column])] # This will return a sliced copy
|
|
|
|
class FoldedLightCurve(LightCurve):
|
|
|
|
@property
|
|
def phase(self):
|
|
"""Alias for `LightCurve.time`."""
|
|
return self.time
|
|
|
|
@property
|
|
def cycle(self):
|
|
cycle_epoch_start = self.epoch_time - self.period / 2
|
|
result = np.asarray(np.floor(((self.time_original - cycle_epoch_start) / self.period).value), dtype=int)
|
|
result = result - result.min()
|
|
return result
|
|
|
|
@property
|
|
def odd_mask(self):
|
|
return self.cycle % 2 == 1
|
|
|
|
@property
|
|
def even_mask(self):
|
|
return ~self.odd_mask
|
|
|
|
def _set_xlabel(self, kwargs):
|
|
if "xlabel" not in kwargs:
|
|
kwargs["xlabel"] = "Phase"
|
|
if isinstance(self.time, TimeDelta):
|
|
kwargs["xlabel"] += f" [{self.time.format.upper()}]"
|
|
return kwargs
|
|
|
|
def plot(self, **kwargs):
|
|
kwargs = self._set_xlabel(kwargs)
|
|
return super(FoldedLightCurve, self).plot(**kwargs)
|
|
|
|
def scatter(self, **kwargs):
|
|
kwargs = self._set_xlabel(kwargs)
|
|
return super(FoldedLightCurve, self).scatter(**kwargs)
|
|
|
|
def errorbar(self, **kwargs):
|
|
kwargs = self._set_xlabel(kwargs)
|
|
return super(FoldedLightCurve, self).errorbar(**kwargs)
|
|
|
|
def plot_river(self, **kwargs):
|
|
ax = super(FoldedLightCurve, self).plot_river(
|
|
period=self.period, epoch_time=self.epoch_time, **kwargs
|
|
)
|
|
return ax
|
|
|
|
class KeplerLightCurve(LightCurve):
|
|
"""Subclass of :class:`LightCurve <lightkurve.lightcurve.LightCurve>`
|
|
to represent data from NASA's Kepler and K2 mission."""
|
|
|
|
_deprecated_keywords = (
|
|
"targetid",
|
|
"label",
|
|
"time_format",
|
|
"time_scale",
|
|
"flux_unit",
|
|
"quality_bitmask",
|
|
"channel",
|
|
"campaign",
|
|
"quarter",
|
|
"mission",
|
|
"ra",
|
|
"dec",
|
|
)
|
|
|
|
_default_time_format = "bkjd"
|
|
|
|
@classmethod
|
|
def read(cls, *args, **kwargs):
|
|
# Default to Kepler file format
|
|
if kwargs.get("format") is None:
|
|
kwargs["format"] = "kepler2"
|
|
return super().read(*args, **kwargs)
|
|
|
|
class TessLightCurve(LightCurve):
|
|
"""Subclass of :class:`LightCurve <lightkurve.lightcurve.LightCurve>`
|
|
to represent data from NASA's TESS mission."""
|
|
|
|
_deprecated_keywords = (
|
|
"targetid",
|
|
"label",
|
|
"time_format",
|
|
"time_scale",
|
|
"flux_unit",
|
|
"quality_bitmask",
|
|
"sector",
|
|
"camera",
|
|
"ccd",
|
|
"mission",
|
|
"ra",
|
|
"dec",
|
|
)
|
|
|
|
_default_time_format = "btjd"
|
|
|
|
@classmethod
|
|
def read(cls, *args, **kwargs):
|
|
# Default to TESS file format
|
|
if kwargs.get("format") is None:
|
|
kwargs["format"] = "tess2"
|
|
return super().read(*args, **kwargs)
|
|
|
|
# Helper functions
|
|
|
|
def _boolean_mask_to_bitmask(aperture_mask):
|
|
# Masks can either be boolean input or Kepler pipeline style
|
|
clean_mask = np.nan_to_num(aperture_mask)
|
|
|
|
contains_bit2 = (clean_mask.astype(np.int_) & 2).any()
|
|
all_zeros_or_ones = (clean_mask.dtype in ["float", "int"]) & (
|
|
(set(np.unique(clean_mask)) - {0, 1}) == set()
|
|
)
|
|
is_bool_mask = (aperture_mask.dtype == "bool") | all_zeros_or_ones
|
|
|
|
if is_bool_mask:
|
|
out_mask = np.ones(aperture_mask.shape, dtype=np.uint8)
|
|
out_mask[aperture_mask == 1] = 3
|
|
out_mask = out_mask.astype(np.uint8)
|
|
elif contains_bit2:
|
|
out_mask = aperture_mask.astype(np.uint8)
|
|
else:
|
|
out_mask = None
|
|
return out_mask
|
|
|
|
|
|
def _make_aperture_extension(hdu_list, aperture_mask):
|
|
"""Returns an `ImageHDU` object containing the 'APERTURE' extension
|
|
of a light curve file."""
|
|
if aperture_mask is not None:
|
|
bitmask = _boolean_mask_to_bitmask(aperture_mask)
|
|
hdu = fits.ImageHDU(bitmask)
|
|
hdu.header["EXTNAME"] = "APERTURE"
|
|
hdu_list.append(hdu)
|
|
return hdu_list
|
|
|
|
def read_generic_lightcurve(
|
|
filename,
|
|
time_column="time",
|
|
flux_column="flux",
|
|
flux_err_column="flux_err",
|
|
quality_column="quality",
|
|
cadenceno_column="cadenceno",
|
|
centroid_col_column="mom_centr1",
|
|
centroid_row_column="mom_centr2",
|
|
time_format=None,
|
|
ext=1,
|
|
):
|
|
"""Generic helper function to convert a Kepler ot TESS light curve file
|
|
into a generic `LightCurve` object.
|
|
"""
|
|
if isinstance(filename, fits.HDUList):
|
|
hdulist = filename # Allow HDUList to be passed
|
|
else:
|
|
with fits.open(filename) as hdulist:
|
|
hdulist = deepcopy(hdulist)
|
|
|
|
# Raise an exception if the requested extension is invalid
|
|
if isinstance(ext, str):
|
|
validate_method(ext, supported_methods=[hdu.name.lower() for hdu in hdulist])
|
|
with warnings.catch_warnings():
|
|
# By default, AstroPy emits noisy warnings about units commonly used
|
|
# in archived TESS data products (e.g., "e-/s" and "pixels").
|
|
# We ignore them here because they don't affect Lightkurve's features.
|
|
# Inconsistencies between TESS data products and the FITS standard
|
|
# out to be addressed at the archive level. (See issue #1216.)
|
|
warnings.simplefilter("ignore", category=UnitsWarning)
|
|
tab = Table.read(hdulist[ext], format="fits")
|
|
|
|
# Make sure the meta data also includes header fields from extension #0
|
|
tab.meta.update(hdulist[0].header)
|
|
|
|
tab.meta = {k: v for k, v in tab.meta.items()}
|
|
|
|
for colname in tab.colnames:
|
|
# Ensure units have the correct astropy format
|
|
# Speed-up: comparing units by their string representation is 1000x
|
|
# faster than performing full-blown unit comparison
|
|
unitstr = str(tab[colname].unit)
|
|
if unitstr == "e-/s":
|
|
tab[colname].unit = "electron/s"
|
|
elif unitstr == "pixels":
|
|
tab[colname].unit = "pixel"
|
|
elif unitstr == "ppm" and repr(tab[colname].unit).startswith("Unrecognized"):
|
|
# Workaround for issue #956
|
|
tab[colname].unit = ppm
|
|
elif unitstr == "ADU":
|
|
tab[colname].unit = "adu"
|
|
elif unitstr.lower() == "unitless":
|
|
tab[colname].unit = ""
|
|
elif unitstr.lower() == "degcelcius":
|
|
# CDIPS has non-astropy units
|
|
tab[colname].unit = "deg_C"
|
|
# Rename columns to lowercase
|
|
tab.rename_column(colname, colname.lower())
|
|
|
|
# Some KEPLER files used to have a T column instead of TIME.
|
|
if time_column == "time" and "time" not in tab.columns and "t" in tab.colnames:
|
|
tab.rename_column("t", "time")
|
|
if time_column != "time":
|
|
tab.rename_column(time_column, "time")
|
|
|
|
# We *have* to remove rows with TIME=NaN because the Astropy Time
|
|
# object does not support the presence of NaNs.
|
|
# Fortunately, such rows are always bad data.
|
|
nans = np.isnan(tab["time"].data)
|
|
if np.any(nans):
|
|
log.debug("Ignoring {} rows with NaN times".format(np.sum(nans)))
|
|
tab = tab[~nans]
|
|
|
|
# Prepare a special time column
|
|
if not time_format:
|
|
if hdulist[ext].header.get("BJDREFI") == 2454833:
|
|
time_format = "bkjd"
|
|
elif hdulist[ext].header.get("BJDREFI") == 2457000:
|
|
time_format = "btjd"
|
|
else:
|
|
raise ValueError(f"Input file has unclear time format: {filename}")
|
|
time = Time(
|
|
tab["time"].data,
|
|
scale=hdulist[ext].header.get("TIMESYS", "tdb").lower(),
|
|
format=time_format,
|
|
)
|
|
tab.remove_column("time")
|
|
|
|
# For backwards compatibility with Lightkurve v1.x,
|
|
# we make sure standard columns and attributes exist.
|
|
if "flux" not in tab.columns:
|
|
tab.add_column(tab[flux_column], name="flux", index=0)
|
|
if "flux_err" not in tab.columns:
|
|
# Try falling back to `{flux_column}_err` if possible
|
|
if flux_err_column not in tab.columns:
|
|
flux_err_column = flux_column + "_err"
|
|
if flux_err_column in tab.columns:
|
|
tab.add_column(tab[flux_err_column], name="flux_err", index=1)
|
|
if "quality" not in tab.columns and quality_column in tab.columns:
|
|
tab.add_column(tab[quality_column], name="quality", index=2)
|
|
if "cadenceno" not in tab.columns and cadenceno_column in tab.columns:
|
|
tab.add_column(tab[cadenceno_column], name="cadenceno", index=3)
|
|
if "centroid_col" not in tab.columns and centroid_col_column in tab.columns:
|
|
tab.add_column(tab[centroid_col_column], name="centroid_col", index=4)
|
|
if "centroid_row" not in tab.columns and centroid_row_column in tab.columns:
|
|
tab.add_column(tab[centroid_row_column], name="centroid_row", index=5)
|
|
|
|
tab.meta["LABEL"] = hdulist[0].header.get("OBJECT")
|
|
tab.meta["MISSION"] = hdulist[0].header.get(
|
|
"MISSION", hdulist[0].header.get("TELESCOP")
|
|
)
|
|
tab.meta["RA"] = hdulist[0].header.get("RA_OBJ")
|
|
tab.meta["DEC"] = hdulist[0].header.get("DEC_OBJ")
|
|
tab.meta["FILENAME"] = filename
|
|
tab.meta["FLUX_ORIGIN"] = flux_column
|
|
|
|
return LightCurve(time=time, data=tab)
|
|
|
|
def read_kepler_lightcurve(
|
|
filename, flux_column="pdcsap_flux", quality_bitmask="default"
|
|
):
|
|
lc = read_generic_lightcurve(
|
|
filename,
|
|
flux_column=flux_column,
|
|
quality_column="sap_quality",
|
|
time_format="bkjd",
|
|
)
|
|
quality_mask = KeplerQualityFlags.create_quality_mask(
|
|
quality_array=lc["sap_quality"], bitmask=quality_bitmask
|
|
)
|
|
lc = lc[quality_mask]
|
|
|
|
lc.meta["AUTHOR"] = "Kepler"
|
|
lc.meta["TARGETID"] = lc.meta.get("KEPLERID")
|
|
lc.meta["QUALITY_BITMASK"] = quality_bitmask
|
|
lc.meta["QUALITY_MASK"] = quality_mask
|
|
return KeplerLightCurve(data=lc)
|
|
|
|
def read_tess_lightcurve(
|
|
filename, flux_column="pdcsap_flux", quality_bitmask="default"
|
|
):
|
|
lc = read_generic_lightcurve(filename, flux_column=flux_column, time_format="btjd")
|
|
|
|
quality_mask = TessQualityFlags.create_quality_mask(
|
|
quality_array=lc["quality"], bitmask=quality_bitmask
|
|
)
|
|
lc = lc[quality_mask]
|
|
|
|
lc.meta["AUTHOR"] = "SPOC"
|
|
lc.meta["TARGETID"] = lc.meta.get("TICID")
|
|
lc.meta["QUALITY_BITMASK"] = quality_bitmask
|
|
lc.meta["QUALITY_MASK"] = quality_mask
|
|
return TessLightCurve(data=lc)
|
|
|
|
def read_qlp_lightcurve(filename, flux_column="sap_flux", flux_err_column="kspsap_flux_err", quality_bitmask="default"):
|
|
lc = read_generic_lightcurve(filename, flux_column=flux_column, flux_err_column=flux_err_column, time_format="btjd")
|
|
|
|
quality_mask = TessQualityFlags.create_quality_mask(
|
|
quality_array=lc["quality"], bitmask=quality_bitmask
|
|
)
|
|
lc = lc[quality_mask]
|
|
|
|
lc.meta["AUTHOR"] = "QLP"
|
|
lc.meta["TARGETID"] = lc.meta.get("TICID")
|
|
lc.meta["QUALITY_BITMASK"] = quality_bitmask
|
|
lc.meta["QUALITY_MASK"] = quality_mask
|
|
|
|
# QLP light curves are normalized by default
|
|
lc.meta["NORMALIZED"] = True
|
|
|
|
return TessLightCurve(data=lc)
|
|
|
|
def read_eleanor_lightcurve(filename,
|
|
flux_column="CORR_FLUX",
|
|
quality_bitmask="default"
|
|
):
|
|
lc = read_generic_lightcurve(
|
|
filename,
|
|
time_column="TIME".lower(),
|
|
flux_column=flux_column.lower(),
|
|
flux_err_column = "FLUX_ERR".lower(),
|
|
time_format="btjd",
|
|
quality_column= "QUALITY".lower(),
|
|
centroid_col_column = "X_CENTROID".lower(),
|
|
centroid_row_column = "Y_CENTROID".lower(),
|
|
cadenceno_column = "FFIINDEX".lower()
|
|
)
|
|
|
|
if quality_bitmask == "hardest":
|
|
# Eleanor has 2 additional bits on top of the 16 TESS SPOC bits
|
|
# they are excluded when hardest is specified.
|
|
quality_bitmask = TessQualityFlags.HARDEST_BITMASK | 2** 17 | 2**18
|
|
|
|
quality_mask = TessQualityFlags.create_quality_mask(
|
|
quality_array=lc["quality"], bitmask=quality_bitmask
|
|
)
|
|
|
|
lc = lc[quality_mask]
|
|
|
|
# Eleanor FITS file do not have units specified. re-add them.
|
|
for colname in ["flux", "flux_err", "raw_flux", "corr_flux", "pca_flux", "psf_flux"]:
|
|
if colname in lc.colnames:
|
|
if lc[colname].unit is not None:
|
|
# for case flux, flux_err, lightkurve has forced it to be u.dimensionless_unscaled
|
|
# can't reset a unit, so we create a new column
|
|
lc[colname] = u.Quantity(lc[colname].value, "electron/s")
|
|
else:
|
|
lc[colname].unit = "electron/s"
|
|
|
|
for colname in ["flux_bkg"]:
|
|
if colname in lc.colnames:
|
|
lc[colname].unit = u.percent
|
|
|
|
for colname in ["centroid_col", "centroid_row", "x_centroid", "y_centroid", "x_com", "y_com"]:
|
|
if colname in lc.colnames:
|
|
lc[colname].unit = u.pix
|
|
|
|
for colname in ["barycorr"]:
|
|
if colname in lc.colnames:
|
|
lc[colname].unit = u.day
|
|
|
|
# In Eleanor fits file, raw_flux's error is in flux_err, which breaks Lightkurve convention.
|
|
# To account for this, the corr_flux error is calculated from corr_flux_err = corr_flux*raw_flux_err/raw_flux. For completeness,
|
|
# the original raw_flux's error is added as a "raw_flux_err" column
|
|
lc["raw_flux_err"] = lc["flux_err"]
|
|
if flux_column.lower() != 'raw_flux':
|
|
lc["flux_err"] = lc[flux_column.lower()]*lc["raw_flux_err"]/lc["raw_flux"]
|
|
|
|
# vanilla eleanor has cadence saved as float,
|
|
# convert to int to ensure we stick with the convention
|
|
for colname in ["ffiindex", "cadenceno"]:
|
|
if colname in lc.colnames:
|
|
if not np.issubdtype(lc[colname].dtype, np.int_):
|
|
lc[colname] = np.asarray(lc[colname].value, dtype=int)
|
|
|
|
if (
|
|
lc.meta.get("TVERSION") is not None
|
|
and lc.meta.get("GITHUB") == "https://github.com/afeinstein20/eleanor"
|
|
):
|
|
# the above headers are GSFC-ELEANOR-LITE-specific, and are not present in vanilla eleanor
|
|
# cf. https://github.com/afeinstein20/eleanor/blob/main/eleanor/targetdata.py
|
|
lc.meta["AUTHOR"] = "GSFC-ELEANOR-LITE"
|
|
else:
|
|
lc.meta["AUTHOR"] = "ELEANOR"
|
|
|
|
# Eleanor light curves are not normalized by default
|
|
lc.meta["NORMALIZED"] = False
|
|
|
|
tic = lc.meta.get("TIC_ID")
|
|
if tic is not None:
|
|
# compatibility with SPOC, QLP, etc.
|
|
lc.meta["TARGETID"] = tic
|
|
lc.meta["TICID"] = tic
|
|
lc.meta["OBJECT"] = f"TIC {tic}"
|
|
# for Lightkurve's plotting methods
|
|
lc.meta["LABEL"] = f"TIC {tic}"
|
|
|
|
return TessLightCurve(data=lc)
|
|
|
|
def read_k2sff_lightcurve(filename, ext="BESTAPER", **kwargs):
|
|
lc = read_generic_lightcurve(
|
|
filename, flux_column="fcor", time_format="bkjd", ext=ext
|
|
)
|
|
|
|
lc.meta["AUTHOR"] = "K2SFF"
|
|
lc.meta["TARGETID"] = lc.meta.get("KEPLERID")
|
|
|
|
return KeplerLightCurve(data=lc, **kwargs)
|
|
|
|
def read_everest_lightcurve(
|
|
filename, flux_column="flux", quality_bitmask="default", **kwargs
|
|
):
|
|
lc = read_generic_lightcurve(
|
|
filename,
|
|
flux_column=flux_column,
|
|
quality_column="quality",
|
|
cadenceno_column="cadn",
|
|
time_format="bkjd",
|
|
)
|
|
|
|
quality_mask = KeplerQualityFlags.create_quality_mask(
|
|
quality_array=lc["quality"], bitmask=quality_bitmask
|
|
)
|
|
lc = lc[quality_mask]
|
|
|
|
lc.meta["AUTHOR"] = "EVEREST"
|
|
lc.meta["TARGETID"] = lc.meta.get("KEPLERID")
|
|
lc.meta["QUALITY_BITMASK"] = quality_bitmask
|
|
lc.meta["QUALITY_MASK"] = quality_mask
|
|
|
|
return KeplerLightCurve(data=lc, **kwargs)
|
|
|
|
def read_pathos_lightcurve(
|
|
filename, flux_column="PSF_FLUX_COR", quality_bitmask="default"
|
|
):
|
|
lc = read_generic_lightcurve(
|
|
filename,
|
|
flux_column=flux_column.lower(),
|
|
time_format="btjd",
|
|
quality_column="DQUALITY",
|
|
)
|
|
|
|
quality_mask = TessQualityFlags.create_quality_mask(
|
|
quality_array=lc["dquality"], bitmask=quality_bitmask
|
|
)
|
|
lc = lc[quality_mask]
|
|
|
|
lc.meta["AUTHOR"] = "PATHOS"
|
|
lc.meta["TARGETID"] = lc.meta.get("TICID")
|
|
lc.meta["QUALITY_BITMASK"] = quality_bitmask
|
|
lc.meta["QUALITY_MASK"] = quality_mask
|
|
|
|
# QLP light curves are normalized by default
|
|
lc.meta["NORMALIZED"] = True
|
|
|
|
return TessLightCurve(data=lc)
|
|
|
|
def read_cdips_lightcurve(filename,
|
|
flux_column="IRM1",
|
|
include_inst_errs=False,
|
|
quality_bitmask=None):
|
|
ap = flux_column[-1]
|
|
|
|
if include_inst_errs:
|
|
# If fluxes are requested, return flux errors
|
|
if flux_column[:-1].lower()=="ifl":
|
|
flux_err_column = f"ife{ap}"
|
|
# Otherwise magnitudes are being requested, return magnitude errors
|
|
else:
|
|
flux_err_column = f"ire{ap}"
|
|
else:
|
|
flux_err_column = ""
|
|
|
|
# Set the appropriate error column for this aperture
|
|
quality_column = f"irq{ap}"
|
|
|
|
lc = read_generic_lightcurve(filename,
|
|
time_column="tmid_bjd",
|
|
flux_column=flux_column.lower(),
|
|
flux_err_column=flux_err_column,
|
|
quality_column=quality_column,
|
|
time_format='btjd')
|
|
|
|
quality_mask = (lc['quality']=="G") | (lc['quality']=="0")
|
|
lc = lc[quality_mask]
|
|
|
|
lc.meta["AUTHOR"] = "CDIPS"
|
|
lc.meta['TARGETID'] = lc.meta.get('TICID')
|
|
lc.meta['QUALITY_BITMASK'] = 36
|
|
lc.meta['QUALITY_MASK'] = quality_mask
|
|
|
|
return TessLightCurve(data=lc)
|
|
|
|
def read_tasoc_lightcurve(filename, flux_column="FLUX_CORR", quality_bitmask=None):
|
|
lc = read_generic_lightcurve(
|
|
filename, flux_column=flux_column.lower(), time_format="btjd"
|
|
)
|
|
|
|
lc.meta["AUTHOR"] = "TASOC"
|
|
lc.meta["TARGETID"] = lc.meta.get("TICID")
|
|
# TASOC light curves are normalized by default
|
|
lc.meta["NORMALIZED"] = True
|
|
return TessLightCurve(data=lc)
|
|
|
|
def read_kepseismic_lightcurve(filename, **kwargs):
|
|
lc = read_generic_lightcurve(
|
|
filename,
|
|
time_format='bkjd')
|
|
|
|
lc.meta["AUTHOR"] = "KEPSEISMIC"
|
|
lc.meta["TARGETID"] = lc.meta.get("KEPLERID")
|
|
|
|
# KEPSEISMIC light curves are normalized by default
|
|
lc.meta["NORMALIZED"] = True
|
|
|
|
return KeplerLightCurve(data=lc, **kwargs)
|
|
|
|
def read_tglc_lightcurve(
|
|
filename, flux_column="cal_psf_flux", quality_bitmask="default"
|
|
):
|
|
lc = read_generic_lightcurve(
|
|
filename,
|
|
time_column="time",
|
|
flux_column=flux_column.lower(),
|
|
quality_column="tess_flags",
|
|
cadenceno_column="cadence_num",
|
|
time_format="btjd",
|
|
)
|
|
|
|
quality_mask = TessQualityFlags.create_quality_mask(
|
|
quality_array=lc["quality"], bitmask=quality_bitmask
|
|
)
|
|
|
|
# TGLC FITS file do not have units specified. re-add them.
|
|
for colname in ["psf_flux", "aperture_flux", "background"]:
|
|
if colname in lc.colnames:
|
|
if lc[colname].unit is not None:
|
|
# for case flux, flux_err, lightkurve has forced it to be u.dimensionless_unscaled
|
|
# can't reset a unit, so we create a new column
|
|
lc[colname] = u.Quantity(
|
|
lc[colname].value, "electron/s", dtype=np.float32
|
|
)
|
|
else:
|
|
lc[colname].unit = "electron/s"
|
|
|
|
# Calibrated columns are normalized, so they are unitless
|
|
for colname in ["cal_psf_flux", "cal_aper_flux"]:
|
|
if colname in lc.colnames:
|
|
if lc[colname].unit is not None:
|
|
# for case flux, flux_err, lightkurve has forced it to be u.dimensionless_unscaled
|
|
# can't reset a unit, so we create a new column
|
|
lc[colname] = u.Quantity(lc[colname].value, "", dtype=np.float32)
|
|
else:
|
|
lc[colname].unit = ""
|
|
|
|
lc = lc[quality_mask]
|
|
lc.meta["AUTHOR"] = "TGLC"
|
|
lc.meta["TARGETID"] = lc.meta.get("OBJECT")
|
|
lc.meta["QUALITY_BITMASK"] = quality_bitmask
|
|
lc.meta["QUALITY_MASK"] = quality_mask
|
|
lc.meta["NORMALIZED"] = True
|
|
tic = lc.meta.get("TICID")
|
|
if tic is not None:
|
|
tic = int(tic)
|
|
# compatibility with SPOC, QLP, etc.
|
|
lc.meta["TARGETID"] = tic
|
|
lc.meta["TICID"] = tic
|
|
lc.meta["OBJECT"] = f"TIC {tic}"
|
|
# for Lightkurve's plotting methods
|
|
lc.meta["LABEL"] = f"TIC {tic}"
|
|
return TessLightCurve(data=lc)
|
|
|
|
try:
|
|
registry.register_reader("kepler2", LightCurve, read_kepler_lightcurve)
|
|
registry.register_reader(data_format="tess2", data_class=LightCurve, function=read_tess_lightcurve)
|
|
registry.register_reader("qlp2", LightCurve, read_qlp_lightcurve)
|
|
registry.register_reader("eleanor2", LightCurve, read_eleanor_lightcurve)
|
|
registry.register_reader("k2sff2", LightCurve, read_k2sff_lightcurve)
|
|
registry.register_reader("everest2", LightCurve, read_everest_lightcurve)
|
|
registry.register_reader("pathos2", LightCurve, read_pathos_lightcurve)
|
|
registry.register_reader("cdips2", LightCurve, read_cdips_lightcurve)
|
|
registry.register_reader("tasoc2", LightCurve, read_tasoc_lightcurve)
|
|
registry.register_reader(
|
|
"kepseismic2", LightCurve, read_kepseismic_lightcurve
|
|
)
|
|
registry.register_reader("tglc2", LightCurve, read_tglc_lightcurve)
|
|
except registry.IORegistryError as exc:
|
|
print(exc)
|
|
pass # necessary to enable autoreload during debugging
|
|
|
|
|