Coverage for tests / test_diaForcedSource.py: 18%
110 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-18 09:05 +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/>.
22import numpy as np
23import pandas as pd
24import unittest
25import unittest.mock
27from lsst.afw.cameraGeom.testUtils import DetectorWrapper
28import lsst.afw.geom as afwGeom
29import lsst.afw.image as afwImage
30import lsst.daf.base as dafBase
31import lsst.meas.algorithms as measAlg
32from lsst.ap.association import DiaForcedSourceTask
33from lsst.meas.base import IdGenerator
34import lsst.utils.tests
37def create_test_dia_objects(n_points, wcs, startPos=100):
38 """Create dummy DIASources or DIAObjects for use in our tests.
40 Parameters
41 ----------
42 n_points : `int`
43 Number of DiaObject test points to create.
44 wcs : `lsst.afw.geom.SkyWcs`
45 Wcs to convert RA/Dec to pixel x/y.
46 startPos : `int`
47 Start position to iterate from when creating test DiaObjects
49 Returns
50 -------
51 test_points : `pandas.DataFrame`
52 Catalog of points to test.
53 """
54 ids = np.arange(n_points, dtype=np.int64)
55 points = [wcs.pixelToSky(startPos + src_idx, startPos + src_idx) for src_idx in ids]
56 ra = np.array([point.getRa().asDegrees() for point in points], dtype=float)
57 dec = np.array([point.getDec().asDegrees() for point in points], dtype=float)
58 rng = np.random.default_rng(1234)
59 maxSources = 10
60 nDiaSources = rng.integers(1, maxSources, size=n_points)
62 objects = pd.DataFrame({
63 "diaObjectId": ids,
64 "ra": ra,
65 "dec": dec,
66 "nDiaSources": nDiaSources,
67 })
69 return objects
72class TestDiaForcedSource(unittest.TestCase):
74 def setUp(self):
75 # metadata taken from CFHT data
76 # v695856-e0/v695856-e0-c000-a00.sci_img.fits
77 self.metadata = dafBase.PropertySet()
79 self.metadata.set("SIMPLE", "T")
80 self.metadata.set("BITPIX", -32)
81 self.metadata.set("NAXIS", 2)
82 self.metadata.set("NAXIS1", 1024)
83 self.metadata.set("NAXIS2", 1153)
84 self.metadata.set("RADECSYS", 'FK5')
85 self.metadata.set("EQUINOX", 2000.)
87 self.metadata.setDouble("CRVAL1", 215.604025685476)
88 self.metadata.setDouble("CRVAL2", 53.1595451514076)
89 self.metadata.setDouble("CRPIX1", 1109.99981456774)
90 self.metadata.setDouble("CRPIX2", 560.018167811613)
91 self.metadata.set("CTYPE1", 'RA---SIN')
92 self.metadata.set("CTYPE2", 'DEC--SIN')
94 self.metadata.setDouble("CD1_1", 5.10808596133527E-05)
95 self.metadata.setDouble("CD1_2", 1.85579539217196E-07)
96 self.metadata.setDouble("CD2_2", -5.10281493481982E-05)
97 self.metadata.setDouble("CD2_1", -8.27440751733828E-07)
99 self.wcs = afwGeom.makeSkyWcs(self.metadata)
101 self.calibration = 10000
102 self.calibrationErr = 100
103 self.exposureId = 1234
104 self.exposureTime = 200.
105 self.imageSize = [1024, 1153]
106 self.dateTime = "2014-05-13T17:00:00.000000000"
108 # Make images with one source in them and distinct values and
109 # variance for each image.
110 # Direct Image
111 source_image = afwImage.MaskedImageF(
112 lsst.geom.ExtentI(self.imageSize[0] + 1, self.imageSize[1] + 1))
113 source_image.image[100, 100, afwImage.LOCAL] = 10
114 source_image.getVariance().set(1)
115 bbox = lsst.geom.BoxI(
116 lsst.geom.PointI(1, 1),
117 lsst.geom.ExtentI(self.imageSize[0],
118 self.imageSize[1]))
119 masked_image = afwImage.MaskedImageF(source_image, bbox, afwImage.LOCAL)
120 self.exposure = afwImage.makeExposure(masked_image, self.wcs)
122 detector = DetectorWrapper(
123 id=23, bbox=self.exposure.getBBox()).detector
124 visit = afwImage.VisitInfo(
125 exposureTime=self.exposureTime,
126 date=dafBase.DateTime(self.dateTime,
127 dafBase.DateTime.Timescale.TAI))
128 self.exposure.info.id = self.exposureId
129 self.exposure.setDetector(detector)
130 self.exposure.info.setVisitInfo(visit)
131 self.exposure.setFilter(afwImage.FilterLabel(band='g', physical='g.MP9401'))
132 self.exposure.setPhotoCalib(afwImage.PhotoCalib(self.calibration, self.calibrationErr))
134 # Difference Image
135 source_image = afwImage.MaskedImageF(
136 lsst.geom.ExtentI(self.imageSize[0] + 1, self.imageSize[1] + 1))
137 source_image.image[100, 100, afwImage.LOCAL] = 20
138 source_image.getVariance().set(2)
139 bbox = lsst.geom.BoxI(
140 lsst.geom.PointI(1, 1),
141 lsst.geom.ExtentI(self.imageSize[0],
142 self.imageSize[1]))
143 masked_image = afwImage.MaskedImageF(source_image, bbox, afwImage.LOCAL)
144 self.diffim = afwImage.makeExposure(masked_image, self.wcs)
145 self.diffim.info.id = self.exposureId
146 self.diffim.setDetector(detector)
147 self.diffim.info.setVisitInfo(visit)
148 self.diffim.setFilter(afwImage.FilterLabel(band='g', physical='g.MP9401'))
149 self.diffim.setPhotoCalib(afwImage.PhotoCalib(self.calibration, self.calibrationErr))
151 FWHM = 5
152 psf = measAlg.DoubleGaussianPsf(15, 15, FWHM/(2*np.sqrt(2*np.log(2))))
153 self.exposure.setPsf(psf)
154 self.diffim.setPsf(psf)
156 self.testDiaObjects = create_test_dia_objects(5, self.wcs)
157 # Add additional diaObjects that are outside of the above difference
158 # and calexp visit images.
159 objects = [
160 (10000000, self.wcs.pixelToSky(-100000, -100000)), # xy outside
161 (10000001, self.wcs.pixelToSky(100, -100000)), # y outside
162 (10000002, self.wcs.pixelToSky(-100000, 100)), # x outside
163 ]
164 rng = np.random.default_rng(1234)
165 extra = pd.DataFrame({
166 "diaObjectId": np.array([id for id, point in objects], dtype=np.int64),
167 "ra": [point.getRa().asDegrees() for id, point in objects],
168 "dec": [point.getDec().asDegrees() for id, point in objects],
169 "nDiaSources": rng.integers(1, 5, size=len(objects)),
170 })
171 self.testDiaObjects = pd.concat([self.testDiaObjects, extra], ignore_index=True)
172 # Ids of objects that were "updated" during "ap_association"
173 # processing.
174 self.updatedTestIds = np.array([1, 2, 3, 4, 10000001], dtype=np.uint64)
175 # Expecdted number of sources is the number of updated ids plus
176 # any that are within the CCD footprint but are not in the
177 # above list of ids.
178 self.expectedDiaForcedSources = 6
180 self.expected_n_columns = 14
182 def testRun(self):
183 """Test that forced source catalogs are successfully created and have
184 sensible values.
185 """
186 test_objects = self.testDiaObjects.copy()
187 test_objects.set_index("diaObjectId", inplace=True, drop=False)
188 dfs = DiaForcedSourceTask()
189 dia_forced_sources = dfs.run(
190 test_objects, self.updatedTestIds, self.exposure, self.diffim, IdGenerator())
192 direct_values = [199854.48417094944, 160097.40719241602,
193 82299.17897267535, 27148.604434624354,
194 5746.988388215507]
195 direct_var = [75214.39252] * 5
196 diff_values = [399708.9683418989, 320194.81438483205,
197 164598.3579453507, 54297.20886924871,
198 11493.976776431015]
199 diff_var = [106369.21398772747] * 5
201 # Should be number of test objects minus one as one object is purposely
202 # outside of the ccd area.
203 self.assertEqual(len(dia_forced_sources), self.expectedDiaForcedSources)
204 self.assertEqual(len(dia_forced_sources.columns),
205 self.expected_n_columns)
207 for (diaFS_id, diaFS), dirVal, diffVal, dirVar, diffVar in zip(dia_forced_sources.iterrows(),
208 direct_values,
209 diff_values,
210 direct_var,
211 diff_var):
212 self.assertAlmostEqual(diaFS["psfFlux"] / diffVal, 1.)
213 self.assertAlmostEqual(diaFS["psfFluxErr"] / diffVar, 1.)
215 self.assertAlmostEqual(diaFS["scienceFlux"] / dirVal, 1.)
216 self.assertAlmostEqual(diaFS["scienceFluxErr"] / dirVar, 1.)
218 self.assertEqual(diaFS["visit"], self.exposure.visitInfo.id)
219 self.assertEqual(diaFS["detector"], self.exposure.detector.getId())
222class MemoryTester(lsst.utils.tests.MemoryTestCase):
223 pass
226def setup_module(module):
227 lsst.utils.tests.init()
230if __name__ == "__main__": 230 ↛ 231line 230 didn't jump to line 231 because the condition on line 230 was never true
231 lsst.utils.tests.init()
232 unittest.main()