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 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 _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 class KeplerLightCurve(LightCurve): """Subclass of :class:`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 ` 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