diff --git a/config/coaddDepthMetricsWholeSkyPlot.py b/config/coaddDepthMetricsWholeSkyPlot.py
new file mode 100644
index 000000000..eadd549ad
--- /dev/null
+++ b/config/coaddDepthMetricsWholeSkyPlot.py
@@ -0,0 +1,30 @@
+# Configuration for instances of lsst.analysis.tools.tasks.WholeSkyAnalysisTask
+# that aggregate metrics from the standard configuration of
+# lsst.analysis.tools.tasks.AssociatedSourcesTractAnalysisTask.
+
+from lsst.analysis.tools.atools import WholeSkyPlotTool
+
+# Keys with band:
+keysWithBand = [
+ "coadd_depth_summary_metrics_depth_above_threshold_1_{band}_mean",
+ "coadd_depth_summary_metrics_depth_above_threshold_3_{band}_mean",
+ "coadd_depth_summary_metrics_depth_above_threshold_5_{band}_mean",
+ "coadd_depth_summary_metrics_depth_above_threshold_12_{band}_mean",
+ ]
+
+for key in keysWithBand:
+ atoolName = key.replace("_{band}", "")
+ setattr(config.atools, atoolName, WholeSkyPlotTool)
+
+ atool = getattr(config.atools, atoolName)
+ setattr(atool, "metric", key)
+
+ plot = getattr(getattr(atool, 'produce'), 'plot')
+ setattr(plot, "showOutliers", False)
+ setattr(plot, "showNaNs", False)
+ setattr(plot, "labelTracts", True)
+ setattr(plot, "colorBarMin", 0.0)
+ setattr(plot, "colorBarMax", 100.0)
+ setattr(plot, "colorBarRange", 1.0)
+
+config.addOutputNamePrefix = True
diff --git a/doc/_static/analysis_tools/coaddDepthPlotExample.png b/doc/_static/analysis_tools/coaddDepthPlotExample.png
new file mode 100644
index 000000000..69f94112d
Binary files /dev/null and b/doc/_static/analysis_tools/coaddDepthPlotExample.png differ
diff --git a/pipelines/coaddDepthMetrics.yaml b/pipelines/coaddDepthMetrics.yaml
new file mode 100644
index 000000000..08d9436cf
--- /dev/null
+++ b/pipelines/coaddDepthMetrics.yaml
@@ -0,0 +1,35 @@
+description: |
+ Tier1 plots and metrics to assess coadd quality.
+tasks:
+ analyzeCoaddDepthCore:
+ class: lsst.analysis.tools.tasks.CoaddDepthSummaryTask
+ config:
+ connections.coaddName: template
+ coaddDepthMetricTract:
+ class: lsst.analysis.tools.tasks.CoaddDepthTableTractAnalysisTask
+ config:
+ connections.coaddName: template
+ connections.outputName: "template_coadd_depth_summary"
+ atools.coadd_depth_summary_metrics: CoaddQualityCheck
+ python: |
+ from lsst.analysis.tools.atools import CoaddQualityCheck
+ coaddDepthPlot:
+ class: lsst.analysis.tools.tasks.CoaddDepthSummaryPlotTask
+ config:
+ connections.coaddName: template
+ connections.outputName: "template"
+ atools.n_image: CoaddQualityPlot
+ python: |
+ from lsst.analysis.tools.atools import CoaddQualityPlot
+ makeCoaddDepthSummaryMetricTable:
+ class: lsst.analysis.tools.tasks.MakeMetricTableTask
+ config:
+ connections.metricBundleName: "coadd_depth_summary_metrics"
+ connections.outputTableName: "coadd_depth_summary_metrics_table"
+ inputDataDimensions: ["skymap", "tract"]
+ coaddDepthSummaryMetricTableWholeSkyPlot:
+ class: lsst.analysis.tools.tasks.WholeSkyAnalysisTask
+ config:
+ connections.inputName: "coadd_depth_summary_metrics_table"
+ connections.outputName: "template"
+ file: $ANALYSIS_TOOLS_DIR/config/coaddDepthMetricsWholeSkyPlot.py
diff --git a/python/lsst/analysis/tools/actions/plot/__init__.py b/python/lsst/analysis/tools/actions/plot/__init__.py
index ad8a7f924..2b505a4a9 100644
--- a/python/lsst/analysis/tools/actions/plot/__init__.py
+++ b/python/lsst/analysis/tools/actions/plot/__init__.py
@@ -1,5 +1,6 @@
from .barPlots import *
from .calculateRange import *
+from .coaddDepthPlot import *
from .colorColorFitPlot import *
from .completenessPlot import *
from .diaSkyPlot import *
diff --git a/python/lsst/analysis/tools/actions/plot/coaddDepthPlot.py b/python/lsst/analysis/tools/actions/plot/coaddDepthPlot.py
new file mode 100644
index 000000000..ff586cb3d
--- /dev/null
+++ b/python/lsst/analysis/tools/actions/plot/coaddDepthPlot.py
@@ -0,0 +1,198 @@
+# This file is part of analysis_tools.
+#
+# Developed for the LSST Data Management System.
+# This product includes software developed by the LSST Project
+# (https://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+
+from __future__ import annotations
+
+__all__ = ("CoaddDepthPlot",)
+
+from typing import TYPE_CHECKING, Any, Mapping
+
+import matplotlib.pyplot as plt
+from lsst.skymap.tractInfo import ExplicitTractInfo
+from lsst.utils.plotting import publication_plots, set_rubin_plotstyle
+from matplotlib.figure import Figure
+from matplotlib.lines import Line2D
+
+from ...interfaces import PlotAction, Vector
+from .plotUtils import addPlotInfo
+
+if TYPE_CHECKING:
+ from ...interfaces import KeyedData, KeyedDataSchema
+
+bands_dict = publication_plots.get_band_dicts()
+
+
+class CoaddDepthPlot(PlotAction):
+ """Make a plot of pixels per coadd depth."""
+
+ def setDefaults(self):
+ super().setDefaults()
+
+ def getInputSchema(self) -> KeyedDataSchema:
+ base: list[tuple[str, type[Vector]]] = []
+ base.append(("patch", Vector))
+ base.append(("band", Vector))
+ base.append(("depth", Vector))
+ base.append(("pixels", Vector))
+ return base
+
+ def __call__(
+ self,
+ data: KeyedData,
+ tractInfo: ExplicitTractInfo,
+ **kwargs,
+ ) -> Figure:
+ self._validateInput(data, tractInfo)
+
+ if not isinstance(tractInfo, ExplicitTractInfo):
+ raise TypeError(f"Input `tractInfo` type must be {ExplicitTractInfo} not {type(tractInfo)}.")
+
+ return self.makePlot(data, tractInfo, **kwargs)
+
+ def _validateInput(
+ self,
+ data: KeyedData,
+ tractInfo: ExplicitTractInfo,
+ ) -> None:
+ needed = set(k[0] for k in self.getInputSchema())
+ if not needed.issubset(data.keys()):
+ raise ValueError(f"Input data does not contain all required keys: {self.getInputSchema()}")
+
+ def makePlot(
+ self,
+ data: KeyedData,
+ tractInfo: ExplicitTractInfo,
+ plotInfo: Mapping[str, str] | None = None,
+ **kwargs: Any,
+ ) -> Figure:
+ """Make the plot.
+
+ Parameters
+ ----------
+ `KeyedData`
+ The catalog to plot the points from.
+ tractInfo : `~lsst.skymap.tractInfo.ExplicitTractInfo`
+ The tract info object.
+ plotInfo : `dict`
+ A dictionary of the plot information.
+
+ Returns
+ -------
+ fig : `~matplotlib.figure.Figure`
+ The resulting figure.
+
+ Examples
+ --------
+ An example coaddDepthPlot may be seen below:
+
+ .. image:: /_static/analysis_tools/coaddDepthPlotExample.png
+ """
+ set_rubin_plotstyle()
+ fig = plt.figure(dpi=300, figsize=(20, 20))
+
+ # Get number of vertical and horizontal patches in the tract.
+ num_patches_x, num_patches_y = tractInfo.getNumPatches()
+
+ max_depth = max(data["depth"])
+ max_pixels = max(data["pixels"])
+
+ plt.subplots_adjust(hspace=0, wspace=0)
+
+ patch_counter = (num_patches_x * num_patches_y) - num_patches_x # The top left corner of a tract
+ m = 0 # subplot index
+ while patch_counter >= 0:
+ for n in range(num_patches_x): # column index
+ ax = plt.subplot(num_patches_x, num_patches_y, m + 1)
+ patch_mask = data["patch"] == patch_counter
+
+ if patch_counter in data["patch"]:
+ uniqueBands = set(data["band"][patch_mask])
+ for band in uniqueBands:
+ color = bands_dict["colors"][band]
+ markerstyle = bands_dict["symbols"][band]
+ band_mask = data["band"] == band
+
+ tot_mask = (patch_mask) & (band_mask)
+
+ ax.plot(
+ data["depth"][tot_mask],
+ data["pixels"][tot_mask],
+ color=color,
+ linewidth=0,
+ ls=None,
+ marker=markerstyle,
+ ms=4,
+ alpha=0.75,
+ label=f"{band}",
+ )
+ ax.grid(alpha=0.5)
+
+ # Chart formatting
+ # Need a solution for ax.set_xscale when max_depth is high,
+ # but not all patches/bands have quality coverage.
+ ax.set_yscale("log")
+ ax.set_xlim(0, max_depth + 5)
+ ax.set_ylim(5, max_pixels)
+ # Can we somehow generalize ax.set_xticks?
+ # ax.set_xticks(np.arange(0, max_depth, 20))
+ ax.tick_params(axis="both", which="minor")
+ ax.tick_params(axis="both", which="both", top=False, right=False)
+
+ # Only label axes of the farmost left and bottom row of plots
+ if n != 0:
+ ax.set_yticklabels([])
+ ax.tick_params(axis="y", which="both", length=0)
+ if patch_counter not in range(num_patches_x):
+ ax.set_xticklabels([])
+ ax.tick_params(axis="x", which="both", length=0)
+
+ ax.set_title(f"patch {patch_counter}", y=0.85)
+ patch_counter += 1
+ m += 1
+ patch_counter -= 2 * (n + 1)
+ fig.supxlabel("Number of input visits (n_image depth)", y=0.075)
+ fig.supylabel("Count (pixels)", x=0.075)
+ legend_elements = [
+ Line2D(
+ [0], [0], color=bands_dict["colors"]["u"], lw=0, marker=bands_dict["symbols"]["u"], label="u"
+ ),
+ Line2D(
+ [0], [0], color=bands_dict["colors"]["g"], lw=0, marker=bands_dict["symbols"]["g"], label="g"
+ ),
+ Line2D(
+ [0], [0], color=bands_dict["colors"]["r"], lw=0, marker=bands_dict["symbols"]["r"], label="r"
+ ),
+ Line2D(
+ [0], [0], color=bands_dict["colors"]["i"], lw=0, marker=bands_dict["symbols"]["i"], label="i"
+ ),
+ Line2D(
+ [0], [0], color=bands_dict["colors"]["z"], lw=0, marker=bands_dict["symbols"]["z"], label="z"
+ ),
+ Line2D(
+ [0], [0], color=bands_dict["colors"]["y"], lw=0, marker=bands_dict["symbols"]["y"], label="y"
+ ),
+ ]
+ plt.legend(handles=legend_elements, loc="upper left", bbox_to_anchor=(1.05, 10))
+
+ if plotInfo is not None:
+ fig = addPlotInfo(fig, plotInfo)
+
+ return fig
diff --git a/python/lsst/analysis/tools/actions/plot/wholeSkyPlot.py b/python/lsst/analysis/tools/actions/plot/wholeSkyPlot.py
index 49d6f3fc3..1ea9957b2 100644
--- a/python/lsst/analysis/tools/actions/plot/wholeSkyPlot.py
+++ b/python/lsst/analysis/tools/actions/plot/wholeSkyPlot.py
@@ -53,11 +53,15 @@ class WholeSkyPlot(PlotAction):
xAxisLabel = Field[str](doc="Label to use for the x axis.", default="RA (degrees)")
yAxisLabel = Field[str](doc="Label to use for the y axis.", default="Dec (degrees)")
- zAxisLabel = Field[str](doc="Label to use for the z axis.", optional=False)
+ zAxisLabel = Field[str](doc="Label to use for the z axis.", default="")
autoAxesLimits = Field[bool](doc="Find axes limits automatically.", default=True)
xLimits = ListField[float](doc="Plotting limits for the x axis.", default=[-5.0, 365.0])
yLimits = ListField[float](doc="Plotting limits for the y axis.", default=[-10.0, 60.0])
+ autoAxesLimits = Field[bool](doc="Find axes limits automatically.", default=True)
+ dpi = Field[int](doc="DPI size of the figure.", default=500)
figureSize = ListField[float](doc="Size of the figure.", default=[9.0, 3.5])
+ colorBarMin = Field[float](doc="The minimum value of the color bar.", optional=True)
+ colorBarMax = Field[float](doc="The minimum value of the color bar.", optional=True)
colorBarRange = Field[float](
doc="The multiplier for the color bar range. The max/min range values are: median +/- N * sigmaMad"
", where N is this config value.",
@@ -72,6 +76,13 @@ class WholeSkyPlot(PlotAction):
doc="List of hexidecimal colors for a user-defined color map.",
optional=True,
)
+ showOutliers = Field[bool](
+ doc="Show the outliers on the plot. "
+ "Outliers are values whose absolute value is > colorBarRange * sigmaMAD.",
+ default=True,
+ )
+ showNaNs = Field[bool](doc="Show the NaNs on the plot.", default=True)
+ labelTracts = Field[bool](doc="Label the tracts.", default=False)
def getInputSchema(self, **kwargs) -> KeyedDataSchema:
base = []
@@ -238,6 +249,8 @@ def makePlot(
tracts = []
ras = []
decs = []
+ mid_ras = []
+ mid_decs = []
for i, tract in enumerate(data["tract"]):
corners = getTractCorners(skymap, tract)
patches.append(Polygon(corners, closed=True))
@@ -245,9 +258,11 @@ def makePlot(
tracts.append(tract)
ras.append(corners[0][0])
decs.append(corners[0][1])
+ mid_ras.append((corners[0][0] + corners[1][0]) / 2)
+ mid_decs.append((corners[0][1] + corners[2][1]) / 2)
# Setup figure.
- fig, ax = plt.subplots(1, 1, figsize=self.figureSize, dpi=500)
+ fig, ax = plt.subplots(1, 1, figsize=self.figureSize, dpi=self.dpi)
if self.autoAxesLimits:
xlim, ylim = self._getAxesLimits(ras, decs)
else:
@@ -265,8 +280,14 @@ def makePlot(
# Define color bar range.
med = np.nanmedian(colBarVals)
sigmaMad = nanSigmaMad(colBarVals)
- vmin = med - self.colorBarRange * sigmaMad
- vmax = med + self.colorBarRange * sigmaMad
+ if self.colorBarMin is not None:
+ vmin = np.float64(self.colorBarMin)
+ else:
+ vmin = med - self.colorBarRange * sigmaMad
+ if self.colorBarMax is not None:
+ vmax = np.float64(self.colorBarMax)
+ else:
+ vmax = med + self.colorBarRange * sigmaMad
# Note tracts with metrics outside (vmin, vmax) as outliers.
outlierInds = np.where((colBarVals < vmin) | (colBarVals > vmax))[0]
@@ -274,63 +295,81 @@ def makePlot(
# Initialize legend handles.
handles = []
- # Plot the outlier patches.
- outlierPatches = []
- if len(outlierInds) > 0:
- for ind in outlierInds:
- outlierPatches.append(patches[ind])
- outlierPatchCollection = PatchCollection(
- outlierPatches,
- cmap=colorMap,
- norm=norm,
- facecolors="none",
- edgecolors="k",
- linewidths=0.5,
- zorder=100,
- )
- ax.add_collection(outlierPatchCollection)
- # Add legend information.
- outlierPatch = Patch(
- facecolor="none",
- edgecolor="k",
- linewidth=0.5,
- label="Outlier",
- )
- handles.append(outlierPatch)
-
- # Plot tracts with NaN metric values.
- nanInds = np.where(~np.isfinite(colBarVals))[0]
- nanPatches = []
- if len(nanInds) > 0:
- for ind in nanInds:
- nanPatches.append(patches[ind])
- nanPatchCollection = PatchCollection(
- nanPatches,
- cmap=None,
- norm=norm,
- facecolors="white",
- edgecolors="grey",
- linestyles="dotted",
- linewidths=0.5,
- zorder=10,
- )
- ax.add_collection(nanPatchCollection)
- # Add legend information.
- nanPatch = Patch(
- facecolor="white",
- edgecolor="grey",
- linestyle="dotted",
- linewidth=0.5,
- label="NaN",
- )
- handles.append(nanPatch)
+ if self.showOutliers:
+ # Plot the outlier patches.
+ outlierPatches = []
+ if len(outlierInds) > 0:
+ for ind in outlierInds:
+ outlierPatches.append(patches[ind])
+ outlierPatchCollection = PatchCollection(
+ outlierPatches,
+ cmap=colorMap,
+ norm=norm,
+ facecolors="none",
+ edgecolors="k",
+ linewidths=0.5,
+ zorder=100,
+ )
+ ax.add_collection(outlierPatchCollection)
+ # Add legend information.
+ outlierPatch = Patch(
+ facecolor="none",
+ edgecolor="k",
+ linewidth=0.5,
+ label="Outlier",
+ )
+ handles.append(outlierPatch)
+
+ if self.showNaNs:
+ # Plot tracts with NaN metric values.
+ nanInds = np.where(~np.isfinite(colBarVals))[0]
+ nanPatches = []
+ if len(nanInds) > 0:
+ for ind in nanInds:
+ nanPatches.append(patches[ind])
+ nanPatchCollection = PatchCollection(
+ nanPatches,
+ cmap=None,
+ norm=norm,
+ facecolors="white",
+ edgecolors="grey",
+ linestyles="dotted",
+ linewidths=0.5,
+ zorder=100,
+ )
+ ax.add_collection(nanPatchCollection)
+ # Add legend information.
+ nanPatch = Patch(
+ facecolor="white",
+ edgecolor="grey",
+ linestyle="dotted",
+ linewidth=0.5,
+ label="NaN",
+ )
+ handles.append(nanPatch)
if len(handles) > 0:
fig.legend(handles=handles)
- # Add text boxes to show the number of tracts, number of NaNs,
- # median, sigma MAD, and the five largest outlier values.
- outlierText = self._getMaxOutlierVals(self.colorBarRange, tracts, colBarVals, outlierInds)
+ if self.labelTracts:
+ # Label the tracts
+ for i, tract in enumerate(tracts):
+ ax.text(
+ mid_ras[i],
+ mid_decs[i],
+ f"{tract}",
+ ha="center",
+ va="center",
+ fontsize=2,
+ alpha=0.7,
+ zorder=100,
+ )
+
+ if self.showOutliers:
+ # Add text boxes to show the number of tracts, number of NaNs,
+ # median, sigma MAD, and the five largest outlier values.
+ outlierText = self._getMaxOutlierVals(self.colorBarRange, tracts, colBarVals, outlierInds)
+
# Make vertical text spacing readable for different figure sizes.
multiplier = 3.5 / self.figureSize[1]
verticalSpacing = 0.028 * multiplier
@@ -342,14 +381,15 @@ def makePlot(
fontsize=8,
alpha=0.7,
)
- fig.text(
- 0.01,
- 0.01 + 2 * verticalSpacing,
- f"Num nans: {len(nanInds)}",
- transform=fig.transFigure,
- fontsize=8,
- alpha=0.7,
- )
+ if self.showNaNs:
+ fig.text(
+ 0.01,
+ 0.01 + 2 * verticalSpacing,
+ f"Num nans: {len(nanInds)}",
+ transform=fig.transFigure,
+ fontsize=8,
+ alpha=0.7,
+ )
fig.text(
0.01,
0.01 + verticalSpacing,
@@ -358,7 +398,8 @@ def makePlot(
fontsize=8,
alpha=0.7,
)
- fig.text(0.01, 0.01, outlierText, transform=fig.transFigure, fontsize=8, alpha=0.7)
+ if self.showOutliers:
+ fig.text(0.01, 0.01, outlierText, transform=fig.transFigure, fontsize=8, alpha=0.7)
# Truncate the color range to (vmin, vmax).
colorBarVals = np.clip(np.array(colBarVals), vmin, vmax)
diff --git a/python/lsst/analysis/tools/atools/coaddInputCount.py b/python/lsst/analysis/tools/atools/coaddInputCount.py
index 18bff6166..9374d598d 100644
--- a/python/lsst/analysis/tools/atools/coaddInputCount.py
+++ b/python/lsst/analysis/tools/atools/coaddInputCount.py
@@ -20,18 +20,24 @@
# along with this program. If not, see .
from __future__ import annotations
-__all__ = ("CoaddInputCount",)
+__all__ = ("CoaddInputCount", "CoaddQualityCheck", "CoaddQualityPlot")
+
+from lsst.pex.config import ListField
from ..actions.plot.calculateRange import MinMax
+from ..actions.plot.coaddDepthPlot import CoaddDepthPlot
from ..actions.plot.skyPlot import SkyPlot
-from ..actions.scalar.scalarActions import MeanAction, MedianAction, SigmaMadAction
-from ..actions.vector import CoaddPlotFlagSelector, LoadVector, SnSelector
+from ..actions.scalar.scalarActions import MeanAction, MedianAction, SigmaMadAction, StdevAction
+from ..actions.vector import BandSelector, CoaddPlotFlagSelector, DownselectVector, LoadVector, SnSelector
from ..interfaces import AnalysisTool
class CoaddInputCount(AnalysisTool):
- """skyPlot and associated metrics indicating the number
- of exposures that have gone into creating a coadd.
+ """Tract-wide metrics pertaining to how many exposures have gone into
+ a deep coadd, per band.
+
+ This AnalysisTool is designed to run on an object table, which is only
+ created for deep coadds, not template coadds.
"""
def setDefaults(self):
@@ -68,7 +74,7 @@ def setDefaults(self):
self.process.calculateActions.sigmaMad = SigmaMadAction()
self.process.calculateActions.sigmaMad.vectorKey = "z"
- # SkyPlot of number of contributing exposures across coadd:
+ # SkyPlot of number of contributing exposures in coad, per tract/band:
self.produce.plot = SkyPlot()
self.produce.plot.plotTypes = ["any"]
self.produce.plot.plotName = "{band}_inputCount"
@@ -79,10 +85,131 @@ def setDefaults(self):
self.produce.plot.showExtremeOutliers = False
self.produce.plot.colorbarRange = MinMax()
- # Summary metrics for the whole coadd.
+ # Summary metrics for the whole coadd, per tract/band.
self.produce.metric.units = {"median": "ct", "sigmaMad": "ct", "mean": "ct"}
self.produce.metric.newNames = {
"median": "{band}_inputCount_median",
"mean": "{band}_inputCount_mean",
"sigmaMad": "{band}_inputCount_sigmaMad",
}
+
+
+class CoaddQualityCheck(AnalysisTool):
+ """Compute the percentage of each coadd that has a number of input
+ exposures exceeding a threshold.
+
+ This AnalysisTool is designed to run on any coadd, provided a
+ coadd_depth_table is created first (via CoaddDepthSummaryAnalysisTask).
+
+ For example, if exactly half of a coadd patch contains 15 overlapping
+ constituent visits and half contains fewer, the value computed for
+ `depth_above_threshold_12` would be 50.
+
+ These values come from the n_image data product, which is an image
+ identical to the coadd but with pixel values of the number of input
+ images instead of flux or counts. n_images are persisted during
+ coadd assembly.
+ """
+
+ band_list = ListField(
+ default=["u", "g", "r", "i", "z", "y"],
+ dtype=str,
+ doc="Bands for colors.",
+ )
+
+ threshold_list = ListField(
+ default=[1, 3, 5, 12],
+ dtype=int,
+ doc="The n_image pixel value thresholds.",
+ )
+
+ quantile_list = ListField(
+ default=[5, 10, 25, 50, 75, 90, 95],
+ dtype=int,
+ doc="The percentiles at which to compute n_image values, in ascending order.",
+ )
+
+ def setDefaults(self):
+ super().setDefaults()
+
+ self.process.buildActions.patch = LoadVector(vectorKey="patch")
+ self.process.buildActions.band = LoadVector(vectorKey="band")
+
+ def finalize(self):
+ for threshold in self.threshold_list:
+ # This gives a RuntimeWarning whenever the tract doesn't have
+ # a particular band. Need to use a band list derived from the
+ # bands found in the "band" column of a given tract.
+ for band in self.band_list:
+ name = f"depth_above_threshold_{threshold}"
+ setattr(self.process.buildActions, name, LoadVector(vectorKey=name))
+ setattr(
+ self.process.filterActions,
+ f"{name}_{band}",
+ DownselectVector(vectorKey=name, selector=BandSelector(bands=[band])),
+ )
+ setattr(
+ self.process.calculateActions,
+ f"{name}_{band}_median",
+ MedianAction(vectorKey=f"{name}_{band}"),
+ )
+ setattr(
+ self.process.calculateActions,
+ f"{name}_{band}_mean",
+ MeanAction(vectorKey=f"{name}_{band}"),
+ )
+ setattr(
+ self.process.calculateActions,
+ f"{name}_{band}_stdev",
+ StdevAction(vectorKey=f"{name}_{band}"),
+ )
+
+ # The units for the quantity are dimensionless (percentage)
+ self.produce.metric.units[f"{name}_{band}_median"] = ""
+ self.produce.metric.units[f"{name}_{band}_mean"] = ""
+ self.produce.metric.units[f"{name}_{band}_stdev"] = ""
+
+ for quantile in self.quantile_list:
+ for band in self.band_list:
+ name = f"depth_{quantile}_percentile"
+ setattr(self.process.buildActions, name, LoadVector(vectorKey=name))
+ setattr(
+ self.process.filterActions,
+ f"{name}_{band}",
+ DownselectVector(vectorKey=name, selector=BandSelector(bands=[band])),
+ )
+ setattr(
+ self.process.calculateActions,
+ f"{name}_{band}_median",
+ MedianAction(vectorKey=f"{name}_{band}"),
+ )
+ setattr(
+ self.process.calculateActions,
+ f"{name}_{band}_mean",
+ MeanAction(vectorKey=f"{name}_{band}"),
+ )
+ setattr(
+ self.process.calculateActions,
+ f"{name}_{band}_stdev",
+ StdevAction(vectorKey=f"{name}_{band}"),
+ )
+
+ # The units for the quantity are dimensionless (percentage)
+ self.produce.metric.units[f"{name}_{band}_median"] = ""
+ self.produce.metric.units[f"{name}_{band}_mean"] = ""
+ self.produce.metric.units[f"{name}_{band}_stdev"] = ""
+
+
+class CoaddQualityPlot(AnalysisTool):
+ """Make a plot of coadd depth."""
+
+ parameterizedBand: bool = False
+
+ def setDefaults(self):
+ super().setDefaults()
+ self.process.buildActions.patch = LoadVector(vectorKey="patch")
+ self.process.buildActions.band = LoadVector(vectorKey="band")
+ self.process.buildActions.depth = LoadVector(vectorKey="depth")
+ self.process.buildActions.pixels = LoadVector(vectorKey="pixels")
+
+ self.produce.plot = CoaddDepthPlot()
diff --git a/python/lsst/analysis/tools/atools/wholeSkyPlotTool.py b/python/lsst/analysis/tools/atools/wholeSkyPlotTool.py
index fb21ba18a..fabbcee1a 100644
--- a/python/lsst/analysis/tools/atools/wholeSkyPlotTool.py
+++ b/python/lsst/analysis/tools/atools/wholeSkyPlotTool.py
@@ -53,13 +53,13 @@ def setDefaults(self):
self.process.buildActions.tract = LoadVector()
self.process.buildActions.tract.vectorKey = "tract"
+ self.produce.plot = WholeSkyPlot()
def finalize(self):
-
self.process.buildActions.z = LoadVector(vectorKey=self.metric)
self.process.buildActions.zUnit = KeyedDataUnitAccessAction(key=self.metric)
- self.produce.plot = WholeSkyPlot(zAxisLabel=self.metric)
- sequentialMetrics = ["count", "num", "igma", "tdev", "Repeat"]
+ self.produce.plot.zAxisLabel = self.metric
+ sequentialMetrics = ["count", "ean", "edian", "num", "igma", "tdev", "Repeat"]
if any(sequentialMetric in self.metric for sequentialMetric in sequentialMetrics):
self.produce.plot.colorMapType = "sequential"
super().finalize()
diff --git a/python/lsst/analysis/tools/tasks/__init__.py b/python/lsst/analysis/tools/tasks/__init__.py
index 9f5b18305..fc98d4deb 100644
--- a/python/lsst/analysis/tools/tasks/__init__.py
+++ b/python/lsst/analysis/tools/tasks/__init__.py
@@ -6,6 +6,9 @@
from .calibrationAnalysis import *
from .catalogMatch import *
from .ccdVisitTableAnalysis import *
+from .coaddDepthSummary import *
+from .coaddDepthSummaryPlot import *
+from .coaddDepthTableTractAnalysis import *
from .diaFakesDetectorVisitAnalysis import *
from .diaFakesVisitAnalysis import *
from .diaObjectDetectorVisitAnalysis import *
diff --git a/python/lsst/analysis/tools/tasks/coaddDepthSummary.py b/python/lsst/analysis/tools/tasks/coaddDepthSummary.py
new file mode 100644
index 000000000..45b93624d
--- /dev/null
+++ b/python/lsst/analysis/tools/tasks/coaddDepthSummary.py
@@ -0,0 +1,131 @@
+# This file is part of analysis_tools.
+#
+# Developed for the LSST Data Management System.
+# This product includes software developed by the LSST Project
+# (https://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+from __future__ import annotations
+
+__all__ = (
+ "CoaddDepthSummaryConfig",
+ "CoaddDepthSummaryTask",
+)
+
+
+import numpy as np
+from astropy.table import Table
+from lsst.pex.config import ListField
+from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct
+from lsst.pipe.base import connectionTypes as cT
+
+
+class CoaddDepthSummaryConnections(
+ PipelineTaskConnections,
+ dimensions=("tract", "skymap"),
+ defaultTemplates={"coaddName": ""}, # set as either deep or template in the pipeline
+):
+ 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",
+ multiple=True,
+ dimensions=("tract", "patch", "band", "skymap"),
+ deferLoad=True,
+ )
+
+ statTable = cT.Output(
+ doc="Table with resulting n_image based depth statistics.",
+ name="{coaddName}_coadd_depth_table",
+ storageClass="ArrowAstropy",
+ dimensions=("tract", "skymap"),
+ )
+
+
+class CoaddDepthSummaryConfig(PipelineTaskConfig, pipelineConnections=CoaddDepthSummaryConnections):
+ threshold_list = ListField(
+ default=[1, 3, 5, 12],
+ dtype=int,
+ doc="The n_image pixel value thresholds, in ascending order.",
+ )
+
+ quantile_list = ListField(
+ default=[5, 10, 25, 50, 75, 90, 95],
+ dtype=int,
+ doc="The percentiles at which to compute n_image values, in ascending order.",
+ )
+
+
+class CoaddDepthSummaryTask(PipelineTask):
+ ConfigClass = CoaddDepthSummaryConfig
+ _DefaultName = "coaddDepthSummary"
+
+ def runQuantum(self, butlerQC, inputRefs, outputRefs):
+ inputs = butlerQC.get(inputRefs)
+ outputs = self.run(inputs)
+ butlerQC.put(outputs, outputRefs)
+
+ def run(self, inputs):
+ t = Table()
+ bands = []
+ patches = []
+ medians = []
+ stdevs = []
+ stats = []
+ quantiles = []
+
+ for n_image_handle in inputs["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)
+ median = np.nanmedian(n_image.array)
+ stdev = np.nanstd(n_image.array)
+
+ bands.append(band)
+ patches.append(patch)
+ medians.append(median)
+ stdevs.append(stdev)
+
+ band_patch_stats = []
+ for threshold in self.config.threshold_list:
+ # Calculate the percentage of the image with an image depth
+ # above the given threshold.
+ stat = np.sum(n_image.array > threshold) * 100 / (n_image.getHeight() * n_image.getWidth())
+ band_patch_stats.append(stat)
+
+ stats.append(band_patch_stats)
+
+ # Calculate the quantiles for image depth
+ # across the whole n_image array.
+ quantile = list(np.percentile(n_image.array, q=self.config.quantile_list))
+ quantiles.append(quantile)
+
+ threshold_col_names = [
+ f"depth_above_threshold_{threshold}" for threshold in self.config.threshold_list
+ ]
+ 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
+ dtype = (
+ ["int", "str", "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)
diff --git a/python/lsst/analysis/tools/tasks/coaddDepthSummaryPlot.py b/python/lsst/analysis/tools/tasks/coaddDepthSummaryPlot.py
new file mode 100644
index 000000000..9c614dbc9
--- /dev/null
+++ b/python/lsst/analysis/tools/tasks/coaddDepthSummaryPlot.py
@@ -0,0 +1,103 @@
+# This file is part of analysis_tools.
+#
+# Developed for the LSST Data Management System.
+# This product includes software developed by the LSST Project
+# (https://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+from __future__ import annotations
+
+__all__ = (
+ "CoaddDepthSummaryPlotConfig",
+ "CoaddDepthSummaryPlotTask",
+)
+
+import numpy as np
+from lsst.pipe.base import Struct
+from lsst.pipe.base import connectionTypes as cT
+from lsst.skymap import BaseSkyMap
+
+from ..interfaces import AnalysisBaseConfig, AnalysisBaseConnections, AnalysisPipelineTask, KeyedData
+
+
+class CoaddDepthSummaryPlotConnections(
+ AnalysisBaseConnections,
+ dimensions=("tract", "skymap"),
+ defaultTemplates={
+ "coaddName": "",
+ },
+):
+ skymap = cT.Input(
+ doc="The skymap that covers the tract that the data is from.",
+ name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
+ storageClass="SkyMap",
+ dimensions=("skymap",),
+ )
+
+ 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",
+ multiple=True,
+ dimensions=("tract", "patch", "band", "skymap"),
+ deferLoad=True,
+ )
+
+
+class CoaddDepthSummaryPlotConfig(AnalysisBaseConfig, pipelineConnections=CoaddDepthSummaryPlotConnections):
+ pass
+
+
+class CoaddDepthSummaryPlotTask(AnalysisPipelineTask):
+ ConfigClass = CoaddDepthSummaryPlotConfig
+ _DefaultName = "coaddDepthSummaryPlot"
+
+ def runQuantum(self, butlerQC, inputRefs, outputRefs):
+ inputs = butlerQC.get(inputRefs)
+ skymap = inputs["skymap"]
+ dataId = butlerQC.quantum.dataId
+ tractInfo = skymap[dataId["tract"]]
+ outputs = self.run(data={"n_image_data": inputs["n_image_data"]}, tractInfo=tractInfo)
+ butlerQC.put(outputs, outputRefs)
+
+ def run(self, *, data: KeyedData | None = None, **kwargs) -> Struct:
+ """Use n_images to make a plot illustrating coadd depth."""
+ bands = []
+ patches = []
+
+ depths = []
+ pixels = []
+
+ for n_image_handle in data["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)
+
+ depth, pixel = np.unique(n_image.array, return_counts=True)
+
+ depths.extend(depth)
+ pixels.extend(pixel)
+
+ for i in range(len(depth)):
+ bands.append(band)
+ patches.append(patch)
+
+ pixel_data = {"patch": patches, "band": bands, "depth": depths, "pixels": pixels}
+
+ outputs = super().run(data=pixel_data, **kwargs) # this creates a struct for the output
+
+ return outputs
diff --git a/python/lsst/analysis/tools/tasks/coaddDepthTableTractAnalysis.py b/python/lsst/analysis/tools/tasks/coaddDepthTableTractAnalysis.py
new file mode 100644
index 000000000..9379b4666
--- /dev/null
+++ b/python/lsst/analysis/tools/tasks/coaddDepthTableTractAnalysis.py
@@ -0,0 +1,57 @@
+# This file is part of analysis_tools.
+#
+# Developed for the LSST Data Management System.
+# This product includes software developed by the LSST Project
+# (https://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+from __future__ import annotations
+
+__all__ = (
+ "CoaddDepthTableTractAnalysisConnections",
+ "CoaddDepthTableTractAnalysisConfig",
+ "CoaddDepthTableTractAnalysisTask",
+)
+
+import lsst.pex.config as pexConfig
+from lsst.pipe.base import connectionTypes as cT
+
+from ..interfaces import AnalysisBaseConfig, AnalysisBaseConnections, AnalysisPipelineTask
+
+
+class CoaddDepthTableTractAnalysisConnections(
+ AnalysisBaseConnections,
+ dimensions=("tract", "skymap"),
+ defaultTemplates={"coaddName": ""},
+):
+ data = cT.Input(
+ doc="Table with coadd depth statistics based on n_image values.",
+ name="{coaddName}_coadd_depth_table",
+ storageClass="ArrowAstropy",
+ dimensions=("tract", "skymap"),
+ deferLoad=True,
+ )
+
+
+class CoaddDepthTableTractAnalysisConfig(
+ AnalysisBaseConfig, pipelineConnections=CoaddDepthTableTractAnalysisConnections
+):
+ load_skymap = pexConfig.Field[bool](doc="Whether to load the skymap.", default=True)
+
+
+class CoaddDepthTableTractAnalysisTask(AnalysisPipelineTask):
+ ConfigClass = CoaddDepthTableTractAnalysisConfig
+ _DefaultName = "coaddDepthTableTractAnalysis"