Repository: lightkurve/lightkurve

Issue:
TESS Quality Flag default bitmask update
#### Problem description

This is a follow up of issue #1163 / PR #1218 , regarding what tweaks should be done for the default TESS quality bitmask.

There is some recommendation in [TESS Data Product documentation](https://outerspace.stsci.edu/display/TESS/2.0+-+Data+Product+Overview#:~:text=As%20a%20starting%20point%2C%20the%20following%20flags%20indicate%20that%20the%20cadence%20will%20likely%20be%20of%20lesser%20quality
).

I tabulate the suggestion in the doc and compare it against the current `default` and `hard` bitmasks

| Bit Description                                      | bin digit # | TESS Doc | `default` | `hard` |
| ---------------------------------------------------- | ----------- | -------- | ------- | ---- |
| AttitudeTweak = 1                                    | 1           | ✓        | ✓       | ✓    |
| SafeMode = 2                                         | 2           | ✓        | ✓       | ✓    |
| CoarsePoint = 4                                      | 3           | ✓        | ✓       | ✓    |
| EarthPoint = 8                                       | 4           | ✓        | ✓       | ✓    |
| Argabrightening = 16                                 | 5           | ✓        |         |      |
| Desat = 32                                           | 6           | ✓        | ✓       | ✓    |
| ApertureCosmic = 64                                  | 7           |          |         | ✓    |
| ManualExclude = 128                                  | 8           | ✓        | ✓       | ✓    |
| Discontinuity = 256                                  | 9           |          |         |      |
| ImpulsiveOutlier = 512                               | 10          | ✓        |         |      |
| CollateralCosmic = 1024                              | 11          |          |         | ✓    |
| Straylight = 2048 #: set by MIT manually             | 12          |          |         | ✓    |
| Straylight2 = 4096 #: set by Ames/SPOC automatically | 13          | ✓        |         | ✓    |
| PlanetSearchExclude = 8192                           | 14          |          |         |      |
| BadCalibrationExclude = 16384                        | 15          | ✓        |         |      |
| InsufficientTargets = 32768                          | 16          |          |         |       |


If we follow the recommendation, the following bits will be added to `default` (and indirectly added to `hard`
- bit 5, Argabrightening = 16
- bit 10, ImpulsiveOutlier = 512	
- bit 13, Straylight2 = 4096 #: set by Ames/SPOC automatically, already in hard
- bit 15, BadCalibrationExclude = 16384


#### Environment

-  lightkurve version (e.g. 1.0b6): 2.2

Hints from issue comments:
@Nschanche: I made this change in the above #1336. However, this causes the test_underfit_metric_neighbors to fail. It looks like the noise properties of the updated 'default' setting, particularly from removing the ImpulsiveOutliers flag, change the underfit metric. If we do want to make this change, I can look into altering the test. @Nschanche: I have updated the PR #1336 to add Argabrightening, Impulsive Outlier, and Bad Calibration Exclude flags to the 'default' quality mask. I think that including the stray light flag by default is overly cautious. While the data may be lower in quality, it can in most cases still be used. @christinahedges do you agree with the changes in the default quality flags?@Nschanche: I have merged #1336, so I'll close this now. This adds Argabrightening, Impulsive Outlier, and Bad Calibration Exclude flags (but not stray light) to the 'default' TESS quality mask. 

Solution patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index ee06f6f23..3bcef4f5b 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -8,7 +8,7 @@
   the new ``flux`` column is different from that of the ``flux_err`` column. [#1467]
 - Fixed numpy v2 compatibility for ``tpf.interact()``. [#1473]
 - Modified RegressionCorrector.correct() to add the model flux error in quadrature. [#1439]
-
+- Updated the default TESS quality mask to include Argabrightening, Impulsive Outlier, and Bad Calibration Exclude flags (see issue #1230). [#1336]
 
 2.5.0 (2024-08-29)
 =====================
diff --git a/src/lightkurve/correctors/metrics.py b/src/lightkurve/correctors/metrics.py
index 17dda4bd6..e93f0dc76 100644
--- a/src/lightkurve/correctors/metrics.py
+++ b/src/lightkurve/correctors/metrics.py
@@ -12,6 +12,7 @@
 from scipy.interpolate import PchipInterpolator
 from memoization import cached
 from astropy import units as u
+from typing import Union
 
 from .. import LightCurve
 
@@ -129,6 +130,7 @@ def underfit_metric_neighbors(
     max_targets: int = 50,
     interpolate: bool = False,
     extrapolate: bool = False,
+    quality_bitmask: Union[int, str] = "default",
 ):
     """This goodness metric measures the degree of under-fitting of the
     CBVs to the light curve. It does so by measuring the mean residual target to
@@ -187,6 +189,7 @@ def underfit_metric_neighbors(
         interpolate=interpolate,
         extrapolate=extrapolate,
         flux_column="sap_flux",
+        quality_bitmask=quality_bitmask,
     )
 
     # Create fluxMatrix. The last entry is the target under study
@@ -199,9 +202,11 @@ def underfit_metric_neighbors(
     # Add in the trimmed target under study
     fluxMatrix[:, -1] = corrected_lc_flux
 
-    # Ignore NaNs
-    mask = ~np.isnan(corrected_lc_flux)
-    fluxMatrix = fluxMatrix[mask, :]
+    # Remove NaNs from any flux column
+    mask = np.zeros(fluxMatrix.shape[0], dtype=bool)
+    for i in range(fluxMatrix.shape[1]):
+        mask |= np.isnan(fluxMatrix[:, i])
+    fluxMatrix = fluxMatrix[~mask, :]
 
     # Determine the target-target correlation between target and
     # neighborhood
@@ -264,10 +269,11 @@ def _unique_key_for_processing_neighbors(
     extrapolate: bool = False,
     author: tuple = ("Kepler", "K2", "SPOC"),
     flux_column: str = "sap_flux",
+    quality_bitmask: Union[str,int] = "default",
 ):
     """Returns a unique key that will determine whether a cached version of a
     call to `_download_and_preprocess_neighbors` can be re-used."""
-    return f"{corrected_lc.ra}{corrected_lc.dec}{corrected_lc.cadenceno}{radius}{min_targets}{max_targets}{author}{flux_column}{interpolate}"
+    return f"{corrected_lc.ra}{corrected_lc.dec}{corrected_lc.cadenceno}{radius}{min_targets}{max_targets}{author}{flux_column}{interpolate}{quality_bitmask}"
 
 
 @cached(custom_key_maker=_unique_key_for_processing_neighbors)
@@ -280,6 +286,7 @@ def _download_and_preprocess_neighbors(
     extrapolate: bool = False,
     author: tuple = ("Kepler", "K2", "SPOC"),
     flux_column: str = "sap_flux",
+    quality_bitmask: Union[str, int] = "default",
 ):
     """Returns a collection of neighboring light curves.
 
@@ -318,7 +325,7 @@ def _download_and_preprocess_neighbors(
         raise Exception('interpolate must be True if extrapolate is True')
 
     search = corrected_lc.search_neighbors(
-        limit=max_targets, radius=radius, author=author
+        limit=max_targets, radius=radius, author=author,
     )
     if len(search) < min_targets:
         raise MinTargetsError(
@@ -327,7 +334,7 @@ def _download_and_preprocess_neighbors(
     log.info(
         f"Downloading {len(search)} neighboring light curves. This might take a while."
     )
-    lcfCol = search.download_all(flux_column=flux_column)
+    lcfCol = search.download_all(flux_column=flux_column, quality_bitmask=quality_bitmask)
 
     # Pre-process the neighboring light curves
     # Align or interpolate to the corrected light curve
@@ -336,7 +343,7 @@ def _download_and_preprocess_neighbors(
     # Extract SAP light curves
     # We want zero-centered median normalized light curves
     for lc in lcfCol:
-        lcSAP = lc.remove_nans().normalize()
+        lcSAP = lc.remove_nans(column=flux_column).normalize()
         lcSAP.flux -= 1.0
         # Align or interpolate the neighboring target with the target under study
         if interpolate:
diff --git a/src/lightkurve/correctors/pldcorrector.py b/src/lightkurve/correctors/pldcorrector.py
index d362d9e73..f418aadd4 100644
--- a/src/lightkurve/correctors/pldcorrector.py
+++ b/src/lightkurve/correctors/pldcorrector.py
@@ -250,7 +250,7 @@ def create_design_matrix(
             # Create a DesignMatrix for each PLD order
             all_pld = []
             for order in range(1, pld_order + 1):
-                reg_n = np.product(list(multichoose(regressors_pld.T, order)), axis=1).T
+                reg_n = np.prod(list(multichoose(regressors_pld.T, order)), axis=1).T
                 pld_n = DesignMatrix(
                     reg_n,
                     prior_sigma=np.ones(reg_n.shape[1]) * prior_sigma / reg_n.shape[1],
diff --git a/src/lightkurve/targetpixelfile.py b/src/lightkurve/targetpixelfile.py
index 54b0814df..bdbc3a8ab 100644
--- a/src/lightkurve/targetpixelfile.py
+++ b/src/lightkurve/targetpixelfile.py
@@ -855,8 +855,8 @@ def _estimate_centroids_via_quadratic(self, aperture_mask):
             col, row = centroid_quadratic(self.flux[idx], mask=aperture_mask)
             col_centr.append(col)
             row_centr.append(row)
-        col_centr = np.asfarray(col_centr) + self.column
-        row_centr = np.asfarray(row_centr) + self.row
+        col_centr = np.asarray(col_centr, dtype=float) + self.column
+        row_centr = np.asarray(row_centr, dtype=float) + self.row
         col_centr = Quantity(col_centr, unit="pixel")
         row_centr = Quantity(row_centr, unit="pixel")
         return col_centr, row_centr
@@ -1554,7 +1554,7 @@ def cutout(self, center=None, size=5):
                     datacol._dims = datacol.array.shape[1:]
                     datacol._format = fits.column._ColumnFormat(
                         "{}{}".format(
-                            np.product(datacol.array.shape[1:]), datacol._format[-1]
+                            np.prod(datacol.array.shape[1:]), datacol._format[-1]
                         )
                     )
                     data_columns.append(datacol)
diff --git a/src/lightkurve/utils.py b/src/lightkurve/utils.py
index fe965d3fd..821638e58 100644
--- a/src/lightkurve/utils.py
+++ b/src/lightkurve/utils.py
@@ -250,8 +250,10 @@ class TessQualityFlags(QualityFlags):
     InsufficientTargets = 32768
 
     #: DEFAULT bitmask identifies all cadences which are definitely useless.
+    # See https://outerspace.stsci.edu/display/TESS/2.0+-+Data+Product+Overview
     DEFAULT_BITMASK = (
-        AttitudeTweak | SafeMode | CoarsePoint | EarthPoint | Desat | ManualExclude
+        AttitudeTweak | SafeMode | CoarsePoint | EarthPoint | Argabrightening | 
+        Desat | ManualExclude | ImpulsiveOutlier | BadCalibrationExclude
     )
     #: HARD bitmask is conservative and may identify cadences which are useful.
     HARD_BITMASK = (


Test patch:
diff --git a/tests/correctors/test_metrics.py b/tests/correctors/test_metrics.py
index b8bb63be7..63d3b6f00 100644
--- a/tests/correctors/test_metrics.py
+++ b/tests/correctors/test_metrics.py
@@ -39,14 +39,13 @@ def test_underfit_metric_neighbors():
     """Sanity checks for `underfit_metric_neighbors`."""
     # PDCSAP_FLUX has a very good score (>0.99) because it has been corrected
     lc_pdcsap = search_lightcurve("Proxima Cen", sector=11, author="SPOC").download(
-        flux_column="pdcsap_flux"
+        flux_column="pdcsap_flux", 
     )
     assert underfit_metric_neighbors(lc_pdcsap, min_targets=3, max_targets=3) > 0.99
-
     # SAP_FLUX has a worse score (<0.95) because it hasn't been corrected
-    lc_sap = lc_pdcsap.copy()
-    lc_sap.flux = lc_pdcsap.sap_flux
-    lc_sap.flux_err = lc_pdcsap.sap_flux_err
+    lc_sap = search_lightcurve("Proxima Cen", sector=11, author="SPOC").download(
+        flux_column="sap_flux", 
+    )
     assert underfit_metric_neighbors(lc_sap, min_targets=3, max_targets=3) < 0.95
 
     # A flat light curve should have a perfect score (1)
diff --git a/tests/test_targetpixelfile.py b/tests/test_targetpixelfile.py
index 79103e82d..67f4a1f7a 100644
--- a/tests/test_targetpixelfile.py
+++ b/tests/test_targetpixelfile.py
@@ -473,7 +473,7 @@ def test_tpf_from_images():
     w.wcs.crval = [0, -90]
     w.wcs.ctype = ["RA---AIR", "DEC--AIR"]
     w.wcs.set_pv([(2, 1, 45.0)])
-    pixcrd = np.array([[0, 0], [24, 38], [45, 98]], np.float_)
+    pixcrd = np.asarray([[0, 0], [24, 38], [45, 98]], dtype=float)
     header = w.to_header()
     header["CRVAL1P"] = 10
     header["CRVAL2P"] = 20
@@ -725,7 +725,7 @@ def test_cutout():
         assert ntpf.flux.shape[2] == 1
         ntpf = tpf.cutout(SkyCoord(tpf.ra, tpf.dec, unit="deg"), size=2)
         ntpf = tpf.cutout(size=2)
-        assert np.product(ntpf.flux.shape[1:]) == 4
+        assert np.prod(ntpf.flux.shape[1:]) == 4
         assert ntpf.targetid == tpf.targetid