Coverage for tests / test_referenceObjectLoader.py: 16%
232 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-28 08:52 +0000
1# This file is part of meas_algorithms.
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 os.path
23import tempfile
24import unittest
25import glob
27import numpy as np
28from smatch.matcher import sphdist
29import astropy.time
31import lsst.daf.butler
32import lsst.afw.geom as afwGeom
33import lsst.afw.table as afwTable
34from lsst.daf.butler import DatasetType, DeferredDatasetHandle
35from lsst.daf.butler.script import ingest_files
36from lsst.meas.algorithms import (ConvertReferenceCatalogTask, ReferenceObjectLoader,
37 getRefFluxField, getRefFluxKeys)
38from lsst.meas.algorithms.testUtils import MockReferenceObjectLoaderFromFiles
39from lsst.meas.algorithms.convertReferenceCatalog import _makeSchema
40import lsst.utils
41import lsst.geom
43import convertReferenceCatalogTestBase
46class ReferenceObjectLoaderGenericTests(lsst.utils.tests.TestCase):
47 """Test parts of the reference loader that don't depend on loading a
48 catalog, for example schema creation, filter maps, units, and metadata.
49 """
50 def testFilterMapVsAnyFilterMapsToThis(self):
51 config = ReferenceObjectLoader.ConfigClass()
52 # check that a filterMap-only config passes validation
53 config.filterMap = {"b": "a"}
54 try:
55 config.validate()
56 except lsst.pex.config.FieldValidationError:
57 self.fail("`filterMap`-only LoadReferenceObjectsConfig should not fail validation.")
59 # anyFilterMapsToThis and filterMap are mutually exclusive
60 config.anyFilterMapsToThis = "c"
61 with self.assertRaises(lsst.pex.config.FieldValidationError):
62 config.validate()
64 # check that a anyFilterMapsToThis-only config passes validation
65 config.filterMap = {}
66 try:
67 config.validate()
68 except lsst.pex.config.FieldValidationError:
69 self.fail("`anyFilterMapsToThis`-only LoadReferenceObjectsConfig should not fail validation.")
71 def testFilterAliasMap(self):
72 """Make a schema with filter aliases."""
73 for filterMap in ({}, {"camr": "r"}):
74 config = ReferenceObjectLoader.ConfigClass()
75 config.filterMap = filterMap
76 loader = ReferenceObjectLoader(None, [], name=None, config=config)
77 refSchema = _makeSchema(filterNameList="r")
78 loader._addFluxAliases(refSchema,
79 anyFilterMapsToThis=config.anyFilterMapsToThis,
80 filterMap=config.filterMap)
82 self.assertIn("r_flux", refSchema)
83 self.assertIn("r_fluxErr", refSchema)
85 # camera filters aliases are named <filter>_camFlux
86 if "camr" in filterMap:
87 self.assertEqual(getRefFluxField(refSchema, "camr"), "camr_camFlux")
88 else:
89 with self.assertRaisesRegex(RuntimeError,
90 r"Could not find flux field\(s\) camr_camFlux, camr_flux"):
91 getRefFluxField(refSchema, "camr")
93 refCat = afwTable.SimpleCatalog(refSchema)
94 refObj = refCat.addNew()
95 refObj["r_flux"] = 1.23
96 self.assertAlmostEqual(refCat[0].get(getRefFluxField(refSchema, "r")), 1.23)
97 if "camr" in filterMap:
98 self.assertAlmostEqual(refCat[0].get(getRefFluxField(refSchema, "camr")), 1.23)
99 refObj["r_fluxErr"] = 0.111
100 if "camr" in filterMap:
101 self.assertEqual(refCat[0].get("camr_camFluxErr"), 0.111)
102 fluxKey, fluxErrKey = getRefFluxKeys(refSchema, "r")
103 self.assertEqual(refCat[0].get(fluxKey), 1.23)
104 self.assertEqual(refCat[0].get(fluxErrKey), 0.111)
105 if "camr" in filterMap:
106 fluxKey, fluxErrKey = getRefFluxKeys(refSchema, "camr")
107 self.assertEqual(refCat[0].get(fluxErrKey), 0.111)
108 else:
109 with self.assertRaises(RuntimeError):
110 getRefFluxKeys(refSchema, "camr")
112 def testAnyFilterMapsToThisAlias(self):
113 # test anyFilterMapsToThis
114 config = ReferenceObjectLoader.ConfigClass()
115 config.anyFilterMapsToThis = "gg"
116 loader = ReferenceObjectLoader(None, [], name=None, config=config)
117 refSchema = _makeSchema(filterNameList=["gg"])
118 loader._addFluxAliases(refSchema,
119 anyFilterMapsToThis=config.anyFilterMapsToThis,
120 filterMap=config.filterMap)
121 self.assertEqual(getRefFluxField(refSchema, "r"), "gg_flux")
122 # raise if "gg" is not in the refcat filter list
123 with self.assertRaises(RuntimeError):
124 refSchema = _makeSchema(filterNameList=["rr"])
125 refSchema = loader._addFluxAliases(refSchema,
126 anyFilterMapsToThis=config.anyFilterMapsToThis,
127 filterMap=config.filterMap)
129 def testGetMetadataCircle(self):
130 center = lsst.geom.SpherePoint(100*lsst.geom.degrees, 45*lsst.geom.degrees)
131 radius = lsst.geom.Angle(1*lsst.geom.degrees)
132 loader = ReferenceObjectLoader(None, [], name=None)
133 metadata = loader.getMetadataCircle(center, radius, "fakeR")
134 self.assertEqual(metadata['RA'], center.getLongitude().asDegrees())
135 self.assertEqual(metadata['DEC'], center.getLatitude().asDegrees())
136 self.assertEqual(metadata['RADIUS'], radius.asDegrees())
137 self.assertEqual(metadata['SMATCHV'], 2)
138 self.assertEqual(metadata['FILTER'], 'fakeR')
139 self.assertEqual(metadata['JEPOCH'], None)
140 self.assertEqual(metadata['TIMESYS'], 'TAI')
142 epoch = astropy.time.Time(2023.0, format="jyear", scale="tai")
143 metadata = loader.getMetadataCircle(center, radius, "fakeR", epoch=epoch)
144 self.assertEqual(metadata['JEPOCH'], epoch.jyear)
147class ReferenceObjectLoaderLoadTests(convertReferenceCatalogTestBase.ConvertReferenceCatalogTestBase,
148 lsst.utils.tests.TestCase):
149 """Tests of loading reference catalogs, using an in-memory generated fake
150 sky catalog that is converted to an LSST refcat.
152 This effectively is a partial integration test of the refcat conversion,
153 ingestion, and loading sequence, focusing mostly on testing the different
154 ways to load a refcat. It significantly overlaps in coverage with
155 ``nopytest_convertReferenceCatalog.py``, but uses a very trivial test
156 refcat and only one core during the conversion.
157 """
158 @classmethod
159 def setUpClass(cls):
160 super().setUpClass()
162 # Generate a catalog, with arbitrary ids
163 inTempDir = tempfile.TemporaryDirectory()
164 inPath = inTempDir.name
165 skyCatalogFile, _, skyCatalog = cls.makeSkyCatalog(inPath, idStart=25, seed=123)
167 cls.skyCatalog = skyCatalog
169 # override some field names.
170 config = convertReferenceCatalogTestBase.makeConvertConfig(withRaDecErr=True, withMagErr=True,
171 withPm=True, withParallax=True,
172 withFullPositionInformation=True)
173 # use a very small HTM pixelization depth
174 depth = 2
175 config.dataset_config.indexer.active.depth = depth
176 # np.savetxt prepends '# ' to the header lines, so use a reader that understands that
177 config.file_reader.format = 'ascii.commented_header'
178 config.n_processes = 1
179 config.id_name = 'id' # Use the ids from the generated catalogs
180 cls.repoTempDir = tempfile.TemporaryDirectory()
181 repoPath = cls.repoTempDir.name
183 # Convert the input data files to our HTM indexed format.
184 dataTempDir = tempfile.TemporaryDirectory()
185 dataPath = dataTempDir.name
186 converter = ConvertReferenceCatalogTask(output_dir=dataPath, config=config)
187 converter.run([skyCatalogFile])
189 # Make a temporary butler to ingest them into.
190 butler = cls.makeTemporaryRepo(repoPath, config.dataset_config.indexer.active.depth)
191 cls.enterClassContext(butler)
192 dimensions = [f"htm{depth}"]
193 datasetType = DatasetType(config.dataset_config.ref_dataset_name,
194 dimensions,
195 "SimpleCatalog",
196 universe=butler.dimensions,
197 isCalibration=False)
198 butler.registry.registerDatasetType(datasetType)
200 # Ingest the files into the new butler.
201 run = "testingRun"
202 htmTableFile = os.path.join(dataPath, "filename_to_htm.ecsv")
203 ingest_files(repoPath,
204 config.dataset_config.ref_dataset_name,
205 run,
206 htmTableFile,
207 transfer="auto")
209 # Test if we can get back the catalogs, with a new butler.
210 butler = lsst.daf.butler.Butler.from_config(repoPath)
211 cls.enterClassContext(butler)
212 datasetRefs = list(butler.registry.queryDatasets(config.dataset_config.ref_dataset_name,
213 collections=[run]).expanded())
214 handles = []
215 for dataRef in datasetRefs:
216 handles.append(DeferredDatasetHandle(butler=butler, ref=dataRef, parameters=None))
218 cls.datasetRefs = datasetRefs
219 cls.handles = handles
221 inTempDir.cleanup()
222 dataTempDir.cleanup()
224 def test_loadSkyCircle(self):
225 """Test the loadSkyCircle routine."""
226 loader = ReferenceObjectLoader([dataRef.dataId for dataRef in self.datasetRefs],
227 self.handles,
228 name="testrefcat")
229 center = lsst.geom.SpherePoint(180.0*lsst.geom.degrees, 0.0*lsst.geom.degrees)
230 cat = loader.loadSkyCircle(
231 center,
232 30.0*lsst.geom.degrees,
233 filterName='a',
234 ).refCat
235 # Check that the max distance is less than the radius
236 dist = sphdist(180.0, 0.0, np.rad2deg(cat['coord_ra']), np.rad2deg(cat['coord_dec']))
237 self.assertLess(np.max(dist), 30.0)
239 # Check that all the objects from the two catalogs are here.
240 dist = sphdist(180.0, 0.0, self.skyCatalog['ra'], self.skyCatalog['dec'])
241 inside, = (dist < 30.0).nonzero()
242 self.assertEqual(len(cat), len(inside))
244 self.assertTrue(cat.isContiguous())
245 self.assertEqual(len(np.unique(cat['id'])), len(cat))
246 # A default-loaded sky circle should not have centroids
247 self.assertNotIn('centroid_x', cat.schema)
248 self.assertNotIn('centroid_y', cat.schema)
249 self.assertNotIn('hasCentroid', cat.schema)
251 def test_loadPixelBox(self):
252 """Test the loadPixelBox routine."""
253 # This will create a box 50 degrees on a side.
254 loaderConfig = ReferenceObjectLoader.ConfigClass()
255 loaderConfig.pixelMargin = 0
256 loader = ReferenceObjectLoader([dataRef.dataId for dataRef in self.datasetRefs],
257 self.handles,
258 name="testrefcat",
259 config=loaderConfig)
260 bbox = lsst.geom.Box2I(corner=lsst.geom.Point2I(0, 0), dimensions=lsst.geom.Extent2I(1000, 1000))
261 crpix = lsst.geom.Point2D(500, 500)
262 crval = lsst.geom.SpherePoint(180.0*lsst.geom.degrees, 0.0*lsst.geom.degrees)
263 cdMatrix = afwGeom.makeCdMatrix(scale=0.05*lsst.geom.degrees)
264 wcs = afwGeom.makeSkyWcs(crpix=crpix, crval=crval, cdMatrix=cdMatrix)
266 cat = loader.loadPixelBox(bbox, wcs, 'a', bboxToSpherePadding=0).refCat
268 # This is a sanity check on the ranges; the exact selection depends
269 # on cos(dec) and the tangent-plane projection.
270 self.assertLess(np.max(np.rad2deg(cat['coord_ra'])), 180.0 + 25.0)
271 self.assertGreater(np.max(np.rad2deg(cat['coord_ra'])), 180.0 - 25.0)
272 self.assertLess(np.max(np.rad2deg(cat['coord_dec'])), 25.0)
273 self.assertGreater(np.min(np.rad2deg(cat['coord_dec'])), -25.0)
275 # The following is to ensure the reference catalog coords are
276 # getting corrected for proper motion when an epoch is provided.
277 # Use an extreme epoch so that differences in corrected coords
278 # will be significant. Note that this simply tests that the coords
279 # do indeed change when the epoch is passed. It makes no attempt
280 # at assessing the correctness of the change. This is left to the
281 # explicit testProperMotion() test below.
282 catWithEpoch = loader.loadPixelBox(
283 bbox,
284 wcs,
285 'a',
286 bboxToSpherePadding=0,
287 epoch=astropy.time.Time(30000, format='mjd', scale='tai')).refCat
289 self.assertFloatsNotEqual(cat['coord_ra'], catWithEpoch['coord_ra'], rtol=1.0e-4)
290 self.assertFloatsNotEqual(cat['coord_dec'], catWithEpoch['coord_dec'], rtol=1.0e-4)
292 def test_filterMap(self):
293 """Test filterMap parameters."""
294 loaderConfig = ReferenceObjectLoader.ConfigClass()
295 loaderConfig.filterMap = {'aprime': 'a'}
296 loader = ReferenceObjectLoader([dataRef.dataId for dataRef in self.datasetRefs],
297 self.handles,
298 name="testrefcat",
299 config=loaderConfig)
300 center = lsst.geom.SpherePoint(180.0*lsst.geom.degrees, 0.0*lsst.geom.degrees)
301 result = loader.loadSkyCircle(
302 center,
303 30.0*lsst.geom.degrees,
304 filterName='aprime',
305 )
306 self.assertEqual(result.fluxField, 'aprime_camFlux')
307 self.assertFloatsEqual(result.refCat['aprime_camFlux'], result.refCat['a_flux'])
309 def test_properMotion(self):
310 """Test proper motion correction."""
311 loaderConfig = ReferenceObjectLoader.ConfigClass()
312 loaderConfig.filterMap = {'aprime': 'a'}
313 loader = ReferenceObjectLoader([dataRef.dataId for dataRef in self.datasetRefs],
314 self.handles,
315 name="testrefcat",
316 config=loaderConfig)
317 center = lsst.geom.SpherePoint(180.0*lsst.geom.degrees, 0.0*lsst.geom.degrees)
318 cat = loader.loadSkyCircle(
319 center,
320 30.0*lsst.geom.degrees,
321 filterName='a'
322 ).refCat
324 # Zero epoch change --> no proper motion correction (except minor numerical effects)
325 cat_pm = loader.loadSkyCircle(
326 center,
327 30.0*lsst.geom.degrees,
328 filterName='a',
329 epoch=self.epoch
330 ).refCat
332 self.assertFloatsAlmostEqual(cat_pm['coord_ra'], cat['coord_ra'], rtol=1.0e-14)
333 self.assertFloatsAlmostEqual(cat_pm['coord_dec'], cat['coord_dec'], rtol=1.0e-14)
334 self.assertFloatsEqual(cat_pm['coord_raErr'], cat['coord_raErr'])
335 self.assertFloatsEqual(cat_pm['coord_decErr'], cat['coord_decErr'])
337 # One year difference
338 cat_pm = loader.loadSkyCircle(
339 center,
340 30.0*lsst.geom.degrees,
341 filterName='a',
342 epoch=self.epoch + 1.0*astropy.units.yr
343 ).refCat
345 self.assertFloatsEqual(cat_pm['pm_raErr'], cat['pm_raErr'])
346 self.assertFloatsEqual(cat_pm['pm_decErr'], cat['pm_decErr'])
348 separations = np.array([cat[i].getCoord().separation(cat_pm[i].getCoord()).asArcseconds()
349 for i in range(len(cat))])
350 bearings = np.array([cat[i].getCoord().bearingTo(cat_pm[i].getCoord()).asArcseconds()
351 for i in range(len(cat))])
352 self.assertFloatsAlmostEqual(separations, self.properMotionAmt.asArcseconds(), rtol=1.0e-10)
353 self.assertFloatsAlmostEqual(bearings, self.properMotionDir.asArcseconds(), rtol=1.0e-10)
355 predictedRaErr = np.hypot(cat["coord_raErr"], cat["pm_raErr"])
356 predictedDecErr = np.hypot(cat["coord_decErr"], cat["pm_decErr"])
357 self.assertFloatsAlmostEqual(cat_pm["coord_raErr"], predictedRaErr)
358 self.assertFloatsAlmostEqual(cat_pm["coord_decErr"], predictedDecErr)
360 # One year negative difference. This demonstrates a fix for DM-38808,
361 # when the refcat epoch is later in time than the data.
362 cat_pm = loader.loadSkyCircle(
363 center,
364 30.0*lsst.geom.degrees,
365 filterName='a',
366 epoch=self.epoch - 1.0*astropy.units.yr
367 ).refCat
369 self.assertFloatsEqual(cat_pm['pm_raErr'], cat['pm_raErr'])
370 self.assertFloatsEqual(cat_pm['pm_decErr'], cat['pm_decErr'])
372 separations = np.array([cat[i].getCoord().separation(cat_pm[i].getCoord()).asArcseconds()
373 for i in range(len(cat))])
374 bearings = np.array([cat[i].getCoord().bearingTo(cat_pm[i].getCoord()).asArcseconds()
375 for i in range(len(cat))])
376 reverse_proper_motion_dir = self.properMotionDir + 180 * lsst.geom.degrees
377 self.assertFloatsAlmostEqual(separations, self.properMotionAmt.asArcseconds(), rtol=1.0e-10)
378 self.assertFloatsAlmostEqual(bearings, reverse_proper_motion_dir.asArcseconds(), rtol=1.0e-10)
380 predictedRaErr = np.hypot(cat["coord_raErr"], cat["pm_raErr"])
381 predictedDecErr = np.hypot(cat["coord_decErr"], cat["pm_decErr"])
382 self.assertFloatsAlmostEqual(cat_pm["coord_raErr"], predictedRaErr)
383 self.assertFloatsAlmostEqual(cat_pm["coord_decErr"], predictedDecErr)
385 def test_requireProperMotion(self):
386 """Tests of the requireProperMotion config field."""
387 loaderConfig = ReferenceObjectLoader.ConfigClass()
388 loaderConfig.requireProperMotion = True
389 loader = ReferenceObjectLoader([dataRef.dataId for dataRef in self.datasetRefs],
390 self.handles,
391 name="testrefcat",
392 config=loaderConfig)
393 center = lsst.geom.SpherePoint(180.0*lsst.geom.degrees, 0.0*lsst.geom.degrees)
395 # Test that we require an epoch set.
396 msg = 'requireProperMotion=True but epoch not provided to loader'
397 with self.assertRaisesRegex(RuntimeError, msg):
398 loader.loadSkyCircle(
399 center,
400 30.0*lsst.geom.degrees,
401 filterName='a'
402 )
405class Version0Version1ReferenceObjectLoaderTestCase(lsst.utils.tests.TestCase):
406 """Test cases for reading version 0 and version 1 catalogs."""
407 def testLoadVersion0(self):
408 """Attempting to read version 0 refcats should raise.
409 """
410 path = os.path.join(
411 os.path.dirname(os.path.abspath(__file__)),
412 'data',
413 'version0',
414 'ref_cats',
415 'cal_ref_cat'
416 )
418 filenames = sorted(glob.glob(os.path.join(path, '????.fits')))
420 loader = MockReferenceObjectLoaderFromFiles(filenames, name='cal_ref_cat', htmLevel=4)
421 with self.assertRaisesRegex(ValueError, "Version 0 refcats are no longer supported"):
422 loader.loadSkyCircle(convertReferenceCatalogTestBase.make_coord(10, 20),
423 5*lsst.geom.degrees,
424 'a')
426 def testLoadVersion1(self):
427 """Test reading a format_version=1 catalog (fluxes unchanged)."""
428 path = os.path.join(
429 os.path.dirname(os.path.abspath(__file__)),
430 'data',
431 'version1',
432 'ref_cats',
433 'cal_ref_cat'
434 )
436 filenames = sorted(glob.glob(os.path.join(path, '????.fits')))
438 loader = MockReferenceObjectLoaderFromFiles(filenames, name='cal_ref_cat', htmLevel=4)
439 result = loader.loadSkyCircle(convertReferenceCatalogTestBase.make_coord(10, 20),
440 5*lsst.geom.degrees,
441 'a')
443 # version>=1 should not change units on read (they're already nJy).
444 catalog = afwTable.SimpleCatalog.readFits(filenames[0])
445 self.assertFloatsEqual(catalog['a_flux'], result.refCat['a_flux'])
446 self.assertFloatsEqual(catalog['a_fluxErr'], result.refCat['a_fluxErr'])
447 self.assertFloatsEqual(catalog['b_flux'], result.refCat['b_flux'])
448 self.assertFloatsEqual(catalog['b_fluxErr'], result.refCat['b_fluxErr'])
451class TestMemory(lsst.utils.tests.MemoryTestCase):
452 pass
455def setup_module(module):
456 lsst.utils.tests.init()
459if __name__ == "__main__": 459 ↛ 460line 459 didn't jump to line 460 because the condition on line 459 was never true
460 lsst.utils.tests.init()
461 unittest.main()