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

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 if written_version >= version.parse("0.6"): 

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

303 else: 

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

305 

306 visit_dict = { 

307 row["visit"]: VisitRecord( 

308 visit=row["visit"], 

309 physical_filter=( 

310 row["physical_filter"] 

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

312 else header["FILTER"] 

313 ), 

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

315 ) 

316 for row in visit_hdu 

317 } 

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

319 

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

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

322 mask_plane_dict = {} 

323 for mask_name, mask_value in mask_definitions: 

324 afwImage.Mask.addMaskPlane(mask_name) 

325 mask_plane_dict[mask_name] = mask_value 

326 else: 

327 logger.warning( 

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

329 "FILE_FORMAT_VERSION %s < 0.6", 

330 header["VERSION"], 

331 ) 

332 mask_plane_dict = None 

333 

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

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

336 for row in visit_summary_hdu.data: 

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

338 obs_id = ObservationIdentifiers( 

339 instrument=header["INSTRUME"], 

340 visit=visit, 

341 detector=row["detector"], 

342 day_obs=visit_dict[visit].day_obs, 

343 physical_filter=visit_dict[visit].physical_filter, 

344 ) 

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

346 num_vertices = row["num_vertices"] 

347 else: 

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

349 visit_polygons[obs_id] = afwGeom.Polygon( 

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

351 ) 

352 

353 for link_row in link_table: 

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

355 visit = link_row["visit"] 

356 obs_id = ObservationIdentifiers( 

357 instrument=header["INSTRUME"], 

358 visit=visit, 

359 detector=link_row["detector"], 

360 day_obs=visit_dict[visit].day_obs, 

361 physical_filter=visit_dict[visit].physical_filter, 

362 ) 

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

364 coadd_input = CoaddInputs( 

365 overlaps_center=link_row["overlaps_center"], 

366 overlap_fraction=link_row["overlap_fraction"], 

367 weight=link_row["weight"], 

368 psf_shape=afwGeom.Quadrupole( 

369 ixx=link_row["psf_shape_ixx"], 

370 iyy=link_row["psf_shape_iyy"], 

371 ixy=link_row["psf_shape_ixy"], 

372 ), 

373 psf_shape_flag=link_row["psf_shape_flag"], 

374 ) 

375 

376 if cell_id in inputs: 

377 inputs[cell_id][obs_id] = coadd_input 

378 else: 

379 inputs[cell_id] = {obs_id: coadd_input} 

380 else: 

381 logger.info( 

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

383 "VERSION = %s", 

384 written_version, 

385 ) 

386 

387 try: 

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

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

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

391 aperture_correction_grid = self._readApCorr(aperture_correction_hdu, grid_shape) 

392 except KeyError: 

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

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

395 # Create an empty grid container regardless. 

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

397 

398 common = CommonComponents( 

399 units=CoaddUnits(coadd_units), 

400 wcs=wcs, 

401 band=header["FILTER"], 

402 identifiers=PatchIdentifiers( 

403 skymap=header["SKYMAP"], 

404 tract=header["TRACT"], 

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

406 band=header["FILTER"], 

407 ), 

408 visit_polygons=visit_polygons, 

409 ) 

410 

411 coadd = MultipleCellCoadd( 

412 ( 

413 self._readSingleCellCoadd( 

414 data=data[row_id], 

415 header=header, 

416 common=common, 

417 inputs=inputs[ 

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

419 ], 

420 outer_cell_size=outer_cell_size, 

421 psf_image_size=psf_image_size, 

422 inner_cell_size=grid_cell_size, 

423 mask_plane_dict=mask_plane_dict, 

424 aperture_correction_map=aperture_correction_grid.get( 

425 Index2D( 

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

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

428 ), 

429 EMPTY_AP_CORR_MAP, 

430 ), 

431 ) 

432 for row_id in range(len(data)) 

433 ), 

434 grid=grid, 

435 outer_cell_size=outer_cell_size, 

436 psf_image_size=psf_image_size, 

437 inner_bbox=inner_bbox, 

438 common=common, 

439 ) 

440 

441 return coadd 

442 

443 @staticmethod 

444 def _readSingleCellCoadd( 

445 data: Mapping[str, Any], 

446 common: CommonComponents, 

447 header: Mapping[str, Any], 

448 *, 

449 inputs: Mapping[ObservationIdentifiers, CoaddInputs], 

450 outer_cell_size: Extent2I, 

451 inner_cell_size: Extent2I, 

452 psf_image_size: Extent2I, 

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

454 aperture_correction_map: SingleCellCoaddApCorrMap = EMPTY_AP_CORR_MAP, 

455 ) -> SingleCellCoadd: 

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

457 

458 Parameters 

459 ---------- 

460 data : `Mapping` 

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

462 table representation. 

463 common : `CommonComponents` 

464 The common components of the coadd. 

465 header : `Mapping` 

466 The header of the FITS file as a dictionary. 

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

468 A mapping of ObservationIdentifiers that contributed to this cell 

469 to their corresponding CoaddInputs. 

470 outer_cell_size : `Extent2I` 

471 The size of the outer cell. 

472 psf_image_size : `Extent2I` 

473 The size of the PSF image. 

474 inner_cell_size : `Extent2I` 

475 The size of the inner cell. 

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

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

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

479 Mapping of algorithm name to aperture correction value. 

480 If None, no aperture correction is applied. 

481 

482 Returns 

483 ------- 

484 coadd : `SingleCellCoadd` 

485 The coadd read from the file. 

486 """ 

487 buffer = (outer_cell_size - inner_cell_size) // 2 

488 

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

490 

491 psf = ImageD( 

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

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

494 ) # use the variable 

495 xy0 = Point2I( 

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

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

498 ) 

499 

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

501 # the definitions. 

502 mask = afwImage.Mask( 

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

504 initialValue=0, 

505 ) 

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

507 if mask_plane_dict is not None: 

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

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

510 else: 

511 mask.array[:, :] = mask_array 

512 

513 try: 

514 maskfrac = data["maskfrac"] 

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

516 except KeyError: 

517 mask_fractions = None 

518 

519 image_planes = OwnedImagePlanes( 

520 image=ImageF( 

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

522 xy0=xy0, 

523 ), 

524 mask=mask, 

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

526 noise_realizations=[ 

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

528 ], 

529 mask_fractions=mask_fractions, 

530 ) 

531 

532 identifiers = CellIdentifiers( 

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

534 skymap=common.identifiers.skymap, 

535 tract=common.identifiers.tract, 

536 patch=common.identifiers.patch, 

537 band=common.identifiers.band, 

538 ) 

539 

540 return SingleCellCoadd( 

541 outer=image_planes, 

542 psf=psf, 

543 inner_bbox=Box2I( 

544 corner=Point2I( 

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

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

547 ), 

548 dimensions=inner_cell_size, 

549 ), 

550 common=common, 

551 identifiers=identifiers, 

552 inputs=inputs, 

553 aperture_correction_map=aperture_correction_map, 

554 ) 

555 

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

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

558 

559 Returns 

560 ------- 

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

562 The WCS information read from the FITS file. 

563 """ 

564 # Read in WCS 

565 ps = PropertySet() 

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

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

568 wcs = afwGeom.makeSkyWcs(ps) 

569 return wcs 

570 

571 def _readApCorr( 

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

573 ) -> GridContainer[SingleCellCoaddApCorrMap]: 

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

575 column_names = aperture_correction_hdu.names 

576 column_names.remove("x") 

577 column_names.remove("y") 

578 

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

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

581 cell_id = Index2D( 

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

583 ) 

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

585 gc[cell_id] = frozendict(ap_corr_map_dict) 

586 

587 return gc 

588 

589 

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

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

592 appropriate data types. 

593 

594 Parameters 

595 ---------- 

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

597 Aperture correction values for each algorithm 

598 xy: `~lsst.skymap.Index2D` 

599 An object identifying the cell 

600 dtypes: `~numpy.typing.DTypeLike` 

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

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

603 

604 Returns 

605 ------- 

606 record : `np.record` 

607 Per-cell aperture correction map as a record. 

608 """ 

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

610 for field_name in ap_corr_map: 

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

612 

613 record["x"] = xy.x 

614 record["y"] = xy.y 

615 

616 return record 

617 

618 

619def writeMultipleCellCoaddAsFits( 

620 multiple_cell_coadd: MultipleCellCoadd, 

621 filename: str, 

622 overwrite: bool = False, 

623 metadata: PropertySet | None = None, 

624) -> fits.HDUList: 

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

626 

627 Parameters 

628 ---------- 

629 multiple_cell_coadd : `MultipleCellCoadd` 

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

631 filename : `str` 

632 The name of the file to write to. 

633 overwrite : `bool`, optional 

634 Whether to overwrite the file if it already exists? 

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

636 Additional metadata to write to the FITS file. 

637 

638 Returns 

639 ------- 

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

641 The FITS file as an HDUList. 

642 

643 Notes 

644 ----- 

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

646 must be accompanied with a change to FILE_FORMAT_VERSION. 

647 """ 

648 # Create metadata tables: 

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

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

651 visit_records: list[Any] = [] 

652 cell_records: list[Any] = [] 

653 instrument_set = set() 

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

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

656 visit_records.append( 

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

658 ) 

659 cell_records.append( 

660 ( 

661 cell_id.x, 

662 cell_id.y, 

663 observation_id.visit, 

664 observation_id.detector, 

665 coadd_input.overlaps_center, 

666 coadd_input.overlap_fraction, 

667 coadd_input.weight, 

668 coadd_input.psf_shape.getIxx(), 

669 coadd_input.psf_shape.getIyy(), 

670 coadd_input.psf_shape.getIxy(), 

671 coadd_input.psf_shape_flag, 

672 ) 

673 ) 

674 instrument_set.add(observation_id.instrument) 

675 

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

677 instrument = instrument_set.pop() 

678 

679 # Create visit_summary equivalent table 

680 visit_column = fits.Column( 

681 name="visit", 

682 format="K", 

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

684 ) 

685 detector_column = fits.Column( 

686 name="detector", 

687 format="I", 

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

689 ) 

690 physical_filter_column = fits.Column( 

691 name="physical_filter", 

692 format="32A", 

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

694 ) 

695 day_obs_column = fits.Column( 

696 name="day_obs", 

697 format="K", 

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

699 ) 

700 number_of_vertices = [] 

701 polygon_vertices_array = [] 

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

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

704 logger.warning( 

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

706 "This geometry should be impossible for two intersecting " 

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

708 "the FITS file.", 

709 obs_id.visit, 

710 obs_id.detector, 

711 num_vertices, 

712 MAX_POLYGON_VERTICES, 

713 ) 

714 number_of_vertices.append(num_vertices) 

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

716 vertices = vertices[:MAX_POLYGON_VERTICES] 

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

718 polygon_column = fits.Column( 

719 name="polygon_vertices", 

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

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

722 array=polygon_vertices_array, 

723 ) 

724 number_of_vertices_column = fits.Column( 

725 name="num_vertices", 

726 format="I", 

727 array=number_of_vertices, 

728 ) 

729 visit_summary_hdu = fits.BinTableHDU.from_columns( 

730 [ 

731 visit_column, 

732 detector_column, 

733 physical_filter_column, 

734 day_obs_column, 

735 number_of_vertices_column, 

736 polygon_column, 

737 ], 

738 name="VISIT_SUMMARY", 

739 ) 

740 visit_summary_hdu.header["POLYVERT"] = MAX_POLYGON_VERTICES 

741 

742 visit_recarray = np.rec.fromrecords( 

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

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

745 names=( 

746 "visit", 

747 "physical_filter", 

748 "day_obs", 

749 ), 

750 ) 

751 cell_recarray = np.rec.fromrecords( 

752 recList=cell_records, 

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

754 names=( 

755 "cell_x", 

756 "cell_y", 

757 "visit", 

758 "detector", 

759 "overlaps_center", 

760 "overlap_fraction", 

761 "weight", 

762 "psf_shape_ixx", 

763 "psf_shape_iyy", 

764 "psf_shape_ixy", 

765 "psf_shape_flag", 

766 ), 

767 ) 

768 

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

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

771 

772 cell_id = fits.Column( 

773 name="cell_id", 

774 format="2I", 

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

776 ) 

777 

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

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

780 image = fits.Column( 

781 name="image", 

782 unit=unit_array[0], 

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

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

785 array=image_array, 

786 ) 

787 

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

789 mask = fits.Column( 

790 name="mask", 

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

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

793 array=mask_array, 

794 ) 

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

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

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

798 mask_definition_table.add_row((name, value)) 

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

800 

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

802 variance = fits.Column( 

803 name="variance", 

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

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

806 array=variance_array, 

807 ) 

808 

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

810 psf = fits.Column( 

811 name="psf", 

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

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

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

815 ) 

816 

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

818 maskfrac_array = [ 

819 cell.outer.mask_fractions.array 

820 for cell in multiple_cell_coadd.cells.values() 

821 if cell.outer.mask_fractions 

822 ] 

823 if maskfrac_array: 

824 maskfrac = fits.Column( 

825 name="maskfrac", 

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

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

828 array=maskfrac_array, 

829 ) 

830 columns.append(maskfrac) 

831 

832 n_noise_realizations = multiple_cell_coadd.n_noise_realizations 

833 for n in range(n_noise_realizations): 

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

835 columns.append( 

836 fits.Column( 

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

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

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