DM-54825: Add template mask columns to coadd_depth_table.

python/lsst/analysis/tools/tasks/coaddDepthSummary.py

@@ -27,7 +27,7 @@
 
 
 import numpy as np
-from astropy.table import Table
+from astropy.table import Table, join
 
 from lsst.pex.config import ListField
 from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct
@@ -39,7 +39,16 @@ class CoaddDepthSummaryConnections(
     dimensions=("tract", "skymap"),
     defaultTemplates={"coaddName": ""},  # set as either deep or template in the pipeline
 ):
-    data = cT.Input(
+    mask_data = cT.Input(
+        doc="Coadd to load from the butler.",
+        name="{coaddName}_coadd.mask",
+        storageClass="MaskX",
+        multiple=True,
+        dimensions=("tract", "patch", "band", "skymap"),
+        deferLoad=True,
+    )
+
+    n_image_data = cT.Input(
         doc="Coadd n_image to load from the butler (pixel values are the number of input images).",
         name="{coaddName}_coadd_n_image",
         storageClass="ImageU",
@@ -80,24 +89,107 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs):
         butlerQC.put(outputs, outputRefs)
 
     def run(self, inputs):
-        t = Table()
-        bands = []
-        patches = []
+        # Calculate coadd metrics.
+        coadd_bands = []
+        coadd_patches = []
+        chip_gap_percents = []
+        no_data_percents = []
+        rejected_percents = []
+        inexact_psf_percents = []
+        nd_intrp_percents = []
+        nd_r_ip_percents = []
+
+        for mask_handle in inputs["mask_data"]:
+            mask = mask_handle.get()
+            data_id = mask_handle.dataId
+            band = str(data_id.band.name)
+            patch = int(data_id.patch.id)
+
+            coadd_bands.append(band)
+            coadd_patches.append(patch)
+
+            no_data = mask.getPlaneBitMask("NO_DATA")
+            rejected = mask.getPlaneBitMask("REJECTED")
+            inexact_psf = mask.getPlaneBitMask("INEXACT_PSF")
+            intrp = mask.getPlaneBitMask("INTRP")
+            mask_array = mask.array
+            tot_pixels = mask.getHeight() * mask.getWidth()
+
+            # Calculate pixels with various mask combinations.
+            no_data_flag = (mask_array & no_data) != 0
+            rejected_flag = (mask_array & rejected) != 0
+            inexact_psf_flag = (mask_array & inexact_psf) != 0
+            intrp_flag = (mask_array & intrp) != 0
+
+            chip_gap_percent = (no_data_flag & ~rejected_flag).sum() * 100 / tot_pixels
+
+            chip_gap_percents.append(chip_gap_percent)
+
+            no_data_percent = no_data_flag.sum() * 100 / tot_pixels
+
+            no_data_percents.append(no_data_percent)
+
+            rejected_percent = rejected_flag.sum() * 100 / tot_pixels
+
+            rejected_percents.append(rejected_percent)
+
+            inexact_psf_percent = inexact_psf_flag.sum() * 100 / tot_pixels
+
+            inexact_psf_percents.append(inexact_psf_percent)
+
+            nd_intrp_percent = (no_data_flag & intrp_flag).sum() * 100 / tot_pixels
+
+            nd_intrp_percents.append(nd_intrp_percent)
+
+            nd_r_ip_percent = (no_data_flag & rejected_flag & inexact_psf_flag).sum() * 100 / tot_pixels
+
+            nd_r_ip_percents.append(nd_r_ip_percent)
+
+        # Construct the Astropy table for coadd information.
+        data = [
+            coadd_patches,
+            coadd_bands,
+            chip_gap_percents,
+            no_data_percents,
+            rejected_percents,
+            inexact_psf_percents,
+            nd_intrp_percents,
+            nd_r_ip_percents,
+        ]
+        names = [
+            "patch",
+            "band",
+            "chip_gap_percent",
+            "no_data_percent",
+            "rejected_percent",
+            "inexact_psf_percent",
+            "nd_intrp_percent",
+            "nd_r_ip_percent",
+        ]
+        dtype = ["int", "str", "float", "float", "float", "float", "float", "float"]
+        coadd_table = Table(data=data, names=names, dtype=dtype)
+
+        # Calculate n_image metrics.
+        n_image_bands = []
+        n_image_patches = []
+        means = []
         medians = []
         stdevs = []
         stats = []
         quantiles = []
 
-        for n_image_handle in inputs["data"]:
+        for n_image_handle in inputs["n_image_data"]:
             n_image = n_image_handle.get()
             data_id = n_image_handle.dataId
             band = str(data_id.band.name)
             patch = int(data_id.patch.id)
+            mean = np.nanmean(n_image.array)
             median = np.nanmedian(n_image.array)
             stdev = np.nanstd(n_image.array)
 
-            bands.append(band)
-            patches.append(patch)
+            n_image_bands.append(band)
+            n_image_patches.append(patch)
+            means.append(mean)
             medians.append(median)
             stdevs.append(stdev)
 
@@ -110,8 +202,7 @@ def run(self, inputs):
 
             stats.append(band_patch_stats)
 
-            # Calculate the quantiles for image depth
-            # across the whole n_image array.
+            # Calculate the quantiles for image depth across n_image array.
             quantile = list(np.percentile(n_image.array, q=self.config.quantile_list))
             quantiles.append(quantile)
 
@@ -120,13 +211,20 @@ def run(self, inputs):
         ]
         quantile_col_names = [f"depth_{q}_percentile" for q in self.config.quantile_list]
 
-        # Construct the Astropy table
-        data = [patches, bands, medians, stdevs] + list(zip(*stats)) + list(zip(*quantiles))
-        names = ["patch", "band", "medians", "stdevs"] + threshold_col_names + quantile_col_names
+        # Construct the Astropy table for n_image information.
+        data = (
+            [n_image_patches, n_image_bands, means, medians, stdevs]
+            + list(zip(*stats))
+            + list(zip(*quantiles))
+        )
+        names = ["patch", "band", "mean", "median", "stdevs"] + threshold_col_names + quantile_col_names
         dtype = (
-            ["int", "str", "float", "float"]
+            ["int", "str", "float", "float", "float"]
             + ["float" for x in range(len(list(zip(*stats))))]
             + ["int" for y in range(len(list(zip(*quantiles))))]
         )
-        t = Table(data=data, names=names, dtype=dtype)
-        return Struct(statTable=t)
+        n_image_table = Table(data=data, names=names, dtype=dtype)
+
+        # Combine tables.
+        combined_table = join(coadd_table, n_image_table, keys=["patch", "band"])
+        return Struct(statTable=combined_table)