Coverage for tests/test_butlerFits.py: 15%
269 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-09 03:42 -0700
1# This file is part of obs_base.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (http://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 <http://www.gnu.org/licenses/>.
22from __future__ import annotations
24import os
25import shutil
26import tempfile
27import unittest
28import unittest.mock
29import uuid
31import astropy.table
33import lsst.afw.cameraGeom.testUtils # for test asserts injected into TestCase
34import lsst.afw.image
35import lsst.pex.config
36import lsst.utils.tests
37from lsst.afw.fits import readMetadata
38from lsst.afw.image import LOCAL, PARENT, ExposureFitsReader, MaskedImageFitsReader
39from lsst.afw.math import flipImage
40from lsst.daf.base import PropertyList, PropertySet
41from lsst.daf.butler import Config, DatasetProvenance, DatasetRef, DatasetType, StorageClassFactory
42from lsst.daf.butler.tests import addDatasetType, makeTestCollection, makeTestRepo
43from lsst.geom import Box2I, Extent2I, Point2I
44from lsst.obs.base.exposureAssembler import ExposureAssembler
45from lsst.obs.base.tests import make_ramp_exposure_trimmed, make_ramp_exposure_untrimmed
47TESTDIR = os.path.dirname(__file__)
49BUTLER_CONFIG = """
50storageClasses:
51 ExposureCompositeF:
52 inheritsFrom: ExposureF
53datastore:
54 # Want to check disassembly so can't use InMemory
55 cls: lsst.daf.butler.datastores.fileDatastore.FileDatastore
56 formatters:
57 ExposureCompositeF: lsst.obs.base.formatters.fitsExposure.FitsExposureFormatter
58 lossless:
59 formatter: lsst.obs.base.formatters.fitsExposure.FitsExposureFormatter
60 parameters:
61 recipe: lossless
62 uncompressed:
63 formatter: lsst.obs.base.formatters.fitsExposure.FitsExposureFormatter
64 parameters:
65 recipe: noCompression
66 lossy:
67 formatter: lsst.obs.base.formatters.fitsExposure.FitsExposureFormatter
68 parameters:
69 recipe: lossyBasic
70 composites:
71 disassembled:
72 ExposureCompositeF: True
73 # Always run with cache disabled. None of these datasets are remote but
74 # we want to check that URI reading does work properly.
75 cached:
76 expiry:
77 default: false
78 mode: disabled
79"""
81# Components present in the test file
82COMPONENTS = {
83 "wcs",
84 "image",
85 "mask",
86 "coaddInputs",
87 "psf",
88 "visitInfo",
89 "variance",
90 "metadata",
91 "photoCalib",
92 "filter",
93 "validPolygon",
94 "transmissionCurve",
95 "detector",
96 "apCorrMap",
97 "summaryStats",
98 "id",
99}
100READ_COMPONENTS = {
101 "bbox",
102 "xy0",
103 "dimensions",
104}
107class ButlerFitsTests(lsst.utils.tests.TestCase):
108 """Tests for butler interaction with FITS files."""
110 @classmethod
111 def setUpClass(cls):
112 """Create a new butler once only."""
113 cls.storageClassFactory = StorageClassFactory()
115 cls.root = tempfile.mkdtemp(dir=TESTDIR)
117 dataIds = {
118 "instrument": ["DummyCam"],
119 "physical_filter": ["d-r"],
120 "visit": [42, 43, 44],
121 }
123 # Ensure that we test in a directory that will include some
124 # metacharacters
125 subdir = "sub?#dir"
126 butlerRoot = os.path.join(cls.root, subdir)
128 cls.creatorButler = makeTestRepo(butlerRoot, dataIds, config=Config.fromYaml(BUTLER_CONFIG))
130 # Create dataset types used by the tests
131 for datasetTypeName, storageClassName in (
132 ("calexp", "ExposureF"),
133 ("noise", "MaskedImageF"),
134 ("unknown", "ExposureCompositeF"),
135 ("testCatalog", "SourceCatalog"),
136 ("lossless", "ExposureF"),
137 ("uncompressed", "ExposureF"),
138 ("lossy", "ExposureF"),
139 ):
140 storageClass = cls.storageClassFactory.getStorageClass(storageClassName)
141 addDatasetType(cls.creatorButler, datasetTypeName, set(dataIds), storageClass)
143 # And some dataset types that have no dimensions for easy testing
144 for datasetTypeName, storageClassName in (
145 ("ps", "PropertySet"),
146 ("pl", "PropertyList"),
147 ("int_exp_trimmed", "ExposureI"),
148 ("int_exp_untrimmed", "ExposureI"),
149 ):
150 storageClass = cls.storageClassFactory.getStorageClass(storageClassName)
151 addDatasetType(cls.creatorButler, datasetTypeName, {}, storageClass)
153 @classmethod
154 def tearDownClass(cls):
155 if cls.root is not None:
156 shutil.rmtree(cls.root, ignore_errors=True)
158 def setUp(self):
159 self.butler = makeTestCollection(self.creatorButler, uniqueId=self.id())
161 def makeExampleCatalog(self) -> lsst.afw.table.SourceCatalog:
162 catalogPath = os.path.join(TESTDIR, "data", "source_catalog.fits")
163 return lsst.afw.table.SourceCatalog.readFits(catalogPath)
165 def assertCatalogEqual(
166 self, inputCatalog: lsst.afw.table.SourceCatalog, outputCatalog: lsst.afw.table.SourceCatalog
167 ) -> None:
168 self.assertIsInstance(outputCatalog, lsst.afw.table.SourceCatalog)
169 inputTable = inputCatalog.getTable()
170 inputRecord = inputCatalog[0]
171 outputTable = outputCatalog.getTable()
172 outputRecord = outputCatalog[0]
173 self.assertEqual(inputRecord.getPsfInstFlux(), outputRecord.getPsfInstFlux())
174 self.assertEqual(inputRecord.getPsfFluxFlag(), outputRecord.getPsfFluxFlag())
175 self.assertEqual(
176 inputTable.getSchema().getAliasMap().get("slot_Centroid"),
177 outputTable.getSchema().getAliasMap().get("slot_Centroid"),
178 )
179 self.assertEqual(inputRecord.getCentroid(), outputRecord.getCentroid())
180 self.assertFloatsAlmostEqual(
181 inputRecord.getCentroidErr()[0, 0], outputRecord.getCentroidErr()[0, 0], rtol=1e-6
182 )
183 self.assertFloatsAlmostEqual(
184 inputRecord.getCentroidErr()[1, 1], outputRecord.getCentroidErr()[1, 1], rtol=1e-6
185 )
186 self.assertEqual(
187 inputTable.getSchema().getAliasMap().get("slot_Shape"),
188 outputTable.getSchema().getAliasMap().get("slot_Shape"),
189 )
190 self.assertFloatsAlmostEqual(
191 inputRecord.getShapeErr()[0, 0], outputRecord.getShapeErr()[0, 0], rtol=1e-6
192 )
193 self.assertFloatsAlmostEqual(
194 inputRecord.getShapeErr()[1, 1], outputRecord.getShapeErr()[1, 1], rtol=1e-6
195 )
196 self.assertFloatsAlmostEqual(
197 inputRecord.getShapeErr()[2, 2], outputRecord.getShapeErr()[2, 2], rtol=1e-6
198 )
200 def runFundamentalTypeTest(self, datasetTypeName, entity):
201 """Put and get the supplied entity and compare."""
202 ref = self.butler.put(entity, datasetTypeName)
203 butler_ps = self.butler.get(ref)
204 self.assertEqual(butler_ps, entity)
206 # Break the contact by ensuring that we are writing YAML
207 uri = self.butler.getURI(ref)
208 self.assertEqual(uri.getExtension(), ".yaml", f"Check extension of {uri}")
210 def testFundamentalTypes(self) -> None:
211 """Ensure that some fundamental stack types round trip."""
212 ps = PropertySet()
213 ps["a.b"] = 5
214 ps["c.d.e"] = "string"
215 self.runFundamentalTypeTest("ps", ps)
217 pl = PropertyList()
218 pl["A"] = 1
219 pl.setComment("A", "An int comment")
220 pl["B"] = "string"
221 pl.setComment("B", "A string comment")
222 self.runFundamentalTypeTest("pl", pl)
224 def _make_provenance(self):
225 """Return provenance information for testing."""
226 prov = DatasetProvenance(quantum_id=uuid.uuid4())
227 ref = DatasetRef(
228 datasetType=self.butler.get_dataset_type("ps"),
229 dataId=[],
230 run="run",
231 )
232 prov.add_input(ref)
233 prov.add_extra_provenance(ref.id, {"EXTRA": 42})
234 return prov
236 def testFitsCatalog(self) -> None:
237 """Test reading of a FITS catalog."""
238 catalog = self.makeExampleCatalog()
239 dataId = {"visit": 42, "instrument": "DummyCam", "physical_filter": "d-r"}
240 prov = self._make_provenance()
241 ref = self.butler.put(catalog, "testCatalog", dataId, provenance=prov)
242 stored = self.butler.get(ref)
243 self.assertCatalogEqual(catalog, stored)
244 self.assertEqual(stored.metadata["LSST BUTLER ID"], str(ref.id))
245 self.assertEqual(stored.metadata["LSST BUTLER QUANTUM"], str(prov.quantum_id))
246 self.assertEqual(stored.metadata["LSST BUTLER INPUT 0 RUN"], "run")
247 self.assertEqual(stored.metadata["LSST BUTLER INPUT 0 EXTRA"], 42)
249 # Override the storage class.
250 astropy_table = self.butler.get(ref, storageClass="AstropyTable")
251 self.assertIsInstance(astropy_table, astropy.table.Table)
252 self.assertEqual(len(astropy_table), len(stored))
254 def testExposureCompositePutGetConcrete(self) -> None:
255 """Test composite with no disassembly."""
256 ref = self.runExposureCompositePutGetTest("calexp")
258 uri = self.butler.getURI(ref)
259 self.assertTrue(uri.exists(), f"Checking URI {uri} existence")
261 def testExposureCompositePutGetConcreteAstropy(self) -> None:
262 """Test composite with no disassembly and astropy component read."""
263 with unittest.mock.patch(
264 "lsst.obs.base.formatters.fitsExposure._ALWAYS_USE_ASTROPY_FOR_COMPONENT_READ", True
265 ):
266 self.runExposureCompositePutGetTest("calexp")
268 def testExposureCompositePutGetVirtual(self) -> None:
269 """Testing composite disassembly."""
270 ref = self.runExposureCompositePutGetTest("unknown")
272 primary, components = self.butler.getURIs(ref)
273 self.assertIsNone(primary)
274 self.assertEqual(set(components), COMPONENTS)
275 for compName, uri in components.items():
276 self.assertTrue(uri.exists(), f"Checking URI {uri} existence for component {compName}")
278 def runExposureCompositePutGetTest(self, datasetTypeName: str) -> DatasetRef:
279 example = os.path.join(TESTDIR, "data", "calexp.fits")
280 exposure = lsst.afw.image.ExposureF(example)
282 prov = self._make_provenance()
283 dataId = {"visit": 42, "instrument": "DummyCam", "physical_filter": "d-r"}
284 ref = self.butler.put(exposure, datasetTypeName, dataId, provenance=prov)
286 # Get the full thing
287 composite = self.butler.get(datasetTypeName, dataId)
289 # Check that provenance has been written.
290 self.assertEqual(composite.metadata["LSST BUTLER ID"], str(ref.id))
291 self.assertEqual(composite.metadata["LSST BUTLER QUANTUM"], str(prov.quantum_id))
292 self.assertEqual(composite.metadata["LSST BUTLER INPUT 0 RUN"], "run")
293 self.assertEqual(composite.metadata["LSST BUTLER INPUT 0 EXTRA"], 42)
295 # There is no assert for Exposure so just look at maskedImage
296 self.assertMaskedImagesEqual(composite.maskedImage, exposure.maskedImage)
298 # Helper for extracting components
299 assembler = ExposureAssembler(ref.datasetType.storageClass)
301 # Check all possible components that can be read
302 allComponents = set()
303 allComponents.update(COMPONENTS, READ_COMPONENTS)
305 # Get each component from butler independently
306 for compName in allComponents:
307 compTypeName = DatasetType.nameWithComponent(datasetTypeName, compName)
308 component = self.butler.get(compTypeName, dataId)
310 reference = assembler.getComponent(exposure, compName)
312 self.assertIsInstance(component, type(reference), f"Checking type of component {compName}")
314 if compName in ("image", "variance"):
315 self.assertImagesEqual(component, reference)
316 elif compName == "mask":
317 self.assertMasksEqual(component, reference)
318 elif compName == "wcs":
319 self.assertWcsAlmostEqualOverBBox(component, reference, exposure.getBBox())
320 elif compName == "coaddInputs":
321 self.assertEqual(
322 len(component.visits), len(reference.visits), f"cf visits {component.visits}"
323 )
324 self.assertEqual(len(component.ccds), len(reference.ccds), f"cf CCDs {component.ccds}")
325 elif compName == "psf":
326 # Equality for PSF does not work
327 pass
328 elif compName == "filter":
329 self.assertEqual(component, reference, compName)
330 elif compName == "id":
331 self.assertEqual(component, reference, compName)
332 elif compName == "visitInfo":
333 self.assertEqual(component, reference, "VisitInfo comparison")
334 elif compName == "metadata":
335 # The component metadata has extra fields in it so cannot
336 # compare directly.
337 for k, v in reference.items():
338 self.assertEqual(component[k], v, f"{compName} key: {k}")
339 elif compName == "photoCalib":
340 # This example has a
341 # "spatially constant with mean: inf error: nan" entry
342 # which does not compare directly.
343 self.assertEqual(str(component), str(reference), compName)
344 self.assertIn("spatially constant with mean: 1.99409", str(component), "Checking photoCalib")
345 elif compName in ("bbox", "xy0", "dimensions", "validPolygon"):
346 self.assertEqual(component, reference, compName)
347 elif compName == "apCorrMap":
348 self.assertEqual(set(component.keys()), set(reference.keys()), compName)
349 elif compName == "transmissionCurve":
350 self.assertEqual(
351 component.getThroughputAtBounds(), reference.getThroughputAtBounds(), compName
352 )
353 elif compName == "detector":
354 c_amps = {a.getName() for a in component.getAmplifiers()}
355 r_amps = {a.getName() for a in reference.getAmplifiers()}
356 self.assertEqual(c_amps, r_amps, compName)
357 elif compName == "summaryStats":
358 self.assertEqual(component.psfSigma, reference.psfSigma, compName)
359 else:
360 raise RuntimeError(f"Unexpected component '{compName}' encountered in test")
362 # Full Exposure with parameters
363 for origin in (LOCAL, PARENT):
364 inBBox = Box2I(minimum=Point2I(5, 3), maximum=Point2I(42, 16))
365 parameters = {"bbox": inBBox, "origin": origin}
366 subset = self.butler.get(datasetTypeName, dataId, parameters=parameters)
367 outBBox = subset.getBBox()
368 self.assertEqual(inBBox, outBBox)
369 self.assertImagesEqual(subset.getImage(), exposure.subset(inBBox, origin=LOCAL).getImage())
371 return ref
373 def putFits(self, exposure, datasetTypeName, visit):
374 """Put different datasetTypes and return information."""
375 dataId = {"visit": visit, "instrument": "DummyCam", "physical_filter": "d-r"}
376 refC = self.butler.put(exposure, datasetTypeName, dataId)
377 uriC = self.butler.getURI(refC)
378 stat = os.stat(uriC.ospath)
379 size = stat.st_size
380 # We can't use butler's Exposure storage class metadata component,
381 # because that intentionally strips keywords that science code
382 # shouldn't ever read (to at least weaken our assumptions that we write
383 # to FITS). Instead use a lower-level method on the URI for this test.
384 meta = readMetadata(uriC.ospath)
385 return meta, size
387 def testCompression(self):
388 """Test that we can write compressed and uncompressed FITS."""
389 example = os.path.join(TESTDIR, "data", "small.fits")
390 exposure = lsst.afw.image.ExposureF(example)
392 # Write a lossless compressed
393 metaC, sizeC = self.putFits(exposure, "lossless", 42)
394 self.assertEqual(metaC["TTYPE1"], "COMPRESSED_DATA")
395 self.assertEqual(metaC["ZCMPTYPE"], "GZIP_2")
397 # Write an uncompressed FITS file
398 metaN, sizeN = self.putFits(exposure, "uncompressed", 43)
399 self.assertNotIn("ZCMPTYPE", metaN)
401 # Write an uncompressed FITS file
402 metaL, sizeL = self.putFits(exposure, "lossy", 44)
403 self.assertEqual(metaL["TTYPE1"], "COMPRESSED_DATA")
404 self.assertEqual(metaL["ZCMPTYPE"], "RICE_1")
406 self.assertNotEqual(sizeC, sizeN)
407 # Data file is so small that Lossy and Compressed are dominated
408 # by the extra compression tables
409 self.assertEqual(sizeL, sizeC)
411 def testExposureFormatterAmpParameter(self):
412 """Test the FitsExposureFormatter implementation of the Exposure
413 StorageClass's 'amp' and 'detector' parameters.
414 """
415 # Our example exposure file has a realistic detector (looks like HSC),
416 # but the image itself doesn't match it. So we just load the detector,
417 # and use it to make our own more useful images and put them.
418 detector = ExposureFitsReader(os.path.join(TESTDIR, "data", "calexp.fits")).readDetector()
419 trimmed_full = make_ramp_exposure_trimmed(detector)
420 untrimmed_full = make_ramp_exposure_untrimmed(detector)
421 trimmed_ref = self.butler.put(trimmed_full, "int_exp_trimmed")
422 untrimmed_ref = self.butler.put(untrimmed_full, "int_exp_untrimmed")
423 for n, amp in enumerate(detector):
424 # Try to read each amp as it is on disk, with a variety of
425 # parameters that should all do the same thing.
426 for amp_parameter in [amp, amp.getName(), n]:
427 for parameters in [{"amp": amp_parameter}, {"amp": amp_parameter, "detector": detector}]:
428 with self.subTest(parameters=parameters):
429 test_trimmed = self.butler.get(trimmed_ref, parameters=parameters)
430 test_untrimmed = self.butler.get(untrimmed_ref, parameters=parameters)
431 self.assertImagesEqual(test_trimmed.image, trimmed_full[amp.getBBox()].image)
432 self.assertImagesEqual(test_untrimmed.image, untrimmed_full[amp.getRawBBox()].image)
433 self.assertEqual(len(test_trimmed.getDetector()), 1)
434 self.assertEqual(len(test_untrimmed.getDetector()), 1)
435 self.assertAmplifiersEqual(test_trimmed.getDetector()[0], amp)
436 self.assertAmplifiersEqual(test_untrimmed.getDetector()[0], amp)
437 # Try to read various transformed versions of the original amp,
438 # to make sure flips and offsets are applied correctly.
439 # First flip X only.
440 amp_t1 = amp.rebuild().transform(outFlipX=True).finish()
441 test_t1_trimmed = self.butler.get(trimmed_ref, parameters={"amp": amp_t1})
442 self.assertImagesEqual(
443 test_t1_trimmed.image, flipImage(trimmed_full[amp.getBBox()].image, flipLR=True, flipTB=False)
444 )
445 self.assertAmplifiersEqual(test_t1_trimmed.getDetector()[0], amp_t1)
446 test_t1_untrimmed = self.butler.get(untrimmed_ref, parameters={"amp": amp_t1})
447 self.assertImagesEqual(
448 test_t1_untrimmed.image,
449 flipImage(untrimmed_full[amp.getRawBBox()].image, flipLR=True, flipTB=False),
450 )
451 self.assertAmplifiersEqual(test_t1_trimmed.getDetector()[0], amp_t1)
452 # Flip Y only.
453 amp_t2 = amp.rebuild().transform(outFlipY=True).finish()
454 test_t2_trimmed = self.butler.get(trimmed_ref, parameters={"amp": amp_t2})
455 self.assertImagesEqual(
456 test_t2_trimmed.image, flipImage(trimmed_full[amp.getBBox()].image, flipLR=False, flipTB=True)
457 )
458 self.assertAmplifiersEqual(test_t2_trimmed.getDetector()[0], amp_t2)
459 test_t2_untrimmed = self.butler.get(untrimmed_ref, parameters={"amp": amp_t2})
460 self.assertImagesEqual(
461 test_t2_untrimmed.image,
462 flipImage(untrimmed_full[amp.getRawBBox()].image, flipLR=False, flipTB=True),
463 )
464 self.assertAmplifiersEqual(test_t2_trimmed.getDetector()[0], amp_t2)
465 # Add an XY offset only.
466 amp_t3 = amp.rebuild().transform(outOffset=Extent2I(5, 4)).finish()
467 test_t3_trimmed = self.butler.get(trimmed_ref, parameters={"amp": amp_t3})
468 self.assertImagesEqual(test_t3_trimmed.image, trimmed_full[amp.getBBox()].image)
469 self.assertAmplifiersEqual(test_t3_trimmed.getDetector()[0], amp_t3)
470 test_t3_untrimmed = self.butler.get(untrimmed_ref, parameters={"amp": amp_t3})
471 self.assertImagesEqual(test_t3_untrimmed.image, untrimmed_full[amp.getRawBBox()].image)
472 self.assertAmplifiersEqual(test_t3_trimmed.getDetector()[0], amp_t3)
474 def testMaskedImageFormatter(self):
475 """Test that a MaskedImage can be persisted and read from a Butler."""
476 # Read in an Exposure as MaskedImage using MaskedImageFitsReader.
477 reader = MaskedImageFitsReader(os.path.join(TESTDIR, "data", "calexp.fits"))
478 mi = reader.read()
480 # Put the MaskedImage into the Butler and get a reference to it.
481 dataId = {"visit": 42, "instrument": "DummyCam", "physical_filter": "d-r"}
482 ref = self.butler.put(mi, "noise", dataId)
484 # Check that the MaskedImage can be retrieved from the butler.
485 maskedImage = self.butler.get(ref)
486 self.assertImagesEqual(maskedImage.image, mi.image)
487 self.assertImagesEqual(maskedImage.mask, mi.mask)
488 self.assertImagesEqual(maskedImage.variance, mi.variance)
490 # Get a DeferredDatasetHandle to load parts of the MaskedImage.
491 handle = self.butler.getDeferred(ref)
493 for parameters in (
494 {},
495 {"bbox": reader.readBBox()},
496 {"bbox": Box2I(minimum=Point2I(3, 3), maximum=Point2I(21, 16))},
497 {"bbox": Box2I(minimum=Point2I(3, 3), maximum=Point2I(21, 16)), "origin": LOCAL},
498 ):
499 bbox = parameters.get("bbox", reader.readBBox())
500 with self.subTest(parameters=parameters):
501 # Check that the reader supports reading sub-regions.
502 subMaskedImage = reader.read(**parameters)
503 self.assertImagesEqual(subMaskedImage.image, mi.image[bbox])
504 self.assertImagesEqual(subMaskedImage.mask, mi.mask[bbox])
505 self.assertImagesEqual(subMaskedImage.variance, mi.variance[bbox])
507 # Get a maskedImage within a bounding box from the butler.
508 subMaskedImage = handle.get(parameters=parameters)
509 self.assertImagesEqual(subMaskedImage.image, mi.image[bbox])
510 self.assertImagesEqual(subMaskedImage.mask, mi.mask[bbox])
511 self.assertImagesEqual(subMaskedImage.variance, mi.variance[bbox])
513 # Get one component at a time from the butler.
514 subImage = handle.get(parameters=parameters, component="image")
515 subMask = handle.get(parameters=parameters, component="mask")
516 subVariance = handle.get(parameters=parameters, component="variance")
517 self.assertImagesEqual(subImage, mi.image[bbox])
518 self.assertImagesEqual(subMask, mi.mask[bbox])
519 self.assertImagesEqual(subVariance, mi.variance[bbox])
522if __name__ == "__main__": 522 ↛ 523line 522 didn't jump to line 523 because the condition on line 522 was never true
523 unittest.main()