diff --git a/docs/tutorials/tutorial1_basics.ipynb b/docs/tutorials/tutorial1_basics.ipynb index c4ba4a9..ea61643 100644 --- a/docs/tutorials/tutorial1_basics.ipynb +++ b/docs/tutorials/tutorial1_basics.ipynb @@ -24724,7 +24724,7 @@ "This is flagging bad quality data with a mask which combines the SPOC quality flags and the `tess_asteroids` lightcurve quality flags. By default it masks the following bits:\n", "\n", "- For SPOC, it masks bits [1, 2, 3, 4, 5, 6, 8, 10, 13, 15] as recommended in the [TESS Archive Manual](https://outerspace.stsci.edu/spaces/TESS/pages/14563420/2.0+-+Data+Product+Overview).\n", - "- For `tess_asteroids`, it masks bits [2, 4, 10] as defined by `_create_lc_quality()`. This corresponds to cadences with at least one non-science pixel inside the mask, at least one saturated pixel inside the mask or a pixel with a negative raw flux (i.e. before background correction).\n", + "- For `tess_asteroids`, it masks bits [2, 4, 10, 14] as defined by `_create_lc_quality()`. This corresponds to cadences with at least one non-science pixel inside the mask, at least one saturated pixel inside the mask or at least one pixel with a negative raw flux (i.e. before background correction) inside the mask, or cadences which are in the region of non-science data during sector 3.\n", "\n", "In the above example, there is no bad quality data when using the default mask. The bits to mask can be updated independently for SPOC, the aperture lightcurve and the PSF lightcurve:" ] @@ -83531,7 +83531,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.18" + "version": "3.10.14" } }, "nbformat": 4, diff --git a/pyproject.toml b/pyproject.toml index a5cadad..2efec11 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "tess-asteroids" -version = "1.5.0" +version = "1.5.1" description = "Create TPFs and LCs for solar system asteroids observed by NASA's TESS mission." license = "MIT" authors = ["Amy Tuson ", diff --git a/src/tess_asteroids/__init__.py b/src/tess_asteroids/__init__.py index 9812753..9e88eba 100644 --- a/src/tess_asteroids/__init__.py +++ b/src/tess_asteroids/__init__.py @@ -41,10 +41,10 @@ # (https://outerspace.stsci.edu/display/TESS/2.0+-+Data+Product+Overview) default_bad_spoc_bits = [1, 2, 3, 4, 5, 6, 8, 10, 13, 15] # Default bad bits for LC quality masking: at least one pixel inside mask is non-science, saturated -# or negative raw flux. -default_bad_lc_bits = [2, 4, 10] +# or negative raw flux, or cadence is in region of non-science data during sector 3. +default_bad_lc_bits = [2, 4, 10, 14] -__version__ = "1.5.0" +__version__ = "1.5.1" __all__ = ["MovingTPF"] from .movingtpf import MovingTPF # noqa: E402 diff --git a/src/tess_asteroids/movingtpf.py b/src/tess_asteroids/movingtpf.py index 859c0ef..9501fff 100644 --- a/src/tess_asteroids/movingtpf.py +++ b/src/tess_asteroids/movingtpf.py @@ -382,6 +382,12 @@ def get_data(self, shape: Tuple[int, int] = (11, 11)): row, column = row[bound_mask], column[bound_mask] pixel_coordinates = np.asarray([row.ravel(), column.ravel()]).T + # Check there are at least two cadences after masking. + if len(self.time_original) <= 1: + raise RuntimeError( + "There are less than two cadences with pixels fully inside bounds of FFI (1<=row<=2078, 1<=col<=2136)." + ) + # Check there are pixels inside FFI bounds. if len(pixel_coordinates) == 0: raise RuntimeError( @@ -841,7 +847,11 @@ def _data_chunks(self, data_chunks: Optional[Union[list, np.ndarray]] = None): self.sector, len(sector_downlinks) ) ) - elif self.sector > 55 and len(sector_downlinks) != 4: + elif ( + self.sector > 55 + and self.sector not in [97, 98] + and len(sector_downlinks) != 4 + ): logger.warning( "For sector {0} there should be 4 data chunks, but there are actually {1} data chunks. Investigate or define your own `data_chunks`.".format( self.sector, len(sector_downlinks) @@ -1252,8 +1262,6 @@ def _create_spoc_quality_mask( """ Creates a boolean mask using SPOC quality flags. - In all cases except `bad_spoc_bits="none"`, non-science data in sector 3 is also masked. - Parameters ---------- spoc_quality : ndarray @@ -1297,13 +1305,6 @@ def _create_spoc_quality_mask( else: spoc_quality_mask = spoc_quality & bad_spoc_bit_value == 0 - # Add non-science data in sector 3 to quality mask, as defined in data release notes. - if self.sector == 3: - t = self.time_original - self.timecorr_original - spoc_quality_mask = np.logical_and( - spoc_quality_mask, np.logical_and(t >= 1385.89663, t <= 1406.29247) - ) - return spoc_quality_mask, bad_spoc_bits def _bg_linear_model( @@ -2618,7 +2619,7 @@ def _psf_photometry( # If all data has been masked, raise warning and return empty arrays. if len(time) == 0: logger.warning( - "During PSF photometry, all times were masked and no PSF light curve was derived." + "During PSF photometry, all times were masked and no PSF lightcurve was derived." ) return { "time": np.array([]), @@ -2895,6 +2896,7 @@ def _psf_photometry( pixel_mask=pixel_masks, pixel_quality=pixel_qualities, prf_nan_mask=prf_nan_mask, + time_sc=(time - time_corr, time_uerr, time_lerr), psf_flux=flux, ), "n_cadences": n_cadences.astype(int), @@ -3078,6 +3080,7 @@ def _aperture_photometry(self, bad_bits: list = [1, 3, 7], **kwargs): pixel_mask=self.aperture_mask, pixel_quality=self.pixel_quality, prf_nan_mask=self.ap_prf_nan_mask, + time_sc=self.time - self.timecorr, bad_bit_value=bad_bit_value, ), "flux_fraction": np.asarray(flux_fraction), @@ -3103,6 +3106,7 @@ def _create_lc_quality( pixel_mask: Union[list, np.ndarray], pixel_quality: Union[list, np.ndarray], prf_nan_mask: np.ndarray, + time_sc: Union[np.ndarray, Tuple[np.ndarray, np.ndarray, np.ndarray]], bad_bit_value: Optional[int] = None, psf_flux: Optional[np.ndarray] = None, ): @@ -3113,12 +3117,14 @@ def _create_lc_quality( If you are creating LC quality for aperture photometry: - `pixel_mask` should correspond to the pixels inside the aperture mask. + - `time_sc` should be an array corresponding to time at the spacecraft. - `bad_bit_value` should correspond to `bad_bits` defined by user in `_aperture_photometry()`. - `psf_flux` should be None. If you are creating LC quality for PSF photometry: - `pixel_mask` should correspond to the pixels that were used to fit the PRF model in `_psf_photometry()`. + - `time_sc` should be a tuple of arrays for time at the spacecraft, upper error on bin and lower error on bin. - `bad_bit_value` should be None. - `psf_flux` should be the PSF flux time-series from `_psf_photometry()`. @@ -3139,6 +3145,7 @@ def _create_lc_quality( - 11 - PSF fit failed due to singular matrix (see np.linalg.LinAlgError) or because all pixels used to fit the PRF model had NaN flux values. Only relevant if `method=psf`. - 12 - at least one pixel inside mask had a poor fitting background star model. Only relevant if `linear_model` background correction was used. - 13 - mask has two or more discrete regions. + - 14 - non-science data from sector 3, as defined in data release notes: https://archive.stsci.edu/missions/tess/doc/tess_drn/tess_sector_03_drn04_v02.pdf Parameters ---------- @@ -3150,6 +3157,10 @@ def _create_lc_quality( prf_nan_mask : ndarray Cadences where the PRF model is NaN, as returned by `_create_target_prf_model()`. This must have length equal to the first dimension of `pixel_mask`. + time_sc : ndarray or tuple + For aperture photometry, this should be an ndarray corresponding to time at the spacecraft. + For PSF photometry, this should be a tuple of arrays for time at the spacecraft, upper error on bin and lower error on bin. + This is used to identify non-science data in sector 3, and each array must have length equal to the first dimension of `pixel_mask`. bad_bit_value : int Value corresponding to `bad_bits`, as defined by user in `_aperture_photometry()`. If creating quality mask for PSF photmetry, this should be None. @@ -3170,6 +3181,20 @@ def _create_lc_quality( "Must define either `bad_bit_value` or `psf_flux`, but not both." ) + # Define sector 3 non-science mask + # https://archive.stsci.edu/missions/tess/doc/tess_drn/tess_sector_03_drn04_v02.pdf + if self.sector == 3: + if isinstance(time_sc, np.ndarray): + time_upper, time_lower = time_sc, time_sc + elif isinstance(time_sc, tuple): + time_upper = time_sc[0] + np.nan_to_num(time_sc[1], nan=0) + time_lower = time_sc[0] - np.nan_to_num(time_sc[2], nan=0) + sector3_mask = np.logical_or( + time_lower < 1385.89663, time_upper > 1406.29247 + ) + else: + sector3_mask = np.zeros(len(pixel_mask), dtype=bool) + # Define bits masks = { # No pixels in mask @@ -3290,6 +3315,11 @@ def _create_lc_quality( for t in range(len(pixel_mask)) ], }, + # Non-science data from sector 3 + "sector3_mask": { + "bit": 14, + "value": sector3_mask, + }, # Add flag for negative pixels (after BG correction) in aperture? } @@ -4508,62 +4538,81 @@ def _make_primary_hdu( after="AP{0}_PAR".format(i if len(ap_masks) > 1 else ""), ) - # Add keywords for object properties - hdu.header.set( - "VMAG", - round( - np.nanmean( - self.ephem.loc[ - np.logical_and( - self.ephem["time"] - >= ( - self.time_original[0] - if self.barycentric - else self.time_original[0] - self.timecorr_original[0] + # Compute average predicted V and H magnitude of target + # Catch warnings that arise if arrays are empty: + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message="Mean of empty slice", + category=RuntimeWarning, + ) + vmag = ( + round( + np.nanmean( + self.ephem.loc[ + np.logical_and( + self.ephem["time"] + >= ( + self.time_original[0] + if self.barycentric + else self.time_original[0] + - self.timecorr_original[0] + ), + self.ephem["time"] + <= ( + self.time_original[-1] + if self.barycentric + else self.time_original[-1] + - self.timecorr_original[-1] + ), ), - self.ephem["time"] - <= ( - self.time_original[-1] - if self.barycentric - else self.time_original[-1] - self.timecorr_original[-1] + "vmag", + ] + ), + 3, + ) + if "vmag" in self.ephem + else np.nan + ) + hmag = ( + round( + np.nanmean( + self.ephem.loc[ + np.logical_and( + self.ephem["time"] + >= ( + self.time_original[0] + if self.barycentric + else self.time_original[0] + - self.timecorr_original[0] + ), + self.ephem["time"] + <= ( + self.time_original[-1] + if self.barycentric + else self.time_original[-1] + - self.timecorr_original[-1] + ), ), - ), - "vmag", - ] - ), - 3, + "hmag", + ] + ), + 3, + ) + if "hmag" in self.ephem and ~np.isnan(self.ephem["hmag"]).all() + else np.nan ) - if "vmag" in self.ephem - else None, + + # Add keywords for object properties + hdu.header.set( + "VMAG", + vmag if ~np.isnan(vmag) else None, comment="[mag] predicted V magnitude", after="EQUINOX", ) hdu.header.set( "HMAG", - round( - np.nanmean( - self.ephem.loc[ - np.logical_and( - self.ephem["time"] - >= ( - self.time_original[0] - if self.barycentric - else self.time_original[0] - self.timecorr_original[0] - ), - self.ephem["time"] - <= ( - self.time_original[-1] - if self.barycentric - else self.time_original[-1] - self.timecorr_original[-1] - ), - ), - "hmag", - ] - ), - 3, - ) - if "hmag" in self.ephem and ~np.isnan(self.ephem["hmag"]).all() - else None, + hmag if ~np.isnan(hmag) else None, comment="[mag] predicted H absolute magnitude", after="VMAG", ) @@ -4886,6 +4935,13 @@ def create_lc_quality_mask( `True` indicates a good quality cadence and `False` indicates a bad quality cadence. """ + # If either `spoc_quality` or `lc_quality` is an empty array, return an empty array and warn user. + if len(spoc_quality) == 0 or len(lc_quality) == 0: + logger.warning( + "`spoc_quality` and/or `lc_quality` are empty arrays. The lightcurve `quality_mask` is returned as an empty array." + ) + return np.asarray([]) + # Create SPOC quality mask using user-defined bad bits. spoc_quality_mask, _ = self._create_spoc_quality_mask( spoc_quality=spoc_quality, bad_spoc_bits=bad_spoc_bits @@ -5016,35 +5072,38 @@ def plot_lc( if key == "aperture": # Run through each available aperture lightcurve for i, lc_ap in enumerate(lc[key]): - # Define aperture lightcurve quality mask. - quality_mask = self.create_lc_quality_mask( - spoc_quality=self.quality, - lc_quality=lc_ap["quality"], - bad_spoc_bits=bad_spoc_bits, - bad_lc_bits=bad_ap_bits, - ) + if len(lc_ap["time"]) > 0: + # Define aperture lightcurve quality mask. + quality_mask = self.create_lc_quality_mask( + spoc_quality=self.quality, + lc_quality=lc_ap["quality"], + bad_spoc_bits=bad_spoc_bits, + bad_lc_bits=bad_ap_bits, + ) - ax[i].errorbar( - lc_ap["time"][quality_mask], - lc_ap["flux"][quality_mask], - yerr=lc_ap["flux_err"][quality_mask] if plot_err else None, - color="deeppink", - marker="o", - ms=2, - ls="", - ) - if plot_bad_quality: ax[i].errorbar( - lc_ap["time"][~quality_mask], - lc_ap["flux"][~quality_mask], - yerr=lc_ap["flux_err"][~quality_mask] if plot_err else None, - color="black", - marker="x", - ms=6, + lc_ap["time"][quality_mask], + lc_ap["flux"][quality_mask], + yerr=lc_ap["flux_err"][quality_mask] if plot_err else None, + color="deeppink", + marker="o", + ms=2, ls="", - label="Bad quality data", ) - ax[i].legend() + if plot_bad_quality: + ax[i].errorbar( + lc_ap["time"][~quality_mask], + lc_ap["flux"][~quality_mask], + yerr=lc_ap["flux_err"][~quality_mask] + if plot_err + else None, + color="black", + marker="x", + ms=6, + ls="", + label="Bad quality data", + ) + ax[i].legend() ax[i].tick_params( axis="x", which="both", labelbottom=True, bottom=True ) @@ -5063,47 +5122,50 @@ def plot_lc( ) elif key == "psf": - # Define PSF lightcurve quality mask. - quality_mask = self.create_lc_quality_mask( - spoc_quality=lc[key]["spoc_quality"], - lc_quality=lc[key]["quality"], - bad_spoc_bits=bad_spoc_bits, - bad_lc_bits=bad_psf_bits, - ) - - ax[-1].errorbar( - lc[key]["time"][quality_mask], - lc[key]["flux"][quality_mask], - yerr=lc[key]["flux_err"][quality_mask] if plot_err else None, - xerr=( - lc[key]["time_lerr"][quality_mask], - lc[key]["time_uerr"][quality_mask], + if len(lc[key]["time"]) > 0: + # Define PSF lightcurve quality mask. + quality_mask = self.create_lc_quality_mask( + spoc_quality=lc[key]["spoc_quality"], + lc_quality=lc[key]["quality"], + bad_spoc_bits=bad_spoc_bits, + bad_lc_bits=bad_psf_bits, ) - if plot_err - else None, - color="mediumorchid", - marker="o", - ms=2, - ls="", - ) - if plot_bad_quality: + ax[-1].errorbar( - lc[key]["time"][~quality_mask], - lc[key]["flux"][~quality_mask], - yerr=lc[key]["flux_err"][~quality_mask] if plot_err else None, + lc[key]["time"][quality_mask], + lc[key]["flux"][quality_mask], + yerr=lc[key]["flux_err"][quality_mask] if plot_err else None, xerr=( - lc[key]["time_lerr"][~quality_mask], - lc[key]["time_uerr"][~quality_mask], + lc[key]["time_lerr"][quality_mask], + lc[key]["time_uerr"][quality_mask], ) if plot_err else None, - color="black", - marker="x", - ms=6, + color="mediumorchid", + marker="o", + ms=2, ls="", - label="Bad quality data", ) - ax[-1].legend() + if plot_bad_quality: + ax[-1].errorbar( + lc[key]["time"][~quality_mask], + lc[key]["flux"][~quality_mask], + yerr=lc[key]["flux_err"][~quality_mask] + if plot_err + else None, + xerr=( + lc[key]["time_lerr"][~quality_mask], + lc[key]["time_uerr"][~quality_mask], + ) + if plot_err + else None, + color="black", + marker="x", + ms=6, + ls="", + label="Bad quality data", + ) + ax[-1].legend() ax[-1].set_ylabel("Flux [e-/s]") ax[-1].grid(ls=":") ax[-1].set_axisbelow(True) diff --git a/tests/test_movingtpf.py b/tests/test_movingtpf.py index 5011f95..3c41a68 100644 --- a/tests/test_movingtpf.py +++ b/tests/test_movingtpf.py @@ -264,9 +264,7 @@ def test_pixel_quality(): """ # Create artificial track for a region that includes a saturated star, straps and non-science pixels. - # Note: the track must have length > 1, but we use two timestamps very close together so resulting - # moving TPF only has one frame. - time = np.linspace(2458328.50, 2458328.52, 2) - 2457000 + time = np.linspace(2458328.5, 2458329.0, 2) - 2457000 track = pd.DataFrame( { "time": time,