Coverage for python/lsst/obs/base/formatters/fitsExposure.py: 83%

1# This file is part of obs_base. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (http://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <http://www.gnu.org/licenses/>. 

21 

22__all__ = ( 

23 "FitsExposureFormatter", 

24 "FitsImageFormatter", 

25 "FitsImageFormatterBase", 

26 "FitsMaskFormatter", 

27 "FitsMaskedImageFormatter", 

28 "StandardFitsImageFormatterBase", 

29 "standardizeAmplifierParameters", 

30) 

31 

32import hashlib 

33import json 

34import logging 

35import threading 

36import uuid 

37import warnings 

38from abc import abstractmethod 

39from collections.abc import Mapping, Set 

40from io import BytesIO 

41from typing import TYPE_CHECKING, Any, ClassVar, NamedTuple, Protocol 

42 

43import astropy.io.fits 

44import numpy as np 

45 

46import lsst.geom 

47from lsst.afw.cameraGeom import AmplifierGeometryComparison, AmplifierIsolator 

48from lsst.afw.fits import CompressionOptions, MemFileManager 

49from lsst.afw.geom.wcsUtils import getImageXY0FromMetadata 

50from lsst.afw.image import ( 

51 ExposureFitsReader, 

52 ExposureInfo, 

53 FilterLabel, 

54 ImageFitsReader, 

55 MaskedImageFitsReader, 

56 MaskFitsReader, 

57) 

58 

59# Needed for ApCorrMap to resolve properly 

60from lsst.afw.math import BoundedField # noqa: F401 

61from lsst.daf.base import PropertyList 

62from lsst.daf.butler import DatasetProvenance, FormatterV2 

63from lsst.resources import ResourcePath 

64from lsst.utils.classes import cached_getter 

65from lsst.utils.introspection import find_outside_stacklevel 

66 

67from ..utils import add_provenance_to_fits_header 

68 

69if TYPE_CHECKING: 

70 import lsst.afw.cameraGeom 

71 

72 

73_LOG = logging.getLogger(__name__) 

74 

75_ALWAYS_USE_ASTROPY_FOR_COMPONENT_READ = False 

76"""If True, the astropy code will always be used to read component and cutouts 

77even if the file is local, the cutout is too large, or the dataset type is 

78wrong. This should mostly be used for testing. 

79""" 

80 

81 

82class _ReaderClassLike(Protocol): 

83 def __init__(self, path: str) -> None: ... 

84 def readBBox(self) -> lsst.geom.Box2I: ... 

85 def read(self, bbox: lsst.geom.Box2I = lsst.geom.Box2I(), dtype: Any = None) -> Any: ... 

86 def readImage(self, bbox: lsst.geom.Box2I = lsst.geom.Box2I(), dtype: Any = None) -> Any: ... 

87 def readMask(self, bbox: lsst.geom.Box2I = lsst.geom.Box2I(), dtype: Any = None) -> Any: ... 

88 def readVariance(self, bbox: lsst.geom.Box2I = lsst.geom.Box2I(), dtype: Any = None) -> Any: ... 

89 def readDetector(self) -> lsst.afw.cameraGeom.Detector: ... 

90 def readComponent(self, component: str) -> Any: ... 

91 def readMetadata(self) -> PropertyList: ... 

92 def readSerializationVersion(self) -> int: ... 

93 

94 

95class FitsImageFormatterBase(FormatterV2): 

96 """Base class formatter for image-like storage classes stored via FITS. 

97 

98 Notes 

99 ----- 

100 This class makes no assumptions about how many HDUs are used to represent 

101 the image on disk, and includes no support for writing. It's really just a 

102 collection of miscellaneous boilerplate common to all FITS image 

103 formatters. 

104 

105 Concrete subclasses must implement `readComponent`, `readFull`, and 

106 `write_local_file` (even if just to disable them by raising an exception). 

107 """ 

108 

109 can_read_from_local_file = True 

110 default_extension = ".fits" 

111 supported_extensions: ClassVar[Set[str]] = frozenset({".fits", ".fits.gz", ".fits.fz", ".fz", ".fit"}) 

112 

113 unsupported_parameters: ClassVar[Set[str]] = frozenset() 

114 """Support all parameters.""" 

115 

116 _reader = None 

117 _reader_path: str | None = None 

118 

119 ReaderClass: type[_ReaderClassLike] # must be set by concrete subclasses 

120 """Class to use for reading FITS files in the expected way. 

121 (e.g., `type` [`lsst.afw.image.ImageFitsReader]) 

122 """ 

123 

124 @property 

125 def reader(self) -> _ReaderClassLike: 

126 """The reader object that backs this formatter's read operations. 

127 

128 This is computed on first use and then cached. It should never be 

129 accessed when writing. Currently assumes a local file. 

130 """ 

131 if self._reader is None: 

132 if self._reader_path is None: 132 ↛ 133line 132 didn't jump to line 133 because the condition on line 132 was never true

133 raise RuntimeError("Internal error in formatter; failing to set path.") 

134 self._reader = self.ReaderClass(self._reader_path) 

135 return self._reader 

136 

137 @property 

138 @cached_getter 

139 def checked_parameters(self) -> dict[str, Any]: 

140 """The parameters passed by the butler user, after checking them 

141 against the storage class and transforming `None` into an empty `dict` 

142 (`dict`). 

143 

144 This is computed on first use and then cached. It should never be 

145 accessed when writing. Subclasses that need additional checking should 

146 delegate to `super` and then check the result before returning it. 

147 """ 

148 parameters = self.file_descriptor.parameters 

149 if parameters is None: 

150 parameters = {} 

151 self.file_descriptor.storageClass.validateParameters(parameters) 

152 return parameters 

153 

154 @property 

155 def storageClass_dtype(self) -> np.dtype | None: 

156 """The numpy data type associated with the storage class.""" 

157 dtype: np.dtype | None = None 

158 try: 

159 # lsst.afw.image.Exposure is generic base class and does not have 

160 # the dtype attribute. 

161 dtype = np.dtype(self.file_descriptor.storageClass.pytype.dtype) # type: ignore[attr-defined] 

162 except AttributeError: 

163 pass 

164 return dtype 

165 

166 def read_from_local_file(self, path: str, component: str | None = None, expected_size: int = -1) -> Any: 

167 # Docstring inherited. 

168 if future_type := _get_future_image_type(self.file_descriptor.readStorageClass.name, component): 168 ↛ 169line 168 didn't jump to line 169 because the condition on line 168 was never true

169 return future_type.read_legacy( 

170 path, 

171 component=component, 

172 preserve_quantization=self.checked_parameters.get("preserve_quantization", False), 

173 ) 

174 elif self.checked_parameters.get("preserve_quantization", False): 174 ↛ 175line 174 didn't jump to line 175 because the condition on line 174 was never true

175 raise NotImplementedError( 

176 "preserve_quantization=True only works when converting to VisitImage on read." 

177 ) 

178 

179 # The methods doing the reading all currently assume local file 

180 # and assume that the file descriptor refers to a local file. 

181 # With FormatterV2 that file descriptor does not refer to a local 

182 # file. 

183 self._reader_path = path 

184 self._reader = None # Ensure the reader class is reset. 

185 try: 

186 if component is not None: 

187 in_memory_dataset = self.readComponent(component) 

188 else: 

189 in_memory_dataset = self.readFull() 

190 finally: 

191 self._reader = None # Release the file handle. 

192 return in_memory_dataset 

193 

194 @abstractmethod 

195 def readComponent(self, component: str) -> Any: 

196 """Read a component dataset. 

197 

198 Parameters 

199 ---------- 

200 component : `str`, optional 

201 Component to read from the file. 

202 

203 Returns 

204 ------- 

205 obj : `typing.Any` 

206 In-memory component object. 

207 

208 Raises 

209 ------ 

210 KeyError 

211 Raised if the requested component cannot be handled. 

212 """ 

213 raise NotImplementedError() 

214 

215 @abstractmethod 

216 def readFull(self) -> Any: 

217 """Read the full dataset (while still accounting for parameters). 

218 

219 Returns 

220 ------- 

221 obj : `typing.Any` 

222 In-memory component object. 

223 

224 """ 

225 raise NotImplementedError() 

226 

227 

228class StandardFitsImageFormatterBase(FitsImageFormatterBase): 

229 """Base class interface for image-like storage stored via FITS, 

230 written using LSST code. 

231 

232 Notes 

233 ----- 

234 Concrete subclasses must provide at least the ``ReaderClass`` attribute. 

235 

236 The provided implementation of `readComponent` handles only the 'bbox', 

237 'dimensions', and 'xy0' components common to all image-like storage 

238 classes. Subclasses with additional components should handle them first, 

239 then delegate to ``super()`` for these (or, if necessary, delegate first 

240 and catch `KeyError`). 

241 

242 The provided implementation of `readFull` handles only parameters that 

243 can be forwarded directly to the reader class (usually ``bbox`` and 

244 ``origin``). Concrete subclasses that need to handle additional parameters 

245 should generally reimplement without delegating (the implementation is 

246 trivial). 

247 

248 This Formatter supports write recipes, and assumes its in-memory type has 

249 ``writeFits`` and (for write recipes) ``writeFitsWithOptions`` methods. 

250 

251 Each ``StandardFitsImageFormatterBase`` recipe for FITS compression should 

252 define ``image``, ``mask`` and ``variance`` entries, each of which may 

253 contain entries supported by 

254 `lsst.afw.fits.CompressionOptions.from_mapping` (``null`` disables 

255 compression). 

256 

257 A very simple example YAML recipe (for the `lsst.afw.image.Exposure` 

258 specialization): 

259 

260 .. code-block:: yaml 

261 

262 lsst.obs.base.fitsExposureFormatter.FitsExposureFormatter: 

263 default: 

264 image: &default 

265 algorithm: GZIP_2 

266 mask: *default 

267 variance: *default 

268 

269 """ 

270 

271 supported_write_parameters = frozenset({"recipe"}) 

272 

273 def readComponent(self, component: str) -> Any: 

274 # Docstring inherited. 

275 if component in ("bbox", "dimensions", "xy0"): 275 ↛ 284line 275 didn't jump to line 284 because the condition on line 275 was always true

276 bbox = self.reader.readBBox() 

277 if component == "dimensions": 

278 return bbox.getDimensions() 

279 elif component == "xy0": 

280 return bbox.getMin() 

281 else: 

282 return bbox 

283 else: 

284 raise KeyError(f"Unknown component requested: {component}") 

285 

286 def readFull(self) -> Any: 

287 # Docstring inherited. 

288 return self.reader.read(**self.checked_parameters, dtype=self.storageClass_dtype) 

289 

290 def write_local_file(self, in_memory_dataset: Any, uri: ResourcePath) -> None: 

291 """Serialize the image to FITS. 

292 

293 Parameters 

294 ---------- 

295 in_memory_dataset : `object` 

296 Image to write. Must support a ``writeFits`` or 

297 ``writeFitsWithOptions`` interface. 

298 uri : `lsst.resources.ResourcePath` 

299 Location to write the local file. 

300 """ 

301 # check to see if we have a recipe requested 

302 recipeName = self.write_parameters.get("recipe") 

303 recipe = self.get_image_compression_settings(recipeName) 

304 if recipe: 

305 in_memory_dataset.writeFitsWithOptions(uri.ospath, options=recipe) 

306 else: 

307 in_memory_dataset.writeFits(uri.ospath) 

308 

309 def get_image_compression_settings(self, recipeName: str | None) -> dict: 

310 """Retrieve the relevant compression settings for this recipe. 

311 

312 Parameters 

313 ---------- 

314 recipeName : `str` or `None` 

315 Label associated with the collection of compression parameters 

316 to select. 

317 

318 Returns 

319 ------- 

320 settings : `dict` 

321 The selected settings. 

322 """ 

323 # if no recipe has been provided and there is no default 

324 # return immediately 

325 if not recipeName: 

326 if "default" not in self.write_recipes: 326 ↛ 327line 326 didn't jump to line 327 because the condition on line 326 was never true

327 return {} 

328 recipeName = "default" 

329 

330 if recipeName not in self.write_recipes: 330 ↛ 331line 330 didn't jump to line 331 because the condition on line 330 was never true

331 raise RuntimeError(f"Unrecognized recipe option given for compression: {recipeName}") 

332 

333 recipe = self.write_recipes[recipeName] 

334 if recipe is None: 

335 return {} 

336 seed: int | None = None 

337 for plane in ("image", "mask", "variance"): 

338 if plane in recipe and (quantization := recipe[plane].get("quantization")) is not None: 

339 if quantization.get("seed", 0) == 0: 339 ↛ 364line 339 didn't jump to line 364 because the condition on line 339 was always true

340 if seed is None: 

341 # Set the seed based on data ID. We can't just use 

342 # 'hash', since like 'set' that's not deterministic. 

343 # And we can't rely on a DimensionPacker because those 

344 # are only defined for certain combinations of 

345 # dimensions. Doing an MD5 of the JSON feels like 

346 # overkill but I don't really see anything much 

347 # simpler. 

348 hash_bytes = hashlib.md5( 

349 json.dumps(list(self.data_id.required_values)).encode(), 

350 usedforsecurity=False, 

351 ).digest() 

352 # And it *really* feels like overkill when we squash 

353 # that into the [1, 10000] range allowed by FITS. 

354 seed = 1 + int.from_bytes(hash_bytes) % 9999 

355 _LOG.debug( 

356 "Setting compression quantization seed for %s %s %s to %s.", 

357 self.data_id, 

358 self.dataset_ref.datasetType.name, 

359 plane, 

360 seed, 

361 ) 

362 quantization["seed"] = seed 

363 else: 

364 _LOG.warning( 

365 "Compression quantization seed for %s %s %s was set explicitly to %s.", 

366 self.dataset_ref.datasetType.name, 

367 self.data_id, 

368 plane, 

369 quantization["seed"], 

370 ) 

371 else: 

372 _LOG.debug( 

373 "No quantization found for %s %s %s.", 

374 self.dataset_ref.datasetType.name, 

375 self.data_id, 

376 plane, 

377 ) 

378 return recipe 

379 

380 @classmethod 

381 def validate_write_recipes(cls, recipes: Mapping[str, Any] | None) -> Mapping[str, Any] | None: 

382 """Validate supplied recipes for this formatter. 

383 

384 The recipes are supplemented with default values where appropriate. 

385 

386 Parameters 

387 ---------- 

388 recipes : `dict` or `None` 

389 Recipes to validate. Can be empty dict or `None`. 

390 

391 Returns 

392 ------- 

393 validated : `dict` 

394 Validated recipes. Returns what was given if there are no 

395 recipes listed. 

396 

397 Raises 

398 ------ 

399 RuntimeError 

400 Raised if validation fails. 

401 """ 

402 if not recipes: 

403 # We can not insist on recipes being specified. 

404 return recipes 

405 

406 validated: dict[str, Any] = {} 

407 for name, recipe in recipes.items(): 

408 if recipe is not None: 

409 validated[name] = {} 

410 for plane in ["image", "mask", "variance"]: 

411 try: 

412 options = CompressionOptions.from_mapping(recipe[plane]) 

413 except Exception as err: 

414 err.add_note(f"Validating write recipe {name!r} ({plane!r} section).") 

415 raise 

416 validated[name][plane] = options.to_dict() 

417 else: 

418 validated[name] = None 

419 return validated 

420 

421 

422class FitsImageFormatter(StandardFitsImageFormatterBase): 

423 """Concrete formatter for reading/writing `~lsst.afw.image.Image` 

424 from/to FITS. 

425 """ 

426 

427 ReaderClass = ImageFitsReader 

428 

429 

430class FitsMaskFormatter(StandardFitsImageFormatterBase): 

431 """Concrete formatter for reading/writing `~lsst.afw.image.Mask` 

432 from/to FITS. 

433 """ 

434 

435 ReaderClass = MaskFitsReader 

436 

437 

438class FitsMaskedImageFormatter(StandardFitsImageFormatterBase): 

439 """Concrete formatter for reading/writing `~lsst.afw.image.MaskedImage` 

440 from/to FITS. 

441 """ 

442 

443 ReaderClass = MaskedImageFitsReader 

444 

445 def readComponent(self, component: str) -> Any: 

446 # Docstring inherited. 

447 if component == "image": 

448 return self.reader.readImage(**self.checked_parameters, dtype=self.storageClass_dtype) 

449 elif component == "mask": 

450 return self.reader.readMask(**self.checked_parameters) 

451 elif component == "variance": 

452 return self.reader.readVariance(**self.checked_parameters, dtype=self.storageClass_dtype) 

453 else: 

454 # Delegate to base for bbox, dimensions, xy0. 

455 return super().readComponent(component) 

456 

457 

458def standardizeAmplifierParameters( 

459 parameters: dict[str, Any], on_disk_detector: lsst.afw.cameraGeom.Detector | None 

460) -> tuple[lsst.afw.cameraGeom.Amplifier, lsst.afw.cameraGeom.Detector, bool]: 

461 """Preprocess the Exposure storage class's "amp" and "detector" parameters. 

462 

463 This checks the given objects for consistency with the on-disk geometry and 

464 converts amplifier IDs/names to Amplifier instances. 

465 

466 Parameters 

467 ---------- 

468 parameters : `dict` 

469 Dictionary of parameters passed to formatter. See the Exposure storage 

470 class definition in daf_butler for allowed keys and values. 

471 on_disk_detector : `lsst.afw.cameraGeom.Detector` or `None` 

472 Detector that represents the on-disk image being loaded, or `None` if 

473 this is unknown (and hence the user must provide one in 

474 ``parameters`` if "amp" is in ``parameters``). 

475 

476 Returns 

477 ------- 

478 amplifier : `lsst.afw.cameraGeom.Amplifier` or `None` 

479 An amplifier object that defines a subimage to load, or `None` if there 

480 was no "amp" parameter. 

481 detector : `lsst.afw.cameraGeom.Detector` or `None` 

482 A detector object whose amplifiers are in the same s/orientation 

483 state as the on-disk image. If there is no "amp" parameter, 

484 ``on_disk_detector`` is simply passed through. 

485 regions_differ : `bool` 

486 `True` if the on-disk detector and the detector given in the parameters 

487 had different bounding boxes for one or more regions. This can happen 

488 if the true overscan region sizes can only be determined when the image 

489 is actually read, but otherwise it should be considered user error. 

490 """ 

491 if (amplifier := parameters.get("amp")) is None: 

492 return None, on_disk_detector, False 

493 if "bbox" in parameters or "origin" in parameters: 493 ↛ 494line 493 didn't jump to line 494 because the condition on line 493 was never true

494 raise ValueError("Cannot pass 'amp' with 'bbox' or 'origin'.") 

495 if isinstance(amplifier, int | str): 

496 amp_key = amplifier 

497 target_amplifier = None 

498 else: 

499 amp_key = amplifier.getName() 

500 target_amplifier = amplifier 

501 if (detector := parameters.get("detector")) is not None: 

502 if on_disk_detector is not None: 502 ↛ 523line 502 didn't jump to line 523 because the condition on line 502 was always true

503 # User passed a detector and we also found one on disk. Check them 

504 # for consistency. Note that we are checking the amps we'd get 

505 # from the two detectors against each other, not the amplifier we 

506 # got directly from the user, as the latter is allowed to differ in 

507 # assembly/orientation state. 

508 comparison = on_disk_detector[amp_key].compareGeometry(detector[amp_key]) 

509 if comparison & comparison.ASSEMBLY_DIFFERS: 509 ↛ 510line 509 didn't jump to line 510 because the condition on line 509 was never true

510 raise ValueError( 

511 "The given 'detector' has a different assembly state and/or orientation from " 

512 f"the on-disk one for amp {amp_key}." 

513 ) 

514 else: 

515 if on_disk_detector is None: 515 ↛ 516line 515 didn't jump to line 516 because the condition on line 515 was never true

516 raise ValueError( 

517 f"No on-disk detector and no detector given; cannot load amplifier from key {amp_key}. " 

518 "Please provide either a 'detector' parameter or an Amplifier instance in the " 

519 "'amp' parameter." 

520 ) 

521 comparison = AmplifierGeometryComparison.EQUAL 

522 detector = on_disk_detector 

523 if target_amplifier is None: 

524 target_amplifier = detector[amp_key] 

525 return target_amplifier, detector, comparison & comparison.REGIONS_DIFFER 

526 

527 

528class _ComponentCache(NamedTuple): 

529 id_: uuid.UUID | None = None 

530 reader: ExposureFitsReader | None = None 

531 bbox: lsst.geom.Box2I | None = None 

532 mem: MemFileManager | None = None 

533 

534 

535class FitsExposureFormatter(FitsMaskedImageFormatter): 

536 """Concrete formatter for reading/writing `~lsst.afw.image.Exposure` 

537 from/to FITS. 

538 

539 Notes 

540 ----- 

541 This class inherits from `FitsMaskedImageFormatter` even though 

542 `lsst.afw.image.Exposure` doesn't inherit from 

543 `lsst.afw.image.MaskedImage`; this is just an easy way to be able to 

544 delegate to `FitsMaskedImageFormatter` for component-handling, and 

545 should be replaced with e.g. both calling a free function if that slight 

546 type covariance violation ever becomes a practical problem. 

547 """ 

548 

549 can_read_from_uri = True 

550 ReaderClass = ExposureFitsReader 

551 # TODO: Remove MemFileManager from cache when DM-49640 is fixed. 

552 _lock = threading.Lock() 

553 _cached_fits: _ComponentCache = _ComponentCache() 

554 

555 def read_from_uri(self, uri: ResourcePath, component: str | None = None, expected_size: int = -1) -> Any: 

556 # For now only support small non-pixel components. In future 

557 # could work with cutouts. 

558 self._reader = None # Guarantee things are reset. 

559 

560 parameters = self.checked_parameters.copy() 

561 preserve_quantization = parameters.pop("preserve_quantization", False) 

562 # Full read, always use local file read. 

563 if preserve_quantization or (not component and not parameters): 

564 return NotImplemented 

565 

566 if not _ALWAYS_USE_ASTROPY_FOR_COMPONENT_READ and uri.isLocal: 

567 # For a local URI allow afw to read it directly. 

568 return NotImplemented 

569 pixel_components = ("mask", "image", "variance") 

570 

571 if component in pixel_components: 

572 # For pixel access currently this can not be cached in memory 

573 # and the performance gains are unclear. Assume local file 

574 # read with file caching for now. 

575 return NotImplemented 

576 

577 # With current file layouts the non-pixel extensions account for 1/3 

578 # of the file size and it is more efficient to download the entire 

579 # file. 

580 if not ( 580 ↛ 584line 580 didn't jump to line 584 because the condition on line 580 was never true

581 _ALWAYS_USE_ASTROPY_FOR_COMPONENT_READ 

582 or self._dataset_ref.dataId.mapping.keys().isdisjoint({"tract", "patch"}) 

583 ): 

584 return NotImplemented 

585 

586 # Cutouts can be optimized. For now only use this optimization 

587 # if bbox is the only parameter and the number of pixels in the 

588 # bounding box is reasonable. 

589 bbox = None 

590 origin = lsst.afw.image.PARENT 

591 if not component: 

592 # Try to support PARENT and LOCAL origin but if there are any 

593 # other parameters do not attempt a cutout. 

594 if parameters.keys() - {"bbox", "origin"}: 594 ↛ 595line 594 didn't jump to line 595 because the condition on line 594 was never true

595 return NotImplemented 

596 bbox = parameters["bbox"] 

597 origin = parameters.get("origin", lsst.afw.image.PARENT) 

598 # For larger cutouts use the full file. 

599 max_cutout_size = 500 * 500 

600 if not _ALWAYS_USE_ASTROPY_FOR_COMPONENT_READ and bbox.width * bbox.height > max_cutout_size: 600 ↛ 601line 600 didn't jump to line 601 because the condition on line 600 was never true

601 return NotImplemented 

602 

603 # We only cache component reads since those are small. 

604 if component: 

605 with self._lock: 

606 cache = type(self)._cached_fits 

607 if self.dataset_ref.id == cache.id_: 

608 if component in {"xy0", "dimensions", "bbox"} and cache.bbox is not None: 

609 match component: 

610 case "xy0": 

611 return cache.bbox.getMin() 

612 case "dimensions": 612 ↛ 614line 612 didn't jump to line 614 because the pattern on line 612 always matched

613 return cache.bbox.getDimensions() 

614 case "bbox": 

615 return cache.bbox 

616 else: 

617 self._reader = cache.reader 

618 return self.readComponent(component) 

619 

620 try: 

621 fs, fspath = uri.to_fsspec() 

622 except Exception: 

623 # fsspec cannot be initialized, fall back to downloading the file. 

624 return NotImplemented 

625 

626 bbox_component = None 

627 try: 

628 hdul = [] 

629 with fs.open(fspath) as f, astropy.io.fits.open(f) as fits_obj: 

630 # Read all non-pixel components and cache. 

631 for hdu in fits_obj: 

632 hdr = hdu.header 

633 extname = hdr.get("EXTNAME") 

634 # Older files have IMAGE in EXTNAME in PRIMARY so check 

635 # for EXTEND=T. 

636 extend = hdr.get("EXTEND") 

637 if not extend and extname and extname.lower() in pixel_components: 

638 # Calculate the dimensional components for later 

639 # caching. Do not derive from cached FITS reader 

640 # because they depend on the dimensionality of the 

641 # pixel data and we do not want to cache the pixel 

642 # data. 

643 if bbox_component is None: 

644 shape = hdu.shape 

645 dimensions = lsst.geom.Extent2I(shape[1], shape[0]) 

646 

647 # XY0 is defined in the A WCS. 

648 pl = PropertyList() 

649 pl.update(hdr) 

650 xy0 = getImageXY0FromMetadata(pl, "A", strip=False) 

651 

652 # This is the PARENT bbox. 

653 bbox_component = lsst.geom.Box2I(xy0, dimensions) 

654 

655 # Handle cutout request. 

656 if bbox: 

657 if origin == lsst.afw.image.PARENT: 

658 full_bbox = bbox_component 

659 else: 

660 full_bbox = lsst.geom.Box2I( 

661 lsst.geom.Point2I(0, 0), bbox_component.getDimensions 

662 ) 

663 minX = bbox.getBeginX() - full_bbox.getBeginX() 

664 maxX = bbox.getEndX() - full_bbox.getBeginX() 

665 minY = bbox.getBeginY() - full_bbox.getBeginY() 

666 maxY = bbox.getEndY() - full_bbox.getBeginY() 

667 data = hdu.section[minY:maxY, minX:maxX] 

668 

669 # Must correct the header WCS to take into 

670 # account the offset. 

671 if (k := "CRPIX1") in hdr: 671 ↛ 673line 671 didn't jump to line 673 because the condition on line 671 was always true

672 hdr[k] -= minX 

673 if (k := "CRPIX2") in hdr: 673 ↛ 675line 673 didn't jump to line 675 because the condition on line 673 was always true

674 hdr[k] -= minY 

675 if (k := "LTV1") in hdr: 675 ↛ 677line 675 didn't jump to line 677 because the condition on line 675 was always true

676 hdr[k] = -bbox.getBeginX() 

677 if (k := "LTV2") in hdr: 677 ↛ 679line 677 didn't jump to line 679 because the condition on line 677 was always true

678 hdr[k] = -bbox.getBeginY() 

679 if (k := "CRVAL1A") in hdr: 679 ↛ 681line 679 didn't jump to line 681 because the condition on line 679 was always true

680 hdr[k] = bbox.getBeginX() 

681 if (k := "CRVAL2A") in hdr: 681 ↛ 687line 681 didn't jump to line 687 because the condition on line 681 was always true

682 hdr[k] = bbox.getBeginY() 

683 else: 

684 data = np.zeros([1, 1], dtype=np.int32) 

685 

686 # Construct a new HDU and copy the header. 

687 stripped_hdu = astropy.io.fits.ImageHDU(data=data, header=hdr) 

688 hdul.append(stripped_hdu) 

689 else: 

690 hdul.append(hdu) 

691 stripped_fits = astropy.io.fits.HDUList(hdus=hdul) 

692 # Write the FITS file to in-memory FITS. 

693 buffer = BytesIO() 

694 stripped_fits.writeto(buffer) 

695 except Exception as e: 

696 # For some reason we can't open the remote file so fall back. 

697 _LOG.debug( 

698 "Attempted remote read of components but encountered an error. " 

699 "Falling back to file download. Error: %s", 

700 str(e), 

701 ) 

702 return NotImplemented 

703 

704 # Pass the new FITS buffer to the reader class without going through 

705 # a temporary file. We can assume this is relatively small for 

706 # components. 

707 fits_data = buffer.getvalue() 

708 mem = MemFileManager(len(fits_data)) 

709 mem.setData(fits_data, len(fits_data)) 

710 self._reader = self.ReaderClass(mem) 

711 

712 if component: 

713 with self._lock: 

714 type(self)._cached_fits = _ComponentCache( 

715 id_=self.dataset_ref.id, 

716 reader=self._reader, 

717 mem=mem, 

718 bbox=bbox_component, 

719 ) 

720 match component: 

721 case "xy0": 721 ↛ 722line 721 didn't jump to line 722 because the pattern on line 721 never matched

722 if bbox_component is None: # For mypy. 

723 return None 

724 return bbox_component.getMin() 

725 case "dimensions": 725 ↛ 726line 725 didn't jump to line 726 because the pattern on line 725 never matched

726 if bbox_component is None: 

727 return None 

728 return bbox_component.getDimensions() 

729 case "bbox": 729 ↛ 731line 729 didn't jump to line 731 because the pattern on line 729 always matched

730 return bbox_component 

731 case _: 

732 return self.readComponent(component) 

733 else: 

734 # Must be a cutout. We have applied the bbox parameter so no 

735 # parameters should be passed here. 

736 cutout = self.reader.read(dtype=self.storageClass_dtype) 

737 cutout.getInfo().setFilter(self._fixFilterLabels(cutout.getInfo().getFilter())) 

738 return cutout 

739 

740 def add_provenance( 

741 self, in_memory_dataset: Any, /, *, provenance: DatasetProvenance | None = None 

742 ) -> Any: 

743 # Add provenance via FITS headers. 

744 add_provenance_to_fits_header(in_memory_dataset.metadata, self.dataset_ref, provenance) 

745 return in_memory_dataset 

746 

747 def readComponent(self, component: str) -> Any: 

748 # Docstring inherited. 

749 # Generic components can be read via a string name; DM-27754 will make 

750 # this mapping larger at the expense of the following one. 

751 genericComponents = { 

752 "summaryStats": ExposureInfo.KEY_SUMMARY_STATS, 

753 } 

754 if (genericComponentName := genericComponents.get(component)) is not None: 

755 return self.reader.readComponent(genericComponentName) 

756 # Other components have hard-coded method names, but don't take 

757 # parameters. 

758 standardComponents = { 

759 "id": "readExposureId", 

760 "metadata": "readMetadata", 

761 "wcs": "readWcs", 

762 "coaddInputs": "readCoaddInputs", 

763 "psf": "readPsf", 

764 "photoCalib": "readPhotoCalib", 

765 "filter": "readFilter", 

766 "validPolygon": "readValidPolygon", 

767 "apCorrMap": "readApCorrMap", 

768 "visitInfo": "readVisitInfo", 

769 "transmissionCurve": "readTransmissionCurve", 

770 "detector": "readDetector", 

771 "exposureInfo": "readExposureInfo", 

772 } 

773 if (methodName := standardComponents.get(component)) is not None: 

774 result = getattr(self.reader, methodName)() 

775 if component == "filter": 

776 return self._fixFilterLabels(result) 

777 return result 

778 # Delegate to MaskedImage and ImageBase implementations for the rest. 

779 return super().readComponent(component) 

780 

781 def readFull(self) -> Any: 

782 # Docstring inherited. 

783 amplifier, detector, _ = standardizeAmplifierParameters( 

784 self.checked_parameters, 

785 self.reader.readDetector(), 

786 ) 

787 if amplifier is not None: 

788 amplifier_isolator = AmplifierIsolator( 

789 amplifier, 

790 self.reader.readBBox(), 

791 detector, 

792 ) 

793 result = amplifier_isolator.transform_subimage( 

794 self.reader.read(bbox=amplifier_isolator.subimage_bbox, dtype=self.storageClass_dtype) 

795 ) 

796 result.setDetector(amplifier_isolator.make_detector()) 

797 else: 

798 result = self.reader.read(**self.checked_parameters, dtype=self.storageClass_dtype) 

799 result.getInfo().setFilter(self._fixFilterLabels(result.getInfo().getFilter())) 

800 return result 

801 

802 def _fixFilterLabels( 

803 self, file_filter_label: lsst.afw.image.FilterLabel, should_be_standardized: bool | None = None 

804 ) -> lsst.afw.image.FilterLabel: 

805 """Compare the filter label read from the file with the one in the 

806 data ID. 

807 

808 Parameters 

809 ---------- 

810 file_filter_label : `lsst.afw.image.FilterLabel` or `None` 

811 Filter label read from the file, if there was one. 

812 should_be_standardized : `bool`, optional 

813 If `True`, expect ``file_filter_label`` to be consistent with the 

814 data ID and warn only if it is not. If `False`, expect it to be 

815 inconsistent and warn only if the data ID is incomplete and hence 

816 the `FilterLabel` cannot be fixed. If `None` (default) guess 

817 whether the file should be standardized by looking at the 

818 serialization version number in file, which requires this method to 

819 have been run after `readFull` or `readComponent`. 

820 

821 Returns 

822 ------- 

823 filter_label : `lsst.afw.image.FilterLabel` or `None` 

824 The preferred filter label; may be the given one or one built from 

825 the data ID. `None` is returned if there should never be any 

826 filters associated with this dataset type. 

827 

828 Notes 

829 ----- 

830 Most test coverage for this method is in ci_hsc_gen3, where we have 

831 much easier access to test data that exhibits the problems it attempts 

832 to solve. 

833 """ 

834 # Remember filter data ID keys that weren't in this particular data ID, 

835 # so we can warn about them later. 

836 missing = [] 

837 band = None 

838 physical_filter = None 

839 if "band" in self.data_id.dimensions.names: 

840 band = self.data_id.get("band") 

841 # band isn't in the data ID; is that just because this data ID 

842 # hasn't been filled in with everything the Registry knows, or 

843 # because this dataset is never associated with a band? 

844 if band is None and not self.data_id.hasFull() and "band" in self.data_id.dimensions.implied: 844 ↛ 845line 844 didn't jump to line 845 because the condition on line 844 was never true

845 missing.append("band") 

846 if "physical_filter" in self.data_id.dimensions.names: 

847 physical_filter = self.data_id.get("physical_filter") 

848 # Same check as above for band, but for physical_filter. 

849 if ( 849 ↛ 854line 849 didn't jump to line 854 because the condition on line 849 was never true

850 physical_filter is None 

851 and not self.data_id.hasFull() 

852 and "physical_filter" in self.data_id.dimensions.implied 

853 ): 

854 missing.append("physical_filter") 

855 if should_be_standardized is None: 855 ↛ 858line 855 didn't jump to line 858 because the condition on line 855 was always true

856 version = self.reader.readSerializationVersion() 

857 should_be_standardized = version >= 2 

858 if missing: 858 ↛ 864line 858 didn't jump to line 864 because the condition on line 858 was never true

859 # Data ID identifies a filter but the actual filter label values 

860 # haven't been fetched from the database; we have no choice but 

861 # to use the one in the file. 

862 # Warn if that's more likely than not to be bad, because the file 

863 # predates filter standardization. 

864 if not should_be_standardized: 

865 warnings.warn( 

866 f"Data ID {self.data_id} is missing (implied) value(s) for {missing}; " 

867 "the correctness of this Exposure's FilterLabel cannot be guaranteed. " 

868 "Call Registry.expandDataId before Butler.get to avoid this.", 

869 # Report the warning from outside of middleware or the 

870 # relevant runQuantum method. 

871 stacklevel=find_outside_stacklevel( 

872 "lsst.obs.base", "lsst.pipe.base", "lsst.daf.butler", allow_methods={"runQuantum"} 

873 ), 

874 ) 

875 return file_filter_label 

876 if band is None and physical_filter is None: 

877 data_id_filter_label = None 

878 else: 

879 data_id_filter_label = FilterLabel(band=band, physical=physical_filter) 

880 if data_id_filter_label != file_filter_label and should_be_standardized: 880 ↛ 885line 880 didn't jump to line 885 because the condition on line 880 was never true

881 # File was written after FilterLabel and standardization, but its 

882 # FilterLabel doesn't agree with the data ID: this indicates a bug 

883 # in whatever code produced the Exposure (though it may be one that 

884 # has been fixed since the file was written). 

885 warnings.warn(