Coverage for tests / test_filterDiaSourceCatalog.py: 12%
200 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-14 00:44 -0700
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/>.
22import unittest
23import numpy as np
25from lsst.ap.association.filterDiaSourceCatalog import (FilterDiaSourceCatalogConfig,
26 FilterDiaSourceCatalogTask,
27 FilterDiaSourceReliabilityConfig,
28 FilterDiaSourceReliabilityTask)
29import lsst.geom as geom
30import lsst.meas.base.tests as measTests
31import lsst.utils.tests
32import lsst.afw.image as afwImage
33import lsst.afw.table as afwTable
34import lsst.daf.base as dafBase
37class TestFilterDiaSourceCatalogTask(unittest.TestCase):
39 def setUp(self):
40 self.config = FilterDiaSourceCatalogConfig()
42 self.nSkySources = 5
43 self.nCrCenterSources = 6
44 self.nFakeFlagSources = 4
45 self.nNegativeSources = 7
46 self.nTrailedSources = 10
47 self.pixelScale = 0.2 # arcseconds/pixel
48 self.nSources = (self.nSkySources + self.nCrCenterSources + self.nFakeFlagSources
49 + self.nNegativeSources + self.nTrailedSources)
50 self.yLoc = 100
51 self.expId = 4321
52 self.bbox = geom.Box2I(geom.Point2I(0, 0),
53 geom.Extent2I(1024, 1153))
54 dataset = measTests.TestDataset(self.bbox)
55 for srcIdx in range(self.nSources):
56 dataset.addSource(10000.0, geom.Point2D(srcIdx, self.yLoc))
57 schema = dataset.makeMinimalSchema()
58 schema.addField("sky_source", type="Flag", doc="Sky objects.")
59 schema.addField("base_PixelFlags_flag_crCenter", type="Flag", doc="A cosmic ray was detected "
60 "and interpolated in this object's center.")
61 schema.addField("base_PixelFlags_flag_high_varianceCenterAll", type="Flag", doc="The object was "
62 "detected in a region with exceptionally high template variance.")
63 schema.addField("fakeBadFlag", type="Flag", doc="A fake flag to test a badFlagList longer "
64 "than one item.")
65 schema.addField("ip_diffim_forced_PsfFlux_instFlux", type="F",
66 doc="Forced photometry flux for a point source model measured on the visit image "
67 "centered at DiaSource position.")
68 schema.addField("ip_diffim_forced_PsfFlux_instFluxErr", type="F",
69 doc="Estimated uncertainty of ip_diffim_forced_PsfFlux_instFlux.")
70 schema.addField('ext_trailedSources_Naive_flag_off_image', type="Flag",
71 doc="Trail extends off image")
72 schema.addField('ext_trailedSources_Naive_flag_suspect_long_trail',
73 type="Flag", doc="Trail length is greater than three times the psf radius")
74 schema.addField('ext_trailedSources_Naive_flag_edge', type="Flag",
75 doc="Trail contains edge pixels")
76 schema.addField('ext_trailedSources_Naive_flag_nan', type="Flag",
77 doc="One or more trail coordinates are missing")
78 schema.addField('ext_trailedSources_Naive_length', type="F",
79 doc="Length of the source trail")
80 schema.addField("reliability", type="F",
81 doc="Reliability of the source")
82 _, self.diaSourceCat = dataset.realize(10.0, schema, randomSeed=1234)
84 # set the sky_source flag for the first set
85 self.diaSourceCat[0:self.nSkySources]["sky_source"] = True
87 # set the pixelFlags_crCenter flag
88 crCenter_offset = self.nSkySources
89 self.diaSourceCat[crCenter_offset:crCenter_offset+self.nCrCenterSources][
90 "base_PixelFlags_flag_crCenter"
91 ] = True
93 # set the fakeBadFlag flag
94 fakeFlag_offset = crCenter_offset + self.nCrCenterSources
95 self.diaSourceCat[fakeFlag_offset:fakeFlag_offset+self.nFakeFlagSources][
96 "fakeBadFlag"
97 ] = True
99 # create increasingly negative ip_diffim_forced_PsfFlux_instFlux/ip_diffim_forced_PsfFlux_instFluxErr
100 self.nRemovedNegativeSources = 0
101 negativeSources_offset = fakeFlag_offset + self.nFakeFlagSources
102 for i, srcIdx in enumerate(range(negativeSources_offset,
103 negativeSources_offset+self.nNegativeSources)):
104 self.diaSourceCat[srcIdx]["ip_diffim_forced_PsfFlux_instFlux"] = -0.5 * i
105 self.diaSourceCat[srcIdx]["ip_diffim_forced_PsfFlux_instFluxErr"] = 1.01
106 if (-0.5 * i)/1.01 < self.config.minAllowedDirectSnr:
107 self.nRemovedNegativeSources += 1
109 # The last 10 sources will all contained trail length measurements,
110 # increasing in size by 1.5 arcseconds. Only the last three will have
111 # lengths which are too long and will be filtered out.
112 self.nFilteredTrailedSources = 0
113 trail_offset = negativeSources_offset + self.nNegativeSources
114 for i, srcIdx in enumerate(range(trail_offset, trail_offset+self.nTrailedSources)):
115 self.diaSourceCat[srcIdx]["ext_trailedSources_Naive_length"] = 1.5*(i+1)/self.pixelScale
116 if 1.5*(i+1) > 36000/3600.0/24.0 * 30.0:
117 self.nFilteredTrailedSources += 1
118 # Setting a combination of flags for filtering in tests
119 self.diaSourceCat[trail_offset+1]["ext_trailedSources_Naive_flag_off_image"] = True
120 self.diaSourceCat[trail_offset+2]["ext_trailedSources_Naive_flag_suspect_long_trail"] = True
121 self.diaSourceCat[trail_offset+2]["ext_trailedSources_Naive_flag_edge"] = True
122 # As only two of these flags are set, the total number of filtered
123 # sources will be self.nFilteredTrailedSources + 2
124 self.nFilteredTrailedSources += 2
126 mjd = 57071.0
127 self.utc_jd = mjd + 2_400_000.5 - 35.0 / (24.0 * 60.0 * 60.0)
129 self.visitInfo = afwImage.VisitInfo(
130 # This incomplete visitInfo is sufficient for testing because the
131 # Python constructor sets all other required values to some
132 # default.
133 exposureTime=30.0,
134 darkTime=3.0,
135 date=dafBase.DateTime(mjd, system=dafBase.DateTime.MJD),
136 boresightRaDec=geom.SpherePoint(0.0, 0.0, geom.degrees),
137 )
139 def test_run_without_filter(self):
140 """Test that when all filters are turned off all sources in the catalog
141 are returned.
142 """
143 self.config.doRemoveSkySources = False
144 self.config.badFlagList = []
145 self.config.doRemoveNegativeDirectImageSources = False
146 self.config.doWriteRejectedSkySources = False
147 self.config.doTrailedSourceFilter = False
148 filterDiaSourceCatalogTask = FilterDiaSourceCatalogTask(config=self.config)
149 result = filterDiaSourceCatalogTask.run(self.diaSourceCat, self.visitInfo)
150 self.assertEqual(len(result.filteredDiaSourceCat), len(self.diaSourceCat))
151 self.assertEqual(len(result.rejectedDiaSources), 0)
152 self.assertEqual(len(self.diaSourceCat), self.nSources)
154 def test_run_with_filter_sky_only(self):
155 """Test that when only the sky filter is turned on the first five
156 sources which are flagged as sky objects are filtered out of the
157 catalog and the rest are returned.
158 """
159 self.config.doRemoveSkySources = True
160 self.config.badFlagList = []
161 self.config.doRemoveNegativeDirectImageSources = False
162 self.config.doWriteRejectedSkySources = True
163 self.config.doTrailedSourceFilter = False
164 filterDiaSourceCatalogTask = FilterDiaSourceCatalogTask(config=self.config)
165 result = filterDiaSourceCatalogTask.run(self.diaSourceCat, self.visitInfo)
166 nExpectedFilteredSources = self.nSources - self.nSkySources
167 self.assertEqual(len(result.filteredDiaSourceCat),
168 len(self.diaSourceCat[~self.diaSourceCat['sky_source']]))
169 self.assertEqual(len(result.filteredDiaSourceCat), nExpectedFilteredSources)
170 self.assertEqual(len(result.rejectedDiaSources), self.nSkySources)
171 self.assertEqual(len(self.diaSourceCat), self.nSources)
173 def test_run_with_filter_defaultBadFlagList_only(self):
174 """Test that when only the CR center filter is turned on the six sources which are flagged
175 as base_PixelFlags_flag_crCenter are filtered out of the catalog and the rest are returned.
176 """
177 self.config.doRemoveSkySources = False
178 self.config.doRemoveNegativeDirectImageSources = False
179 self.config.doWriteRejectedSkySources = False
180 self.config.doTrailedSourceFilter = False
181 filterDiaSourceCatalogTask = FilterDiaSourceCatalogTask(config=self.config)
182 result = filterDiaSourceCatalogTask.run(self.diaSourceCat, self.visitInfo)
183 nExpectedFilteredSources = self.nSources - self.nCrCenterSources
184 self.assertEqual(len(result.filteredDiaSourceCat),
185 len(self.diaSourceCat[~self.diaSourceCat['base_PixelFlags_flag_crCenter']]))
186 self.assertEqual(len(result.filteredDiaSourceCat), nExpectedFilteredSources)
187 self.assertEqual(len(result.rejectedDiaSources), self.nCrCenterSources)
188 self.assertEqual(len(self.diaSourceCat), self.nSources)
190 def test_run_with_filter_nonDefaultBadFlagList_only(self):
191 """Test that the badFlagList filters appropriately when it is not when the default configuration.
192 The six sources flagged base_PixelFlags_flag_crCenter and the four sources flagged fakeBadFlag
193 should be filtered out of the catalog and the rest are returned.
194 """
195 self.config.doRemoveSkySources = False
196 self.config.badFlagList = ["base_PixelFlags_flag_crCenter", "fakeBadFlag"]
197 self.config.doRemoveNegativeDirectImageSources = False
198 self.config.doWriteRejectedSkySources = False
199 self.config.doTrailedSourceFilter = False
200 filterDiaSourceCatalogTask = FilterDiaSourceCatalogTask(config=self.config)
201 result = filterDiaSourceCatalogTask.run(self.diaSourceCat, self.visitInfo)
202 nExpectedFilteredSources = self.nSources - self.nCrCenterSources - self.nFakeFlagSources
203 self.assertEqual(len(result.filteredDiaSourceCat),
204 len(self.diaSourceCat[~self.diaSourceCat['base_PixelFlags_flag_crCenter']
205 & ~self.diaSourceCat['fakeBadFlag']]))
206 self.assertEqual(len(result.filteredDiaSourceCat), nExpectedFilteredSources)
207 self.assertEqual(len(result.rejectedDiaSources), self.nCrCenterSources + self.nFakeFlagSources)
208 self.assertEqual(len(self.diaSourceCat), self.nSources)
210 def test_run_with_filter_negative_only(self):
211 """Test that when only the negative filter is turned on then
212 sources which below the negtive snr cut are filtered out of the
213 catalog and the rest are returned.
214 """
215 self.config.doRemoveSkySources = False
216 self.config.badFlagList = []
217 self.config.doRemoveNegativeDirectImageSources = True
218 self.config.doWriteRejectedSkySources = True
219 self.config.doTrailedSourceFilter = False
220 filterDiaSourceCatalogTask = FilterDiaSourceCatalogTask(config=self.config)
221 result = filterDiaSourceCatalogTask.run(self.diaSourceCat, self.visitInfo)
222 nExpectedFilteredSources = self.nSources - self.nRemovedNegativeSources
223 self.assertEqual(len(result.filteredDiaSourceCat), nExpectedFilteredSources)
224 self.assertEqual(len(result.rejectedDiaSources), self.nRemovedNegativeSources)
225 self.assertEqual(len(self.diaSourceCat), self.nSources)
226 self.assertEqual(np.sum(result.filteredDiaSourceCat['ip_diffim_forced_PsfFlux_instFlux']
227 / result.filteredDiaSourceCat['ip_diffim_forced_PsfFlux_instFluxErr']
228 < self.config.minAllowedDirectSnr), 0)
229 self.assertEqual(np.sum(result.rejectedDiaSources['ip_diffim_forced_PsfFlux_instFlux']
230 / result.rejectedDiaSources['ip_diffim_forced_PsfFlux_instFluxErr']
231 < self.config.minAllowedDirectSnr),
232 self.nRemovedNegativeSources)
234 def test_run_with_filter_negative_and_sky(self):
235 """Test concatenating rejects when both sky and negative filtering
236 are on.
237 """
238 self.config.doRemoveSkySources = True
239 self.config.badFlagList = []
240 self.config.doRemoveNegativeDirectImageSources = True
241 self.config.doWriteRejectedSkySources = True
242 self.config.doTrailedSourceFilter = False
243 filterDiaSourceCatalogTask = FilterDiaSourceCatalogTask(config=self.config)
244 result = filterDiaSourceCatalogTask.run(self.diaSourceCat, self.visitInfo)
245 nExpectedFilteredSources = self.nSources - self.nSkySources - self.nRemovedNegativeSources
246 self.assertEqual(len(result.filteredDiaSourceCat), nExpectedFilteredSources)
247 self.assertEqual(len(result.rejectedDiaSources), self.nSkySources + self.nRemovedNegativeSources)
248 self.assertEqual(len(self.diaSourceCat), self.nSources)
250 def test_run_with_filter_reliability_only(self):
251 """Test that when only the reliability filter is turned on,
252 sources below the reliability threshold are filtered out."""
254 reliability_threshold = 0.7
255 config = FilterDiaSourceReliabilityConfig()
256 config.minReliability = reliability_threshold
258 schema = afwTable.SourceTable.makeMinimalSchema()
259 schema.addField("score", type="F",
260 doc="Reliability of the source")
261 reliabilityCat = afwTable.SourceCatalog(schema)
262 reliabilityCat.reserve(self.nSources)
264 # Set reliability: first half below threshold, second half above
265 for srcIdx in range(self.nSources):
266 reliabilityCat.addNew()
267 if srcIdx < self.nSources // 2:
268 reliabilityCat[srcIdx]["score"] = 0.25
269 else:
270 reliabilityCat[srcIdx]["score"] = 0.95
271 reliabilityCat['id'] = self.diaSourceCat['id']
272 nLowReliability = np.sum(reliabilityCat["score"] < 0.5)
273 filterDiaSourceCatalogTask = FilterDiaSourceReliabilityTask(config=config)
274 result = filterDiaSourceCatalogTask.run(self.diaSourceCat, reliabilityCat)
275 self.assertEqual(len(result.filteredDiaSources), self.nSources - nLowReliability)
276 self.assertEqual(len(result.rejectedDiaSources), nLowReliability)
277 self.assertTrue(np.all(result.filteredDiaSources["reliability"] >= reliability_threshold))
278 self.assertTrue(np.all(result.rejectedDiaSources["reliability"] < reliability_threshold))
280 def test_run_with_filter_trailed_sources_only(self):
281 """Test that when only the trail filter is turned on the correct number
282 of sources are filtered out. The filtered sources should be the last
283 three sources which have long trails, one source where both the suspect
284 trail and edge trail flag are set, and one source where off_image is
285 set. All sky objects should remain in the catalog.
286 """
287 self.config.doRemoveSkySources = False
288 self.config.badFlagList = []
289 self.config.doRemoveNegativeDirectImageSources = False
290 self.config.doWriteRejectedSkySources = False
291 self.config.doTrailedSourceFilter = True
292 filterDiaSourceCatalogTask = FilterDiaSourceCatalogTask(config=self.config)
293 result = filterDiaSourceCatalogTask.run(self.diaSourceCat, self.visitInfo)
294 nExpectedFilteredSources = self.nSources - self.nFilteredTrailedSources
295 self.assertEqual(len(result.filteredDiaSourceCat), nExpectedFilteredSources)
296 self.assertEqual(len(self.diaSourceCat), self.nSources)
298 def test_run_with_all_filters(self):
299 """Test that all sources are filtered out correctly. Only 15 sources
300 should remain in the catalog after filtering.
301 """
302 self.config.doRemoveSkySources = True
303 self.config.doRemoveNegativeDirectImageSources = True
304 self.config.doWriteRejectedSkySources = True
305 self.config.doTrailedSourceFilter = True
306 filterDiaSourceCatalogTask = FilterDiaSourceCatalogTask(config=self.config)
307 result = filterDiaSourceCatalogTask.run(self.diaSourceCat, self.visitInfo)
308 nExpectedFilteredSources = (self.nSources - self.nSkySources
309 - self.nCrCenterSources
310 - self.nFilteredTrailedSources
311 - self.nRemovedNegativeSources)
312 nExpectedRejectedSources = (self.nSkySources
313 + self.nCrCenterSources
314 + self.nRemovedNegativeSources)
315 # 32 total sources
316 # 5 filtered out sky sources
317 # 6 filtered out sources with cosmic ray detections
318 # 2 filtered out negative sources
319 # 4 filtered out trailed sources, 2 with long trails 2 with flags
320 # 15 sources left
321 self.assertEqual(len(result.filteredDiaSourceCat), nExpectedFilteredSources)
322 # 17 sources rejected, 4 trailed sources not included, 13 rejected sources in catalog.
323 self.assertEqual(len(result.rejectedDiaSources), nExpectedRejectedSources)
324 self.assertEqual(len(self.diaSourceCat), self.nSources)
326 def test_pixelScale_calculation(self):
327 """Check the calculation of the pixel scale from the input catalog.
328 """
329 self.config.doTrailedSourceFilter = True
330 filterDiaSourceCatalogTask = FilterDiaSourceCatalogTask(config=self.config)
331 scale = filterDiaSourceCatalogTask._estimate_pixel_scale(self.diaSourceCat)
332 # Should be almost but not actually equal
333 self.assertNotEqual(self.config.estimatedPixelScale, scale)
334 self.assertAlmostEqual(self.config.estimatedPixelScale, scale, places=6)
336 # If the estimatedPixelScale is very different, that value should be
337 # used exactly and it should not raise an error.
338 self.config.estimatedPixelScale = 1.2
339 filterDiaSourceCatalogTask = FilterDiaSourceCatalogTask(config=self.config)
340 scale = filterDiaSourceCatalogTask._estimate_pixel_scale(self.diaSourceCat)
341 self.assertEqual(self.config.estimatedPixelScale, scale)
344class MemoryTester(lsst.utils.tests.MemoryTestCase):
345 pass
348def setup_module(module):
349 lsst.utils.tests.init()
352if __name__ == "__main__": 352 ↛ 353line 352 didn't jump to line 353 because the condition on line 352 was never true
353 lsst.utils.tests.init()
354 unittest.main()