Coverage for python / lsst / analysis / tools / tasks / limitingSurfaceBrightnessAnalysis.py: 0%
84 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-23 09:00 +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 numpy as np
34from astropy.table import Table
36from lsst.afw.math import StatisticsControl, makeStatistics, stringToStatisticsProperty
37from lsst.pex.config import Field, ListField
38from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct
39from lsst.pipe.base.connectionTypes import Input, Output
40from lsst.skymap import BaseSkyMap
42from ..interfaces import AnalysisBaseConfig, AnalysisBaseConnections, AnalysisPipelineTask
44log = logging.getLogger(__name__)
47class LimitingSurfaceBrightnessConnections(
48 PipelineTaskConnections,
49 dimensions=(),
50 defaultTemplates={
51 "detectionTableName": "object_all",
52 "photoCalibName": "deep_coadd.photoCalib",
53 "wcsName": "deep_coadd.wcs",
54 },
55):
56 """Connections class for LimitingSurfaceBrightnessTask."""
58 data = Input(
59 doc="Visit- or coadd-level object table",
60 name="{detectionTableName}",
61 storageClass="ArrowAstropy",
62 multiple=True,
63 dimensions=(),
64 deferLoad=True,
65 )
67 skymap = Input(
68 doc="The skymap that covers the originating data's tract.",
69 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
70 storageClass="SkyMap",
71 dimensions=("skymap",),
72 )
74 photoCalib = Input(
75 name="{photoCalibName}",
76 storageClass="PhotoCalib",
77 doc="Photometric calibration associated with the originating data image.",
78 multiple=True,
79 dimensions=(),
80 deferLoad=True,
81 )
83 wcs = Input(
84 name="{wcsName}",
85 storageClass="Wcs",
86 doc="WCS header associated with the originating data image.",
87 multiple=True,
88 dimensions=(),
89 deferLoad=True,
90 )
92 limiting_surface_brightness_table = Output(
93 name="limiting_surface_brightness_table",
94 storageClass="ArrowAstropy",
95 doc="A table containing two columns: the detector or patch IDs and the values of limiting surface "
96 "brightness derived for those detectors or patches.",
97 dimensions=(),
98 )
100 def __init__(self, *, config=None):
101 super().__init__(config=config)
102 # Update output table name for configurable dimensions
103 dimen = "_visit" if "visit" in config.inputTableDimensions else "_tract"
104 self.data = Input(
105 doc=self.data.doc,
106 name=self.data.name,
107 storageClass=self.data.storageClass,
108 deferLoad=self.data.deferLoad,
109 dimensions=frozenset(sorted(config.inputTableDimensions)),
110 multiple=self.data.multiple,
111 )
112 self.photoCalib = Input(
113 doc=self.photoCalib.doc,
114 name=self.photoCalib.name,
115 storageClass=self.photoCalib.storageClass,
116 deferLoad=self.photoCalib.deferLoad,
117 dimensions=frozenset(sorted(config.inputCalibDimensions)),
118 multiple=self.photoCalib.multiple,
119 )
120 self.wcs = Input(
121 doc=self.wcs.doc,
122 name=self.wcs.name,
123 storageClass=self.wcs.storageClass,
124 deferLoad=self.wcs.deferLoad,
125 dimensions=frozenset(sorted(config.inputCalibDimensions)),
126 multiple=self.wcs.multiple,
127 )
128 self.limiting_surface_brightness_table = Output(
129 doc=self.limiting_surface_brightness_table.doc,
130 name=self.limiting_surface_brightness_table.name + dimen,
131 storageClass=self.limiting_surface_brightness_table.storageClass,
132 dimensions=frozenset(sorted(config.outputDataDimensions)),
133 )
135 assert config is not None, "Missing required config object."
137 if "tract" not in config.inputTableDimensions:
138 del self.skymap
140 self.dimensions.update(frozenset(sorted(config.outputDataDimensions)))
143class LimitingSurfaceBrightnessConfig(
144 PipelineTaskConfig, pipelineConnections=LimitingSurfaceBrightnessConnections
145):
146 """Config class for LimitingSurfaceBrightnessTask."""
148 inputTableDimensions = ListField[str](
149 doc="Dimensions of the input object table data.",
150 default=("skymap", "tract"),
151 optional=False,
152 )
153 inputCalibDimensions = ListField[str](
154 doc="Dimensions of the input calibration data.",
155 default=("tract", "band"),
156 optional=False,
157 )
158 outputDataDimensions = ListField[str](
159 doc="Dimensions of the output table data.",
160 default=("tract", "band"),
161 optional=False,
162 )
163 apertureSize = Field[int](
164 doc="The size of the sky objects photometry aperture.",
165 default=9,
166 )
169class LimitingSurfaceBrightnessTask(PipelineTask):
170 """A task for measuring the 3 sigma limiting surface brightness on 10
171 arcsecond scales for a given image. This is currently a widely accepted
172 metric for depth within the low surface brightness community.
174 Reference: Roman, J., Trujillo, I., & Montes, M., 2020, A & A, 644, 42
175 """
177 ConfigClass = LimitingSurfaceBrightnessConfig
178 _DefaultName = "limitingSurfaceBrightness"
180 def __init__(self, initInputs=None, *args, **kwargs):
181 super().__init__(*args, **kwargs)
183 def runQuantum(self, butlerQC, inputRefs, outputRefs):
184 inputs = butlerQC.get(inputRefs)
185 limiting_surface_brightness_struct = self.run(**{k: v for k, v in inputs.items() if k != "skymap"})
186 butlerQC.put(limiting_surface_brightness_struct, outputRefs)
188 def run(self, data, photoCalib, wcs):
189 """Output a table of per-detector or per-patch limiting surface
190 brightnesses for the input visit or tract coadd
192 Parameters
193 ----------
194 data : `list`
195 List of dicts with information respecting the extracted source
196 detection catalogues
197 photoCalib : `list`
198 List of dicts with information respecting the extracted image
199 photometric calibrations
200 wcs : `list`
201 List of dicts with information respecting the extracted image
202 world coordinate systems
204 Returns
205 -------
206 `pipe.base.Struct` containing `astropy.table.Table`
207 Table containing per-detector or per-patch limiting surface
208 brightness values
209 """
210 # Generate lookup tables
211 lookup_photoCalib = {x.dataId: x for x in photoCalib}
212 lookup_wcs = {x.dataId: x for x in wcs}
214 # Retrieve the source catalogue and photometric band if required
215 source_catalogue = data[0].get()
216 if "band" in self.config.inputCalibDimensions:
217 band = photoCalib[0].dataId["band"] + "_"
218 else:
219 band = ""
221 # Set up the global table
222 limiting_surface_brightness_table = Table(
223 names=["imageID", "limiting_surface_brightness"],
224 dtype=[np.int64, np.float64],
225 )
227 # Loop over the detector or patch calibration data
228 for dataId in lookup_photoCalib.keys():
229 # And the photometric calibration
230 instFluxToMagnitude = lookup_photoCalib[dataId].get().instFluxToMagnitude(1)
231 # And the pixel scale
232 pxScale = lookup_wcs[dataId].get().getPixelScale().asArcseconds()
234 # Isolate the sky sources per image, 9px radius apertures
235 idKey = "detector" if "detector" in dataId else "patch"
236 skyKey = "sky_source" if "detector" in dataId else "sky_object"
237 isImage = source_catalogue[idKey] == dataId[idKey]
238 isSky = source_catalogue[skyKey] > 0
239 skySources = source_catalogue[isImage & isSky][f"{band}ap{self.config.apertureSize:02d}Flux"]
240 # Some patches contain no detections
241 if len(skySources) == 0:
242 muLim = np.nan
243 limiting_surface_brightness_table.add_row([dataId[idKey], muLim])
244 continue
246 # Derive the clipped standard deviation of sky sources in nJy
247 nPix = np.pi * self.config.apertureSize**2 # Number of pixels within the circular aperture
248 ctrl = StatisticsControl(3, 3)
249 ctrl.setNanSafe(True)
250 statistic = stringToStatisticsProperty("STDEVCLIP")
251 sigSkySources = makeStatistics(skySources / nPix, statistic, ctrl).getValue()
253 # Derive limiting surface brightness. 3sigma, on 10"x10" scales
254 pixScaleRatio = np.sqrt(pxScale**2 / (nPix * pxScale**2))
255 sigma = sigSkySources / pixScaleRatio
256 muLim = -2.5 * np.log10((3 * sigma) / (pxScale * 10)) + instFluxToMagnitude
258 # Append a new row to the table
259 limiting_surface_brightness_table.add_row([dataId[idKey], muLim])
261 # Return results.
262 return Struct(limiting_surface_brightness_table=limiting_surface_brightness_table)
265class LimitingSurfaceBrightnessAnalysisConnections(
266 AnalysisBaseConnections,
267 dimensions=(),
268 defaultTemplates={"outputName": "limitingSurfaceBrightness"},
269):
270 data = Input(
271 name="limiting_surface_brightness_table",
272 storageClass="ArrowAstropy",
273 doc="A table containing two columns: the detector or patch IDs and the values of limiting surface "
274 "brightness derived for those detectors or patches.",
275 deferLoad=True,
276 dimensions=(),
277 )
279 def __init__(self, *, config=None):
280 super().__init__(config=config)
281 # Update input table name for configurable dimensions
282 dimen = "_visit" if "visit" in config.inputDataDimensions else "_tract"
283 self.data = Input(
284 name=self.data.name + dimen,
285 storageClass=self.data.storageClass,
286 doc=self.data.doc,
287 deferLoad=self.data.deferLoad,
288 dimensions=frozenset(sorted(config.inputDataDimensions)),
289 )
290 self.dimensions.update(frozenset(sorted(config.inputDataDimensions)))
293class LimitingSurfaceBrightnessAnalysisConfig(
294 AnalysisBaseConfig, pipelineConnections=LimitingSurfaceBrightnessAnalysisConnections
295):
296 inputDataDimensions = ListField[str](
297 doc="Dimensions of the input table data.",
298 default=(),
299 optional=False,
300 )
303class LimitingSurfaceBrightnessAnalysisTask(AnalysisPipelineTask):
304 ConfigClass = LimitingSurfaceBrightnessAnalysisConfig
305 _DefaultName = "limitingSurfaceBrightnessAnalysis"