Coverage for python/lsst/cell_coadds/_fits.py: 16%

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/>. 

21 

22"""Module to handle FITS serialization and de-serialization. 

23 

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. 

28 

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: 

34 

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. 

38 

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. 

45 

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. 

55 

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. 

61 

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. 

66 

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. 

72 

73Note that major version 0 is considered unstable and experimental and none of 

74the guarantee above applies. 

75""" 

76 

77from __future__ import annotations 

78 

79__all__ = ( 

80 "CellCoaddFitsFormatter", 

81 "CellCoaddFitsReader", 

82 "IncompatibleVersionError", 

83 "writeMultipleCellCoaddAsFits", 

84) 

85 

86import logging 

87import os 

88from collections.abc import Mapping 

89from dataclasses import dataclass 

90from typing import Any 

91 

92import numpy as np 

93from astropy.io import fits 

94from astropy.table import Table 

95from frozendict import frozendict 

96from numpy.typing import DTypeLike 

97from packaging import version 

98 

99import lsst.afw.geom as afwGeom 

100import lsst.afw.image as afwImage 

101from lsst.afw.image import ImageD, ImageF 

102from lsst.daf.base import PropertySet 

103from lsst.geom import Box2I, Extent2I, Point2D, Point2I 

104from lsst.obs.base.formatters.fitsGeneric import FitsGenericFormatter 

105from lsst.skymap import Index2D 

106 

107from ._coadd_ap_corr_map import EMPTY_AP_CORR_MAP 

108from ._common_components import CoaddUnits, CommonComponents 

109from ._grid_container import GridContainer 

110from ._identifiers import CellIdentifiers, ObservationIdentifiers, PatchIdentifiers 

111from ._image_planes import OwnedImagePlanes 

112from ._multiple_cell_coadd import MultipleCellCoadd 

113from ._single_cell_coadd import CoaddInputs, SingleCellCoadd 

114from ._uniform_grid import UniformGrid 

115from .typing_helpers import SingleCellCoaddApCorrMap 

116 

117FILE_FORMAT_VERSION = "0.7" 

118"""Version number for the file format as persisted, presented as a string of 

119the form M.m, where M is the major version, m is the minor version. 

120""" 

121 

122MAX_POLYGON_VERTICES = 6 

123"""Maximum number of vertices the overlap polygon region between a per-detector 

124warp and the patch bounding box can have.""" 

125# 3 vertices from the detector, 3 vertices from the patch bounding box. 

126 

127logger = logging.getLogger(__name__) 

128 

129 

130class IncompatibleVersionError(RuntimeError): 

131 """Exception raised when the CellCoaddFitsReader version is not compatible 

132 with the FITS file attempted to read. 

133 """ 

134 

135 

136@dataclass 

137class VisitRecord: 

138 """A dataclass to hold relevant info about a visit. 

139 

140 This is intended for use with this module. 

141 """ 

142 

143 visit: int 

144 day_obs: int 

145 physical_filter: str 

146 

147 

148class CellCoaddFitsFormatter(FitsGenericFormatter): 

149 """Interface for writing and reading cell coadds to/from FITS files. 

150 

151 This assumes the existence of readFits and writeFits methods (for now). 

152 """ 

153 

154 

155class CellCoaddFitsReader: 

156 """A reader class to read from a FITS file and produce cell-based coadds. 

157 

158 This reader class has read methods that can either return a single 

159 component without reading the entire file (e.g., readBBox, readWcs) 

160 and read methods that return a full coadd (e.g., 

161 readAsMultipleCellCoadd, readAsExplodedCellCoadd, readAsStitchedCoadd). 

162 

163 Parameters 

164 ---------- 

165 filename : `str` 

166 The name of the FITS file to read. 

167 """ 

168 

169 # Minimum and maximum compatible file format versions are listed as 

170 # iterables so as to allow for discontiguous intervals. 

171 MINIMUM_FILE_FORMAT_VERSIONS = ("0.1",) 

172 MAXIMUM_FILE_FORMAT_VERSIONS = ("1.0",) 

173 

174 def __init__(self, filename: str) -> None: 

175 if not os.path.exists(filename): 

176 raise FileNotFoundError(f"File {filename} not found") 

177 

178 self.filename = filename 

179 

180 @classmethod 

181 def isCompatibleWith(cls, written_version: str, /) -> bool: 

182 """Check if the serialization version is compatible with the reader. 

183 

184 This is a convenience method to ask if the current version of this 

185 class can read a file, based on the VERSION in its header. 

186 

187 Parameters 

188 ---------- 

189 written_version: `str` 

190 The VERSION of the file to be read. 

191 

192 Returns 

193 ------- 

194 compatible : `bool` 

195 Whether the reader can read a file whose VERSION is 

196 ``written_version``. 

197 

198 Notes 

199 ----- 

200 This accepts the other version as a positional argument only. 

201 """ 

202 written_version_object = version.parse(written_version) 

203 for min_version, max_version in zip( 

204 cls.MINIMUM_FILE_FORMAT_VERSIONS, 

205 cls.MAXIMUM_FILE_FORMAT_VERSIONS, 

206 strict=True, 

207 ): 

208 if version.parse(min_version) <= written_version_object < version.parse(max_version): 

209 return True 

210 

211 return False 

212 

213 def readAsMultipleCellCoadd(self) -> MultipleCellCoadd: 

214 """Read the FITS file as a MultipleCellCoadd object. 

215 

216 Raises 

217 ------ 

218 IncompatibleError 

219 Raised if the version of this module that wrote the file is 

220 incompatible with this module that is reading it in. 

221 ValueError 

222 Raised if the outer cell size is not padded equally in either 

223 direction or by an even number of pixels. 

224 """ 

225 with fits.open(self.filename) as hdu_list: 

226 header = hdu_list[1].header 

227 written_version = header.get("VERSION", "0.1") 

228 if not self.isCompatibleWith(written_version): 

229 raise IncompatibleVersionError( 

230 f"{self.filename} was written with version {written_version}" 

231 f"but attempting to read it with a reader designed for {FILE_FORMAT_VERSION}" 

232 ) 

233 if written_version != FILE_FORMAT_VERSION: 

234 logger.info( 

235 "Reading %s having version %s with reader designed for %s", 

236 self.filename, 

237 written_version, 

238 FILE_FORMAT_VERSION, 

239 ) 

240 

241 written_version = version.parse(written_version) 

242 

243 # TODO: Remove this when FILE_FORMAT_VERSION is bumped to 1.0 

244 if written_version < version.parse("0.3"): 

245 header.rename_keyword("BAND", "FILTER") 

246 

247 data = hdu_list[1].data 

248 

249 # Read in WCS 

250 ps = PropertySet() 

251 ps.update(hdu_list[0].header) 

252 wcs = afwGeom.makeSkyWcs(ps) 

253 

254 # Build the quantities needed to construct a MultipleCellCoadd. 

255 if written_version >= version.parse("0.4"): 

256 coadd_units = header["TUNIT2"] 

257 else: 

258 logger.info("Attemping to set pixel units from TUNIT1.") 

259 coadd_units = header.get("TUNIT1", CoaddUnits.legacy.name) 

260 

261 outer_cell_size = Extent2I(header["OCELL1"], header["OCELL2"]) 

262 psf_image_size = Extent2I(header["PSFSIZE1"], header["PSFSIZE2"]) 

263 

264 grid_cell_size = Extent2I(header["GRCELL1"], header["GRCELL2"]) # Inner size of a single cell. 

265 grid_shape = Extent2I(header["GRSHAPE1"], header["GRSHAPE2"]) 

266 grid_min = Point2I(header["GRMIN1"], header["GRMIN2"]) 

267 grid_padding_extent = outer_cell_size - grid_cell_size 

268 # In hindsight, it would have been easier to store the padding as a 

269 # separate keyword in the header instead of OCELL1 and OCELL2. 

270 # This can be done when the file format is bumped to 1.0. 

271 if grid_padding_extent.x != grid_padding_extent.y: 

272 raise ValueError( 

273 "Outer cell size is not padded equally in either directions. " 

274 f"Got {outer_cell_size} and {grid_cell_size}." 

275 ) 

276 if grid_padding_extent.x % 2 != 0: 

277 raise ValueError( 

278 "Outer cell size is not padded by an even number of pixels. " 

279 f"Got {outer_cell_size} and {grid_cell_size}." 

280 ) 

281 grid_padding = grid_padding_extent.x // 2 

282 grid = UniformGrid(cell_size=grid_cell_size, shape=grid_shape, padding=grid_padding, min=grid_min) 

283 

284 # This is the inner bounding box for the multiple cell coadd 

285 inner_bbox = Box2I( 

286 Point2I(header["INBBOX11"], header["INBBOX12"]), 

287 Point2I(header["INBBOX21"], header["INBBOX22"]), 

288 ) 

289 

290 # Attempt to get inputs for each cell. 

291 inputs = GridContainer[dict[ObservationIdentifiers, CoaddInputs]](shape=grid.shape) 

292 coadd_input = CoaddInputs( # default version for written_version <= 0.5. 

293 overlaps_center=True, 

294 overlap_fraction=1.0, 

295 weight=0.0, 

296 psf_shape=afwGeom.Quadrupole(), 

297 psf_shape_flag=True, 

298 ) 

299 visit_polygons = {} 

300 if written_version >= version.parse("0.3"): 

301 visit_dict = { 

302 row["visit"]: VisitRecord( 

303 visit=row["visit"], 

304 physical_filter=row["physical_filter"], 

305 day_obs=row["day_obs"], 

306 ) 

307 for row in hdu_list[hdu_list.index_of("VISIT")].data 

308 } 

309 link_table = hdu_list[hdu_list.index_of("CELL")].data 

310 

311 if written_version >= version.parse("0.6"): 

312 mask_definitions = hdu_list[hdu_list.index_of("MASKDEF")].data 

313 mask_plane_dict = {} 

314 for mask_name, mask_value in mask_definitions: 

315 afwImage.Mask.addMaskPlane(mask_name) 

316 mask_plane_dict[mask_name] = mask_value 

317 else: 

318 logger.warning( 

319 "Mask definitions can be incorrect since this was written with " 

320 "FILE_FORMAT_VERSION %s < 0.6", 

321 header["VERSION"], 

322 ) 

323 mask_plane_dict = None 

324 

325 if written_version >= version.parse("0.6"): 

326 visit_summary_hdu = hdu_list[hdu_list.index_of("VISIT_SUMMARY")] 

327 for row in visit_summary_hdu.data: 

328 visit = int(row["visit"]) 

329 obs_id = ObservationIdentifiers( 

330 instrument=header["INSTRUME"], 

331 visit=visit, 

332 detector=row["detector"], 

333 day_obs=visit_dict[visit].day_obs, 

334 physical_filter=visit_dict[visit].physical_filter, 

335 ) 

336 if written_version >= version.parse("0.7"): 

337 num_vertices = row["num_vertices"] 

338 else: 

339 num_vertices = None # li[:None] returns the entire list. 

340 visit_polygons[obs_id] = afwGeom.Polygon( 

341 [Point2D(vertex) for vertex in row["polygon_vertices"][:num_vertices]] 

342 ) 

343 

344 for link_row in link_table: 

345 cell_id = Index2D(link_row["cell_x"], link_row["cell_y"]) 

346 visit = link_row["visit"] 

347 obs_id = ObservationIdentifiers( 

348 instrument=header["INSTRUME"], 

349 visit=visit, 

350 detector=link_row["detector"], 

351 day_obs=visit_dict[visit].day_obs, 

352 physical_filter=visit_dict[visit].physical_filter, 

353 ) 

354 if written_version >= version.parse("0.6"): 

355 coadd_input = CoaddInputs( 

356 overlaps_center=link_row["overlaps_center"], 

357 overlap_fraction=link_row["overlap_fraction"], 

358 weight=link_row["weight"], 

359 psf_shape=afwGeom.Quadrupole( 

360 ixx=link_row["psf_shape_ixx"], 

361 iyy=link_row["psf_shape_iyy"], 

362 ixy=link_row["psf_shape_ixy"], 

363 ), 

364 psf_shape_flag=link_row["psf_shape_flag"], 

365 ) 

366 

367 if cell_id in inputs: 

368 inputs[cell_id][obs_id] = coadd_input 

369 else: 

370 inputs[cell_id] = {obs_id: coadd_input} 

371 else: 

372 logger.info( 

373 "Cell inputs are available for VERSION=0.3 or later. The file provided has ", 

374 "VERSION = %s", 

375 written_version, 

376 ) 

377 

378 try: 

379 aperture_correction_hdu = hdu_list[hdu_list.index_of("APCORR")].data 

380 if len(aperture_correction_hdu) < len(data): 

381 logger.warning("Aperture correction map is not available for all cells.") 

382 aperture_correction_grid = self._readApCorr(aperture_correction_hdu, grid_shape) 

383 except KeyError: 

384 if written_version >= version.parse("0.4"): 

385 logger.warning("Unable to read aperture correction map from the file.") 

386 # Create an empty grid container regardless. 

387 aperture_correction_grid = GridContainer[SingleCellCoaddApCorrMap](shape=grid_shape) 

388 

389 common = CommonComponents( 

390 units=CoaddUnits(coadd_units), 

391 wcs=wcs, 

392 band=header["FILTER"], 

393 identifiers=PatchIdentifiers( 

394 skymap=header["SKYMAP"], 

395 tract=header["TRACT"], 

396 patch=Index2D(x=header["PATCH_X"], y=header["PATCH_Y"]), 

397 band=header["FILTER"], 

398 ), 

399 visit_polygons=visit_polygons, 

400 ) 

401 

402 coadd = MultipleCellCoadd( 

403 ( 

404 self._readSingleCellCoadd( 

405 data=data[row_id], 

406 header=header, 

407 common=common, 

408 inputs=inputs[ 

409 Index2D(int(data[row_id]["cell_id"][0]), int(data[row_id]["cell_id"][1])) 

410 ], 

411 outer_cell_size=outer_cell_size, 

412 psf_image_size=psf_image_size, 

413 inner_cell_size=grid_cell_size, 

414 mask_plane_dict=mask_plane_dict, 

415 aperture_correction_map=aperture_correction_grid.get( 

416 Index2D( 

417 data[row_id]["cell_id"][0].tolist(), 

418 data[row_id]["cell_id"][1].tolist(), 

419 ), 

420 EMPTY_AP_CORR_MAP, 

421 ), 

422 ) 

423 for row_id in range(len(data)) 

424 ), 

425 grid=grid, 

426 outer_cell_size=outer_cell_size, 

427 psf_image_size=psf_image_size, 

428 inner_bbox=inner_bbox, 

429 common=common, 

430 ) 

431 

432 return coadd 

433 

434 @staticmethod 

435 def _readSingleCellCoadd( 

436 data: Mapping[str, Any], 

437 common: CommonComponents, 

438 header: Mapping[str, Any], 

439 *, 

440 inputs: Mapping[ObservationIdentifiers, CoaddInputs], 

441 outer_cell_size: Extent2I, 

442 inner_cell_size: Extent2I, 

443 psf_image_size: Extent2I, 

444 mask_plane_dict: Mapping[str, int] | None, 

445 aperture_correction_map: SingleCellCoaddApCorrMap = EMPTY_AP_CORR_MAP, 

446 ) -> SingleCellCoadd: 

447 """Read a coadd from a FITS file. 

448 

449 Parameters 

450 ---------- 

451 data : `Mapping` 

452 The data from the FITS file. Usually, a single row from the binary 

453 table representation. 

454 common : `CommonComponents` 

455 The common components of the coadd. 

456 header : `Mapping` 

457 The header of the FITS file as a dictionary. 

458 inputs : `Mapping` [`ObservationIdentifiers`, `CoaddInputs`] 

459 A mapping of ObservationIdentifiers that contributed to this cell 

460 to their corresponding CoaddInputs. 

461 outer_cell_size : `Extent2I` 

462 The size of the outer cell. 

463 psf_image_size : `Extent2I` 

464 The size of the PSF image. 

465 inner_cell_size : `Extent2I` 

466 The size of the inner cell. 

467 mask_plane_dict: `Mapping[str, int]` | None 

468 A mapping defining the mask planes to their bit value. 

469 aperture_correction_map : `Mapping` [`str`, `float`], optional 

470 Mapping of algorithm name to aperture correction value. 

471 If None, no aperture correction is applied. 

472 

473 Returns 

474 ------- 

475 coadd : `SingleCellCoadd` 

476 The coadd read from the file. 

477 """ 

478 buffer = (outer_cell_size - inner_cell_size) // 2 

479 

480 n_noise_realizations = header.get("NNOISE", 0) 

481 

482 psf = ImageD( 

483 array=data["psf"].astype(np.float64), 

484 xy0=(-(psf_image_size // 2)).asPoint(), # integer division and negation do not commute. 

485 ) # use the variable 

486 xy0 = Point2I( 

487 inner_cell_size.x * data["cell_id"][0] - buffer.x + header["GRMIN1"], 

488 inner_cell_size.y * data["cell_id"][1] - buffer.y + header["GRMIN2"], 

489 ) 

490 

491 # Parse the mask array carefully, accounting for a change in 

492 # the definitions. 

493 mask = afwImage.Mask( 

494 bbox=Box2I(corner=xy0, dimensions=Extent2I(outer_cell_size.x, outer_cell_size.y)), 

495 initialValue=0, 

496 ) 

497 mask_array = data["mask"].astype(np.int32) 

498 if mask_plane_dict is not None: 

499 for mask_name, bit_value in mask_plane_dict.items(): 

500 mask.array[(mask_array & 2**bit_value) > 0] |= afwImage.Mask.getPlaneBitMask(mask_name) 

501 else: 

502 mask.array[:, :] = mask_array 

503 

504 try: 

505 maskfrac = data["maskfrac"] 

506 mask_fractions = ImageF(maskfrac.astype(np.float32), xy0=xy0) 

507 except KeyError: 

508 mask_fractions = None 

509 

510 image_planes = OwnedImagePlanes( 

511 image=ImageF( 

512 data["image"].astype(np.float32), 

513 xy0=xy0, 

514 ), 

515 mask=mask, 

516 variance=ImageF(data["variance"].astype(np.float32), xy0=xy0), 

517 noise_realizations=[ 

518 ImageF(data[f"noise_{n:02}"].astype(np.float32), xy0=xy0) for n in range(n_noise_realizations) 

519 ], 

520 mask_fractions=mask_fractions, 

521 ) 

522 

523 identifiers = CellIdentifiers( 

524 cell=Index2D(data["cell_id"][0], data["cell_id"][1]), 

525 skymap=common.identifiers.skymap, 

526 tract=common.identifiers.tract, 

527 patch=common.identifiers.patch, 

528 band=common.identifiers.band, 

529 ) 

530 

531 return SingleCellCoadd( 

532 outer=image_planes, 

533 psf=psf, 

534 inner_bbox=Box2I( 

535 corner=Point2I( 

536 inner_cell_size.x * data["cell_id"][0] + header["GRMIN1"], 

537 inner_cell_size.y * data["cell_id"][1] + header["GRMIN2"], 

538 ), 

539 dimensions=inner_cell_size, 

540 ), 

541 common=common, 

542 identifiers=identifiers, 

543 inputs=inputs, 

544 aperture_correction_map=aperture_correction_map, 

545 ) 

546 

547 def readWcs(self) -> afwGeom.SkyWcs: 

548 """Read the WCS information from the FITS file. 

549 

550 Returns 

551 ------- 

552 wcs : `~lsst.afw.geom.SkyWcs` 

553 The WCS information read from the FITS file. 

554 """ 

555 # Read in WCS 

556 ps = PropertySet() 

557 with fits.open(self.filename) as hdu_list: 

558 ps.update(hdu_list[0].header) 

559 wcs = afwGeom.makeSkyWcs(ps) 

560 return wcs 

561 

562 def _readApCorr( 

563 self, aperture_correction_hdu: fits.fitsrec.FITS_rec, grid_shape: Index2D 

564 ) -> GridContainer[SingleCellCoaddApCorrMap]: 

565 """Read the aperture correction map from the FITS file.""" 

566 column_names = aperture_correction_hdu.names 

567 column_names.remove("x") 

568 column_names.remove("y") 

569 

570 gc = GridContainer[SingleCellCoaddApCorrMap](shape=grid_shape) 

571 for row_idx in range(len(aperture_correction_hdu)): 

572 cell_id = Index2D( 

573 x=aperture_correction_hdu["x"][row_idx], y=aperture_correction_hdu["y"][row_idx] 

574 ) 

575 ap_corr_map_dict = {name: float(aperture_correction_hdu[name][row_idx]) for name in column_names} 

576 gc[cell_id] = frozendict(ap_corr_map_dict) 

577 

578 return gc 

579 

580 

581def to_numpy_record(ap_corr_map: Mapping[str, float], xy: Index2D, dtypes: DTypeLike) -> np.record: 

582 """Convert the aperture correction map to a numpy record with 

583 appropriate data types. 

584 

585 Parameters 

586 ---------- 

587 ap_corr_map: `Mapping` [`str`, `float`] 

588 Aperture correction values for each algorithm 

589 xy: `~lsst.skymap.Index2D` 

590 An object identifying the cell 

591 dtypes: `~numpy.typing.DTypeLike` 

592 An iterable of tuples describing the name and data type. It must 

593 contain ("x", int) and ("y", int) to accommodate ``xy``. 

594 

595 Returns 

596 ------- 

597 record : `np.record` 

598 Per-cell aperture correction map as a record. 

599 """ 

600 record = np.recarray((1,), dtype=dtypes)[0] 

601 for field_name in ap_corr_map: 

602 record[field_name] = ap_corr_map.get(field_name, np.nan) 

603 

604 record["x"] = xy.x 

605 record["y"] = xy.y 

606 

607 return record 

608 

609 

610def writeMultipleCellCoaddAsFits( 

611 multiple_cell_coadd: MultipleCellCoadd, 

612 filename: str, 

613 overwrite: bool = False, 

614 metadata: PropertySet | None = None, 

615) -> fits.HDUList: 

616 """Write a MultipleCellCoadd object to a FITS file. 

617 

618 Parameters 

619 ---------- 

620 multiple_cell_coadd : `MultipleCellCoadd` 

621 The multiple cell coadd to write to a FITS file. 

622 filename : `str` 

623 The name of the file to write to. 

624 overwrite : `bool`, optional 

625 Whether to overwrite the file if it already exists? 

626 metadata : `~lsst.daf.base.PropertySet`, optional 

627 Additional metadata to write to the FITS file. 

628 

629 Returns 

630 ------- 

631 hdu_list : `~astropy.io.fits.HDUList` 

632 The FITS file as an HDUList. 

633 

634 Notes 

635 ----- 

636 Changes to this function that modify the way the file is written to disk 

637 must be accompanied with a change to FILE_FORMAT_VERSION. 

638 """ 

639 # Create metadata tables: 

640 # 1. Visit table containing information about the visits. 

641 # 2. Cell table containing info about the visit+detector for each cell. 

642 visit_records: list[Any] = [] 

643 cell_records: list[Any] = [] 

644 instrument_set = set() 

645 for cell_id, single_cell_coadd in multiple_cell_coadd.cells.items(): 

646 for observation_id, coadd_input in single_cell_coadd.inputs.items(): 

647 visit_records.append( 

648 (observation_id.visit, observation_id.physical_filter, observation_id.day_obs) 

649 ) 

650 cell_records.append( 

651 ( 

652 cell_id.x, 

653 cell_id.y, 

654 observation_id.visit, 

655 observation_id.detector, 

656 coadd_input.overlaps_center, 

657 coadd_input.overlap_fraction, 

658 coadd_input.weight, 

659 coadd_input.psf_shape.getIxx(), 

660 coadd_input.psf_shape.getIyy(), 

661 coadd_input.psf_shape.getIxy(), 

662 coadd_input.psf_shape_flag, 

663 ) 

664 ) 

665 instrument_set.add(observation_id.instrument) 

666 

667 assert len(instrument_set) == 1, "All cells must have the same instrument." 

668 instrument = instrument_set.pop() 

669 

670 # Create visit_summary equivalent table 

671 visit_column = fits.Column( 

672 name="visit", 

673 format="K", 

674 array=[obs_id.visit for obs_id in multiple_cell_coadd.common.visit_polygons], 

675 ) 

676 detector_column = fits.Column( 

677 name="detector", 

678 format="I", 

679 array=[obs_id.detector for obs_id in multiple_cell_coadd.common.visit_polygons], 

680 ) 

681 number_of_vertices = [] 

682 polygon_vertices_array = [] 

683 for obs_id, poly in multiple_cell_coadd.common.visit_polygons.items(): 

684 if num_vertices := len(poly.getVertices()) > MAX_POLYGON_VERTICES: 

685 logger.warning( 

686 "Visit %d, detector %d has a polygon with %d vertices. " 

687 "This geometry should be impossible for two intersecting " 

688 "convex quadrilaterals. Only the first %d will be stored in " 

689 "the FITS file.", 

690 obs_id.visit, 

691 obs_id.detector, 

692 num_vertices, 

693 MAX_POLYGON_VERTICES, 

694 ) 

695 number_of_vertices.append(num_vertices) 

696 vertices = poly.getVertices() + poly.getVertices() 

697 vertices = vertices[:MAX_POLYGON_VERTICES] 

698 polygon_vertices_array.append(np.array(vertices)) 

699 polygon_column = fits.Column( 

700 name="polygon_vertices", 

701 format=f"{2 * MAX_POLYGON_VERTICES}E", 

702 dim=f"(2,{MAX_POLYGON_VERTICES})", 

703 array=polygon_vertices_array, 

704 ) 

705 number_of_vertices_column = fits.Column( 

706 name="num_vertices", 

707 format="I", 

708 array=number_of_vertices, 

709 ) 

710 visit_summary_hdu = fits.BinTableHDU.from_columns( 

711 [visit_column, detector_column, number_of_vertices_column, polygon_column], 

712 name="VISIT_SUMMARY", 

713 ) 

714 visit_summary_hdu.header["POLYVERT"] = MAX_POLYGON_VERTICES 

715 

716 visit_recarray = np.rec.fromrecords( 

717 recList=sorted(set(visit_records), key=lambda x: x[0]), # Sort by visit. 

718 formats=None, # formats has specified to please mypy. See numpy#26376. 

719 names=( 

720 "visit", 

721 "physical_filter", 

722 "day_obs", 

723 ), 

724 ) 

725 cell_recarray = np.rec.fromrecords( 

726 recList=cell_records, 

727 formats=None, # formats has specified to please mypy. See numpy#26376. 

728 names=( 

729 "cell_x", 

730 "cell_y", 

731 "visit", 

732 "detector", 

733 "overlaps_center", 

734 "overlap_fraction", 

735 "weight", 

736 "psf_shape_ixx", 

737 "psf_shape_iyy", 

738 "psf_shape_ixy", 

739 "psf_shape_flag", 

740 ), 

741 ) 

742 

743 visit_hdu = fits.BinTableHDU.from_columns(visit_recarray, name="VISIT") 

744 cell_hdu = fits.BinTableHDU.from_columns(cell_recarray, name="CELL") 

745 

746 cell_id = fits.Column( 

747 name="cell_id", 

748 format="2I", 

749 array=[cell.identifiers.cell for cell in multiple_cell_coadd.cells.values()], 

750 ) 

751 

752 image_array = [cell.outer.image.array for cell in multiple_cell_coadd.cells.values()] 

753 unit_array = [cell.common.units.name for cell in multiple_cell_coadd.cells.values()] 

754 image = fits.Column( 

755 name="image", 

756 unit=unit_array[0], 

757 format=f"{image_array[0].size}E", 

758 dim=f"({image_array[0].shape[1]}, {image_array[0].shape[0]})", 

759 array=image_array, 

760 ) 

761 

762 mask_array = [cell.outer.mask.array for cell in multiple_cell_coadd.cells.values()] 

763 mask = fits.Column( 

764 name="mask", 

765 format=f"{mask_array[0].size}I", 

766 dim=f"({mask_array[0].shape[1]}, {mask_array[0].shape[0]})", 

767 array=mask_array, 

768 ) 

769 # Add in mask bit definitions as its own table. 

770 mask_definition_table = Table(names=("MASK_NAME", "BIT_VALUE"), dtype=(str, np.uint32)) 

771 for name, value in multiple_cell_coadd.cells.arbitrary.outer.mask.getMaskPlaneDict().items(): 

772 mask_definition_table.add_row((name, value)) 

773 mask_definition_hdu = fits.BinTableHDU.from_columns(mask_definition_table.as_array(), name="MASKDEF") 

774 

775 variance_array = [cell.outer.variance.array for cell in multiple_cell_coadd.cells.values()] 

776 variance = fits.Column( 

777 name="variance", 

778 format=f"{variance_array[0].size}E", 

779 dim=f"({variance_array[0].shape[1]}, {variance_array[0].shape[0]})", 

780 array=variance_array, 

781 ) 

782 

783 psf_array = [cell.psf_image.array for cell in multiple_cell_coadd.cells.values()] 

784 psf = fits.Column( 

785 name="psf", 

786 format=f"{psf_array[0].size}D", 

787 dim=f"({psf_array[0].shape[1]}, {psf_array[0].shape[0]})", 

788 array=[cell.psf_image.array for cell in multiple_cell_coadd.cells.values()], 

789 ) 

790 

791 columns = [cell_id, image, mask, variance, psf] 

792 maskfrac_array = [ 

793 cell.outer.mask_fractions.array 

794 for cell in multiple_cell_coadd.cells.values() 

795 if cell.outer.mask_fractions 

796 ] 

797 if maskfrac_array: 

798 maskfrac = fits.Column( 

799 name="maskfrac", 

800 format=f"{maskfrac_array[0].size}E", 

801 dim=f"({maskfrac_array[0].shape[1]}, {maskfrac_array[0].shape[0]})", 

802 array=maskfrac_array, 

803 ) 

804 columns.append(maskfrac) 

805 

806 n_noise_realizations = multiple_cell_coadd.n_noise_realizations 

807 for n in range(n_noise_realizations): 

808 noise_array = [cell.outer.noise_realizations[n].array for cell in multiple_cell_coadd.cells.values()] 

809 columns.append( 

810 fits.Column( 

811 name=f"noise_{n:02}", 

812 format=f"{noise_array[0].size}E", 

813 dim=f"({noise_array[0].shape[1]}, {noise_array[0].shape[0]})", 

814 array=noise_array, 

815 ) 

816 ) 

817 

818 if multiple_cell_coadd.cells.arbitrary.aperture_correction_map: 

819 apcorr_fields: set[str] = set() 

820 for cell in multiple_cell_coadd.cells.values(): 

821 if cell.aperture_correction_map: 

822 apcorr_fields.update(cell.aperture_correction_map) 

823 dtypes = [("x", int), ("y", int)] + [(key, float) for key in apcorr_fields] 

824 aperture_correction_recarray = np.rec.fromrecords( 

825 recList=[ 

826 to_numpy_record(cell.aperture_correction_map, cell.identifiers.cell, dtypes) 

827 for cell in multiple_cell_coadd.cells.values() 

828 if cell.aperture_correction_map 

829 ], 

830 dtype=dtypes, 

831 ) 

832 aperture_correction_hdu = fits.BinTableHDU.from_columns(aperture_correction_recarray, name="APCORR")