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

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 future_type, read_kwargs = _get_future_image_type_and_kwargs( 

169 self.file_descriptor.readStorageClass.name, 

170 component=component, 

171 parameters=self.checked_parameters, 

172 ) 

173 if future_type is not None: 173 ↛ 174line 173 didn't jump to line 174 because the condition on line 173 was never true

174 return future_type.read_legacy(path, component=component, **read_kwargs) 

175 

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

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

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

179 # file. 

180 self._reader_path = path 

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

182 try: 

183 if component is not None: 

184 in_memory_dataset = self.readComponent(component) 

185 else: 

186 in_memory_dataset = self.readFull() 

187 finally: 

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

189 return in_memory_dataset 

190 

191 @abstractmethod 

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

193 """Read a component dataset. 

194 

195 Parameters 

196 ---------- 

197 component : `str`, optional 

198 Component to read from the file. 

199 

200 Returns 

201 ------- 

202 obj : `typing.Any` 

203 In-memory component object. 

204 

205 Raises 

206 ------ 

207 KeyError 

208 Raised if the requested component cannot be handled. 

209 """ 

210 raise NotImplementedError() 

211 

212 @abstractmethod 

213 def readFull(self) -> Any: 

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

215 

216 Returns 

217 ------- 

218 obj : `typing.Any` 

219 In-memory component object. 

220 

221 """ 

222 raise NotImplementedError() 

223 

224 

225class StandardFitsImageFormatterBase(FitsImageFormatterBase): 

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

227 written using LSST code. 

228 

229 Notes 

230 ----- 

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

232 

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

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

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

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

237 and catch `KeyError`). 

238 

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

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

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

242 should generally reimplement without delegating (the implementation is 

243 trivial). 

244 

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

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

247 

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

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

250 contain entries supported by 

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

252 compression). 

253 

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

255 specialization): 

256 

257 .. code-block:: yaml 

258 

259 lsst.obs.base.fitsExposureFormatter.FitsExposureFormatter: 

260 default: 

261 image: &default 

262 algorithm: GZIP_2 

263 mask: *default 

264 variance: *default 

265 

266 """ 

267 

268 supported_write_parameters = frozenset({"recipe"}) 

269 

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

271 # Docstring inherited. 

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

273 bbox = self.reader.readBBox() 

274 if component == "dimensions": 

275 return bbox.getDimensions() 

276 elif component == "xy0": 

277 return bbox.getMin() 

278 else: 

279 return bbox 

280 else: 

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

282 

283 def readFull(self) -> Any: 

284 # Docstring inherited. 

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

286 

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

288 """Serialize the image to FITS. 

289 

290 Parameters 

291 ---------- 

292 in_memory_dataset : `object` 

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

294 ``writeFitsWithOptions`` interface. 

295 uri : `lsst.resources.ResourcePath` 

296 Location to write the local file. 

297 """ 

298 # check to see if we have a recipe requested 

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

300 recipe = self.get_image_compression_settings(recipeName) 

301 if recipe: 

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

303 else: 

304 in_memory_dataset.writeFits(uri.ospath) 

305 

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

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

308 

309 Parameters 

310 ---------- 

311 recipeName : `str` or `None` 

312 Label associated with the collection of compression parameters 

313 to select. 

314 

315 Returns 

316 ------- 

317 settings : `dict` 

318 The selected settings. 

319 """ 

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

321 # return immediately 

322 if not recipeName: 

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

324 return {} 

325 recipeName = "default" 

326 

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

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

329 

330 recipe = self.write_recipes[recipeName] 

331 if recipe is None: 

332 return {} 

333 seed: int | None = None 

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

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

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

337 if seed is None: 

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

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

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

341 # are only defined for certain combinations of 

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

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

344 # simpler. 

345 hash_bytes = hashlib.md5( 

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

347 usedforsecurity=False, 

348 ).digest() 

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

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

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

352 _LOG.debug( 

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

354 self.data_id, 

355 self.dataset_ref.datasetType.name, 

356 plane, 

357 seed, 

358 ) 

359 quantization["seed"] = seed 

360 else: 

361 _LOG.warning( 

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

363 self.dataset_ref.datasetType.name, 

364 self.data_id, 

365 plane, 

366 quantization["seed"], 

367 ) 

368 else: 

369 _LOG.debug( 

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

371 self.dataset_ref.datasetType.name, 

372 self.data_id, 

373 plane, 

374 ) 

375 return recipe 

376 

377 @classmethod 

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

379 """Validate supplied recipes for this formatter. 

380 

381 The recipes are supplemented with default values where appropriate. 

382 

383 Parameters 

384 ---------- 

385 recipes : `dict` or `None` 

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

387 

388 Returns 

389 ------- 

390 validated : `dict` 

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

392 recipes listed. 

393 

394 Raises 

395 ------ 

396 RuntimeError 

397 Raised if validation fails. 

398 """ 

399 if not recipes: 

400 # We can not insist on recipes being specified. 

401 return recipes 

402 

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

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

405 if recipe is not None: 

406 validated[name] = {} 

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

408 try: 

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

410 except Exception as err: 

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

412 raise 

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

414 else: 

415 validated[name] = None 

416 return validated 

417 

418 

419class FitsImageFormatter(StandardFitsImageFormatterBase): 

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

421 from/to FITS. 

422 """ 

423 

424 ReaderClass = ImageFitsReader 

425 

426 

427class FitsMaskFormatter(StandardFitsImageFormatterBase): 

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

429 from/to FITS. 

430 """ 

431 

432 ReaderClass = MaskFitsReader 

433 

434 

435class FitsMaskedImageFormatter(StandardFitsImageFormatterBase): 

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

437 from/to FITS. 

438 """ 

439 

440 ReaderClass = MaskedImageFitsReader 

441 

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

443 # Docstring inherited. 

444 if component == "image": 

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

446 elif component == "mask": 

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

448 elif component == "variance": 

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

450 else: 

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

452 return super().readComponent(component) 

453 

454 

455def standardizeAmplifierParameters( 

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

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

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

459 

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

461 converts amplifier IDs/names to Amplifier instances. 

462 

463 Parameters 

464 ---------- 

465 parameters : `dict` 

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

467 class definition in daf_butler for allowed keys and values. 

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

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

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

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

472 

473 Returns 

474 ------- 

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

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

477 was no "amp" parameter. 

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

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

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

481 ``on_disk_detector`` is simply passed through. 

482 regions_differ : `bool` 

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

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

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

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

487 """ 

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

489 return None, on_disk_detector, False 

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

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

492 if isinstance(amplifier, int | str): 

493 amp_key = amplifier 

494 target_amplifier = None 

495 else: 

496 amp_key = amplifier.getName() 

497 target_amplifier = amplifier 

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

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

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

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

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

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

504 # assembly/orientation state. 

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

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

507 raise ValueError( 

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

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

510 ) 

511 else: 

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

513 raise ValueError( 

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

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

516 "'amp' parameter." 

517 ) 

518 comparison = AmplifierGeometryComparison.EQUAL 

519 detector = on_disk_detector 

520 if target_amplifier is None: 

521 target_amplifier = detector[amp_key] 

522 return target_amplifier, detector, comparison & comparison.REGIONS_DIFFER 

523 

524 

525class _ComponentCache(NamedTuple): 

526 id_: uuid.UUID | None = None 

527 reader: ExposureFitsReader | None = None 

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

529 mem: MemFileManager | None = None 

530 

531 

532class FitsExposureFormatter(FitsMaskedImageFormatter): 

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

534 from/to FITS. 

535 

536 Notes 

537 ----- 

538 This class inherits from `FitsMaskedImageFormatter` even though 

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

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

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

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

543 type covariance violation ever becomes a practical problem. 

544 """ 

545 

546 can_read_from_uri = True 

547 ReaderClass = ExposureFitsReader 

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

549 _lock = threading.Lock() 

550 _cached_fits: _ComponentCache = _ComponentCache() 

551 

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

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

554 # could work with cutouts. 

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

556 

557 parameters = self.checked_parameters.copy() 

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

559 # Full read, always use local file read. 

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

561 return NotImplemented 

562 

563 if not _ALWAYS_USE_ASTROPY_FOR_COMPONENT_READ and uri.isLocal: 

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

565 return NotImplemented 

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

567 

568 if component in pixel_components: 

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

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

571 # read with file caching for now. 

572 return NotImplemented 

573 

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

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

576 # file. 

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

578 _ALWAYS_USE_ASTROPY_FOR_COMPONENT_READ 

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

580 ): 

581 return NotImplemented 

582 

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

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

585 # bounding box is reasonable. 

586 bbox = None 

587 origin = lsst.afw.image.PARENT 

588 if not component: 

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

590 # other parameters do not attempt a cutout. 

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

592 return NotImplemented 

593 bbox = parameters["bbox"] 

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

595 # For larger cutouts use the full file. 

596 max_cutout_size = 500 * 500 

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

598 return NotImplemented 

599 

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

601 if component: 

602 with self._lock: 

603 cache = type(self)._cached_fits 

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

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

606 match component: 

607 case "xy0": 

608 return cache.bbox.getMin() 

609 case "dimensions": 

610 return cache.bbox.getDimensions() 

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

612 return cache.bbox 

613 else: 

614 self._reader = cache.reader 

615 return self.readComponent(component) 

616 

617 try: 

618 fs, fspath = uri.to_fsspec() 

619 except Exception: 

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

621 return NotImplemented 

622 

623 bbox_component = None 

624 try: 

625 hdul = [] 

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

627 # Read all non-pixel components and cache. 

628 for hdu in fits_obj: 

629 hdr = hdu.header 

630 extname = hdr.get("EXTNAME") 

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

632 # for EXTEND=T. 

633 extend = hdr.get("EXTEND") 

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

635 # Calculate the dimensional components for later 

636 # caching. Do not derive from cached FITS reader 

637 # because they depend on the dimensionality of the 

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

639 # data. 

640 if bbox_component is None: 

641 shape = hdu.shape 

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

643 

644 # XY0 is defined in the A WCS. 

645 pl = PropertyList() 

646 pl.update(hdr) 

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

648 

649 # This is the PARENT bbox. 

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

651 

652 # Handle cutout request. 

653 if bbox: 

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

655 full_bbox = bbox_component 

656 else: 

657 full_bbox = lsst.geom.Box2I( 

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

659 ) 

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

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

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

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

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

665 

666 # Must correct the header WCS to take into 

667 # account the offset. 

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

669 hdr[k] -= minX 

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

671 hdr[k] -= minY 

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

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

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

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

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

677 hdr[k] = bbox.getBeginX() 

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

679 hdr[k] = bbox.getBeginY() 

680 else: 

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

682 

683 # Construct a new HDU and copy the header. 

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

685 hdul.append(stripped_hdu) 

686 else: 

687 hdul.append(hdu) 

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

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

690 buffer = BytesIO() 

691 stripped_fits.writeto(buffer) 

692 except Exception as e: 

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

694 _LOG.debug( 

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

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

697 str(e), 

698 ) 

699 return NotImplemented 

700 

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

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

703 # components. 

704 fits_data = buffer.getvalue() 

705 mem = MemFileManager(len(fits_data)) 

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

707 self._reader = self.ReaderClass(mem) 

708 

709 if component: 

710 with self._lock: 

711 type(self)._cached_fits = _ComponentCache( 

712 id_=self.dataset_ref.id, 

713 reader=self._reader, 

714 mem=mem, 

715 bbox=bbox_component, 

716 ) 

717 match component: 

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

719 if bbox_component is None: # For mypy. 

720 return None 

721 return bbox_component.getMin() 

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

723 if bbox_component is None: 

724 return None 

725 return bbox_component.getDimensions() 

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

727 return bbox_component 

728 case _: 

729 return self.readComponent(component) 

730 else: 

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

732 # parameters should be passed here. 

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

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

735 return cutout 

736 

737 def add_provenance( 

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

739 ) -> Any: 

740 # Add provenance via FITS headers. 

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

742 return in_memory_dataset 

743 

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

745 # Docstring inherited. 

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

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

748 genericComponents = { 

749 "summaryStats": ExposureInfo.KEY_SUMMARY_STATS, 

750 } 

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

752 return self.reader.readComponent(genericComponentName) 

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

754 # parameters. 

755 standardComponents = { 

756 "id": "readExposureId", 

757 "metadata": "readMetadata", 

758 "wcs": "readWcs", 

759 "coaddInputs": "readCoaddInputs", 

760 "psf": "readPsf", 

761 "photoCalib": "readPhotoCalib", 

762 "filter": "readFilter", 

763 "validPolygon": "readValidPolygon", 

764 "apCorrMap": "readApCorrMap", 

765 "visitInfo": "readVisitInfo", 

766 "transmissionCurve": "readTransmissionCurve", 

767 "detector": "readDetector", 

768 "exposureInfo": "readExposureInfo", 

769 } 

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

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

772 if component == "filter": 

773 return self._fixFilterLabels(result) 

774 return result 

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

776 return super().readComponent(component) 

777 

778 def readFull(self) -> Any: 

779 # Docstring inherited. 

780 amplifier, detector, _ = standardizeAmplifierParameters( 

781 self.checked_parameters, 

782 self.reader.readDetector(), 

783 ) 

784 if amplifier is not None: 

785 amplifier_isolator = AmplifierIsolator( 

786 amplifier, 

787 self.reader.readBBox(), 

788 detector, 

789 ) 

790 result = amplifier_isolator.transform_subimage( 

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

792 ) 

793 result.setDetector(amplifier_isolator.make_detector()) 

794 else: 

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

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

797 return result 

798 

799 def _fixFilterLabels( 

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

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

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

803 data ID. 

804 

805 Parameters 

806 ---------- 

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

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

809 should_be_standardized : `bool`, optional 

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

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

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

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

814 whether the file should be standardized by looking at the 

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

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

817 

818 Returns 

819 ------- 

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

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

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

823 filters associated with this dataset type. 

824 

825 Notes 

826 ----- 

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

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

829 to solve. 

830 """ 

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

832 # so we can warn about them later. 

833 missing = [] 

834 band = None 

835 physical_filter = None 

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

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

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

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

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

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

842 missing.append("band") 

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

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

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

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

847 physical_filter is None 

848 and not self.data_id.hasFull() 

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

850 ): 

851 missing.append("physical_filter") 

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

853 version = self.reader.readSerializationVersion() 

854 should_be_standardized = version >= 2 

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

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

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

858 # to use the one in the file. 

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

860 # predates filter standardization. 

861 if not should_be_standardized: 

862 warnings.warn( 

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

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

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

866 # Report the warning from outside of middleware or the 

867 # relevant runQuantum method. 

868 stacklevel=find_outside_stacklevel( 

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

870 ), 

871 ) 

872 return file_filter_label 

873 if band is None and physical_filter is None: 

874 data_id_filter_label = None 

875 else: 

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

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

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

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

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

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

882 warnings.warn( 

883 f"Reading {self.file_descriptor.location} with data ID {self.data_id}: " 

884 f"filter label mismatch (file is {file_filter_label}, data ID is "