Coverage for python / lsst / pipe / tasks / matchDiffimSourceInjected.py: 29%
119 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-17 09:21 +0000
1# This file is part of ap_pipe.
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__ = ["MatchInjectedToDiaSourceTask",
23 "MatchInjectedToDiaSourceConfig",
24 "MatchInjectedToAssocDiaSourceTask",
25 "MatchInjectedToAssocDiaSourceConfig"]
27import astropy.units as u
28import numpy as np
29from scipy.spatial import cKDTree
31from lsst.afw import table as afwTable
32from lsst import geom as lsstGeom
33import lsst.pex.config as pexConfig
34from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct
35import lsst.pipe.base.connectionTypes as connTypes
36from lsst.meas.base import ForcedMeasurementTask, ForcedMeasurementConfig
39class MatchInjectedToDiaSourceConnections(
40 PipelineTaskConnections,
41 defaultTemplates={"coaddName": "deep",
42 "fakesType": "fakes_"},
43 dimensions=("instrument",
44 "visit",
45 "detector")):
46 injectedCat = connTypes.Input(
47 doc="Catalog of sources injected in the images.",
48 name="{fakesType}_pvi_catalog",
49 storageClass="ArrowAstropy",
50 dimensions=("instrument", "visit", "detector"),
51 )
52 diffIm = connTypes.Input(
53 doc="Difference image on which the DiaSources were detected.",
54 name="{fakesType}{coaddName}Diff_differenceExp",
55 storageClass="ExposureF",
56 dimensions=("instrument", "visit", "detector"),
57 )
58 diaSources = connTypes.Input(
59 doc="A DiaSource catalog to match against fakeCat.",
60 name="{fakesType}{coaddName}Diff_diaSrc",
61 storageClass="SourceCatalog",
62 dimensions=("instrument", "visit", "detector"),
63 )
64 matchDiaSources = connTypes.Output(
65 doc="A catalog of those fakeCat sources that have a match in "
66 "diaSrc. The schema is the union of the schemas for "
67 "``fakeCat`` and ``diaSrc``.",
68 name="{fakesType}{coaddName}Diff_matchDiaSrc",
69 storageClass="DataFrame",
70 dimensions=("instrument", "visit", "detector"),
71 )
74class MatchInjectedToDiaSourceConfig(
75 PipelineTaskConfig,
76 pipelineConnections=MatchInjectedToDiaSourceConnections):
77 """Config for MatchFakesTask.
78 """
79 matchDistanceArcseconds = pexConfig.RangeField(
80 doc="Distance in arcseconds to match within.",
81 dtype=float,
82 default=0.5,
83 min=0,
84 max=10,
85 )
86 doMatchVisit = pexConfig.Field(
87 dtype=bool,
88 default=True,
89 doc="Match visit to trim the fakeCat"
90 )
91 trimBuffer = pexConfig.Field(
92 doc="Size of the pixel buffer surrounding the image."
93 "Only those fake sources with a centroid"
94 "falling within the image+buffer region will be considered matches.",
95 dtype=int,
96 default=50,
97 )
98 doForcedMeasurement = pexConfig.Field(
99 dtype=bool,
100 default=True,
101 doc="Force measurement of the fakes at the injection locations."
102 )
103 forcedMeasurement = pexConfig.ConfigurableField(
104 target=ForcedMeasurementTask,
105 doc="Task to force photometer difference image at injection locations.",
106 )
109class MatchInjectedToDiaSourceTask(PipelineTask):
111 _DefaultName = "matchInjectedToDiaSource"
112 ConfigClass = MatchInjectedToDiaSourceConfig
114 def run(self, injectedCat, diffIm, diaSources):
115 """Match injected sources to detected diaSources within a difference image bound.
117 Parameters
118 ----------
119 injectedCat : `astropy.table.table.Table`
120 Table of catalog of synthetic sources to match to detected diaSources.
121 diffIm : `lsst.afw.image.Exposure`
122 Difference image where ``diaSources`` were detected.
123 diaSources : `afw.table.SourceCatalog`
124 Catalog of difference image sources detected in ``diffIm``.
125 Returns
126 -------
127 result : `lsst.pipe.base.Struct`
128 Results struct with components.
130 - ``matchedDiaSources`` : Fakes matched to input diaSources. Has
131 length of ``injectedCalexpCat``. (`pandas.DataFrame`)
132 """
134 if self.config.doMatchVisit:
135 fakeCat = self._trimFakeCat(injectedCat, diffIm)
136 else:
137 fakeCat = injectedCat
138 if self.config.doForcedMeasurement:
139 self._estimateFakesSNR(fakeCat, diffIm)
141 return self._processFakes(fakeCat, diaSources)
143 def _estimateFakesSNR(self, injectedCat, diffIm):
144 """Estimate the signal-to-noise ratio of the fakes in the given catalog.
146 Parameters
147 ----------
148 injectedCat : `astropy.table.Table`
149 Catalog of synthetic sources to estimate the S/N of. **This table
150 will be modified in place**.
151 diffIm : `lsst.afw.image.Exposure`
152 Difference image where the sources were detected.
153 """
154 # Create a schema for the forced measurement task
155 schema = afwTable.SourceTable.makeMinimalSchema()
156 schema.addField("x", "D", "x position in image.", units="pixel")
157 schema.addField("y", "D", "y position in image.", units="pixel")
158 schema.addField("deblend_nChild", "I", "Need for minimal forced phot schema")
160 pluginList = [
161 "base_PixelFlags",
162 "base_SdssCentroid",
163 "base_CircularApertureFlux",
164 "base_PsfFlux",
165 "base_LocalBackground"
166 ]
167 forcedMeasConfig = ForcedMeasurementConfig(plugins=pluginList)
168 forcedMeasConfig.slots.centroid = 'base_SdssCentroid'
169 forcedMeasConfig.slots.shape = None
171 # Create the forced measurement task
172 forcedMeas = ForcedMeasurementTask(schema, config=forcedMeasConfig)
174 # Specify the columns to copy from the input catalog to the output catalog
175 forcedMeas.copyColumns = {"coord_ra": "ra", "coord_dec": "dec"}
177 # Create an afw table from the input catalog
178 outputCatalog = afwTable.SourceCatalog(schema)
179 outputCatalog.reserve(len(injectedCat))
180 for row in injectedCat:
181 outputRecord = outputCatalog.addNew()
182 outputRecord.setId(row['injection_id'])
183 outputRecord.setCoord(lsstGeom.SpherePoint(row["ra"], row["dec"], lsstGeom.degrees))
184 outputRecord.set("x", row["x"])
185 outputRecord.set("y", row["y"])
187 # Generate the forced measurement catalog
188 forcedSources = forcedMeas.generateMeasCat(diffIm, outputCatalog, diffIm.getWcs())
189 # Attach the PSF shape footprints to the forced measurement catalog
190 forcedMeas.attachPsfShapeFootprints(forcedSources, diffIm)
192 # Copy the x and y positions from the forced measurement catalog back
193 # to the input catalog
194 for src, tgt in zip(forcedSources, outputCatalog):
195 src.set('base_SdssCentroid_x', tgt['x'])
196 src.set('base_SdssCentroid_y', tgt['y'])
198 # Define the centroid for the forced measurement catalog
199 forcedSources.defineCentroid('base_SdssCentroid')
200 # Run the forced measurement task
201 forcedMeas.run(forcedSources, diffIm, outputCatalog, diffIm.getWcs())
202 # Convert the forced measurement catalog to an astropy table
203 forcedSources_table = forcedSources.asAstropy()
205 # Add the forced measurement columns to the input catalog
206 for column in forcedSources_table.columns:
207 if "Flux" in column or "flag" in column:
208 injectedCat["forced_"+column] = forcedSources_table[column]
210 # Add the SNR columns to the input catalog
211 for column in injectedCat.colnames:
212 if column.endswith("instFlux"):
213 flux = injectedCat[column]
214 fluxErr = injectedCat[column+"Err"].copy()
215 fluxErr = np.where(
216 (fluxErr <= 0) | (np.isnan(fluxErr)), np.nanmax(fluxErr), fluxErr)
218 injectedCat[column+"_SNR"] = flux / fluxErr
220 def _processFakes(self, injectedCat, diaSources):
221 """Match fakes to detected diaSources within a difference image bound.
223 Parameters
224 ----------
225 injectedCat : `astropy.table.table.Table`
226 Catalog of injected sources to match to detected diaSources.
227 diaSources : `afw.table.SourceCatalog`
228 Catalog of difference image sources detected in ``diffIm``.
229 associatedDiaSources : `pandas.DataFrame`
230 Catalog of associated difference image sources detected in ``diffIm``.
232 Returns
233 -------
234 result : `lsst.pipe.base.Struct`
235 Results struct with components.
237 - ``matchedDiaSources`` : Fakes matched to input diaSources. Has
238 length of ``fakeCat``. (`pandas.DataFrame`)
239 """
240 # First match the diaSrc to the injected fakes
241 injectedCat = injectedCat.to_pandas()
242 nPossibleFakes = len(injectedCat)
244 fakeVects = self._getVectors(
245 np.radians(injectedCat.ra),
246 np.radians(injectedCat.dec))
247 diaSrcVects = self._getVectors(
248 diaSources['coord_ra'],
249 diaSources['coord_dec'])
251 diaSrcTree = cKDTree(diaSrcVects)
252 dist, idxs = diaSrcTree.query(
253 fakeVects,
254 distance_upper_bound=np.radians(self.config.matchDistanceArcseconds / 3600))
255 nFakesFound = np.isfinite(dist).sum()
257 self.log.info("Found %d out of %d possible in diaSources.", nFakesFound, nPossibleFakes)
259 # assign diaSourceId to the matched fakes
260 diaSrcIds = diaSources['id'][np.where(np.isfinite(dist), idxs, 0)]
261 matchedFakes = injectedCat.assign(diaSourceId=np.where(np.isfinite(dist), diaSrcIds, 0))
262 matchedFakes['dist_diaSrc'] = np.where(np.isfinite(dist), 3600*np.rad2deg(dist), -1)
264 return Struct(matchDiaSources=matchedFakes)
266 def _getVectors(self, ras, decs):
267 """Convert ra dec to unit vectors on the sphere.
269 Parameters
270 ----------
271 ras : `numpy.ndarray`, (N,)
272 RA coordinates in radians.
273 decs : `numpy.ndarray`, (N,)
274 Dec coordinates in radians.
276 Returns
277 -------
278 vectors : `numpy.ndarray`, (N, 3)
279 Vectors on the unit sphere for the given RA/DEC values.
280 """
281 vectors = np.empty((len(ras), 3))
283 vectors[:, 2] = np.sin(decs)
284 vectors[:, 0] = np.cos(decs) * np.cos(ras)
285 vectors[:, 1] = np.cos(decs) * np.sin(ras)
287 return vectors
289 def _addPixCoords(self, fakeCat, image):
290 """Add pixel coordinates to the catalog of fakes.
292 Parameters
293 ----------
294 fakeCat : `astropy.table.table.Table`
295 The catalog of fake sources to be input
296 image : `lsst.afw.image.exposure.exposure.ExposureF`
297 The image into which the fake sources should be added
298 Returns
299 -------
300 fakeCat : `astropy.table.table.Table`
301 """
303 wcs = image.getWcs()
305 # Get x/y pixel coordinates for injected sources.
306 xs, ys = wcs.skyToPixelArray(
307 fakeCat["ra"],
308 fakeCat["dec"],
309 degrees=True
310 )
311 fakeCat["x"] = xs
312 fakeCat["y"] = ys
314 return fakeCat
316 def _trimFakeCat(self, fakeCat, image):
317 """Trim the fake cat to the exact size of the input image.
319 Parameters
320 ----------
321 fakeCat : `astropy.table.table.Table`
322 The catalog of fake sources that was input
323 image : `lsst.afw.image.exposure.exposure.ExposureF`
324 The image into which the fake sources were added
325 Returns
326 -------
327 fakeCat : `astropy.table.table.Table`
328 The original fakeCat trimmed to the area of the image
329 """
331 # fakeCat must be processed with _addPixCoords before trimming
332 fakeCat = self._addPixCoords(fakeCat, image)
334 # Prefilter in ra/dec to avoid cases where the wcs incorrectly maps
335 # input fakes which are really off the chip onto it.
336 ras = fakeCat["ra"] * u.deg
337 decs = fakeCat["dec"] * u.deg
339 isContainedRaDec = image.containsSkyCoords(ras, decs, padding=0)
341 # now use the exact pixel BBox to filter to only fakes that were inserted
342 xs = fakeCat["x"]
343 ys = fakeCat["y"]
345 bbox = lsstGeom.Box2D(image.getBBox())
346 isContainedXy = xs - self.config.trimBuffer >= bbox.minX
347 isContainedXy &= xs + self.config.trimBuffer <= bbox.maxX
348 isContainedXy &= ys - self.config.trimBuffer >= bbox.minY
349 isContainedXy &= ys + self.config.trimBuffer <= bbox.maxY
351 return fakeCat[isContainedRaDec & isContainedXy]
354class MatchInjectedToAssocDiaSourceConnections(
355 PipelineTaskConnections,
356 defaultTemplates={"coaddName": "deep",
357 "fakesType": "fakes_"},
358 dimensions=("instrument",
359 "visit",
360 "detector")):
362 assocDiaSources = connTypes.Input(
363 doc="An assocDiaSource catalog to match against fakeCat from the"
364 "diaPipe run. Assumed to be SDMified.",
365 name="{fakesType}{coaddName}Diff_assocDiaSrc",
366 storageClass="DataFrame",
367 dimensions=("instrument", "visit", "detector"),
368 )
369 matchDiaSources = connTypes.Input(
370 doc="A catalog of those fakeCat sources that have a match in "
371 "diaSrc. The schema is the union of the schemas for "
372 "``fakeCat`` and ``diaSrc``.",
373 name="{fakesType}{coaddName}Diff_matchDiaSrc",
374 storageClass="DataFrame",
375 dimensions=("instrument", "visit", "detector"),
376 )
377 matchAssocDiaSources = connTypes.Output(
378 doc="A catalog of those fakeCat sources that have a match in "
379 "associatedDiaSources. The schema is the union of the schemas for "
380 "``fakeCat`` and ``associatedDiaSources``.",
381 name="{fakesType}{coaddName}Diff_matchAssocDiaSrc",
382 storageClass="DataFrame",
383 dimensions=("instrument", "visit", "detector"),
384 )
387class MatchInjectedToAssocDiaSourceConfig(
388 PipelineTaskConfig,
389 pipelineConnections=MatchInjectedToAssocDiaSourceConnections):
390 """Config for MatchFakesTask.
391 """
394class MatchInjectedToAssocDiaSourceTask(PipelineTask):
396 _DefaultName = "matchInjectedToAssocDiaSource"
397 ConfigClass = MatchInjectedToAssocDiaSourceConfig
399 def run(self, assocDiaSources, matchDiaSources):
400 """Tag matched injected sources to associated diaSources.
402 Parameters
403 ----------
404 matchDiaSources : `pandas.DataFrame`
405 Catalog of matched diaSrc to injected sources
406 assocDiaSources : `pandas.DataFrame`
407 Catalog of associated difference image sources detected in ``diffIm``.
408 Returns
409 -------
410 result : `lsst.pipe.base.Struct`
411 Results struct with components.
413 - ``matchAssocDiaSources`` : Fakes matched to associated diaSources. Has
414 length of ``matchDiaSources``. (`pandas.DataFrame`)
415 """
416 # Match the fakes to the associated sources. For this we don't use the coordinates
417 # but instead check for the diaSources. Since they were present in the table already
418 nPossibleFakes = len(matchDiaSources)
419 matchDiaSources['isAssocDiaSource'] = matchDiaSources.diaSourceId.isin(assocDiaSources.diaSourceId)
420 assocNFakesFound = matchDiaSources.isAssocDiaSource.sum()
421 self.log.info("Found %d out of %d possible in assocDiaSources."%(assocNFakesFound, nPossibleFakes))
423 return Struct(
424 matchAssocDiaSources=matchDiaSources.merge(
425 assocDiaSources.reset_index(drop=True),
426 on="diaSourceId",
427 how="left",
428 suffixes=('_ssi', '_diaSrc')
429 )
430 )