Coverage for python / lsst / ap / association / filterDiaSourceCatalog.py: 30%
119 statements
« prev ^ index » next coverage.py v7.13.2, created at 2026-02-03 02:01 -0800
1# This file is part of ap_association
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 "FilterDiaSourceCatalogConfig",
24 "FilterDiaSourceCatalogTask",
25 "FilterDiaSourceReliabilityConfig",
26 "FilterDiaSourceReliabilityTask"
27)
29import numpy as np
31import lsst.pex.config as pexConfig
32import lsst.pipe.base as pipeBase
33import lsst.pipe.base.connectionTypes as connTypes
34from lsst.utils.timer import timeMethod
37class FilterDiaSourceCatalogConnections(
38 pipeBase.PipelineTaskConnections,
39 dimensions=("instrument", "visit", "detector"),
40 defaultTemplates={"coaddName": "deep", "fakesType": ""},
41):
42 """Connections class for FilterDiaSourceCatalogTask."""
44 diaSourceCat = connTypes.Input(
45 doc="Catalog of DiaSources produced during image differencing.",
46 name="{fakesType}{coaddName}Diff_diaSrc",
47 storageClass="SourceCatalog",
48 dimensions=("instrument", "visit", "detector"),
49 )
51 diffImVisitInfo = connTypes.Input(
52 doc="VisitInfo of diffIm.",
53 name="{fakesType}{coaddName}Diff_differenceExp.visitInfo",
54 storageClass="VisitInfo",
55 dimensions=("instrument", "visit", "detector"),
56 )
58 filteredDiaSourceCat = connTypes.Output(
59 doc="Output catalog of DiaSources after filtering.",
60 name="{fakesType}{coaddName}Diff_candidateDiaSrc",
61 storageClass="SourceCatalog",
62 dimensions=("instrument", "visit", "detector"),
63 )
65 rejectedDiaSources = connTypes.Output(
66 doc="Optional output storing all the rejected DiaSources.",
67 name="{fakesType}{coaddName}Diff_rejectedDiaSrc",
68 storageClass="SourceCatalog",
69 dimensions={"instrument", "visit", "detector"},
70 )
72 longTrailedSources = connTypes.Output(
73 doc="Optional output temporarily storing long trailed diaSources.",
74 dimensions=("instrument", "visit", "detector"),
75 storageClass="ArrowAstropy",
76 name="{fakesType}{coaddName}Diff_longTrailedSrc",
77 )
79 def __init__(self, *, config=None):
80 super().__init__(config=config)
81 if not self.config.doWriteRejectedSkySources:
82 self.outputs.remove("rejectedDiaSources")
83 if not self.config.doTrailedSourceFilter:
84 self.outputs.remove("longTrailedSources")
85 if not self.config.doWriteTrailedSources:
86 self.outputs.remove("longTrailedSources")
89class FilterDiaSourceCatalogConfig(
90 pipeBase.PipelineTaskConfig, pipelineConnections=FilterDiaSourceCatalogConnections
91):
92 """Config class for FilterDiaSourceCatalogTask."""
94 doRemoveSkySources = pexConfig.Field(
95 dtype=bool,
96 default=False,
97 doc="Input DiaSource catalog contains SkySources that should be "
98 "removed before storing the output DiaSource catalog.",
99 )
101 doWriteRejectedSkySources = pexConfig.Field(
102 dtype=bool,
103 default=True,
104 doc="Store the output DiaSource catalog containing all the rejected "
105 "sky sources."
106 )
108 badFlagList = pexConfig.ListField(
109 dtype=str,
110 default=[
111 "base_PixelFlags_flag_crCenter",
112 "base_PixelFlags_flag_high_varianceCenterAll"
113 ],
114 doc="List of flags which cause a source to be removed.",
115 )
117 doRemoveNegativeDirectImageSources = pexConfig.Field(
118 dtype=bool,
119 default=True,
120 doc="Remove DIASources with negative scienceFlux/scienceFluxErr "
121 "according to a configurable threshold.",
122 )
124 minAllowedDirectSnr = pexConfig.Field(
125 dtype=float,
126 doc="Minimum allowed ratio of scienceFlux/scienceFluxErr.",
127 default=-2.0,
128 )
130 doTrailedSourceFilter = pexConfig.Field(
131 doc="Run trailedSourceFilter to remove long trailed sources from the"
132 "diaSource output catalog.",
133 dtype=bool,
134 default=True,
135 )
137 doWriteTrailedSources = pexConfig.Field(
138 doc="Write trailed diaSources sources to a table.",
139 dtype=bool,
140 default=True,
141 deprecated="Trailed sources will not be written out during production."
142 )
144 max_trail_length = pexConfig.Field(
145 dtype=float,
146 doc="Length of long trailed sources to remove from the input catalog, "
147 "in arcseconds per second. Default comes from DMTN-199, which "
148 "requires removal of sources with trails longer than 10 "
149 "degrees/day, which is 36000/3600/24 arcsec/second, or roughly"
150 "0.416 arcseconds per second.",
151 default=36000/3600.0/24.0,
152 )
153 estimatedPixelScale = pexConfig.Field(
154 dtype=float,
155 doc="Approximate plate scale, in arcseconds/pixel."
156 "Used to convert trail length if the catalog calculation fails.",
157 default=0.2,
158 )
161class FilterDiaSourceCatalogTask(pipeBase.PipelineTask):
162 """Filter sources from a DiaSource catalog."""
164 ConfigClass = FilterDiaSourceCatalogConfig
165 _DefaultName = "filterDiaSourceCatalog"
167 @timeMethod
168 def run(self, diaSourceCat, diffImVisitInfo):
169 """Filter sources from the supplied DiaSource catalog.
171 Parameters
172 ----------
173 diaSourceCat : `lsst.afw.table.SourceCatalog`
174 Catalog of sources measured on the difference image.
175 diffImVisitInfo: `lsst.afw.image.VisitInfo`
176 VisitInfo for the difference image corresponding to diaSourceCat.
178 Returns
179 -------
180 filterResults : `lsst.pipe.base.Struct`
182 ``filteredDiaSourceCat`` : `lsst.afw.table.SourceCatalog`
183 The catalog of filtered sources.
184 ``rejectedDiaSources`` : `lsst.afw.table.SourceCatalog`
185 The catalog of rejected sources.
186 ``longTrailedDiaSources`` : `astropy.table.Table`
187 DiaSources which have trail lengths greater than
188 max_trail_length*exposure_time.
189 """
190 rejectedSources = None
191 exposure_time = diffImVisitInfo.exposureTime
192 rejected_mask = np.zeros(len(diaSourceCat), dtype=bool)
193 if self.config.doRemoveSkySources:
194 sky_source_column = diaSourceCat["sky_source"]
195 self.log.info(f"Filtered {np.sum(sky_source_column & ~rejected_mask)} sky sources.")
196 rejected_mask |= sky_source_column
198 for flag in self.config.badFlagList:
199 flag_mask = diaSourceCat[flag]
200 self.log.info(f"Filtered {np.sum(flag_mask & ~rejected_mask)} sources with flag {flag}.")
201 rejected_mask |= flag_mask
203 if self.config.doRemoveNegativeDirectImageSources:
204 direct_snr = (diaSourceCat["ip_diffim_forced_PsfFlux_instFlux"]
205 / diaSourceCat["ip_diffim_forced_PsfFlux_instFluxErr"])
206 too_negative = direct_snr < self.config.minAllowedDirectSnr
207 self.log.info(f"Filtered {np.sum(too_negative & ~rejected_mask)} negative direct sources.")
208 rejected_mask |= too_negative
210 if self.config.doTrailedSourceFilter:
211 trail_mask = self._check_dia_source_trail(diaSourceCat, exposure_time)
212 longTrailedDiaSources = diaSourceCat[trail_mask].copy(deep=True)
213 rejectedSources = diaSourceCat[rejected_mask].copy(deep=True)
214 rejected_mask |= trail_mask
215 diaSourceCat = diaSourceCat[~rejected_mask].copy(deep=True)
217 self.log.info("%i DiaSources exceed max_trail_length %f arcseconds per second, "
218 "dropping from source catalog."
219 % (len(longTrailedDiaSources), self.config.max_trail_length))
220 self.metadata.add("num_filtered", len(longTrailedDiaSources))
222 if self.config.doWriteTrailedSources:
223 filterResults = pipeBase.Struct(filteredDiaSourceCat=diaSourceCat,
224 rejectedDiaSources=rejectedSources,
225 longTrailedSources=longTrailedDiaSources.asAstropy())
226 else:
227 filterResults = pipeBase.Struct(filteredDiaSourceCat=diaSourceCat,
228 rejectedDiaSources=rejectedSources)
229 else:
230 rejectedSources = diaSourceCat[rejected_mask].copy(deep=True)
231 diaSourceCat = diaSourceCat[~rejected_mask].copy(deep=True)
232 filterResults = pipeBase.Struct(filteredDiaSourceCat=diaSourceCat,
233 rejectedDiaSources=rejectedSources)
234 return filterResults
236 def _check_dia_source_trail(self, dia_sources, exposure_time):
237 """Find DiaSources that have long trails or trails with indeterminant
238 end points.
240 Return a mask of sources with lengths greater than
241 (``config.max_trail_length`` multiplied by the exposure time)
242 arcseconds.
243 Additionally, set mask if
244 ``ext_trailedSources_Naive_flag_off_image`` is set or if
245 ``ext_trailedSources_Naive_flag_suspect_long_trail`` and
246 ``ext_trailedSources_Naive_flag_edge`` are both set.
248 Parameters
249 ----------
250 dia_sources : `pandas.DataFrame`
251 Input diaSources to check for trail lengths.
252 exposure_time : `float`
253 Exposure time from difference image.
255 Returns
256 -------
257 trail_mask : `pandas.DataFrame`
258 Boolean mask for diaSources which are greater than the
259 Boolean mask for diaSources which are greater than the
260 cutoff length or have trails which extend beyond the edge of the
261 detector (off_image set). Also checks if both
262 suspect_long_trail and edge are set and masks those sources out.
263 """
264 pixelScale = self._estimate_pixel_scale(dia_sources)
265 trail_mask = (dia_sources["ext_trailedSources_Naive_length"]
266 >= (self.config.max_trail_length*exposure_time/pixelScale))
267 trail_mask |= dia_sources['ext_trailedSources_Naive_flag_off_image']
268 trail_mask |= (dia_sources['ext_trailedSources_Naive_flag_suspect_long_trail']
269 & dia_sources['ext_trailedSources_Naive_flag_edge'])
271 return trail_mask
273 def _estimate_pixel_scale(self, catalog):
274 """Quickly calculate the pixel scale from catalog values
276 Will return a fallback value from the task config if there is any error
278 Parameters
279 ----------
280 catalog : `lsst.afw.table.SourceCatalog`
281 Catalog of sources measured on the difference image.
283 Returns
284 -------
285 scale : `float`
286 Pixel scale of the catalog, in arcseconds/pixel
287 """
288 nSrc = len(catalog)
289 if nSrc < 2:
290 return self.config.estimatedPixelScale
291 try:
292 coordKey = catalog.getCoordKey()
293 decVals = catalog[coordKey.getDec()] # in radians
294 raVals = catalog[coordKey.getRa()] # in radians
295 xVals = catalog.getX()
296 yVals = catalog.getY()
297 # Find two points that are well separated for the calculation
298 # Start with a point near one edge, and find the furthest point
299 # from there
300 iMin = np.argmin(xVals)
301 dist = np.sqrt((xVals[iMin] - xVals)**2 + (yVals[iMin] - yVals)**2)
302 iMax = np.argmax(dist)
303 # Use the spherical law of cosines:
304 t1 = np.sin(decVals[iMin])*np.sin(decVals[iMax])
305 t2 = np.cos(decVals[iMin])*np.cos(decVals[iMax])*np.cos(raVals[iMin] - raVals[iMax])
306 separation = np.arccos(t1 + t2) # in radians
307 scale = separation/max(dist)*3600*180/np.pi # convert to arcseconds/pixel
308 except Exception as e:
309 self.log.warning("Error encountered estimating the pixel scale from the catalog: %s", e)
310 return self.config.estimatedPixelScale
311 else:
312 if abs(scale - self.config.estimatedPixelScale)/self.config.estimatedPixelScale > 0.1:
313 self.log.warning(f"Calculated pixel scale of {scale} too different from estimated value "
314 f"{self.config.estimatedPixelScale}. Falling back on estimate.")
315 return self.config.estimatedPixelScale
316 else:
317 return scale
320class FilterDiaSourceReliabilityConnections(
321 pipeBase.PipelineTaskConnections,
322 dimensions=("instrument", "visit", "detector"),
323 defaultTemplates={"coaddName": "deep", "fakesType": ""}
324):
325 diaSourceCat = connTypes.Input(
326 doc="Catalog of DiaSources produced during image differencing.",
327 name="{fakesType}{coaddName}Diff_candidateDiaSrc",
328 storageClass="SourceCatalog",
329 dimensions=("instrument", "visit", "detector"),
330 )
331 reliability = connTypes.Input(
332 doc="Reliability (e.g. real/bogus) classificiation of diaSourceCat sources.",
333 name="{fakesType}{coaddName}RealBogusSources",
334 storageClass="Catalog",
335 dimensions=("instrument", "visit", "detector"),
336 )
337 filteredDiaSources = connTypes.Output(
338 doc="Accepted diaSource catalog filtered by reliability score.",
339 name="dia_source_high_reliability",
340 storageClass="SourceCatalog",
341 dimensions=("instrument", "visit", "detector"),
342 )
343 rejectedDiaSources = connTypes.Output(
344 doc="Rejected diaSource catalog with low reliability scores.",
345 name="dia_source_low_reliability",
346 storageClass="SourceCatalog",
347 dimensions=("instrument", "visit", "detector"),
348 )
351class FilterDiaSourceReliabilityConfig(
352 pipeBase.PipelineTaskConfig, pipelineConnections=FilterDiaSourceReliabilityConnections
353):
354 """Configuration for the FilterDiaSourceReliabilityTask."""
355 minReliability = pexConfig.Field(
356 doc="Minimum reliability score to keep a source in the DiaSource catalog.",
357 dtype=float,
358 default=0.0,
359 )
362class FilterDiaSourceReliabilityTask(pipeBase.PipelineTask):
363 """Filter DiaSource catalog by reliability score.
365 Parameters
366 ----------
367 diaSourceSchema: `lsst.afw.table.Schema`
368 Schema for the input DiaSource catalog.
369 diaSourceCat : `lsst.afw.table.SourceCatalog`
370 Catalog of DiaSources produced during image differencing.
371 reliability : `lsst.afw.table.SourceCatalog`, optional
372 Reliability (e.g. real/bogus) classification of the sources in `diaSourceCat`.
374 Returns
375 -------
376 filteredResults : `lsst.pipe.base.Struct`
378 ``filteredDiaSources`` : `lsst.afw.table.SourceCatalog`
379 Catalog of unstandardized DiaSources filtered by reliability score.
380 ``rejectedDiaSources`` : `lsst.afw.table.SourceCatalog`
381 Catalog of unstandardized DiaSources that were rejected due to low
382 reliability scores.
383 """
385 ConfigClass = FilterDiaSourceReliabilityConfig
386 _DefaultName = "filterDiaSourceReliability"
388 def runQuantum(self, butlerQC, inputRefs, outputRefs):
389 inputs = butlerQC.get(inputRefs)
390 outputs = self.run(**inputs)
391 butlerQC.put(outputs, outputRefs)
393 def run(self, diaSourceCat, reliability):
394 """Run the task to filter DiaSources by reliability."""
396 # Copy the scores in the output catalog
397 if np.all(diaSourceCat['id'] == reliability['id']):
398 diaSourceCat['reliability'] = reliability['score']
399 else:
400 # If the identifiers do not match, we cannot filter reliably.
401 raise ValueError(
402 "Reliability ids do not match DiaSource ids.")
404 # Filter the DiaSource catalog by reliability score
405 low_reliability = diaSourceCat["reliability"] < self.config.minReliability
406 rejectedDiaSources = diaSourceCat[low_reliability].copy(deep=True)
407 filteredDiaSources = diaSourceCat[~low_reliability].copy(deep=True)
409 self.log.info(f"Filtered {np.sum(low_reliability)} sources with low reliability.")
411 return pipeBase.Struct(
412 filteredDiaSources=filteredDiaSources,
413 rejectedDiaSources=rejectedDiaSources
414 )