diff --git a/main/astrodatagui/CalcAllFlaresThread.py b/main/astrodatagui/CalcAllFlaresThread.py index 967877a..5f10ac1 100644 --- a/main/astrodatagui/CalcAllFlaresThread.py +++ b/main/astrodatagui/CalcAllFlaresThread.py @@ -1,18 +1,18 @@ from PyQt5.QtCore import pyqtSignal, QThread import multiprocessing import concurrent.futures -import lightkurve as lk +from ..util.MinimalLightCurve import read import pandas as pd from ..flaredetector.flaredetector import calculateFlareFitsForLightcurve def getFlareCount(filesDict): - lc = lk.read(filesDict["FilePath"]) + lc = read(filesDict["FilePath"]) lc.flux = lc["sap_flux"] lc.flux_err = lc["sap_flux_err"] - sapPeaks, sapFits = calculateFlareFitsForLightcurve(lc.flatten()) + sapPeaks, sapFits = calculateFlareFitsForLightcurve(lc.flatten(), minimalLC=True) lc.flux = lc["pdcsap_flux"] lc.flux_err = lc["pdcsap_flux_err"] - pdcsapPeaks, pdcsapFits = calculateFlareFitsForLightcurve(lc.flatten()) + pdcsapPeaks, pdcsapFits = calculateFlareFitsForLightcurve(lc.flatten(), minimalLC=True) filesDict["sapPeaks"] = sapPeaks filesDict["sapPeaksCount"] = len(sapPeaks) @@ -21,7 +21,8 @@ def getFlareCount(filesDict): filesDict["pdcsapPeaksCount"] = len(pdcsapPeaks) filesDict["pdcsapFits"] = pdcsapFits del lc - + if(filesDict["StarName"] == "CD-51 13128" and filesDict["Sequence"] == 1): + print(filesDict) return filesDict class CalcAllFlaresThread(QThread): @@ -35,7 +36,6 @@ class CalcAllFlaresThread(QThread): #@pyqtSlot() def run(self): cpuCount = multiprocessing.cpu_count() - splitList = [self.allFlaresDictList[i:i + cpuCount] for i in range(0, len(self.allFlaresDictList), cpuCount)] executor = concurrent.futures.ProcessPoolExecutor(cpuCount) resFrame = pd.DataFrame(executor.map(getFlareCount, self.allFlaresDictList.to_dict(orient="records"))) diff --git a/main/flaredetector/flaredetector.py b/main/flaredetector/flaredetector.py index f75dc39..754ef7f 100644 --- a/main/flaredetector/flaredetector.py +++ b/main/flaredetector/flaredetector.py @@ -40,8 +40,9 @@ def calculateFitForFlare(lightcurve, peakIndex): flareFitDataPoints = flareFitDataPoints + peakIndex return np.array(Peak+1), flareFitDataPoints, np.array(fit+1) -def calculateFlareFitsForLightcurve(lightcurve, num=100, distance=1, height=0.05): - maxIndices = findMaxIndices(lightcurve, num, distance=distance, height=height, sortByHighest=True) +def calculateFlareFitsForLightcurve(lightcurve, num=100, distance=1, height=0.05, minimalLC=False): + maxIndices = findMaxIndices(lightcurve, num, distance=distance, height=height, + sortByHighest=True, minimalLC=minimalLC) peaks = [] fits = [] for index in maxIndices: @@ -57,7 +58,6 @@ def calculateFlareFitsForLightcurve(lightcurve, num=100, distance=1, height=0.05 "Quality": lightcurve.quality[flareFitDataPoints]}) except Exception as error: pass - #print("Flare fit exception: ", error) return peaks, fits diff --git a/main/flaredetector/util.py b/main/flaredetector/util.py index d4c8f8f..fa7f56b 100644 --- a/main/flaredetector/util.py +++ b/main/flaredetector/util.py @@ -1,16 +1,19 @@ import numpy as np from scipy.signal import find_peaks -from lightkurve.periodogram import Periodogram -from lightkurve.lightcurve import LightCurve - -def findMaxIndices(lc, num=3, distance=100, height=(None, None), sortByHighest=False): - if(isinstance(lc, Periodogram)): - data = lc.power - elif(isinstance(lc, LightCurve)): +def findMaxIndices(lc, num=3, distance=100, height=(None, None), sortByHighest=False, + minimalLC=False): + if(minimalLC): data = lc.flux else: - return np.zeros(num)-1 + from lightkurve.periodogram import Periodogram as OrigPG + from lightkurve.lightcurve import LightCurve as OrigLC + if(isinstance(lc, OrigPG)): + data = lc.power + elif(isinstance(lc, OrigLC)): + data = lc.flux + else: + return np.zeros(num)-1 peak_indices, peak_dict = find_peaks(data, height=height, distance=distance) peak_heights = peak_dict["peak_heights"] diff --git a/main/util/MinimalLightCurve.py b/main/util/MinimalLightCurve.py new file mode 100644 index 0000000..9df8ff5 --- /dev/null +++ b/main/util/MinimalLightCurve.py @@ -0,0 +1,1420 @@ +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