Coverage for tests/test_astrometryTask.py: 18%
179 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-07 02:18 -0800
1#
2# LSST Data Management System
3# Copyright 2008-2017 LSST Corporation.
4#
5# This product includes software developed by the
6# LSST Project (http://www.lsst.org/).
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the LSST License Statement and
19# the GNU General Public License along with this program. If not,
20# see <http://www.lsstcorp.org/LegalNotices/>.
21#
23import os.path
24import math
25import unittest
26import glob
28import astropy.units as u
29import scipy.stats
30import numpy as np
32import lsst.utils.tests
33import lsst.geom
34import lsst.afw.geom as afwGeom
35import lsst.afw.table as afwTable
36import lsst.afw.image as afwImage
37import lsst.meas.base as measBase
38from lsst.meas.algorithms.testUtils import MockReferenceObjectLoaderFromFiles
39from lsst.meas.astrom import AstrometryTask, exceptions
40import lsst.pipe.base as pipeBase
43class TestAstrometricSolver(lsst.utils.tests.TestCase):
45 def setUp(self):
46 refCatDir = os.path.join(os.path.dirname(__file__), "data", "sdssrefcat")
48 self.bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(3001, 3001))
49 crpix = lsst.geom.Box2D(self.bbox).getCenter()
50 self.tanWcs = afwGeom.makeSkyWcs(crpix=crpix,
51 crval=lsst.geom.SpherePoint(215.5, 53.0, lsst.geom.degrees),
52 cdMatrix=afwGeom.makeCdMatrix(scale=5.1e-5*lsst.geom.degrees))
53 self.exposure = afwImage.ExposureF(self.bbox)
54 self.exposure.setWcs(self.tanWcs)
55 self.exposure.info.setVisitInfo(afwImage.VisitInfo(date=lsst.daf.base.DateTime(60000)))
56 self.exposure.setFilter(afwImage.FilterLabel(band="r", physical="rTest"))
57 filenames = sorted(glob.glob(os.path.join(refCatDir, 'ref_cats', 'cal_ref_cat', '??????.fits')))
58 self.refObjLoader = MockReferenceObjectLoaderFromFiles(filenames, htmLevel=8)
60 def testTrivial(self):
61 """Test fit with no distortion
62 """
63 self.doTest(afwGeom.makeIdentityTransform())
65 def testRadial(self):
66 """Test fit with radial distortion
68 The offset comes from the fact that the CCD is not centered
69 """
70 self.doTest(afwGeom.makeRadialTransform([0, 1.01, 1e-7]))
72 def doTest(self, pixelsToTanPixels):
73 """Test using pixelsToTanPixels to distort the source positions.
74 """
75 schema = self._makeSourceCatalogSchema()
76 config = AstrometryTask.ConfigClass()
77 config.wcsFitter.order = 3
78 config.wcsFitter.numRejIter = 0
79 task = AstrometryTask(config=config, refObjLoader=self.refObjLoader, schema=schema)
80 distortedWcs = afwGeom.makeModifiedWcs(pixelTransform=pixelsToTanPixels, wcs=self.tanWcs,
81 modifyActualPixels=False)
82 # Make the source catalog at the distorted positions, but keep the
83 # initial TAN WCS on the exposure, to check that the fitted WCS
84 # is close to the distorted one and different from the input.
85 sourceCat = self.makeSourceCat(distortedWcs, schema)
86 # This test is from before rough magnitude rejection was implemented.
87 config.doMagnitudeOutlierRejection = False
88 results = task.run(sourceCat=sourceCat, exposure=self.exposure)
90 self.assertWcsAlmostEqualOverBBox(distortedWcs, self.exposure.wcs, self.bbox,
91 maxDiffSky=0.002*lsst.geom.arcseconds, maxDiffPix=0.02)
92 # Test that the sources used in the fit are flagged in the catalog.
93 self.assertEqual(sum(sourceCat["calib_astrometry_used"]), len(results.matches))
95 # The following tests that the measured mean astrometry offset is not
96 # equal to its standard deviation (within the default tolerances of
97 # double-precision epsilon). While this condition is not strictly
98 # required to be true, in reality it would be highly unlikely and is
99 # certainly not the case for this test as configured. The motivation
100 # for adding this check is that it would have caught a refactor that
101 # mistakenly set them both as the mean.
102 self.assertFloatsNotEqual(
103 self.exposure.metadata["SFM_ASTROM_OFFSET_MEAN"], self.exposure.metadata["SFM_ASTROM_OFFSET_STD"]
104 )
106 srcCoordKey = afwTable.CoordKey(sourceCat.schema["coord"])
107 refCoordKey = afwTable.CoordKey(results.refCat.schema["coord"])
108 refCentroidKey = afwTable.Point2DKey(results.refCat.schema["centroid"])
109 maxAngSep = 0*lsst.geom.radians
110 maxPixSep = 0
111 for refObj, src, d in results.matches:
112 refCoord = refObj.get(refCoordKey)
113 refPixPos = refObj.get(refCentroidKey)
114 srcCoord = src.get(srcCoordKey)
115 srcPixPos = src.getCentroid()
117 angSep = refCoord.separation(srcCoord)
118 maxAngSep = max(maxAngSep, angSep)
120 pixSep = math.hypot(*(srcPixPos-refPixPos))
121 maxPixSep = max(maxPixSep, pixSep)
122 print("max angular separation = %0.4f arcsec" % (maxAngSep.asArcseconds(),))
123 print("max pixel separation = %0.3f" % (maxPixSep,))
124 self.assertLess(maxAngSep.asArcseconds(), 0.0038)
125 self.assertLess(maxPixSep, 0.021)
127 # try again, invoking the reference selector
128 config.referenceSelector.doUnresolved = True
129 config.referenceSelector.unresolved.name = 'resolved'
130 solverRefSelect = AstrometryTask(config=config, refObjLoader=self.refObjLoader)
131 self.exposure.setWcs(distortedWcs)
132 resultsRefSelect = solverRefSelect.run(
133 sourceCat=sourceCat,
134 exposure=self.exposure,
135 )
136 self.assertLess(len(resultsRefSelect.matches), len(results.matches))
138 # try again, allowing magnitude outlier rejection.
139 config.doMagnitudeOutlierRejection = True
140 solverMagOutlierRejection = AstrometryTask(config=config, refObjLoader=self.refObjLoader)
141 self.exposure.setWcs(distortedWcs)
142 resultsMagOutlierRejection = solverMagOutlierRejection.run(
143 sourceCat=sourceCat,
144 exposure=self.exposure,
145 )
146 self.assertLess(len(resultsMagOutlierRejection.matches), len(resultsRefSelect.matches))
147 config.doMagnitudeOutlierRejection = False
149 # try again, but without fitting the WCS, no reference selector
150 config.referenceSelector.doUnresolved = False
151 config.forceKnownWcs = True
152 solverNoFit = AstrometryTask(config=config, refObjLoader=self.refObjLoader)
153 self.exposure.setWcs(distortedWcs)
154 resultsNoFit = solverNoFit.run(
155 sourceCat=sourceCat,
156 exposure=self.exposure,
157 )
158 self.assertIsNone(resultsNoFit.scatterOnSky)
160 # fitting should result in matches that are at least as good
161 # (strictly speaking fitting might result in a larger match list with
162 # some outliers, but in practice this test passes)
163 meanFitDist = np.mean([match.distance for match in results.matches])
164 meanNoFitDist = np.mean([match.distance for match in resultsNoFit.matches])
165 self.assertLessEqual(meanFitDist, meanNoFitDist)
167 # try once again, without fitting the WCS, with the reference selector
168 # (this goes through a different code path)
169 config.referenceSelector.doUnresolved = True
170 solverNoFitRefSelect = AstrometryTask(config=config, refObjLoader=self.refObjLoader)
171 resultsNoFitRefSelect = solverNoFitRefSelect.run(
172 sourceCat=sourceCat,
173 exposure=self.exposure,
174 )
175 self.assertLess(len(resultsNoFitRefSelect.matches), len(resultsNoFit.matches))
177 @staticmethod
178 def _makeSourceCatalogSchema():
179 """Return a catalog schema with all necessary fields added.
180 """
181 schema = afwTable.SourceTable.makeMinimalSchema()
182 measBase.SingleFrameMeasurementTask(schema=schema) # expand the schema
183 schema.addField("deblend_nChild", type=np.int32,
184 doc="Number of children this object has (defaults to 0)")
185 schema.addField("detect_isPrimary", type=np.int32,
186 doc="true if source has no children and is not a sky source")
187 return schema
189 def makeSourceCat(self, wcs, schema, doScatterCentroids=False):
190 """Make a source catalog by reading the position reference stars using
191 the proviced WCS.
193 Optionally, via doScatterCentroids, add some scatter to the centroids
194 assiged to the source catalog (otherwise they will be identical to
195 those of the reference catalog).
196 """
197 loadRes = self.refObjLoader.loadPixelBox(bbox=self.bbox, wcs=wcs, filterName="r")
198 refCat = loadRes.refCat
200 sourceCat = afwTable.SourceCatalog(schema)
202 sourceCat.resize(len(refCat))
203 scatterFactor = 1.0
204 if doScatterCentroids:
205 np.random.seed(12345)
206 scatterFactor = np.random.uniform(0.999, 1.001, len(sourceCat))
207 sourceCat["slot_Centroid_x"] = scatterFactor*refCat["centroid_x"]
208 sourceCat["slot_Centroid_y"] = scatterFactor*refCat["centroid_y"]
209 sourceCat["slot_PsfFlux_instFlux"] = refCat["r_flux"]
210 sourceCat["slot_PsfFlux_instFluxErr"] = refCat["r_flux"]/100
211 # All of these sources are primary.
212 sourceCat['detect_isPrimary'] = 1
214 # Deliberately add some outliers to check that the magnitude
215 # outlier rejection code is being run.
216 sourceCat["slot_PsfFlux_instFlux"][0: 4] *= 1000.0
218 return sourceCat
220 def testBadAstrometry(self):
221 """Test that an appropriately informative exception is raised for a
222 bad quality fit.
223 """
224 catalog = self.makeSourceCat(self.tanWcs, self._makeSourceCatalogSchema())
225 # Fake match list with 10" match distance for every source to force
226 # a bad quality fit result.
227 matches = [afwTable.ReferenceMatch(r, r, (10*u.arcsecond).to_value(u.radian)) for r in catalog]
228 result = pipeBase.Struct(
229 matches=matches,
230 wcs=self.exposure.wcs,
231 scatterOnSky=20.0*lsst.geom.arcseconds,
232 matchTolerance=None
233 )
234 with unittest.mock.patch("lsst.meas.astrom.AstrometryTask._matchAndFitWcs",
235 return_value=result, autospec=True):
236 with self.assertRaises(exceptions.BadAstrometryFit):
237 task = AstrometryTask(refObjLoader=self.refObjLoader)
238 task.run(catalog, self.exposure)
239 task.check(self.exposure, catalog, len(matches))
240 self.assertIsNone(self.exposure.wcs)
241 self.assertTrue(np.all(np.isnan(catalog["coord_ra"])))
242 self.assertTrue(np.all(np.isnan(catalog["coord_dec"])))
244 def testMatcherFails(self):
245 """Test that a matcher exception has additional metadata attached.
246 """
247 catalog = self.makeSourceCat(self.tanWcs, self._makeSourceCatalogSchema())
248 with unittest.mock.patch("lsst.meas.astrom.AstrometryTask._matchAndFitWcs", autospec=True,
249 side_effect=exceptions.MatcherFailure("some matcher problem")):
250 with self.assertRaises(exceptions.MatcherFailure) as cm:
251 task = AstrometryTask(refObjLoader=self.refObjLoader)
252 task.run(catalog, self.exposure)
253 self.assertEqual(cm.exception.metadata["iterations"], 1)
254 self.assertIsNone(self.exposure.wcs)
255 self.assertTrue(np.all(np.isnan(catalog["coord_ra"])))
256 self.assertTrue(np.all(np.isnan(catalog["coord_dec"])))
258 def testExceptions(self):
259 """Test that the custom astrometry exceptions are well behaved.
260 """
261 error = exceptions.AstrometryError("something", blah=10)
262 self.assertEqual(error.metadata["blah"], 10)
263 self.assertIn("something", str(error))
264 self.assertIn("'blah': 10", str(error))
266 # Metadata cannot contain an astropy unit.
267 error = exceptions.AstrometryError("something", blah=10*u.arcsecond)
268 with self.assertRaisesRegex(TypeError, "blah is of type <class 'astropy.units.quantity.Quantity'>"):
269 error.metadata
272class TestMagnitudeOutliers(lsst.utils.tests.TestCase):
273 def testMagnitudeOutlierRejection(self):
274 """Test rejection of magnitude outliers.
276 This test only tests the outlier rejection, and not any other
277 part of the matching or astrometry fitter.
278 """
279 config = AstrometryTask.ConfigClass()
280 config.doMagnitudeOutlierRejection = True
281 config.magnitudeOutlierRejectionNSigma = 4.0
282 task = AstrometryTask(config=config, refObjLoader=None)
284 nTest = 100
286 refSchema = lsst.afw.table.SimpleTable.makeMinimalSchema()
287 refSchema.addField('refFlux', 'F')
288 refCat = lsst.afw.table.SimpleCatalog(refSchema)
289 refCat.resize(nTest)
291 srcSchema = lsst.afw.table.SourceTable.makeMinimalSchema()
292 srcSchema.addField('srcFlux', 'F')
293 srcCat = lsst.afw.table.SourceCatalog(srcSchema)
294 srcCat.resize(nTest)
296 np.random.seed(12345)
297 refMag = np.full(nTest, 20.0)
298 srcMag = np.random.normal(size=nTest, loc=0.0, scale=1.0)
300 # Determine the sigma of the random sample
301 zp = np.median(refMag[: -4] - srcMag[: -4])
302 sigma = scipy.stats.median_abs_deviation(srcMag[: -4], scale='normal')
304 # Deliberately alter some magnitudes to be outliers.
305 srcMag[-3] = (config.magnitudeOutlierRejectionNSigma + 0.1)*sigma + (20.0 - zp)
306 srcMag[-4] = -(config.magnitudeOutlierRejectionNSigma + 0.1)*sigma + (20.0 - zp)
308 refCat['refFlux'] = (refMag*u.ABmag).to_value(u.nJy)
309 srcCat['srcFlux'] = 10.0**(srcMag/(-2.5))
311 # Deliberately poison some reference fluxes.
312 refCat['refFlux'][-1] = np.inf
313 refCat['refFlux'][-2] = np.nan
315 matchesIn = []
316 for ref, src in zip(refCat, srcCat):
317 matchesIn.append(lsst.afw.table.ReferenceMatch(first=ref, second=src, distance=0.0))
319 matchesOut = task._removeMagnitudeOutliers('srcFlux', 'refFlux', matchesIn)
321 # We should lose the 4 outliers we created.
322 self.assertEqual(len(matchesOut), len(matchesIn) - 4)
325class MemoryTester(lsst.utils.tests.MemoryTestCase):
326 pass
329def setup_module(module):
330 lsst.utils.tests.init()
333if __name__ == "__main__": 333 ↛ 334line 333 didn't jump to line 334 because the condition on line 333 was never true
334 lsst.utils.tests.init()
335 unittest.main()