Coverage for python/lsst/ap/association/filterDiaSourceCatalog.py: 85%
143 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-07-02 08:59 +0000
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 "slot_Centroid_flag",
112 "base_PixelFlags_flag_crCenter",
113 "base_PixelFlags_flag_high_varianceCenterAll"
114 ],
115 doc="List of flags which cause a source to be removed.",
116 )
118 doRemoveNegativeDirectImageSources = pexConfig.Field(
119 dtype=bool,
120 default=True,
121 doc="Remove DIASources with negative scienceFlux/scienceFluxErr "
122 "according to a configurable threshold.",
123 )
125 minAllowedDirectSnr = pexConfig.Field(
126 dtype=float,
127 doc="Minimum allowed ratio of scienceFlux/scienceFluxErr.",
128 default=-2.0,
129 )
131 doTrailedSourceFilter = pexConfig.Field(
132 doc="Run trailedSourceFilter to remove long trailed sources from the"
133 "diaSource output catalog.",
134 dtype=bool,
135 default=True,
136 )
138 doWriteTrailedSources = pexConfig.Field(
139 doc="Write trailed diaSources sources to a table.",
140 dtype=bool,
141 default=True,
142 deprecated="Trailed sources will not be written out during production."
143 )
145 max_trail_length = pexConfig.Field(
146 dtype=float,
147 doc="Length of long trailed sources to remove from the input catalog, "
148 "in arcseconds per second. Default comes from DMTN-199, which "
149 "requires removal of sources with trails longer than 10 "
150 "degrees/day, which is 36000/3600/24 arcsec/second, or roughly"
151 "0.416 arcseconds per second.",
152 default=36000/3600.0/24.0,
153 )
154 estimatedPixelScale = pexConfig.Field(
155 dtype=float,
156 doc="Approximate plate scale, in arcseconds/pixel."
157 "Used to convert trail length if the catalog calculation fails.",
158 default=0.2,
159 )
162class FilterDiaSourceCatalogTask(pipeBase.PipelineTask):
163 """Filter sources from a DiaSource catalog based on their trail length,
164 sky sources, and bad flags.
165 """
167 ConfigClass = FilterDiaSourceCatalogConfig
168 _DefaultName = "filterDiaSourceCatalog"
170 @timeMethod
171 def run(self, diaSourceCat, diffImVisitInfo):
172 """Filter sources from the supplied DiaSource catalog.
174 DiaSources are filtered based on criteria including trail length,
175 whether they are a sky source, and whether they have one or more bad
176 flags. Any source which meets a filtering criteria is removed from the
177 main output catalog and optionally saved in one or more separate output
178 catalogs.
180 Parameters
181 ----------
182 diaSourceCat : `lsst.afw.table.SourceCatalog`
183 Catalog of sources measured on the difference image.
184 diffImVisitInfo: `lsst.afw.image.VisitInfo`
185 VisitInfo for the difference image corresponding to diaSourceCat.
187 Returns
188 -------
189 filterResults : `lsst.pipe.base.Struct`
191 ``filteredDiaSourceCat`` : `lsst.afw.table.SourceCatalog`
192 The catalog of filtered sources.
193 ``rejectedDiaSources`` : `lsst.afw.table.SourceCatalog`
194 The catalog of rejected sources.
195 ``longTrailedDiaSources`` : `astropy.table.Table`
196 DiaSources which have trail lengths greater than
197 max_trail_length*exposure_time.
198 """
199 rejectedSources = None
200 exposure_time = diffImVisitInfo.exposureTime
201 rejected_mask = np.zeros(len(diaSourceCat), dtype=bool)
202 if self.config.doRemoveSkySources:
203 sky_source_column = diaSourceCat["sky_source"]
204 self.log.info(f"Filtered {np.sum(sky_source_column & ~rejected_mask)} sky sources.")
205 rejected_mask |= sky_source_column
207 for flag in self.config.badFlagList:
208 flag_mask = diaSourceCat[flag]
209 self.log.info(f"Filtered {np.sum(flag_mask & ~rejected_mask)} sources with flag {flag}.")
210 rejected_mask |= flag_mask
212 if self.config.doRemoveNegativeDirectImageSources:
213 direct_snr = (diaSourceCat["ip_diffim_forced_PsfFlux_instFlux"]
214 / diaSourceCat["ip_diffim_forced_PsfFlux_instFluxErr"])
215 too_negative = direct_snr < self.config.minAllowedDirectSnr
216 self.log.info(f"Filtered {np.sum(too_negative & ~rejected_mask)} negative direct sources.")
217 rejected_mask |= too_negative
219 if self.config.doTrailedSourceFilter:
220 trail_mask = self._check_dia_source_trail(diaSourceCat, exposure_time)
221 bbox_mask = self._check_dia_source_trail_bbox(diaSourceCat, exposure_time)
222 num_trail_filtered = np.sum(trail_mask & ~bbox_mask)
223 num_bbox_filtered = np.sum(bbox_mask)
224 trail_mask |= bbox_mask
225 longTrailedDiaSources = diaSourceCat[trail_mask].copy(deep=True)
226 rejectedSources = diaSourceCat[rejected_mask].copy(deep=True)
227 rejected_mask |= trail_mask
228 diaSourceCat = diaSourceCat[~rejected_mask].copy(deep=True)
230 if num_trail_filtered > 0: 230 ↛ 237line 230 didn't jump to line 237 because the condition on line 230 was always true
232 self.log.info("%i DiaSources exceed max_trail_length "
233 "(%f arcseconds per second), dropping from "
234 "source catalog. "
235 % (num_trail_filtered, self.config.max_trail_length,))
237 if num_bbox_filtered > 0:
238 self.log.info(
239 " %i DiaSources had no trail calculation and their bounding"
240 " box exceeded max_trail_length (%f arcseconds per second) "
241 "in either direction or across the diagonal, dropping from "
242 "source catalog."
243 % (num_bbox_filtered, self.config.max_trail_length))
245 self.metadata.add("num_filtered", len(longTrailedDiaSources))
247 if self.config.doWriteTrailedSources: 247 ↛ 252line 247 didn't jump to line 252 because the condition on line 247 was always true
248 filterResults = pipeBase.Struct(filteredDiaSourceCat=diaSourceCat,
249 rejectedDiaSources=rejectedSources,
250 longTrailedSources=longTrailedDiaSources.asAstropy())
251 else:
252 filterResults = pipeBase.Struct(filteredDiaSourceCat=diaSourceCat,
253 rejectedDiaSources=rejectedSources)
254 else:
255 rejectedSources = diaSourceCat[rejected_mask].copy(deep=True)
256 diaSourceCat = diaSourceCat[~rejected_mask].copy(deep=True)
257 filterResults = pipeBase.Struct(filteredDiaSourceCat=diaSourceCat,
258 rejectedDiaSources=rejectedSources)
259 return filterResults
261 def _check_dia_source_trail(self, dia_sources, exposure_time):
262 """Find DiaSources that have long trails or trails with indeterminant
263 end points.
265 Return a mask of sources with lengths greater than
266 (``config.max_trail_length`` multiplied by the exposure time)
267 arcseconds.
268 Additionally, set mask if
269 ``ext_trailedSources_Naive_flag_off_image`` is set or if
270 ``ext_trailedSources_Naive_flag_suspect_long_trail`` and
271 ``ext_trailedSources_Naive_flag_edge`` are both set.
273 Parameters
274 ----------
275 dia_sources : `pandas.DataFrame`
276 Input diaSources to check for trail lengths.
277 exposure_time : `float`
278 Exposure time from difference image.
280 Returns
281 -------
282 trail_mask : `pandas.DataFrame`
283 Boolean mask for diaSources which are greater than the
284 Boolean mask for diaSources which are greater than the
285 cutoff length or have trails which extend beyond the edge of the
286 detector (off_image set). Also checks if both
287 suspect_long_trail and edge are set and masks those sources out.
288 """
289 pixelScale = self._estimate_pixel_scale(dia_sources)
290 trail_mask = (dia_sources["ext_trailedSources_Naive_length"]
291 >= (self.config.max_trail_length*exposure_time/pixelScale))
292 trail_mask |= dia_sources['ext_trailedSources_Naive_flag_off_image']
293 trail_mask |= (dia_sources['ext_trailedSources_Naive_flag_suspect_long_trail']
294 & dia_sources['ext_trailedSources_Naive_flag_edge'])
296 return trail_mask
298 def _check_dia_source_trail_bbox(self, dia_sources, exposure_time):
299 """Check bounding boxes of DiaSources with failed trail measurements.
301 For sources where the trail measurement flag is set (indicating the
302 trail length is unable to be measured), fall back to checking the footprint
303 bounding box dimensions against the max trail length threshold.
305 Parameters
306 ----------
307 dia_sources : `lsst.afw.table.SourceCatalog`
308 Input diaSources to check.
309 exposure_time : `float`
310 Exposure time from difference image.
312 Returns
313 -------
314 bbox_mask : `numpy.ndarray`
315 Boolean mask for diaSources whose footprint bounding box width,
316 height, or diagonal exceeds the max trail length threshold.
317 """
318 pixelScale = self._estimate_pixel_scale(dia_sources)
319 max_length_pixels = self.config.max_trail_length * exposure_time / pixelScale
321 trail_flag = dia_sources["ext_trailedSources_Naive_flag"]
322 bbox_mask = np.zeros(len(dia_sources), dtype=bool)
324 for i in np.where(trail_flag)[0]:
325 bbox = dia_sources[i].getFootprint().getBBox()
326 width = bbox.getWidth()
327 height = bbox.getHeight()
328 diagonal = np.sqrt(width ** 2 + height ** 2)
329 if diagonal >= max_length_pixels:
330 bbox_mask[i] = True
332 return bbox_mask
334 def _estimate_pixel_scale(self, catalog):
335 """Quickly calculate the pixel scale from catalog values
337 Will return a fallback value from the task config if there is any error
339 Parameters
340 ----------
341 catalog : `lsst.afw.table.SourceCatalog`
342 Catalog of sources measured on the difference image.
344 Returns
345 -------
346 scale : `float`
347 Pixel scale of the catalog, in arcseconds/pixel
348 """
349 nSrc = len(catalog)
350 if nSrc < 2: 350 ↛ 351line 350 didn't jump to line 351 because the condition on line 350 was never true
351 return self.config.estimatedPixelScale
352 try:
353 coordKey = catalog.getCoordKey()
354 decVals = catalog[coordKey.getDec()] # in radians
355 raVals = catalog[coordKey.getRa()] # in radians
356 xVals = catalog.getX()
357 yVals = catalog.getY()
358 # Find two points that are well separated for the calculation
359 # Start with a point near one edge, and find the furthest point
360 # from there
361 iMin = np.argmin(xVals)
362 dist = np.sqrt((xVals[iMin] - xVals)**2 + (yVals[iMin] - yVals)**2)
363 iMax = np.argmax(dist)
364 # Use the spherical law of cosines:
365 t1 = np.sin(decVals[iMin])*np.sin(decVals[iMax])
366 t2 = np.cos(decVals[iMin])*np.cos(decVals[iMax])*np.cos(raVals[iMin] - raVals[iMax])
367 separation = np.arccos(t1 + t2) # in radians
368 scale = separation/max(dist)*3600*180/np.pi # convert to arcseconds/pixel
369 except Exception as e:
370 self.log.warning("Error encountered estimating the pixel scale from the catalog: %s", e)
371 return self.config.estimatedPixelScale
372 else:
373 if abs(scale - self.config.estimatedPixelScale)/self.config.estimatedPixelScale > 0.1:
374 self.log.warning(f"Calculated pixel scale of {scale} too different from estimated value "
375 f"{self.config.estimatedPixelScale}. Falling back on estimate.")
376 return self.config.estimatedPixelScale
377 else:
378 return scale
381class FilterDiaSourceReliabilityConnections(
382 pipeBase.PipelineTaskConnections,
383 dimensions=("instrument", "visit", "detector"),
384 defaultTemplates={"coaddName": "deep", "fakesType": ""}
385):
386 diaSourceCat = connTypes.Input(
387 doc="Catalog of DiaSources produced during image differencing.",
388 name="{fakesType}{coaddName}Diff_candidateDiaSrc",
389 storageClass="SourceCatalog",
390 dimensions=("instrument", "visit", "detector"),
391 )
392 reliability = connTypes.Input(
393 doc="Reliability (e.g. real/bogus) classificiation of diaSourceCat sources.",
394 name="{fakesType}{coaddName}RealBogusSources",
395 storageClass="Catalog",
396 dimensions=("instrument", "visit", "detector"),
397 )
398 filteredDiaSources = connTypes.Output(
399 doc="Accepted diaSource catalog filtered by reliability score.",
400 name="dia_source_high_reliability",
401 storageClass="SourceCatalog",
402 dimensions=("instrument", "visit", "detector"),
403 )
404 rejectedDiaSources = connTypes.Output(
405 doc="Rejected diaSource catalog with low reliability scores.",
406 name="dia_source_low_reliability",
407 storageClass="SourceCatalog",
408 dimensions=("instrument", "visit", "detector"),
409 )
412class FilterDiaSourceReliabilityConfig(
413 pipeBase.PipelineTaskConfig, pipelineConnections=FilterDiaSourceReliabilityConnections
414):
415 """Configuration for the FilterDiaSourceReliabilityTask."""
416 minReliability = pexConfig.Field(
417 doc="Minimum reliability score to keep a source in the DiaSource catalog.",
418 dtype=float,
419 default=0.0,
420 )
423class FilterDiaSourceReliabilityTask(pipeBase.PipelineTask):
424 """Filter DiaSource catalog by reliability score.
426 Parameters
427 ----------
428 diaSourceSchema: `lsst.afw.table.Schema`
429 Schema for the input DiaSource catalog.
430 diaSourceCat : `lsst.afw.table.SourceCatalog`
431 Catalog of DiaSources produced during image differencing.
432 reliability : `lsst.afw.table.SourceCatalog`, optional
433 Reliability (e.g. real/bogus) classification of the sources in `diaSourceCat`.
435 Returns
436 -------
437 filteredResults : `lsst.pipe.base.Struct`
439 ``filteredDiaSources`` : `lsst.afw.table.SourceCatalog`
440 Catalog of unstandardized DiaSources filtered by reliability score.
441 ``rejectedDiaSources`` : `lsst.afw.table.SourceCatalog`
442 Catalog of unstandardized DiaSources that were rejected due to low
443 reliability scores.
444 """
446 ConfigClass = FilterDiaSourceReliabilityConfig
447 _DefaultName = "filterDiaSourceReliability"
449 def runQuantum(self, butlerQC, inputRefs, outputRefs):
450 inputs = butlerQC.get(inputRefs)
451 outputs = self.run(**inputs)
452 butlerQC.put(outputs, outputRefs)
454 def run(self, diaSourceCat, reliability):
455 """Run the task to filter DiaSources by reliability."""
457 # Copy the scores in the output catalog
458 if np.all(diaSourceCat['id'] == reliability['id']): 458 ↛ 466line 458 didn't jump to line 466 because the condition on line 458 was always true
459 diaSourceCat['reliability'] = reliability['score']
460 src_key = reliability.schema.find('version').key
461 dst_key = diaSourceCat.schema.find('reliabilityVersion').key
462 for i in range(len(diaSourceCat)):
463 diaSourceCat[i].set(dst_key, reliability[i].get(src_key))
464 else:
465 # If the identifiers do not match, we cannot filter reliably.
466 raise ValueError(
467 "Reliability ids do not match DiaSource ids.")
469 # Filter the DiaSource catalog by reliability score
470 low_reliability = diaSourceCat["reliability"] < self.config.minReliability
471 rejectedDiaSources = diaSourceCat[low_reliability].copy(deep=True)
472 filteredDiaSources = diaSourceCat[~low_reliability].copy(deep=True)
474 self.log.info(f"Filtered {np.sum(low_reliability)} sources with low reliability.")
476 return pipeBase.Struct(
477 filteredDiaSources=filteredDiaSources,
478 rejectedDiaSources=rejectedDiaSources
479 )