Coverage for python/lsst/images/_transforms/_sky_projection.py: 58%

1# This file is part of lsst-images. 

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# Use of this source code is governed by a 3-clause BSD-style 

10# license that can be found in the LICENSE file. 

11 

12from __future__ import annotations 

13 

14__all__ = ("SkyProjection", "SkyProjectionAstropyView", "SkyProjectionSerializationModel") 

15 

16import functools 

17import itertools 

18import statistics 

19from collections.abc import Collection 

20from typing import TYPE_CHECKING, Any, ClassVar, Self, TypeVar, assert_type, cast, final, overload 

21 

22import astropy.units as u 

23import astropy.wcs 

24import numpy as np 

25import numpy.typing as npt 

26import pydantic 

27from astropy.coordinates import ICRS, Latitude, Longitude, SkyCoord 

28from astropy.wcs.wcsapi import BaseLowLevelWCS, HighLevelWCSMixin 

29 

30from .._geom import XY, YX, Bounds, Box 

31from ..describe import DescribableMixin, DescribeOptions, FieldRole, Report, ReportField, ReportTable 

32from ..serialization import ArchiveTree, InputArchive, InvalidParameterError, OutputArchive 

33from ..utils import is_none 

34from . import _ast as astshim 

35from ._frames import Frame, SkyFrame 

36from ._transform import Transform, TransformSerializationModel, _ast_apply 

37 

38if TYPE_CHECKING: 

39 try: 

40 from lsst.afw.geom import SkyWcs as LegacySkyWcs 

41 except ImportError: 

42 type LegacySkyWcs = Any # type: ignore[no-redef] 

43 

44 

45# This pre-python-3.12 declaration is needed by Sphinx (probably the 

46# autodoc-typehints plugin. 

47F = TypeVar("F", bound=Frame) 

48P = TypeVar("P", bound=pydantic.BaseModel) 

49 

50 

51def _ast_skyframe_axis_labels(frame: astshim.SkyFrame) -> tuple[str, str]: 

52 """Return the ``(longitude, latitude)`` axis labels of an AST SkyFrame. 

53 

54 The AST ``LonAxis``/``LatAxis`` attributes say which frame axis carries 

55 longitude and which carries latitude, so the labels come back ordered 

56 longitude-first regardless of the underlying axis order. For an ICRS 

57 frame these are ``("Right ascension", "Declination")``; for a Galactic 

58 frame ``("Galactic longitude", "Galactic latitude")``. 

59 """ 

60 return frame.getLabel(frame.lonAxis), frame.getLabel(frame.latAxis) 

61 

62 

63def _sky_parts(sky: SkyCoord) -> tuple[str, float, str, float]: 

64 """Return sexagesimal and decimal-degree pieces of a scalar sky position. 

65 

66 Returns ``(ra_sexagesimal, ra_deg, dec_sexagesimal, dec_deg)``. Right 

67 ascension is rendered in sexagesimal hours (``hms``) to a tenth of a 

68 second and declination in sexagesimal degrees (``dms``) to the nearest 

69 arcsecond, so the units are unambiguous. 

70 """ 

71 return ( 

72 sky.ra.to_string(unit=u.hour, sep="hms", pad=True, precision=1), 

73 sky.ra.deg, 

74 sky.dec.to_string(sep="dms", pad=True, alwayssign=True, precision=0), 

75 sky.dec.deg, 

76 ) 

77 

78 

79def _format_pixel_sky(sky: SkyCoord) -> str: 

80 """Return the sexagesimal plus labeled decimal-degree form of a sky 

81 position, e.g. ``"00h00m29.3s +02d07m51s (RA 0.122076°, Dec 2.130862°)"``. 

82 """ 

83 ra_sex, ra_deg, dec_sex, dec_deg = _sky_parts(sky) 

84 return f"{ra_sex} {dec_sex} (RA {ra_deg:.6f}°, Dec {dec_deg:+.6f}°)" 

85 

86 

87def _angular_unit(size: u.Quantity) -> u.UnitBase: 

88 """Return the unit that states an angular size at a readable scale. 

89 

90 Parameters 

91 ---------- 

92 size 

93 Angular size to choose a unit for. 

94 """ 

95 if size < 1.0 * u.arcmin: 

96 return u.arcsec 

97 if size < 1.0 * u.deg: 

98 return u.arcmin 

99 return u.deg 

100 

101 

102def _format_extent(width: u.Quantity, height: u.Quantity, position_angle: u.Quantity) -> str: 

103 """Return the angular size and orientation of a box on the sky. 

104 

105 Parameters 

106 ---------- 

107 width, height 

108 Great-circle extent along the pixel x and y axes. 

109 position_angle 

110 Orientation of the pixel y axis, measured from North toward East. 

111 

112 Returns 

113 ------- 

114 extent : `str` 

115 For example ``"55 x 54 arcsec @ 40 deg E of N"``. 

116 

117 Notes 

118 ----- 

119 Each extent gets the unit that suits it, named once when the two agree. 

120 A long, thin box has no unit that reads well for both of its sides. 

121 """ 

122 width_unit = _angular_unit(width) 

123 height_unit = _angular_unit(height) 

124 if width_unit == height_unit: 

125 size = f"{width.to_value(width_unit):.3g} x {height.to_value(height_unit):.3g} {width_unit!s}" 

126 else: 

127 size = ( 

128 f"{width.to_value(width_unit):.3g} {width_unit!s}" 

129 f" x {height.to_value(height_unit):.3g} {height_unit!s}" 

130 ) 

131 # Whole degrees are as fine as an orientation needs to be read to, and 

132 # rounding takes an angle just short of a full turn back to zero. 

133 angle = round(position_angle.to_value(u.deg)) % 360 

134 return f"{size} @ {angle} deg E of N" 

135 

136 

137@final 

138class SkyProjection[F: Frame](DescribableMixin): 

139 """A transform from pixel coordinates to sky coordinates. 

140 

141 Parameters 

142 ---------- 

143 pixel_to_sky 

144 A low-level transform that maps pixel coordinates to sky coordinates. 

145 fits_approximation 

146 An approximation to ``pixel_to_sky`` that is guaranteed to have a 

147 `~Transform.as_fits_wcs` method that does not return `None`. This 

148 should not be provided if ``pixel_to_sky`` is itself representable 

149 as a FITS WCS. 

150 

151 Notes 

152 ----- 

153 `Transform` is conceptually immutable (the internal AST Mapping should 

154 never be modified in-place after construction), and hence does not need to 

155 be copied when any object that holds it is copied. 

156 """ 

157 

158 def __init__( 

159 self, pixel_to_sky: Transform[F, SkyFrame], fits_approximation: Transform[F, SkyFrame] | None = None 

160 ) -> None: 

161 self._pixel_to_sky = pixel_to_sky 

162 if pixel_to_sky.in_frame.unit != u.pix: 

163 raise ValueError("Transform is not a mapping from pixel coordinates.") 

164 if pixel_to_sky.out_frame != SkyFrame.ICRS: 

165 raise ValueError("Transform is not a mapping to ICRS.") 

166 self._fits_approximation = fits_approximation 

167 

168 def __eq__(self, other: Any) -> bool: 

169 if self is other: 

170 return True 

171 if not isinstance(other, SkyProjection): 

172 return NotImplemented 

173 # Even though two approximations could be different and yet consistent 

174 # with the primary mapping (for example using different tolerances 

175 # on construction) we require them to be equal to declare that the 

176 # two objects are equal. 

177 if self._fits_approximation != other._fits_approximation: 

178 return False 

179 return self._pixel_to_sky == other._pixel_to_sky 

180 

181 @staticmethod 

182 def from_fits_wcs( 

183 fits_wcs: astropy.wcs.WCS, 

184 pixel_frame: F, 

185 pixel_bounds: Bounds | None = None, 

186 x0: int = 0, 

187 y0: int = 0, 

188 ) -> SkyProjection[F]: 

189 """Construct a transform from a FITS WCS. 

190 

191 Parameters 

192 ---------- 

193 fits_wcs 

194 FITS WCS to convert. 

195 pixel_frame 

196 Coordinate frame for the pixel grid. 

197 pixel_bounds 

198 The region that bounds valid pixels for this transform. 

199 x0 

200 Logical coordinate of the first column in the array this WCS 

201 relates to world coordinates. 

202 y0 

203 Logical coordinate of the first column in the array this WCS 

204 relates to world coordinates. 

205 

206 Notes 

207 ----- 

208 The ``x0`` and ``y0`` parameters reflect the fact that for FITS, the 

209 first row and column are always labeled ``(1, 1)``, while in Astropy 

210 and most other Python libraries, they are ``(0, 0)``. The `types` in 

211 this package (e.g. `Image`, `Mask`) allow them to be any pair of 

212 integers. 

213 

214 See Also 

215 -------- 

216 Transform.from_fits_wcs 

217 """ 

218 return SkyProjection( 

219 Transform.from_fits_wcs( 

220 fits_wcs, pixel_frame, SkyFrame.ICRS, in_bounds=pixel_bounds, x0=x0, y0=y0 

221 ) 

222 ) 

223 

224 @staticmethod 

225 def from_ast_frame_set( 

226 ast_frame_set: astshim.FrameSet, 

227 pixel_frame: F, 

228 pixel_bounds: Bounds | None = None, 

229 ) -> SkyProjection[F]: 

230 """Construct a sky projection from an AST FrameSet. 

231 

232 The current frame of the FrameSet must be an AST SkyFrame. Its 

233 coordinate system is forced to ICRS (AST adjusts the mapping 

234 automatically) so the resulting Projection is always in ICRS 

235 regardless of the original sky system. 

236 

237 Parameters 

238 ---------- 

239 ast_frame_set 

240 An AST FrameSet whose base frame is pixel coordinates and 

241 whose current frame is a SkyFrame (in any supported sky 

242 coordinate system). 

243 pixel_frame 

244 Coordinate frame for the pixel grid. 

245 pixel_bounds 

246 The region that bounds valid pixels for this transform. 

247 

248 Raises 

249 ------ 

250 ValueError 

251 If the current frame of the FrameSet is not a SkyFrame. 

252 """ 

253 current_frame = ast_frame_set.getFrame(ast_frame_set.current, copy=False) 

254 if not isinstance(current_frame, astshim.SkyFrame): 254 ↛ 255line 254 didn't jump to line 255 because the condition on line 254 was never true

255 raise ValueError( 

256 "The current frame of the AST FrameSet is not a SkyFrame " 

257 f"(got {type(current_frame).__name__})." 

258 ) 

259 current_frame.system = "ICRS" 

260 return SkyProjection(Transform(pixel_frame, SkyFrame.ICRS, ast_frame_set, in_bounds=pixel_bounds)) 

261 

262 @property 

263 def pixel_frame(self) -> F: 

264 """Coordinate frame for the pixel grid.""" 

265 return self._pixel_to_sky.in_frame 

266 

267 @property 

268 def sky_frame(self) -> SkyFrame: 

269 """Coordinate frame for the sky.""" 

270 return self._pixel_to_sky.out_frame 

271 

272 def _sky_axis_labels(self) -> tuple[str, str]: 

273 """Return the ``(longitude, latitude)`` sky-axis labels. 

274 

275 The labels are read from the AST SkyFrame when the underlying mapping 

276 exposes one, so a non-ICRS frame (e.g. Galactic) reports its own 

277 labels. A projection is always ICRS in practice, so this falls back 

278 to ``("Right ascension", "Declination")`` when the mapping is not a 

279 frame set. 

280 """ 

281 ast_mapping = self._pixel_to_sky._ast_mapping 

282 if isinstance(ast_mapping, astshim.FrameSet): 282 ↛ 286line 282 didn't jump to line 286 because the condition on line 282 was always true

283 frame = ast_mapping.getFrame(ast_mapping.current, copy=False) 

284 if isinstance(frame, astshim.SkyFrame): 284 ↛ 286line 284 didn't jump to line 286 because the condition on line 284 was always true

285 return _ast_skyframe_axis_labels(frame) 

286 return "Right ascension", "Declination" 

287 

288 @property 

289 def pixel_bounds(self) -> Bounds | None: 

290 """The region that bounds valid pixel points (`Bounds` | `None`).""" 

291 return self._pixel_to_sky.in_bounds 

292 

293 @property 

294 def pixel_to_sky_transform(self) -> Transform[F, SkyFrame]: 

295 """Low-level transform from pixel to sky coordinates (`Transform`).""" 

296 return self._pixel_to_sky 

297 

298 @property 

299 def sky_to_pixel_transform(self) -> Transform[SkyFrame, F]: 

300 """Low-level transform from sky to pixel coordinates (`Transform`).""" 

301 return self._pixel_to_sky.inverted() 

302 

303 @property 

304 def fits_approximation(self) -> SkyProjection[F] | None: 

305 """An approximation to this projection that is guaranteed to have an 

306 `as_fits_wcs` method that does not return `None`. 

307 """ 

308 return SkyProjection(self._fits_approximation) if self._fits_approximation is not None else None 

309 

310 def show(self, simplified: bool = False, comments: bool = False) -> str: 

311 """Return the AST native representation of the transform. 

312 

313 Parameters 

314 ---------- 

315 simplified 

316 Whether to ask AST to simplify the mapping before showing it. 

317 This will make it much more likely that two equivalent transforms 

318 have the same `show` result. 

319 comments 

320 Whether to include descriptive comments. 

321 """ 

322 return self._pixel_to_sky.show(simplified=simplified, comments=comments) 

323 

324 @overload 

325 def pixel_to_sky(self, point: XY[int | float] | YX[int | float], /) -> SkyCoord: ... 

326 

327 @overload 

328 def pixel_to_sky(self, point: XY[npt.ArrayLike] | YX[npt.ArrayLike], /) -> SkyCoord: ... 

329 

330 @overload 

331 def pixel_to_sky( 

332 self, /, *, x: int | float | npt.ArrayLike, y: int | float | npt.ArrayLike 

333 ) -> SkyCoord: ... 

334 

335 def pixel_to_sky( 

336 self, 

337 point: XY[Any] | YX[Any] | None = None, 

338 /, 

339 *, 

340 x: Any = None, 

341 y: Any = None, 

342 ) -> SkyCoord: 

343 """Transform one or more pixel points to sky coordinates. 

344 

345 Parameters 

346 ---------- 

347 point 

348 An `XY` or `YX` coordinate pair of pixel positions to transform. 

349 Mutually exclusive with ``x`` and ``y``. 

350 x : `float` | array-like 

351 ``x`` values of the pixel points to transform, as a scalar or 

352 any array-like. Results are broadcast against ``y``. 

353 Mutually exclusive with ``point``. 

354 y : `float` | array-like 

355 ``y`` values of the pixel points to transform, as a scalar or 

356 any array-like. Results are broadcast against ``x``. 

357 Mutually exclusive with ``point``. 

358 

359 Returns 

360 ------- 

361 astropy.coordinates.SkyCoord 

362 Transformed sky coordinates, with the broadcast shape of ``x`` 

363 and ``y``. 

364 """ 

365 match point: 

366 case None: 

367 if x is None or y is None: 

368 raise TypeError("Pass either a point or both x= and y= to 'pixel_to_sky'.") 

369 case XY() | YX(): 369 ↛ 373line 369 didn't jump to line 373 because the pattern on line 369 always matched

370 if x is not None or y is not None: 

371 raise TypeError("'pixel_to_sky' point argument is mutually exclusive with x= and y=.") 

372 x, y = point.x, point.y 

373 case _: 

374 raise TypeError(f"Unexpected positional argument type: {type(point)!r}.") 

375 sky_rad = self._pixel_to_sky.apply_forward(x=x, y=y) 

376 return SkyCoord(ra=sky_rad.x, dec=sky_rad.y, unit=u.rad) 

377 

378 def sky_to_pixel(self, sky: SkyCoord) -> XY[Any]: 

379 """Transform one or more sky coordinates to pixels. 

380 

381 Parameters 

382 ---------- 

383 sky 

384 Sky coordinates to transform. Any shape is supported; the 

385 result has the same shape as ``sky``. 

386 

387 Returns 

388 ------- 

389 `XY` [`numpy.ndarray` | `float`] 

390 Transformed pixel coordinates with the same shape as ``sky``. 

391 """ 

392 if sky.frame.name != "icrs": 392 ↛ 393line 392 didn't jump to line 393 because the condition on line 392 was never true

393 sky = sky.transform_to("icrs") 

394 ra: Longitude = sky.ra 

395 dec: Latitude = sky.dec 

396 # cast works around a mypy false positive specific to generic 

397 # NamedTuple classes: returning XY[Any] from an overloaded function 

398 # when called with Any-typed arguments produces "Incompatible return 

399 # value type (got XY[Any], expected XY[Any])". 

400 return cast( 

401 XY[Any], 

402 self._pixel_to_sky.apply_inverse( 

403 x=ra.to_value(u.rad), 

404 y=dec.to_value(u.rad), 

405 ), 

406 ) 

407 

408 def as_astropy(self, bbox: Box | None = None) -> SkyProjectionAstropyView: 

409 """Return an `astropy.wcs` view of this `SkyProjection`. 

410 

411 Parameters 

412 ---------- 

413 bbox 

414 Bounding box of the array the view will describe. This 

415 projection object is assumed to work on the same coordinate system 

416 in which ``bbox`` is defined, while the Astropy view will consider 

417 the first row and column in that box to be ``(0, 0)``. 

418 

419 Notes 

420 ----- 

421 This returns an object that satisfies the 

422 `astropy.wcs.wcsapi.BaseHighLevelWCS` and 

423 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces while evaluating the 

424 underlying `SkyProjection` itself. It is *not* an `astropy.wcs.WCS` 

425 instance, which is a type that also satisfies those interfaces but 

426 only supports FITS WCS representations (see `as_fits_wcs`). 

427 """ 

428 return SkyProjectionAstropyView(self._pixel_to_sky._ast_mapping, bbox) 

429 

430 def as_fits_wcs(self, bbox: Box, allow_approximation: bool = False) -> astropy.wcs.WCS | None: 

431 """Return a FITS WCS representation of this projection, if possible. 

432 

433 Parameters 

434 ---------- 

435 bbox 

436 Bounding box of the array the FITS WCS will describe. This 

437 transform object is assumed to work on the same coordinate system 

438 in which ``bbox`` is defined, while the FITS WCS will consider the 

439 first row and column in that box to be ``(0, 0)`` (in Astropy 

440 interfaces) or ``(1, 1)`` (in the FITS representation itself). 

441 allow_approximation 

442 If `True` and this `SkyProjection` holds a FITS approximation to 

443 itself, return that approximation. 

444 """ 

445 if allow_approximation and self._fits_approximation: 445 ↛ 446line 445 didn't jump to line 446 because the condition on line 445 was never true

446 return self._fits_approximation.as_fits_wcs(bbox) 

447 return self._pixel_to_sky.as_fits_wcs(bbox) 

448 

449 def _nominal_pixel_scale(self, bbox: Box) -> list[float]: 

450 """Return the nominal pixel scale in arcsec for each sky axis. 

451 

452 Parameters 

453 ---------- 

454 bbox : `Box` 

455 Pixel bounding box over which the scale is characterized. 

456 

457 Returns 

458 ------- 

459 `list` [`float`] 

460 Nominal pixel scale in arcsec/pixel for the longitude and 

461 latitude axes, in that order. 

462 

463 Notes 

464 ----- 

465 This is a port of the Starlink KAPPA ``KPG1_SCALE``/``KPG1_PXSCL`` 

466 routines. At each of a 3x3 grid of test points it perturbs the pixel 

467 position by unit offsets along both axes, finds the neighbour that 

468 moves farthest along each sky axis, and takes the ratio of the 

469 great-circle sky distance to the pixel-space distance; the per-axis 

470 result is the median over the grid. Great-circle distances make the 

471 result correct near the poles and under coordinate rotation. The 

472 scale attaches to the sky axis, so a 90 degree rotation swaps the two 

473 returned values. 

474 """ 

475 offsets = [o for o in itertools.product((0.0, 1.0, -1.0), repeat=2) if o != (0.0, 0.0)] 

476 step_x = 0.3 * bbox.x.size 

477 step_y = 0.3 * bbox.y.size 

478 lon_scales: list[float] = [] 

479 lat_scales: list[float] = [] 

480 for dx, dy in itertools.product((-step_x, 0.0, step_x), (-step_y, 0.0, step_y)): 

481 cx = bbox.x.center + dx 

482 cy = bbox.y.center + dy 

483 center = self.pixel_to_sky(x=cx, y=cy) 

484 neighbours = self.pixel_to_sky( 

485 x=np.array([cx + o[0] for o in offsets]), 

486 y=np.array([cy + o[1] for o in offsets]), 

487 ) 

488 grid0 = self.sky_to_pixel(center) 

489 gx0, gy0 = float(grid0.x), float(grid0.y) 

490 lon0 = center.ra.wrap_at(180 * u.deg) 

491 lat0 = center.dec 

492 # Longitude axis: neighbour with the largest change in RA. 

493 dlon = (neighbours.ra.wrap_at(180 * u.deg) - lon0).wrap_at(180 * u.deg) 

494 probe = SkyCoord(ra=neighbours.ra[int(np.argmax(np.abs(dlon.rad)))], dec=lat0) 

495 grid = self.sky_to_pixel(probe) 

496 dpix = np.hypot(float(grid.x) - gx0, float(grid.y) - gy0) 

497 lon_scales.append(center.separation(probe).to_value(u.arcsec) / dpix) 

498 # Latitude axis: neighbour with the largest change in Dec. 

499 dlat = neighbours.dec - lat0 

500 probe = SkyCoord(ra=lon0, dec=neighbours.dec[int(np.argmax(np.abs(dlat.rad)))]) 

501 grid = self.sky_to_pixel(probe) 

502 dpix = np.hypot(float(grid.x) - gx0, float(grid.y) - gy0) 

503 lat_scales.append(center.separation(probe).to_value(u.arcsec) / dpix) 

504 return [statistics.median(lon_scales), statistics.median(lat_scales)] 

505 

506 def _describe( 

507 self, options: DescribeOptions = DescribeOptions(), /, *, bbox: Box | None = None 

508 ) -> Report: 

509 """Return a `Report` describing this sky projection. 

510 

511 Parameters 

512 ---------- 

513 options : `DescribeOptions`, optional 

514 Rendering options. `DescribeOptions.brief` returns only the type, 

515 title, and summary, skipping the pixel and WCS characterization. 

516 bbox : `Box`, optional 

517 Pixel bounding box to characterize the projection over. Defaults 

518 to the projection's own `pixel_bounds` when it has them. Given 

519 either, the report gains the sky coordinates of the box center and 

520 of the corners of the area the box covers, and the nominal pixel 

521 scale is characterized over the box; given neither, it falls back 

522 to the pixel origin. 

523 """ 

524 if options.brief: 

525 return Report( 

526 type_name="SkyProjection", 

527 title=f"{self.sky_frame.value} coordinates", 

528 summary=f"{type(self.pixel_frame).__name__} → {self.sky_frame.value}", 

529 ) 

530 fields: list[ReportField] = [] 

531 # Characterize over the box the caller supplied, falling back to the 

532 # region the projection itself declares valid. 

533 box = bbox 

534 if box is None and self.pixel_bounds is not None: 

535 box = self.pixel_bounds.bbox 

536 

537 corners_table: list[ReportTable] = [] 

538 if box is not None: 

539 cx, cy = box.x.center, box.y.center 

540 center_sky = self.pixel_to_sky(x=cx, y=cy) 

541 fields.append( 

542 ReportField( 

543 label="center pixel", 

544 value=f"(x={cx:g}, y={cy:g}) → {_format_pixel_sky(center_sky)}", 

545 role=FieldRole.DERIVED, 

546 ) 

547 ) 

548 # Measure the area the box covers rather than the span between 

549 # outermost pixel centers, so let the polygon own the half-pixel 

550 # expansion; its vertex order is not part of its interface, so 

551 # take the bounds from the vertices themselves. 

552 corners = box.to_polygon() 

553 x0, x1 = corners.x_vertices.min(), corners.x_vertices.max() 

554 y0, y1 = corners.y_vertices.min(), corners.y_vertices.max() 

555 origin = self.pixel_to_sky(x=x0, y=y0) 

556 width = origin.separation(self.pixel_to_sky(x=x1, y=y0)) 

557 height = origin.separation(self.pixel_to_sky(x=x0, y=y1)) 

558 # Orientation of the pixel y axis at the box center, which is what 

559 # "E of N" conventionally describes for an image. 

560 up_sky = self.pixel_to_sky(x=cx, y=cy + 1.0) 

561 fields.append( 

562 ReportField( 

563 label="Image extent", 

564 value=_format_extent(width, height, center_sky.position_angle(up_sky)), 

565 role=FieldRole.DERIVED, 

566 ) 

567 ) 

568 else: 

569 # No box is available from either source, so the pixel origin is 

570 # the only place left to sample. It carries no special meaning 

571 # for the projection -- the array this describes may lie far from 

572 # it -- so name the pixel explicitly rather than implying one. 

573 origin_sky = self.pixel_to_sky(x=0, y=0) 

574 fields.append( 

575 ReportField( 

576 label="origin pixel", 

577 value=f"(x=0, y=0) → {_format_pixel_sky(origin_sky)}", 

578 role=FieldRole.DERIVED, 

579 ) 

580 ) 

581 

582 # The nominal pixel scale is characterized over the box when there is 

583 # one, otherwise over a single pixel at the origin. 

584 lon_label, lat_label = self._sky_axis_labels() 

585 lon_scale, lat_scale = self._nominal_pixel_scale(box if box is not None else Box.factory[0:1, 0:1]) 

586 # Report a single value when the two axes round to the same 0.01 arcsec 

587 # figure; otherwise name each sky axis with its own scale. 

588 if round(lon_scale, 2) == round(lat_scale, 2): 

589 fields.append( 

590 ReportField( 

591 label="Nominal pixel scale", value=f"{lon_scale:.2f} arcsec", role=FieldRole.DERIVED 

592 ) 

593 ) 

594 else: 

595 fields.append( 

596 ReportField( 

597 label="Nominal pixel scales", 

598 value=f"{lon_label} {lon_scale:.2f} arcsec; {lat_label} {lat_scale:.2f} arcsec", 

599 role=FieldRole.DERIVED, 

600 ) 

601 ) 

602 

603 if box is not None: 

604 # Walk the same polygon the extent was measured over, whose 

605 # vertices are expanded by half a pixel to cover the full area the 

606 # image occupies on the sky rather than its pixel centers. 

607 rows = [] 

608 for x, y in zip(corners.x_vertices, corners.y_vertices, strict=True): 

609 ra_sex, ra_deg, dec_sex, dec_deg = _sky_parts(self.pixel_to_sky(x=x, y=y)) 

610 rows.append([f"{x:g}", f"{y:g}", ra_sex, dec_sex, f"{ra_deg:.6f}", f"{dec_deg:+.6f}"]) 

611 corners_table.append( 

612 ReportTable( 

613 title="Corners", 

614 columns=["x", "y", "RA", "Dec", "RA (°)", "Dec (°)"], 

615 rows=rows, 

616 ) 

617 ) 

618 

619 # Representability as a FITS WCS depends on the box; without one the 

620 # answer cannot be determined. 

621 if self._fits_approximation is not None: 621 ↛ 622line 621 didn't jump to line 622 because the condition on line 621 was never true

622 fits_wcs = "approximate" 

623 elif box is None: 

624 fits_wcs = "unknown" 

625 else: 

626 fits_wcs = "available" if self.as_fits_wcs(box) is not None else "none" 

627 fields.append(ReportField(label="fits_wcs", value=fits_wcs, role=FieldRole.DERIVED)) 

628 

629 return Report( 

630 type_name="SkyProjection", 

631 title=f"{self.sky_frame.value} coordinates", 

632 summary=f"{type(self.pixel_frame).__name__} → {self.sky_frame.value}", 

633 fields=fields, 

634 tables=corners_table, 

635 ) 

636 

637 def describe( 

638 self, 

639 *, 

640 brief: bool = False, 

641 detail: bool = False, 

642 exclude: Collection[str] = (), 

643 bbox: Box | None = None, 

644 ) -> Report: 

645 """Return a `~lsst.images.Report` describing this sky projection. 

646 

647 Parameters 

648 ---------- 

649 brief : `bool`, optional 

650 Whether to build only what ``repr`` and ``str`` read. 

651 detail : `bool`, optional 

652 Whether to include extras that are too expensive for a default 

653 report. 

654 exclude : `~collections.abc.Collection` [`str`], optional 

655 Names of report elements to omit. 

656 bbox : `~lsst.images.Box`, optional 

657 Pixel bounding box to characterize the projection over. 

658 

659 Returns 

660 ------- 

661 report : `~lsst.images.Report` 

662 Report describing this sky projection. 

663 """ 

664 return self._describe( 

665 DescribeOptions(brief=brief, detail=detail, exclude=frozenset(exclude)), bbox=bbox 

666 ) 

667 

668 def serialize[P: pydantic.BaseModel]( 

669 self, archive: OutputArchive[P], *, use_frame_sets: bool = False 

670 ) -> SkyProjectionSerializationModel[P]: 

671 """Serialize a projection to an archive. 

672 

673 Parameters 

674 ---------- 

675 archive 

676 Archive to serialize to. 

677 use_frame_sets 

678 If `True`, decompose the underlying transform and try to reference 

679 component mappings that were already serialized into a `FrameSet` 

680 in the archive. The FITS approximation transform is never 

681 decomposed. 

682 

683 Returns 

684 ------- 

685 `SkyProjectionSerializationModel` 

686 Serialized form of the projection. 

687 """ 

688 pixel_to_sky = archive.serialize_direct( 

689 "pixel_to_sky", functools.partial(self._pixel_to_sky.serialize, use_frame_sets=use_frame_sets) 

690 ) 

691 fits_approximation = ( 

692 archive.serialize_direct("fits_approximation", self._fits_approximation.serialize) 

693 if self._fits_approximation is not None 

694 else None 

695 ) 

696 return SkyProjectionSerializationModel( 

697 pixel_to_sky=pixel_to_sky, fits_approximation=fits_approximation 

698 ) 

699 

700 @staticmethod 

701 def _get_archive_tree_type[P: pydantic.BaseModel]( 

702 pointer_type: type[P], 

703 ) -> type[SkyProjectionSerializationModel[P]]: 

704 """Return the serialization model type for this object for an archive 

705 type that uses the given pointer type. 

706 """ 

707 return SkyProjectionSerializationModel[pointer_type] # type: ignore 

708 

709 @staticmethod 

710 def from_legacy( 

711 sky_wcs: LegacySkyWcs, pixel_frame: F, pixel_bounds: Bounds | None = None 

712 ) -> SkyProjection[F]: 

713 """Construct a transform from a legacy `lsst.afw.geom.SkyWcs`. 

714 

715 Parameters 

716 ---------- 

717 sky_wcs : `lsst.afw.geom.SkyWcs` 

718 Legacy WCS object. 

719 pixel_frame 

720 Coordinate frame for the pixel grid. 

721 pixel_bounds 

722 The region that bounds valid pixels for this transform. 

723 """ 

724 fits_approximation: Transform[F, SkyFrame] | None = None 

725 if (legacy_fits_approximation := sky_wcs.getFitsApproximation()) is not None: 

726 fits_approximation = Transform( 

727 pixel_frame, 

728 SkyFrame.ICRS, 

729 legacy_fits_approximation.getFrameDict(), 

730 pixel_bounds, 

731 ) 

732 return SkyProjection( 

733 Transform(pixel_frame, SkyFrame.ICRS, sky_wcs.getFrameDict(), pixel_bounds), 

734 fits_approximation=fits_approximation, 

735 ) 

736 

737 def to_legacy(self) -> LegacySkyWcs: 

738 """Convert to a legacy `lsst.afw.geom.SkyWcs` instance.""" 

739 from lsst.afw.geom import SkyWcs as LegacySkyWcs 

740 

741 try: 

742 ast_mapping = astshim.FrameDict(self._pixel_to_sky._ast_mapping) 

743 except TypeError as err: 

744 err.add_note( 

745 "Only Projections created by from_legacy and from_fits_wcs " 

746 "are guaranteed to be convertible to SkyWcs." 

747 ) 

748 raise 

749 # SkyWcs requires the pixel frame's domain to be PIXELS, while AST 

750 # (and hence NDF and from_fits_wcs) uses PIXEL, so rename it in the 

751 # FrameDict (a deep copy, so this projection itself is unchanged). 

752 if ast_mapping.hasDomain("PIXEL") and not ast_mapping.hasDomain("PIXELS"): 752 ↛ 757line 752 didn't jump to line 757 because the condition on line 752 was always true

753 saved_current = ast_mapping.current 

754 ast_mapping.setCurrent("PIXEL") 

755 ast_mapping.setDomain("PIXELS") 

756 ast_mapping.current = saved_current 

757 legacy_wcs = LegacySkyWcs(ast_mapping) 

758 if self.fits_approximation is not None: 758 ↛ 759line 758 didn't jump to line 759 because the condition on line 758 was never true

759 legacy_wcs = legacy_wcs.copyWithFitsApproximation(self.fits_approximation.to_legacy()) 

760 return legacy_wcs 

761 

762 

763class SkyProjectionAstropyView(BaseLowLevelWCS, HighLevelWCSMixin): 

764 """An Astropy-interface view of a `SkyProjection`. 

765 

766 Parameters 

767 ---------- 

768 ast_pixel_to_sky 

769 AST mapping from pixel coordinates to sky coordinates. 

770 bbox 

771 Bounding box of the projection, or `None` if unbounded. 

772 

773 Notes 

774 ----- 

775 The constructor of this classe is considered a private implementation 

776 detail; use `SkyProjection.as_astropy` instead. 

777 

778 This object satisfies the `astropy.wcs.wcsapi.BaseHighLevelWCS` and 

779 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces while evaluating the 

780 underlying `SkyProjection` itself. It is *not* an `astropy.wcs.WCS` 

781 subclass, which is a type that also satisfies those interfaces but 

782 only supports FITS WCS representations (see `SkyProjection.as_fits_wcs`). 

783 """ 

784 

785 def __init__(self, ast_pixel_to_sky: astshim.Mapping, bbox: Box | None) -> None: 

786 self._bbox = bbox 

787 if bbox is not None: 

788 ast_pixel_to_sky = astshim.ShiftMap(list(bbox.start.xy)).then(ast_pixel_to_sky) 

789 self._ast_pixel_to_sky = ast_pixel_to_sky 

790 

791 def __str__(self) -> str: 

792 region = str(self._bbox) if self._bbox is not None else "unbounded" 

793 return f"SkyProjectionAstropyView({region} → ICRS)" 

794 

795 def __repr__(self) -> str: 

796 # The bbox shift is baked into the mapping, so pixel (0, 0) is the 

797 # first array pixel in the view's own (Astropy) convention. 

798 ra, dec = self.pixel_to_world_values(np.zeros(()), np.zeros(())) 

799 reference = _format_pixel_sky(SkyCoord(ra=float(ra) * u.rad, dec=float(dec) * u.rad, frame=ICRS)) 

800 lines = ["SkyProjectionAstropyView", " world axes : ICRS (ra, dec)"] 

801 if self._bbox is not None: 

802 lines.append(f" array shape : {self._bbox.shape}") 

803 lines.append(f" pixel (0, 0): {reference}") 

804 return "\n".join(lines) 

805 

806 @property 

807 def low_level_wcs(self) -> Self: 

808 return self 

809 

810 @property 

811 def array_shape(self) -> YX[int] | None: 

812 return self._bbox.shape if self._bbox is not None else None 

813 

814 @property 

815 def axis_correlation_matrix(self) -> np.ndarray: 

816 return np.array([[True, True], [True, True]]) 

817 

818 @property 

819 def pixel_axis_names(self) -> XY[str]: 

820 return XY("x", "y") 

821 

822 @property 

823 def pixel_bounds(self) -> XY[tuple[int, int]] | None: 

824 if self._bbox is None: 

825 return None 

826 return XY((self._bbox.x.min, self._bbox.x.max), (self._bbox.y.min, self._bbox.y.max)) 

827 

828 @property 

829 def pixel_n_dim(self) -> int: 

830 return 2 

831 

832 @property 

833 def pixel_shape(self) -> XY[int] | None: 

834 array_shape = self.array_shape