Coverage for python / lsst / analysis / tools / tasks / limitingSurfaceBrightnessAnalysis.py: 0%
79 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-05 18:53 +0000
1# This file is part of analysis_tools.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22__all__ = (
23 "LimitingSurfaceBrightnessTask",
24 "LimitingSurfaceBrightnessConfig",
25 "LimitingSurfaceBrightnessConnections",
26 "LimitingSurfaceBrightnessAnalysisConnections",
27 "LimitingSurfaceBrightnessAnalysisConfig",
28 "LimitingSurfaceBrightnessAnalysisTask",
29)
31import logging
33import lsst.afw.math as afwMath
34import numpy as np
35from astropy.table import Table
36from lsst.pex.config import Field, ListField
37from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct
38from lsst.pipe.base.connectionTypes import Input, Output
39from lsst.skymap import BaseSkyMap
41from ..interfaces import AnalysisBaseConfig, AnalysisBaseConnections, AnalysisPipelineTask
43log = logging.getLogger(__name__)
46class LimitingSurfaceBrightnessConnections(
47 PipelineTaskConnections,
48 dimensions=(),
49 defaultTemplates={
50 "detectionTableName": "object_all",
51 "photoCalibName": "deep_coadd.photoCalib",
52 "wcsName": "deep_coadd.wcs",
53 },
54):
55 """Connections class for LimitingSurfaceBrightnessTask."""
57 data = Input(
58 doc="Visit- or coadd-level object table",
59 name="{detectionTableName}",
60 storageClass="ArrowAstropy",
61 multiple=True,
62 dimensions=(),
63 deferLoad=True,
64 )
66 skymap = Input(
67 doc="The skymap that covers the originating data's tract.",
68 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
69 storageClass="SkyMap",
70 dimensions=("skymap",),
71 )
73 photoCalib = Input(
74 name="{photoCalibName}",
75 storageClass="PhotoCalib",
76 doc="Photometric calibration associated with the originating data image.",
77 multiple=True,
78 dimensions=(),
79 deferLoad=True,
80 )
82 wcs = Input(
83 name="{wcsName}",
84 storageClass="Wcs",
85 doc="WCS header associated with the originating data image.",
86 multiple=True,
87 dimensions=(),
88 deferLoad=True,
89 )
91 limiting_surface_brightness_table = Output(
92 name="limiting_surface_brightness_table",
93 storageClass="ArrowAstropy",
94 doc="A table containing two columns: the detector or patch IDs and the values of limiting surface "
95 "brightness derived for those detectors or patches.",
96 dimensions=(),
97 )
99 def __init__(self, *, config=None):
100 super().__init__(config=config)
101 # Update output table name for configurable dimensions
102 dimen = "_visit" if "visit" in config.inputTableDimensions else "_tract"
103 self.data = Input(
104 doc=self.data.doc,
105 name=self.data.name,
106 storageClass=self.data.storageClass,
107 deferLoad=self.data.deferLoad,
108 dimensions=frozenset(sorted(config.inputTableDimensions)),
109 multiple=self.data.multiple,
110 )
111 self.photoCalib = Input(
112 doc=self.photoCalib.doc,
113 name=self.photoCalib.name,
114 storageClass=self.photoCalib.storageClass,
115 deferLoad=self.photoCalib.deferLoad,
116 dimensions=frozenset(sorted(config.inputCalibDimensions)),
117 multiple=self.photoCalib.multiple,
118 )
119 self.wcs = Input(
120 doc=self.wcs.doc,
121 name=self.wcs.name,
122 storageClass=self.wcs.storageClass,
123 deferLoad=self.wcs.deferLoad,
124 dimensions=frozenset(sorted(config.inputCalibDimensions)),
125 multiple=self.wcs.multiple,
126 )
127 self.limiting_surface_brightness_table = Output(
128 doc=self.limiting_surface_brightness_table.doc,
129 name=self.limiting_surface_brightness_table.name + dimen,
130 storageClass=self.limiting_surface_brightness_table.storageClass,
131 dimensions=frozenset(sorted(config.outputDataDimensions)),
132 )
134 assert config is not None, "Missing required config object."
136 if "tract" not in config.inputTableDimensions:
137 del self.skymap
139 self.dimensions.update(frozenset(sorted(config.outputDataDimensions)))
142class LimitingSurfaceBrightnessConfig(
143 PipelineTaskConfig, pipelineConnections=LimitingSurfaceBrightnessConnections
144):
145 """Config class for LimitingSurfaceBrightnessTask."""
147 inputTableDimensions = ListField[str](
148 doc="Dimensions of the input object table data.",
149 default=("skymap", "tract"),
150 optional=False,
151 )
152 inputCalibDimensions = ListField[str](
153 doc="Dimensions of the input calibration data.",
154 default=("tract", "band"),
155 optional=False,
156 )
157 outputDataDimensions = ListField[str](
158 doc="Dimensions of the output table data.",
159 default=("tract", "band"),
160 optional=False,
161 )
162 apertureSize = Field[int](
163 doc="The size of the sky objects photometry aperture.",
164 default=9,
165 )
168class LimitingSurfaceBrightnessTask(PipelineTask):
169 """A task for measuring the 3 sigma limiting surface brightness on 10
170 arcsecond scales for a given image. This is currently a widely accepted
171 metric for depth within the low surface brightness community.
173 Reference: Roman, J., Trujillo, I., & Montes, M., 2020, A & A, 644, 42
174 """
176 ConfigClass = LimitingSurfaceBrightnessConfig
177 _DefaultName = "limitingSurfaceBrightness"
179 def __init__(self, initInputs=None, *args, **kwargs):
180 super().__init__(*args, **kwargs)
182 def runQuantum(self, butlerQC, inputRefs, outputRefs):
183 inputs = butlerQC.get(inputRefs)
184 limiting_surface_brightness_struct = self.run(**{k: v for k, v in inputs.items() if k != "skymap"})
185 butlerQC.put(limiting_surface_brightness_struct, outputRefs)
187 def run(self, data, photoCalib, wcs):
188 """Output a table of per-detector or per-patch limiting surface
189 brightnesses for the input visit or tract coadd
191 Parameters
192 ----------
193 data : `list`
194 List of dicts with information respecting the extracted source
195 detection catalogues
196 photoCalib : `list`
197 List of dicts with information respecting the extracted image
198 photometric calibrations
199 wcs : `list`
200 List of dicts with information respecting the extracted image
201 world coordinate systems
203 Returns
204 -------
205 `pipe.base.Struct` containing `astropy.table.Table`
206 Table containing per-detector or per-patch limiting surface
207 brightness values
208 """
209 # Generate lookup tables
210 lookup_photoCalib = {x.dataId: x for x in photoCalib}
211 lookup_wcs = {x.dataId: x for x in wcs}
213 # Retrieve the source catalogue and photometric band if required
214 source_catalogue = data[0].get()
215 if "band" in self.config.inputCalibDimensions:
216 band = photoCalib[0].dataId["band"] + "_"
217 else:
218 band = ""
220 # Set up the global table
221 limiting_surface_brightness_table = Table(
222 names=["imageID", "limiting_surface_brightness"],
223 dtype=[np.int64, np.float64],
224 )
226 # Loop over the detector or patch calibration data
227 for dataId in lookup_photoCalib.keys():
228 # And the photometric calibration
229 instFluxToMagnitude = lookup_photoCalib[dataId].get().instFluxToMagnitude(1)
230 # And the pixel scale
231 pxScale = lookup_wcs[dataId].get().getPixelScale().asArcseconds()
233 # Isolate the sky sources per image, 9px radius apertures
234 idKey = "detector" if "detector" in dataId else "patch"
235 skyKey = "sky_source" if "detector" in dataId else "sky_object"
236 isImage = source_catalogue[idKey] == dataId[idKey]
237 isSky = source_catalogue[skyKey] > 0
238 skySources = source_catalogue[isImage & isSky][band + "ap%02dFlux" % (self.config.apertureSize)]
240 # Derive the clipped standard deviation of sky sources in nJy
241 nPix = np.pi * self.config.apertureSize**2 # Number of pixels within the circular aperture
242 ctrl = afwMath.StatisticsControl(3, 3)
243 ctrl.setNanSafe(True)
244 statistic = afwMath.stringToStatisticsProperty("STDEVCLIP")
245 sigSkySources = afwMath.makeStatistics(skySources / nPix, statistic, ctrl).getValue()
247 # Derive limiting surface brightness. 3sigma, on 10"x10" scales
248 pixScaleRatio = np.sqrt(pxScale**2 / (nPix * pxScale**2))
249 sigma = sigSkySources / pixScaleRatio
250 muLim = -2.5 * np.log10((3 * sigma) / (pxScale * 10)) + instFluxToMagnitude
252 # Append a new row to the table
253 limiting_surface_brightness_table.add_row([dataId[idKey], muLim])
255 # Return results.
256 return Struct(limiting_surface_brightness_table=limiting_surface_brightness_table)
259class LimitingSurfaceBrightnessAnalysisConnections(
260 AnalysisBaseConnections,
261 dimensions=(),
262 defaultTemplates={"outputName": "limitingSurfaceBrightness"},
263):
264 data = Input(
265 name="limiting_surface_brightness_table",
266 storageClass="ArrowAstropy",
267 doc="A table containing two columns: the detector or patch IDs and the values of limiting surface "
268 "brightness derived for those detectors or patches.",
269 deferLoad=True,
270 dimensions=(),
271 )
273 def __init__(self, *, config=None):
274 super().__init__(config=config)
275 # Update input table name for configurable dimensions
276 dimen = "_visit" if "visit" in config.inputDataDimensions else "_tract"
277 self.data = Input(
278 name=self.data.name + dimen,
279 storageClass=self.data.storageClass,
280 doc=self.data.doc,
281 deferLoad=self.data.deferLoad,
282 dimensions=frozenset(sorted(config.inputDataDimensions)),
283 )
286class LimitingSurfaceBrightnessAnalysisConfig(
287 AnalysisBaseConfig, pipelineConnections=LimitingSurfaceBrightnessAnalysisConnections
288):
289 inputDataDimensions = ListField[str](
290 doc="Dimensions of the input table data.",
291 default=(),
292 optional=False,
293 )
296class LimitingSurfaceBrightnessAnalysisTask(AnalysisPipelineTask):
297 ConfigClass = LimitingSurfaceBrightnessAnalysisConfig
298 _DefaultName = "limitingSurfaceBrightnessAnalysis"