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

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.9" 

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 = 8 

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

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

125# -------------- 

126# | / \ | 

127# | / \ | 

128# |/ \| 

129# | | 

130# | | 

131# |\ | 

132# | \ /| 

133# | \ / | 

134# -------------- 

135 

136logger = logging.getLogger(__name__) 

137 

138 

139class IncompatibleVersionError(RuntimeError): 

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

141 with the FITS file attempted to read. 

142 """ 

143 

144 

145@dataclass 

146class VisitRecord: 

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

148 

149 This is intended for use with this module. 

150 """ 

151 

152 visit: int 

153 day_obs: int 

154 physical_filter: str 

155 

156 

157class CellCoaddFitsFormatter(FitsGenericFormatter): 

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

159 

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

161 """ 

162 

163 

164class CellCoaddFitsReader: 

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

166 

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

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

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

170 readAsMultipleCellCoadd, readAsExplodedCellCoadd, readAsStitchedCoadd). 

171 

172 Parameters 

173 ---------- 

174 filename : `str` 

175 The name of the FITS file to read. 

176 """ 

177 

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

179 # iterables so as to allow for discontiguous intervals. 

180 MINIMUM_FILE_FORMAT_VERSIONS = ("0.1",) 

181 MAXIMUM_FILE_FORMAT_VERSIONS = ("1.0",) 

182 

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

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

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

186 

187 self.filename = filename 

188 

189 @classmethod 

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

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

192 

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

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

195 

196 Parameters 

197 ---------- 

198 written_version: `str` 

199 The VERSION of the file to be read. 

200 

201 Returns 

202 ------- 

203 compatible : `bool` 

204 Whether the reader can read a file whose VERSION is 

205 ``written_version``. 

206 

207 Notes 

208 ----- 

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

210 """ 

211 written_version_object = version.parse(written_version) 

212 for min_version, max_version in zip( 

213 cls.MINIMUM_FILE_FORMAT_VERSIONS, 

214 cls.MAXIMUM_FILE_FORMAT_VERSIONS, 

215 strict=True, 

216 ): 

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

218 return True 

219 

220 return False 

221 

222 def readAsMultipleCellCoadd(self) -> MultipleCellCoadd: 

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

224 

225 Raises 

226 ------ 

227 IncompatibleError 

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

229 incompatible with this module that is reading it in. 

230 ValueError 

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

232 direction or by an even number of pixels. 

233 """ 

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

235 header = hdu_list[1].header 

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

237 if not self.isCompatibleWith(written_version): 

238 raise IncompatibleVersionError( 

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

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

241 ) 

242 if written_version != FILE_FORMAT_VERSION: 

243 logger.info( 

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

245 self.filename, 

246 written_version, 

247 FILE_FORMAT_VERSION, 

248 ) 

249 

250 written_version = version.parse(written_version) 

251 

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

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

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

255 

256 data = hdu_list[1].data 

257 

258 # Read in WCS 

259 ps = PropertySet() 

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

261 wcs = afwGeom.makeSkyWcs(ps) 

262 

263 # Build the quantities needed to construct a MultipleCellCoadd. 

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

265 coadd_units = header["TUNIT2"] 

266 else: 

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

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

269 

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

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

272 

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

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

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

276 grid_padding_extent = outer_cell_size - grid_cell_size 

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

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

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

280 if grid_padding_extent.x != grid_padding_extent.y: 

281 raise ValueError( 

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

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

284 ) 

285 if grid_padding_extent.x % 2 != 0: 

286 raise ValueError( 

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

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

289 ) 

290 grid_padding = grid_padding_extent.x // 2 

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

292 

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

294 inner_bbox = Box2I( 

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

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

297 ) 

298 

299 # Attempt to get inputs for each cell. 

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

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

302 overlaps_center=True, 

303 overlap_fraction=1.0, 

304 unmasked_overlap_fraction=1.0, 

305 weight=0.0, 

306 psf_shape=afwGeom.Quadrupole(), 

307 psf_shape_flag=True, 

308 ) 

309 visit_polygons = {} 

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

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

312 visit_hdu = hdu_list[hdu_list.index_of("VISIT_SUMMARY")].data 

313 else: 

314 visit_hdu = hdu_list[hdu_list.index_of("VISIT")].data 

315 

316 visit_dict = { 

317 row["visit"]: VisitRecord( 

318 visit=row["visit"], 

319 physical_filter=( 

320 row["physical_filter"] 

321 if written_version >= version.parse("0.8") 

322 else header["FILTER"] 

323 ), 

324 day_obs=row["day_obs"] if written_version >= version.parse("0.8") else 00000000, 

325 ) 

326 for row in visit_hdu 

327 } 

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

329 

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

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

332 mask_plane_dict = {} 

333 for mask_name, mask_value in mask_definitions: 

334 afwImage.Mask.addMaskPlane(mask_name) 

335 mask_plane_dict[mask_name] = mask_value 

336 else: 

337 logger.warning( 

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

339 "FILE_FORMAT_VERSION %s < 0.6", 

340 header["VERSION"], 

341 ) 

342 mask_plane_dict = None 

343 

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

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

346 for row in visit_summary_hdu.data: 

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

348 obs_id = ObservationIdentifiers( 

349 instrument=header["INSTRUME"], 

350 visit=visit, 

351 detector=row["detector"], 

352 day_obs=visit_dict[visit].day_obs, 

353 physical_filter=visit_dict[visit].physical_filter, 

354 ) 

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

356 num_vertices = row["num_vertices"] 

357 else: 

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

359 visit_polygons[obs_id] = afwGeom.Polygon( 

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

361 ) 

362 

363 for link_row in link_table: 

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

365 visit = link_row["visit"] 

366 obs_id = ObservationIdentifiers( 

367 instrument=header["INSTRUME"], 

368 visit=visit, 

369 detector=link_row["detector"], 

370 day_obs=visit_dict[visit].day_obs, 

371 physical_filter=visit_dict[visit].physical_filter, 

372 ) 

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

374 coadd_input = CoaddInputs( 

375 overlaps_center=link_row["overlaps_center"], 

376 overlap_fraction=link_row["overlap_fraction"], 

377 unmasked_overlap_fraction=( 

378 link_row["unmasked_overlap_fraction"] 

379 if written_version >= version.parse("0.9") 

380 else link_row["overlap_fraction"] 

381 ), 

382 weight=link_row["weight"], 

383 psf_shape=afwGeom.Quadrupole( 

384 ixx=link_row["psf_shape_ixx"], 

385 iyy=link_row["psf_shape_iyy"], 

386 ixy=link_row["psf_shape_ixy"], 

387 ), 

388 psf_shape_flag=link_row["psf_shape_flag"], 

389 ) 

390 

391 if cell_id in inputs: 

392 inputs[cell_id][obs_id] = coadd_input 

393 else: 

394 inputs[cell_id] = {obs_id: coadd_input} 

395 else: 

396 logger.info( 

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

398 "VERSION = %s", 

399 written_version, 

400 ) 

401 

402 try: 

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

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

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

406 aperture_correction_grid = self._readApCorr(aperture_correction_hdu, grid_shape) 

407 except KeyError: 

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

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

410 # Create an empty grid container regardless. 

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

412 

413 common = CommonComponents( 

414 units=CoaddUnits(coadd_units), 

415 wcs=wcs, 

416 band=header["FILTER"], 

417 identifiers=PatchIdentifiers( 

418 skymap=header["SKYMAP"], 

419 tract=header["TRACT"], 

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

421 band=header["FILTER"], 

422 ), 

423 visit_polygons=visit_polygons, 

424 ) 

425 

426 coadd = MultipleCellCoadd( 

427 ( 

428 self._readSingleCellCoadd( 

429 data=data[row_id], 

430 header=header, 

431 common=common, 

432 inputs=inputs.get( 

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

434 {}, 

435 ), 

436 outer_cell_size=outer_cell_size, 

437 psf_image_size=psf_image_size, 

438 inner_cell_size=grid_cell_size, 

439 mask_plane_dict=mask_plane_dict, 

440 aperture_correction_map=aperture_correction_grid.get( 

441 Index2D( 

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

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

444 ), 

445 EMPTY_AP_CORR_MAP, 

446 ), 

447 ) 

448 for row_id in range(len(data)) 

449 ), 

450 grid=grid, 

451 outer_cell_size=outer_cell_size, 

452 psf_image_size=psf_image_size, 

453 inner_bbox=inner_bbox, 

454 common=common, 

455 ) 

456 

457 return coadd 

458 

459 @staticmethod 

460 def _readSingleCellCoadd( 

461 data: Mapping[str, Any], 

462 common: CommonComponents, 

463 header: Mapping[str, Any], 

464 *, 

465 inputs: Mapping[ObservationIdentifiers, CoaddInputs], 

466 outer_cell_size: Extent2I, 

467 inner_cell_size: Extent2I, 

468 psf_image_size: Extent2I, 

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

470 aperture_correction_map: SingleCellCoaddApCorrMap = EMPTY_AP_CORR_MAP, 

471 ) -> SingleCellCoadd: 

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

473 

474 Parameters 

475 ---------- 

476 data : `Mapping` 

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

478 table representation. 

479 common : `CommonComponents` 

480 The common components of the coadd. 

481 header : `Mapping` 

482 The header of the FITS file as a dictionary. 

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

484 A mapping of ObservationIdentifiers that contributed to this cell 

485 to their corresponding CoaddInputs. 

486 outer_cell_size : `Extent2I` 

487 The size of the outer cell. 

488 psf_image_size : `Extent2I` 

489 The size of the PSF image. 

490 inner_cell_size : `Extent2I` 

491 The size of the inner cell. 

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

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

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

495 Mapping of algorithm name to aperture correction value. 

496 If None, no aperture correction is applied. 

497 

498 Returns 

499 ------- 

500 coadd : `SingleCellCoadd` 

501 The coadd read from the file. 

502 """ 

503 buffer = (outer_cell_size - inner_cell_size) // 2 

504 

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

506 

507 psf = ImageD( 

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

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

510 ) # use the variable 

511 xy0 = Point2I( 

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

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

514 ) 

515 

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

517 # the definitions. 

518 mask = afwImage.Mask( 

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

520 initialValue=0, 

521 ) 

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

523 if mask_plane_dict is not None: 

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

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

526 else: 

527 mask.array[:, :] = mask_array 

528 

529 try: 

530 maskfrac = data["maskfrac"] 

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

532 except KeyError: 

533 mask_fractions = None 

534 

535 image_planes = OwnedImagePlanes( 

536 image=ImageF( 

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

538 xy0=xy0, 

539 ), 

540 mask=mask, 

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

542 noise_realizations=[ 

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

544 ], 

545 mask_fractions=mask_fractions, 

546 ) 

547 

548 identifiers = CellIdentifiers( 

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

550 skymap=common.identifiers.skymap, 

551 tract=common.identifiers.tract, 

552 patch=common.identifiers.patch, 

553 band=common.identifiers.band, 

554 ) 

555 

556 return SingleCellCoadd( 

557 outer=image_planes, 

558 psf=psf, 

559 inner_bbox=Box2I( 

560 corner=Point2I( 

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

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

563 ), 

564 dimensions=inner_cell_size, 

565 ), 

566 common=common, 

567 identifiers=identifiers, 

568 inputs=inputs, 

569 aperture_correction_map=aperture_correction_map, 

570 ) 

571 

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

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

574 

575 Returns 

576 ------- 

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

578 The WCS information read from the FITS file. 

579 """ 

580 # Read in WCS 

581 ps = PropertySet() 

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

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

584 wcs = afwGeom.makeSkyWcs(ps) 

585 return wcs 

586 

587 def _readApCorr( 

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

589 ) -> GridContainer[SingleCellCoaddApCorrMap]: 

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

591 column_names = aperture_correction_hdu.names 

592 column_names.remove("x") 

593 column_names.remove("y") 

594 

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

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

597 cell_id = Index2D( 

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

599 ) 

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

601 gc[cell_id] = frozendict(ap_corr_map_dict) 

602 

603 return gc 

604 

605 

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

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

608 appropriate data types. 

609 

610 Parameters 

611 ---------- 

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

613 Aperture correction values for each algorithm 

614 xy: `~lsst.skymap.Index2D` 

615 An object identifying the cell 

616 dtypes: `~numpy.typing.DTypeLike` 

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

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

619 

620 Returns 

621 ------- 

622 record : `np.record` 

623 Per-cell aperture correction map as a record. 

624 """ 

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

626 for field_name in ap_corr_map: 

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

628 

629 record["x"] = xy.x 

630 record["y"] = xy.y 

631 

632 return record 

633 

634 

635def writeMultipleCellCoaddAsFits( 

636 multiple_cell_coadd: MultipleCellCoadd, 

637 filename: str, 

638 overwrite: bool = False, 

639 metadata: PropertySet | None = None, 

640) -> fits.HDUList: 

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

642 

643 Parameters 

644 ---------- 

645 multiple_cell_coadd : `MultipleCellCoadd` 

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

647 filename : `str` 

648 The name of the file to write to. 

649 overwrite : `bool`, optional 

650 Whether to overwrite the file if it already exists? 

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

652 Additional metadata to write to the FITS file. 

653 

654 Returns 

655 ------- 

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

657 The FITS file as an HDUList. 

658 

659 Notes 

660 ----- 

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

662 must be accompanied with a change to FILE_FORMAT_VERSION. 

663 """ 

664 # Create metadata tables: 

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

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

667 visit_records: list[Any] = [] 

668 cell_records: list[Any] = [] 

669 instrument_set = set() 

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

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

672 visit_records.append( 

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

674 ) 

675 cell_records.append( 

676 ( 

677 cell_id.x, 

678 cell_id.y, 

679 observation_id.visit, 

680 observation_id.detector, 

681 coadd_input.overlaps_center, 

682 coadd_input.overlap_fraction, 

683 coadd_input.unmasked_overlap_fraction, 

684 coadd_input.weight, 

685 coadd_input.psf_shape.getIxx(), 

686 coadd_input.psf_shape.getIyy(), 

687 coadd_input.psf_shape.getIxy(), 

688 coadd_input.psf_shape_flag, 

689 ) 

690 ) 

691 instrument_set.add(observation_id.instrument) 

692 

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

694 instrument = instrument_set.pop() 

695 

696 # Create visit_summary equivalent table 

697 visit_column = fits.Column( 

698 name="visit", 

699 format="K", 

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

701 ) 

702 detector_column = fits.Column( 

703 name="detector", 

704 format="I", 

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

706 ) 

707 physical_filter_column = fits.Column( 

708 name="physical_filter", 

709 format="32A", 

710 array=[obs_id.physical_filter for obs_id in multiple_cell_coadd.common.visit_polygons], 

711 ) 

712 day_obs_column = fits.Column( 

713 name="day_obs", 

714 format="K", 

715 array=[obs_id.day_obs for obs_id in multiple_cell_coadd.common.visit_polygons], 

716 ) 

717 number_of_vertices = [] 

718 polygon_vertices_array = [] 

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

720 # Polygons are closed, so the first and last vertices are the same. 

721 if (num_vertices := len(poly.getVertices())) > MAX_POLYGON_VERTICES + 1: 

722 logger.warning( 

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

724 "This geometry should be impossible for two intersecting " 

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

726 "the FITS file.", 

727 obs_id.visit, 

728 obs_id.detector, 

729 num_vertices, 

730 MAX_POLYGON_VERTICES, 

731 ) 

732 number_of_vertices.append(num_vertices) 

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

734 vertices = vertices[:MAX_POLYGON_VERTICES] 

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

736 polygon_column = fits.Column( 

737 name="polygon_vertices", 

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

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

740 array=polygon_vertices_array, 

741 ) 

742 number_of_vertices_column = fits.Column( 

743 name="num_vertices", 

744 format="I", 

745 array=number_of_vertices, 

746 ) 

747 visit_summary_hdu = fits.BinTableHDU.from_columns( 

748 [ 

749 visit_column, 

750 detector_column, 

751 physical_filter_column, 

752 day_obs_column, 

753 number_of_vertices_column, 

754 polygon_column, 

755 ], 

756 name="VISIT_SUMMARY", 

757 ) 

758 visit_summary_hdu.header["POLYVERT"] = MAX_POLYGON_VERTICES 

759 

760 visit_recarray = np.rec.fromrecords( 

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

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

763 names=( 

764 "visit", 

765 "physical_filter", 

766 "day_obs", 

767 ), 

768 ) 

769 cell_recarray = np.rec.fromrecords( 

770 recList=cell_records, 

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

772 names=( 

773 "cell_x", 

774 "cell_y", 

775 "visit", 

776 "detector", 

777 "overlaps_center", 

778 "overlap_fraction", 

779 "unmasked_overlap_fraction", 

780 "weight", 

781 "psf_shape_ixx", 

782 "psf_shape_iyy", 

783 "psf_shape_ixy", 

784 "psf_shape_flag", 

785 ), 

786 ) 

787 

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

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

790 

791 cell_id = fits.Column( 

792 name="cell_id", 

793 format="2I", 

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

795 ) 

796 

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

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

799 image = fits.Column( 

800 name="image", 

801 unit=unit_array[0], 

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

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

804 array=image_array, 

805 ) 

806 

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

808 mask = fits.Column( 

809 name="mask", 

810 format=f"{mask_array[0].size}J", 

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

812 array=mask_array, 

813 ) 

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

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

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

817 mask_definition_table.add_row((name, value)) 

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

819 

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

821 variance = fits.Column( 

822 name="variance", 

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

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

825 array=variance_array, 

826 ) 

827 

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

829 psf = fits.Column( 

830 name="psf", 

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

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

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

834 ) 

835 

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

837 maskfrac_array = [ 

838 cell.outer.mask_fractions.array 

839 for cell in multiple_cell_coadd.cells.values() 

840 if cell.outer.mask_fractions 

841 ] 

842 if maskfrac_array: 

843 maskfrac = fits.Column( 

844 name="maskfrac", 

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

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