Coverage for python / lsst / ip / diffim / computeSpatiallySampledMetrics.py: 21%
126 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-30 09:08 +0000
1# This file is part of ip_diffim.
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/>.
22import numpy as np
24import lsst.geom
26import lsst.afw.table as afwTable
27import lsst.pipe.base as pipeBase
28import lsst.pex.config as pexConfig
30from lsst.ip.diffim.utils import getPsfFwhm, angleMean, evaluateMaskFraction, getKernelCenterDisplacement
31from lsst.meas.algorithms import SkyObjectsTask
32from lsst.pex.exceptions import InvalidParameterError, RangeError
33from lsst.utils.timer import timeMethod
35import lsst.utils
37__all__ = ["SpatiallySampledMetricsConfig", "SpatiallySampledMetricsTask"]
40class SpatiallySampledMetricsConnections(pipeBase.PipelineTaskConnections,
41 dimensions=("instrument", "visit", "detector"),
42 defaultTemplates={"coaddName": "deep",
43 "warpTypeSuffix": "",
44 "fakesType": ""}):
45 science = pipeBase.connectionTypes.Input(
46 doc="Input science exposure.",
47 dimensions=("instrument", "visit", "detector"),
48 storageClass="ExposureF",
49 name="{fakesType}calexp"
50 )
51 template = pipeBase.connectionTypes.Input(
52 doc="Warped and not PSF-matched template used to create the difference image.",
53 dimensions=("instrument", "visit", "detector"),
54 storageClass="ExposureF",
55 name="{fakesType}{coaddName}Diff_templateExp",
56 )
57 difference = pipeBase.connectionTypes.Input(
58 doc="Difference image with detection mask plane filled in.",
59 dimensions=("instrument", "visit", "detector"),
60 storageClass="ExposureF",
61 name="{fakesType}{coaddName}Diff_differenceExp",
62 )
63 diaSources = pipeBase.connectionTypes.Input(
64 doc="Filtered diaSources on the difference image.",
65 dimensions=("instrument", "visit", "detector"),
66 storageClass="ArrowAstropy",
67 name="{fakesType}dia_source_detector",
68 )
69 psfMatchingKernel = pipeBase.connectionTypes.Input(
70 doc="Kernel used to PSF match the science and template images.",
71 dimensions=("instrument", "visit", "detector"),
72 storageClass="MatchingKernel",
73 name="{fakesType}{coaddName}Diff_psfMatchKernel",
74 )
75 spatiallySampledMetrics = pipeBase.connectionTypes.Output(
76 doc="Summary metrics computed at randomized locations.",
77 dimensions=("instrument", "visit", "detector"),
78 storageClass="ArrowAstropy",
79 name="{fakesType}{coaddName}Diff_spatiallySampledMetrics",
80 )
83class SpatiallySampledMetricsConfig(pipeBase.PipelineTaskConfig,
84 pipelineConnections=SpatiallySampledMetricsConnections):
85 """Config for SpatiallySampledMetricsTask
86 """
87 metricsMaskPlanes = lsst.pex.config.ListField(
88 dtype=str,
89 doc="List of mask planes to include in metrics",
90 default=('BAD', 'CLIPPED', 'CR', 'DETECTED', 'DETECTED_NEGATIVE', 'EDGE',
91 'INEXACT_PSF', 'INJECTED', 'INJECTED_TEMPLATE', 'INTRP', 'NOT_DEBLENDED',
92 'NO_DATA', 'REJECTED', 'SAT', 'SAT_TEMPLATE', 'SENSOR_EDGE', 'STREAK', 'SUSPECT',
93 'UNMASKEDNAN',
94 ),
95 )
96 metricSources = pexConfig.ConfigurableField(
97 target=SkyObjectsTask,
98 doc="Generate QA metric sources",
99 )
101 def setDefaults(self):
102 self.metricSources.avoidMask = ["NO_DATA", "EDGE"]
105class SpatiallySampledMetricsTask(lsst.pipe.base.PipelineTask):
106 """Detect and measure sources on a difference image.
107 """
108 ConfigClass = SpatiallySampledMetricsConfig
109 _DefaultName = "spatiallySampledMetrics"
111 def __init__(self, **kwargs):
112 super().__init__(**kwargs)
114 self.makeSubtask("metricSources")
115 self.schema = afwTable.SourceTable.makeMinimalSchema()
116 self.schema.addField(
117 "x", "F",
118 "X location of the metric evaluation.",
119 units="pixel")
120 self.schema.addField(
121 "y", "F",
122 "Y location of the metric evaluation.",
123 units="pixel")
124 self.metricSources.skySourceKey = self.schema.addField("sky_source", type="Flag",
125 doc="Metric evaluation objects.")
126 self.schema.addField(
127 "source_density", "F",
128 "Density of diaSources at location.",
129 units="count/degree^2")
130 self.schema.addField(
131 "dipole_density", "F",
132 "Density of dipoles at location.",
133 units="count/degree^2")
134 self.schema.addField(
135 "dipole_direction", "F",
136 "Mean dipole orientation.",
137 units="radian")
138 self.schema.addField(
139 "dipole_separation", "F",
140 "Mean dipole separation.",
141 units="pixel")
142 self.schema.addField(
143 "template_value", "F",
144 "Median of template at location.",
145 units="nJy")
146 self.schema.addField(
147 "template_variance", "F",
148 "Median of template variance at location.",
149 units="nJy^2")
150 self.schema.addField(
151 "science_value", "F",
152 "Median of science at location.",
153 units="nJy")
154 self.schema.addField(
155 "science_variance", "F",
156 "Median of science variance at location.",
157 units="nJy^2")
158 self.schema.addField(
159 "diffim_value", "F",
160 "Median of diffim at location.",
161 units="nJy")
162 self.schema.addField(
163 "diffim_variance", "F",
164 "Median of diffim variance at location.",
165 units="nJy^2")
166 self.schema.addField(
167 "science_psfSize", "F",
168 "Width of the science image PSF at location.",
169 units="pixel")
170 self.schema.addField(
171 "template_psfSize", "F",
172 "Width of the template image PSF at location.",
173 units="pixel")
174 for maskPlane in self.config.metricsMaskPlanes:
175 self.schema.addField(
176 "%s_mask_fraction"%maskPlane.lower(), "F",
177 "Fraction of pixels with %s mask"%maskPlane
178 )
179 self.schema.addField(
180 "psfMatchingKernel_sum", "F",
181 "PSF matching kernel sum at location.")
182 self.schema.addField(
183 "psfMatchingKernel_dx", "F",
184 "PSF matching kernel centroid offset in x at location.",
185 units="pixel")
186 self.schema.addField(
187 "psfMatchingKernel_dy", "F",
188 "PSF matching kernel centroid offset in y at location.",
189 units="pixel")
190 self.schema.addField(
191 "psfMatchingKernel_length", "F",
192 "PSF matching kernel centroid offset module.",
193 units="arcsecond")
194 self.schema.addField(
195 "psfMatchingKernel_position_angle", "F",
196 "PSF matching kernel centroid offset position angle.",
197 units="radian")
198 self.schema.addField(
199 "psfMatchingKernel_direction", "F",
200 "PSF matching kernel centroid offset direction in detector plane.",
201 units="radian")
203 @timeMethod
204 def run(self, science, template, difference, diaSources, psfMatchingKernel):
205 """Calculate difference image metrics on specific locations across the images
207 Parameters
208 ----------
209 science : `lsst.afw.image.ExposureF`
210 Science exposure that the template was subtracted from.
211 template : `lsst.afw.image.ExposureF`
212 Warped and non PSF-matched template that was used produce
213 the difference image.
214 difference : `lsst.afw.image.ExposureF`
215 Result of subtracting template from the science image.
216 diaSources : `lsst.afw.table.SourceCatalog`
217 The catalog of detected sources.
218 psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel`
219 The PSF matching kernel of the subtraction to evaluate.
221 Returns
222 -------
223 results : `lsst.pipe.base.Struct`
224 ``spatiallySampledMetrics`` : `astropy.table.Table`
225 Image quality metrics spatially sampled locations.
226 """
228 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
230 spatiallySampledMetrics = afwTable.SourceCatalog(self.schema)
231 spatiallySampledMetrics.getTable().setIdFactory(idFactory)
233 self.metricSources.run(mask=science.mask, seed=difference.info.id, catalog=spatiallySampledMetrics)
235 metricsMaskPlanes = []
236 for maskPlane in self.config.metricsMaskPlanes:
237 try:
238 metricsMaskPlanes.append(maskPlane)
239 except InvalidParameterError:
240 self.log.info("Unable to calculate metrics for mask plane %s: not in image"%maskPlane)
242 for src in spatiallySampledMetrics:
243 self._evaluateLocalMetric(src, science, template, difference, diaSources,
244 metricsMaskPlanes=metricsMaskPlanes,
245 psfMatchingKernel=psfMatchingKernel)
246 spatiallySampledMetrics = spatiallySampledMetrics.copy(deep=True).asAstropy()
247 return pipeBase.Struct(spatiallySampledMetrics=spatiallySampledMetrics)
249 def _evaluateLocalMetric(self, src, science, template, difference, diaSources,
250 metricsMaskPlanes, psfMatchingKernel):
251 """Calculate image quality metrics at spatially sampled locations.
253 Parameters
254 ----------
255 src : `lsst.afw.table.SourceRecord`
256 The source record to be updated with metric calculations.
257 diaSources : `lsst.afw.table.SourceCatalog`
258 The catalog of detected sources.
259 science : `lsst.afw.image.Exposure`
260 The science image.
261 difference : `lsst.afw.image.Exposure`
262 Result of subtracting template from the science image.
263 metricsMaskPlanes : `list` of `str`
264 Mask planes to calculate metrics from.
265 psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel`
266 The PSF matching kernel of the subtraction to evaluate.
267 """
268 bbox = src.getFootprint().getBBox()
269 pix = bbox.getCenter()
270 src.set('science_psfSize', getPsfFwhm(science.psf, position=pix))
271 try:
272 src.set('template_psfSize', getPsfFwhm(template.psf, position=pix))
273 except (InvalidParameterError, RangeError):
274 src.set('template_psfSize', np.nan)
276 metricRegionSize = 100
277 bbox.grow(metricRegionSize)
278 bbox = bbox.clippedTo(science.getBBox())
279 nPix = bbox.getArea()
280 pixScale = science.wcs.getPixelScale(bbox.getCenter())
281 area = nPix*pixScale.asDegrees()**2
282 peak = src.getFootprint().getPeaks()[0]
283 src.set('x', peak['i_x'])
284 src.set('y', peak['i_y'])
285 src.setCoord(science.wcs.pixelToSky(peak['i_x'], peak['i_y']))
286 selectSources = diaSources[bbox.contains(diaSources['x'], diaSources['y'])]
287 sourceDensity = len(selectSources)/area
288 dipoleSources = selectSources[selectSources["isDipole"]]
289 dipoleDensity = len(dipoleSources)/area
291 if dipoleSources:
292 meanDipoleOrientation = angleMean(dipoleSources["dipoleAngle"])
293 src.set('dipole_direction', meanDipoleOrientation)
294 meanDipoleSeparation = np.mean(dipoleSources["dipoleLength"])
295 src.set('dipole_separation', meanDipoleSeparation)
297 templateVal = np.median(template[bbox].image.array)
298 templateVar = np.median(template[bbox].variance.array)
299 scienceVal = np.median(science[bbox].image.array)
300 scienceVar = np.median(science[bbox].variance.array)
301 diffimVal = np.median(difference[bbox].image.array)
302 diffimVar = np.median(difference[bbox].variance.array)
303 src.set('source_density', sourceDensity)
304 src.set('dipole_density', dipoleDensity)
305 src.set('template_value', templateVal)
306 src.set('template_variance', templateVar)
307 src.set('science_value', scienceVal)
308 src.set('science_variance', scienceVar)
309 src.set('diffim_value', diffimVal)
310 src.set('diffim_variance', diffimVar)
311 for maskPlane in metricsMaskPlanes:
312 src.set("%s_mask_fraction"%maskPlane.lower(),
313 evaluateMaskFraction(difference.mask[bbox], maskPlane)
314 )
316 krnlSum, dx, dy, direction, length = getKernelCenterDisplacement(
317 psfMatchingKernel, src.get('x'), src.get('y'))
319 point1 = lsst.geom.SpherePoint(
320 src.get('coord_ra'), src.get('coord_dec'),
321 lsst.geom.radians)
322 point2 = science.wcs.pixelToSky(src.get('x') + dx, src.get('y') + dy)
323 bearing = point1.bearingTo(point2)
324 pa_ref_angle = lsst.geom.Angle(np.pi/2, lsst.geom.radians)
325 pa = pa_ref_angle - bearing
326 # Wrap around to get Delta_RA from -pi to +pi
327 pa = pa.wrapCtr()
328 position_angle = pa.asRadians()
330 src.set('psfMatchingKernel_sum', krnlSum)
331 src.set('psfMatchingKernel_dx', dx)
332 src.set('psfMatchingKernel_dy', dy)
333 src.set('psfMatchingKernel_length', length*pixScale.asArcseconds())
334 src.set('psfMatchingKernel_position_angle', position_angle) # in E of N position angle
335 src.set('psfMatchingKernel_direction', direction) # direction offset in detector