Coverage for python/lsst/cell_coadds/_fits.py: 19%
220 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-10 02:10 -0700
1# This file is part of cell_coadds.
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/>.
22"""Module to handle FITS serialization and de-serialization.
24The routines to write and read the files are in the same module, as a change to
25one is typically accompanied by a corresponding change to another. Code changes
26relating to writing the file must bump to the version number denoted by the
27module constant FILE_FORMAT_VERSION.
29Although the typical use case is for newer versions of the code to read files
30written by an older version, for the purposes of deciding the newer version
31string, it is helpful to think about an older version of the reader attempting
32to read a newer version of the file on disk. The policy for bumping the version
33is as follows:
351. When the on-disk file format written by this module changes such that the
36previous version of the reader can still read files written by the newer
37version, then there should be a minor bump.
392. When the on-disk format written by this module changes in a way that will
40prevent the previous version of the reader from reading a file produced by the
41current version of the module, then there should be a major bump. This usually
42means that the new version of the reader cannot read older file either,
43save the temporary support with deprecation warnings, possibly until a new
44release of the Science Pipelines is made.
46Examples
47--------
481. A file with VERSION=1.3 should still be readable by the reader in
49this module when the module-level constant FILE_FORMAT_VERSION=1.4. A file
50written with VERSION=1.4 will typically be readable by a reader when the
51module-level FILE_FORMAT_VERSION=1.3, although such a use case is not expected.
52A concrete example of change
53that requires only a minor bump is adding another BinTable that keeps track of
54the input visits.
562. An example of major change would be migrating from using
57BinTableHDU to ImageHDU to save data. Even if the reader supports reading
58either of this formats based on the value of VERSION from the header, it should
59be a major change because the previous version of the reader cannot read data
60from ImageHDUs.
62Unit tests only check that a file written can be read by the concurrent version
63of the module, but not by any of the previous ones. Hence, bumping
64FILE_FORMAT_VERSION to the appropriate value is ultimately at the discretion of
65the developers.
67A major bump must also be recorded in the `isCompatibleWith` method.
68It is plausible that different (non-consequent) major format versions can be
69read by the same reader (due to reverting back to an earlier format, or to
70something very similar). `isCompatibleWith` method offers the convenience of
71checking if a particular format version can be read by the current reader.
73Note that major version 0 is considered unstable and experimental and none of
74the guarantee above applies.
75"""
77from __future__ import annotations
79__all__ = (
80 "CellCoaddFitsFormatter",
81 "CellCoaddFitsReader",
82 "IncompatibleVersionError",
83 "writeMultipleCellCoaddAsFits",
84)
86import logging
87import os
88from collections.abc import Iterable, Mapping
89from dataclasses import dataclass
90from typing import Any
92import numpy as np
93from astropy.io import fits
94from frozendict import frozendict
95from numpy.typing import DTypeLike
96from packaging import version
98import lsst.afw.geom as afwGeom
99import lsst.afw.image as afwImage
100from lsst.afw.image import ImageD, ImageF
101from lsst.daf.base import PropertySet
102from lsst.geom import Box2I, Extent2I, Point2I
103from lsst.obs.base.formatters.fitsGeneric import FitsGenericFormatter
104from lsst.skymap import Index2D
106from ._coadd_ap_corr_map import EMPTY_AP_CORR_MAP
107from ._common_components import CoaddUnits, CommonComponents
108from ._grid_container import GridContainer
109from ._identifiers import CellIdentifiers, ObservationIdentifiers, PatchIdentifiers
110from ._image_planes import OwnedImagePlanes
111from ._multiple_cell_coadd import MultipleCellCoadd, SingleCellCoadd
112from ._uniform_grid import UniformGrid
113from .typing_helpers import SingleCellCoaddApCorrMap
115FILE_FORMAT_VERSION = "0.5"
116"""Version number for the file format as persisted, presented as a string of
117the form M.m, where M is the major version, m is the minor version.
118"""
120logger = logging.getLogger(__name__)
123class IncompatibleVersionError(RuntimeError):
124 """Exception raised when the CellCoaddFitsReader version is not compatible
125 with the FITS file attempted to read.
126 """
129@dataclass
130class VisitRecord:
131 """A dataclass to hold relevant info about a visit.
133 This is intended for use with this module.
134 """
136 visit: int
137 day_obs: int
138 physical_filter: str
141class CellCoaddFitsFormatter(FitsGenericFormatter):
142 """Interface for writing and reading cell coadds to/from FITS files.
144 This assumes the existence of readFits and writeFits methods (for now).
145 """
148class CellCoaddFitsReader:
149 """A reader class to read from a FITS file and produce cell-based coadds.
151 This reader class has read methods that can either return a single
152 component without reading the entire file (e.g., readBBox, readWcs)
153 and read methods that return a full coadd (e.g.,
154 readAsMultipleCellCoadd, readAsExplodedCellCoadd, readAsStitchedCoadd).
156 Parameters
157 ----------
158 filename : `str`
159 The name of the FITS file to read.
160 """
162 # Minimum and maximum compatible file format versions are listed as
163 # iterables so as to allow for discontiguous intervals.
164 MINIMUM_FILE_FORMAT_VERSIONS = ("0.1",)
165 MAXIMUM_FILE_FORMAT_VERSIONS = ("1.0",)
167 def __init__(self, filename: str) -> None:
168 if not os.path.exists(filename):
169 raise FileNotFoundError(f"File {filename} not found")
171 self.filename = filename
173 @classmethod
174 def isCompatibleWith(cls, written_version: str, /) -> bool:
175 """Check if the serialization version is compatible with the reader.
177 This is a convenience method to ask if the current version of this
178 class can read a file, based on the VERSION in its header.
180 Parameters
181 ----------
182 written_version: `str`
183 The VERSION of the file to be read.
185 Returns
186 -------
187 compatible : `bool`
188 Whether the reader can read a file whose VERSION is
189 ``written_version``.
191 Notes
192 -----
193 This accepts the other version as a positional argument only.
194 """
195 written_version_object = version.parse(written_version)
196 for min_version, max_version in zip(
197 cls.MINIMUM_FILE_FORMAT_VERSIONS,
198 cls.MAXIMUM_FILE_FORMAT_VERSIONS,
199 strict=True,
200 ):
201 if version.parse(min_version) <= written_version_object < version.parse(max_version):
202 return True
204 return False
206 def readAsMultipleCellCoadd(self) -> MultipleCellCoadd:
207 """Read the FITS file as a MultipleCellCoadd object.
209 Raises
210 ------
211 IncompatibleError
212 Raised if the version of this module that wrote the file is
213 incompatible with this module that is reading it in.
214 ValueError
215 Raised if the outer cell size is not padded equally in either
216 direction or by an even number of pixels.
217 """
218 with fits.open(self.filename) as hdu_list:
219 header = hdu_list[1].header
220 written_version = header.get("VERSION", "0.1")
221 if not self.isCompatibleWith(written_version):
222 raise IncompatibleVersionError(
223 f"{self.filename} was written with version {written_version}"
224 f"but attempting to read it with a reader designed for {FILE_FORMAT_VERSION}"
225 )
226 if written_version != FILE_FORMAT_VERSION:
227 logger.info(
228 "Reading %s having version %s with reader designed for %s",
229 self.filename,
230 written_version,
231 FILE_FORMAT_VERSION,
232 )
234 written_version = version.parse(written_version)
236 # TODO: Remove this when FILE_FORMAT_VERSION is bumped to 1.0
237 if written_version < version.parse("0.3"):
238 header.rename_keyword("BAND", "FILTER")
240 data = hdu_list[1].data
242 # Read in WCS
243 ps = PropertySet()
244 ps.update(hdu_list[0].header)
245 wcs = afwGeom.makeSkyWcs(ps)
247 # Build the quantities needed to construct a MultipleCellCoadd.
248 if written_version >= version.parse("0.4"):
249 coadd_units = header["TUNIT2"]
250 else:
251 logger.info("Attemping to set pixel units from TUNIT1.")
252 coadd_units = header.get("TUNIT1", CoaddUnits.legacy.name)
254 common = CommonComponents(
255 units=CoaddUnits(coadd_units),
256 wcs=wcs,
257 band=header["FILTER"],
258 identifiers=PatchIdentifiers(
259 skymap=header["SKYMAP"],
260 tract=header["TRACT"],
261 patch=Index2D(x=header["PATCH_X"], y=header["PATCH_Y"]),
262 band=header["FILTER"],
263 ),
264 )
266 outer_cell_size = Extent2I(header["OCELL1"], header["OCELL2"])
267 psf_image_size = Extent2I(header["PSFSIZE1"], header["PSFSIZE2"])
269 grid_cell_size = Extent2I(header["GRCELL1"], header["GRCELL2"]) # Inner size of a single cell.
270 grid_shape = Extent2I(header["GRSHAPE1"], header["GRSHAPE2"])
271 grid_min = Point2I(header["GRMIN1"], header["GRMIN2"])
272 grid_padding_extent = outer_cell_size - grid_cell_size
273 # In hindsight, it would have been easier to store the padding as a
274 # separate keyword in the header instead of OCELL1 and OCELL2.
275 # This can be done when the file format is bumped to 1.0.
276 if grid_padding_extent.x != grid_padding_extent.y:
277 raise ValueError(
278 "Outer cell size is not padded equally in either directions. "
279 f"Got {outer_cell_size} and {grid_cell_size}."
280 )
281 if grid_padding_extent.x % 2 != 0:
282 raise ValueError(
283 "Outer cell size is not padded by an even number of pixels. "
284 f"Got {outer_cell_size} and {grid_cell_size}."
285 )
286 grid_padding = grid_padding_extent.x // 2
287 grid = UniformGrid(cell_size=grid_cell_size, shape=grid_shape, padding=grid_padding, min=grid_min)
289 # This is the inner bounding box for the multiple cell coadd
290 inner_bbox = Box2I(
291 Point2I(header["INBBOX11"], header["INBBOX12"]),
292 Point2I(header["INBBOX21"], header["INBBOX22"]),
293 )
295 # Attempt to get inputs for each cell.
296 inputs = GridContainer[list[ObservationIdentifiers]](shape=grid.shape)
297 if written_version >= version.parse("0.3"):
298 visit_dict = {
299 row["visit"]: VisitRecord(
300 visit=row["visit"],
301 physical_filter=row["physical_filter"],
302 day_obs=row["day_obs"],
303 )
304 for row in hdu_list[hdu_list.index_of("VISIT")].data
305 }
306 link_table = hdu_list[hdu_list.index_of("CELL")].data
307 for link_row in link_table:
308 cell_id = Index2D(link_row["cell_x"], link_row["cell_y"])
309 visit = link_row["visit"]
310 obs_id = ObservationIdentifiers(
311 instrument=header["INSTRUME"],
312 visit=visit,
313 detector=link_row["detector"],
314 day_obs=visit_dict[visit].day_obs,
315 physical_filter=visit_dict[visit].physical_filter,
316 )
317 if cell_id in inputs:
318 inputs[cell_id] += [obs_id]
319 else:
320 inputs[cell_id] = [obs_id]
321 else:
322 logger.info(
323 "Cell inputs are available for VERSION=0.3 or later. The file provided has ",
324 "VERSION = %s",
325 written_version,
326 )
328 try:
329 aperture_correction_hdu = hdu_list[hdu_list.index_of("APCORR")].data
330 assert len(aperture_correction_hdu) == 0 or len(aperture_correction_hdu) == len(
331 data
332 ), "Aperture correction map is not available for all cells."
333 aperture_correction_grid = self._readApCorr(aperture_correction_hdu, grid_shape)
334 except KeyError:
335 if written_version >= version.parse("0.4"):
336 logger.warning("Unable to read aperture correction map from the file.")
337 # Create an empty grid container regardless.
338 aperture_correction_grid = GridContainer[SingleCellCoaddApCorrMap](shape=grid_shape)
340 coadd = MultipleCellCoadd(
341 (
342 self._readSingleCellCoadd(
343 data=data[row_id],
344 header=header,
345 common=common,
346 inputs=inputs[
347 Index2D(int(data[row_id]["cell_id"][0]), int(data[row_id]["cell_id"][1]))
348 ],
349 outer_cell_size=outer_cell_size,
350 psf_image_size=psf_image_size,
351 inner_cell_size=grid_cell_size,
352 aperture_correction_map=aperture_correction_grid.get(
353 Index2D(
354 data[row_id]["cell_id"][0].tolist(),
355 data[row_id]["cell_id"][1].tolist(),
356 ),
357 EMPTY_AP_CORR_MAP,
358 ),
359 )
360 for row_id in range(len(data))
361 ),
362 grid=grid,
363 outer_cell_size=outer_cell_size,
364 psf_image_size=psf_image_size,
365 inner_bbox=inner_bbox,
366 common=common,
367 )
369 return coadd
371 @staticmethod
372 def _readSingleCellCoadd(
373 data: Mapping[str, Any],
374 common: CommonComponents,
375 header: Mapping[str, Any],
376 *,
377 inputs: Iterable[ObservationIdentifiers],
378 outer_cell_size: Extent2I,
379 inner_cell_size: Extent2I,
380 psf_image_size: Extent2I,
381 aperture_correction_map: SingleCellCoaddApCorrMap = EMPTY_AP_CORR_MAP,
382 ) -> SingleCellCoadd:
383 """Read a coadd from a FITS file.
385 Parameters
386 ----------
387 data : `Mapping`
388 The data from the FITS file. Usually, a single row from the binary
389 table representation.
390 common : `CommonComponents`
391 The common components of the coadd.
392 header : `Mapping`
393 The header of the FITS file as a dictionary.
394 inputs : `Iterable` [`ObservationIdentifiers`]
395 Any iterable of ObservationIdentifiers instances that contributed
396 to this cell.
397 outer_cell_size : `Extent2I`
398 The size of the outer cell.
399 psf_image_size : `Extent2I`
400 The size of the PSF image.
401 inner_cell_size : `Extent2I`
402 The size of the inner cell.
403 aperture_correction_map : `Mapping` [`str`, `float`], optional
404 Mapping of algorithm name to aperture correction value.
405 If None, no aperture correction is applied.
407 Returns
408 -------
409 coadd : `SingleCellCoadd`
410 The coadd read from the file.
411 """
412 buffer = (outer_cell_size - inner_cell_size) // 2
414 n_noise_realizations = header.get("NNOISE", 0)
416 psf = ImageD(
417 array=data["psf"].astype(np.float64),
418 xy0=(-(psf_image_size // 2)).asPoint(), # integer division and negation do not commute.
419 ) # use the variable
420 xy0 = Point2I(
421 inner_cell_size.x * data["cell_id"][0] - buffer.x + header["GRMIN1"],
422 inner_cell_size.y * data["cell_id"][1] - buffer.y + header["GRMIN2"],
423 )
424 mask = afwImage.Mask(data["mask"].astype(np.int32), xy0=xy0)
425 try:
426 maskfrac = data["maskfrac"]
427 mask_fractions = ImageF(maskfrac.astype(np.float32), xy0=xy0)
428 except KeyError:
429 mask_fractions = None
431 image_planes = OwnedImagePlanes(
432 image=ImageF(
433 data["image"].astype(np.float32),
434 xy0=xy0,
435 ),
436 mask=mask,
437 variance=ImageF(data["variance"].astype(np.float32), xy0=xy0),
438 noise_realizations=[
439 ImageF(data[f"noise_{n:02}"].astype(np.float32), xy0=xy0) for n in range(n_noise_realizations)
440 ],
441 mask_fractions=mask_fractions,
442 )
444 identifiers = CellIdentifiers(
445 cell=Index2D(data["cell_id"][0], data["cell_id"][1]),
446 skymap=common.identifiers.skymap,
447 tract=common.identifiers.tract,
448 patch=common.identifiers.patch,
449 band=common.identifiers.band,
450 )
452 return SingleCellCoadd(
453 outer=image_planes,
454 psf=psf,
455 inner_bbox=Box2I(
456 corner=Point2I(
457 inner_cell_size.x * data["cell_id"][0] + header["GRMIN1"],
458 inner_cell_size.y * data["cell_id"][1] + header["GRMIN2"],
459 ),
460 dimensions=inner_cell_size,
461 ),
462 common=common,
463 identifiers=identifiers,
464 inputs=inputs,
465 aperture_correction_map=aperture_correction_map,
466 )
468 def readWcs(self) -> afwGeom.SkyWcs:
469 """Read the WCS information from the FITS file.
471 Returns
472 -------
473 wcs : `~lsst.afw.geom.SkyWcs`
474 The WCS information read from the FITS file.
475 """
476 # Read in WCS
477 ps = PropertySet()
478 with fits.open(self.filename) as hdu_list:
479 ps.update(hdu_list[0].header)
480 wcs = afwGeom.makeSkyWcs(ps)
481 return wcs
483 def _readApCorr(
484 self, aperture_correction_hdu: fits.fitsrec.FITS_rec, grid_shape: Index2D
485 ) -> GridContainer[SingleCellCoaddApCorrMap]:
486 """Read the aperture correction map from the FITS file."""
487 column_names = aperture_correction_hdu.names
488 column_names.remove("x")
489 column_names.remove("y")
491 gc = GridContainer[SingleCellCoaddApCorrMap](shape=grid_shape)
492 for row_idx in range(len(aperture_correction_hdu)):
493 cell_id = Index2D(
494 x=aperture_correction_hdu["x"][row_idx], y=aperture_correction_hdu["y"][row_idx]
495 )
496 ap_corr_map_dict = {name: float(aperture_correction_hdu[name][row_idx]) for name in column_names}
497 gc[cell_id] = frozendict(ap_corr_map_dict)
499 return gc
502def to_numpy_record(ap_corr_map: Mapping[str, float], xy: Index2D, dtypes: DTypeLike) -> np.record:
503 """Convert the aperture correction map to a numpy record with
504 appropriate data types.
506 Parameters
507 ----------
508 ap_corr_map: `Mapping` [`str`, `float`]
509 Aperture correction values for each algorithm
510 xy: `~lsst.skymap.Index2D`
511 An object identifying the cell
512 dtypes: `~numpy.typing.DTypeLike`
513 An iterable of tuples describing the name and data type. It must
514 contain ("x", int) and ("y", int) to accommodate ``xy``.
516 Returns
517 -------
518 record : `np.record`
519 Per-cell aperture correction map as a record.
520 """
521 record = np.recarray((1,), dtype=dtypes)[0]
522 for field_name in ap_corr_map:
523 record[field_name] = ap_corr_map.get(field_name, np.nan)
525 record["x"] = xy.x
526 record["y"] = xy.y
528 return record
531def writeMultipleCellCoaddAsFits(
532 multiple_cell_coadd: MultipleCellCoadd,
533 filename: str,
534 overwrite: bool = False,
535 metadata: PropertySet | None = None,
536) -> fits.HDUList:
537 """Write a MultipleCellCoadd object to a FITS file.
539 Parameters
540 ----------
541 multiple_cell_coadd : `MultipleCellCoadd`
542 The multiple cell coadd to write to a FITS file.
543 filename : `str`
544 The name of the file to write to.
545 overwrite : `bool`, optional
546 Whether to overwrite the file if it already exists?
547 metadata : `~lsst.daf.base.PropertySet`, optional
548 Additional metadata to write to the FITS file.
550 Returns
551 -------
552 hdu_list : `~astropy.io.fits.HDUList`
553 The FITS file as an HDUList.
555 Notes
556 -----
557 Changes to this function that modify the way the file is written to disk
558 must be accompanied with a change to FILE_FORMAT_VERSION.
559 """
560 # Create metadata tables:
561 # 1. Visit table containing information about the visits.
562 # 2. Cell table containing info about the visit+detector for each cell.
563 visit_records: list[Any] = []
564 cell_records: list[Any] = []
565 instrument_set = set()
566 for cell_id, single_cell_coadd in multiple_cell_coadd.cells.items():
567 for observation_id in single_cell_coadd.inputs:
568 visit_records.append(
569 (observation_id.visit, observation_id.physical_filter, observation_id.day_obs)
570 )
571 cell_records.append((cell_id.x, cell_id.y, observation_id.visit, observation_id.detector))
572 instrument_set.add(observation_id.instrument)
574 assert len(instrument_set) == 1, "All cells must have the same instrument."
575 instrument = instrument_set.pop()
577 visit_recarray = np.rec.fromrecords(
578 recList=sorted(set(visit_records), key=lambda x: x[0]), # Sort by visit.
579 formats=None, # formats has specified to please mypy. See numpy#26376.
580 names=(
581 "visit",
582 "physical_filter",
583 "day_obs",
584 ),
585 )
586 cell_recarray = np.rec.fromrecords(
587 recList=cell_records,
588 formats=None, # formats has specified to please mypy. See numpy#26376.
589 names=(
590 "cell_x",
591 "cell_y",
592 "visit",
593 "detector",
594 ),
595 )
597 visit_hdu = fits.BinTableHDU.from_columns(visit_recarray, name="VISIT")
598 cell_hdu = fits.BinTableHDU.from_columns(cell_recarray, name="CELL")
600 cell_id = fits.Column(
601 name="cell_id",
602 format="2I",
603 array=[cell.identifiers.cell for cell in multiple_cell_coadd.cells.values()],
604 )
606 image_array = [cell.outer.image.array for cell in multiple_cell_coadd.cells.values()]
607 unit_array = [cell.common.units.name for cell in multiple_cell_coadd.cells.values()]
608 image = fits.Column(
609 name="image",
610 unit=unit_array[0],
611 format=f"{image_array[0].size}E",
612 dim=f"({image_array[0].shape[1]}, {image_array[0].shape[0]})",
613 array=image_array,
614 )
616 mask_array = [cell.outer.mask.array for cell in multiple_cell_coadd.cells.values()]
617 mask = fits.Column(
618 name="mask",
619 format=f"{mask_array[0].size}I",
620 dim=f"({mask_array[0].shape[1]}, {mask_array[0].shape[0]})",
621 array=mask_array,
622 )
624 variance_array = [cell.outer.variance.array for cell in multiple_cell_coadd.cells.values()]
625 variance = fits.Column(
626 name="variance",
627 format=f"{variance_array[0].size}E",
628 dim=f"({variance_array[0].shape[1]}, {variance_array[0].shape[0]})",
629 array=variance_array,
630 )
632 psf_array = [cell.psf_image.array for cell in multiple_cell_coadd.cells.values()]
633 psf = fits.Column(
634 name="psf",
635 format=f"{psf_array[0].size}D",
636 dim=f"({psf_array[0].shape[1]}, {psf_array[0].shape[0]})",
637 array=[cell.psf_image.array for cell in multiple_cell_coadd.cells.values()],
638 )
640 columns = [cell_id, image, mask, variance, psf]
641 maskfrac_array = [
642 cell.outer.mask_fractions.array
643 for cell in multiple_cell_coadd.cells.values()
644 if cell.outer.mask_fractions
645 ]
646 if maskfrac_array:
647 maskfrac = fits.Column(
648 name="maskfrac",
649 format=f"{maskfrac_array[0].size}E",
650 dim=f"({maskfrac_array[0].shape[1]}, {maskfrac_array[0].shape[0]})",
651 array=maskfrac_array,
652 )
653 columns.append(maskfrac)
655 n_noise_realizations = multiple_cell_coadd.n_noise_realizations
656 for n in range(n_noise_realizations):
657 noise_array = [cell.outer.noise_realizations[n].array for cell in multiple_cell_coadd.cells.values()]
658 columns.append(
659 fits.Column(
660 name=f"noise_{n:02}",
661 format=f"{noise_array[0].size}E",
662 dim=f"({noise_array[0].shape[1]}, {noise_array[0].shape[0]})",
663 array=noise_array,
664 )
665 )
667 if multiple_cell_coadd.cells.arbitrary.aperture_correction_map:
668 apcorr_fields: set[str] = set()
669 for cell in multiple_cell_coadd.cells.values():
670 apcorr_fields.update(cell.aperture_correction_map)
671 dtypes = [("x", int), ("y", int)] + [(key, float) for key in apcorr_fields]
672 aperture_correction_recarray = np.rec.fromrecords(
673 recList=[
674 to_numpy_record(cell.aperture_correction_map, cell.identifiers.cell, dtypes)
675 for cell in multiple_cell_coadd.cells.values()
676 if cell.aperture_correction_map
677 ],
678 dtype=dtypes,
679 )
680 aperture_correction_hdu = fits.BinTableHDU.from_columns(aperture_correction_recarray, name="APCORR")
681 else:
682 aperture_correction_hdu = None
684 col_defs = fits.ColDefs(columns)
685 hdu = fits.BinTableHDU.from_columns(col_defs)
687 grid_cell_size = multiple_cell_coadd.grid.cell_size
688 grid_shape = multiple_cell_coadd.grid.shape
689 grid_min = multiple_cell_coadd.grid.bbox.getMin()
690 grid_cards = {
691 "GRCELL1": grid_cell_size.x,
692 "GRCELL2": grid_cell_size.y,
693 "GRSHAPE1": grid_shape.x,
694 "GRSHAPE2": grid_shape.y,
695 "GRMIN1": grid_min.x,
696 "GRMIN2": grid_min.y,
697 }
698 hdu.header.extend(grid_cards)
700 outer_cell_size_cards = {
701 "OCELL1": multiple_cell_coadd.outer_cell_size.x,
702 "OCELL2": multiple_cell_coadd.outer_cell_size.y,
703 }
704 hdu.header.extend(outer_cell_size_cards)
706 psf_image_size_cards = {
707 "PSFSIZE1": multiple_cell_coadd.psf_image_size.x,
708 "PSFSIZE2": multiple_cell_coadd.psf_image_size.y,
709 }
710 hdu.header.extend(psf_image_size_cards)
712 inner_bbox_cards = {
713 "INBBOX11": multiple_cell_coadd.inner_bbox.minX,
714 "INBBOX12": multiple_cell_coadd.inner_bbox.minY,
715 "INBBOX21": multiple_cell_coadd.inner_bbox.maxX,
716 "INBBOX22": multiple_cell_coadd.inner_bbox.maxY,
717 }
718 hdu.header.extend(inner_bbox_cards)
720 wcs = multiple_cell_coadd.common.wcs
721 wcs_cards = wcs.getFitsMetadata().toDict()
722 primary_hdu = fits.PrimaryHDU()
723 primary_hdu.header.extend(wcs_cards)
725 hdu.header["VERSION"] = FILE_FORMAT_VERSION
726 hdu.header["TUNIT2"] = multiple_cell_coadd.common.units.name
727 hdu.header["TUNIT4"] = multiple_cell_coadd.common.units.name + "**2"
728 # This assumed to be the same as multiple_cell_coadd.common.identifers.band
729 # See DM-38843.
730 hdu.header["INSTRUME"] = instrument
731 hdu.header["FILTER"] = multiple_cell_coadd.common.band
732 hdu.header["SKYMAP"] = multiple_cell_coadd.common.identifiers.skymap
733 hdu.header["TRACT"] = multiple_cell_coadd.common.identifiers.tract
734 hdu.header["PATCH_X"] = multiple_cell_coadd.common.identifiers.patch.x
735 hdu.header["PATCH_Y"] = multiple_cell_coadd.common.identifiers.patch.y
736 hdu.header["NNOISE"] = n_noise_realizations
738 if metadata is not None:
739 hdu.header.extend(metadata.toDict())
741 hdu_list = fits.HDUList([primary_hdu, hdu, cell_hdu, visit_hdu])
742 if aperture_correction_hdu:
743 hdu_list.append(aperture_correction_hdu)
745 hdu_list.writeto(filename, overwrite=overwrite)
747 return hdu_list