Coverage for python/lsst/images/_transforms/

1# This file is part of lsst-images.

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.

9# Use of this source code is governed by a 3-clause BSD-style

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

12from __future__ import annotations

14__all__ = ("Projection", "ProjectionAstropyView", "ProjectionSerializationModel")

16import functools

17from typing import Any, Self, TypeVar, final

19import astropy.units as u

20import astropy.wcs

21import numpy as np

22import pydantic

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

24from astropy.wcs.wcsapi import BaseLowLevelWCS, HighLevelWCSMixin

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

27from ..serialization import ArchiveTree, InputArchive, OutputArchive

28from ..utils import is_none

29from ._frames import Frame, SkyFrame

30from ._transform import Transform, TransformSerializationModel, _ast_apply

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

33# autodoc-typehints plugin.

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

37@final

38class Projection[F: Frame]:

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

41 Parameters

42 ----------

43 pixel_to_sky

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

45 fits_approximation

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

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

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

49 as a FITS WCS.

51 Notes

52 -----

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

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

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

56 """

58 def __init__(

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

60 ):

61 self._pixel_to_sky = pixel_to_sky

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

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

64 if pixel_to_sky.out_frame != SkyFrame.ICRS:

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

66 self._fits_approximation = fits_approximation

68 @staticmethod

69 def from_fits_wcs(

70 fits_wcs: astropy.wcs.WCS,

71 pixel_frame: F,

72 pixel_bounds: Bounds | None = None,

73 x0: int = 0,

74 y0: int = 0,

75 ) -> Projection[F]:

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

78 Parameters

79 ----------

80 fits_wcs

81 FITS WCS to convert.

82 pixel_frame

83 Coordinate frame for the pixel grid.

84 pixel_bounds

85 The region that bounds valid pixels for this transform.

86 x0

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

88 relates to world coordinates.

89 y0

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

91 relates to world coordinates.

93 Notes

94 -----

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

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

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

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

99 integers.

100

101 See Also

102 --------

103 Transform.from_fits_wcs

104 """

105 return Projection(

106 Transform.from_fits_wcs(

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

108 )

109 )

110

111 @property

112 def pixel_frame(self) -> F:

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

114 return self._pixel_to_sky.in_frame

115

116 @property

117 def sky_frame(self) -> SkyFrame:

118 """Coordinate frame for the sky."""

119 return self._pixel_to_sky.out_frame

120

121 @property

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

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

124 return self._pixel_to_sky

125

126 @property

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

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

129 return self._pixel_to_sky.inverted()

130

131 @property

132 def fits_approximation(self) -> Projection[F] | None:

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

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

135 """

136 return Projection(self._fits_approximation) if self._fits_approximation is not None else None

137

138 def pixel_to_sky[T: np.ndarray | float](self, *, x: T, y: T) -> SkyCoord:

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

140

141 Parameters

142 ----------

143 x : `numpy.ndarray` | `float`

144 ``x`` values of the pixel points to transform.

145 y : `numpy.ndarray` | `float`

146 ``y`` values of the pixel points to transform.

147

148 Returns

149 -------

150 astropy.coordinates.SkyCoord

151 Transformed sky coordinates.

152 """

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

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

155

156 def sky_to_pixel(self, sky: SkyCoord) -> XY[np.ndarray | float]:

157 """Transform one or more sky coordinates to pixels

158

159 Parameters

160 ----------

161 sky

162 Sky coordinates to transform.

163

164 Returns

165 -------

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

167 Transformed pixel coordinates.

168 """

169 if sky.frame.name != "icrs":

170 sky = sky.transform_to("icrs")

171 ra: Longitude = sky.ra

172 dec: Latitude = sky.dec

173 return self._pixel_to_sky.apply_inverse(

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

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

176 )

177

178 def as_astropy(self, bbox: Box | None = None) -> ProjectionAstropyView:

179 """Return an `astropy.wcs` view of this `Projection`.

180

181 Parameters

182 ----------

183 bbox

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

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

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

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

188

189 Notes

190 -----

191 This returns an object that satisfies the

192 `astropy.wcs.wcsapi.BaseHighLevelWCS` and

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

194 underlying `Projection` itself. It is *not* an `astropy.wcs.WCS`

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

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

197 """

198 return ProjectionAstropyView(self._pixel_to_sky._ast_mapping, bbox)

199

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

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

202

203 Parameters

204 ----------

205 bbox

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

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

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

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

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

211 allow_approximation

212 If `True` and this `Projection` holds a FITS approximation to

213 itself, return that approximation.

214 """

215 if allow_approximation and self._fits_approximation:

216 return self._fits_approximation.as_fits_wcs(bbox)

217 return self._pixel_to_sky.as_fits_wcs(bbox)

218

219 def serialize[P: pydantic.BaseModel](

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

221 ) -> ProjectionSerializationModel[P]:

222 """Serialize a projection to an archive.

223

224 Parameters

225 ----------

226 archive

227 Archive to serialize to.

228 use_frame_sets

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

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

231 in the archive. The FITS approximation transform is never

232 decomposed.

233

234 Returns

235 -------

236 `ProjectionSerializationModel`

237 Serialized form of the projection.

238 """

239 pixel_to_sky = archive.serialize_direct(

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

241 )

242 fits_approximation = (

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

244 if self._fits_approximation is not None

245 else None

246 )

247 return ProjectionSerializationModel(pixel_to_sky=pixel_to_sky, fits_approximation=fits_approximation)

248

249 @staticmethod

250 def deserialize[P: pydantic.BaseModel](

251 model: ProjectionSerializationModel[P], archive: InputArchive[P]

252 ) -> Projection[Any]:

253 """Deserialize a projection from an archive.

254

255 Parameters

256 ----------

257 model

258 Seralized form of the projection.

259 archive

260 Archive to read from.

261 """

262 pixel_to_sky = Transform.deserialize(model.pixel_to_sky, archive)

263 fits_approximation = (

264 Transform.deserialize(model.fits_approximation, archive)

265 if model.fits_approximation is not None

266 else None

267 )

268 return Projection(pixel_to_sky, fits_approximation=fits_approximation)

269

270 @staticmethod

271 def _get_archive_tree_type[P: pydantic.BaseModel](

272 pointer_type: type[P],

273 ) -> type[ProjectionSerializationModel[P]]:

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

275 type that uses the given pointer type.

276 """

277 return ProjectionSerializationModel[pointer_type] # type: ignore

278

279 @staticmethod

280 def from_legacy(sky_wcs: Any, pixel_frame: F, pixel_bounds: Bounds | None = None) -> Projection[F]:

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

282

283 Parameters

284 ----------

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

286 Legacy WCS object.

287 pixel_frame

288 Coordinate frame for the pixel grid.

289 pixel_bounds

290 The region that bounds valid pixels for this transform.

291 """

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

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

294 fits_approximation = Transform.from_legacy(

295 legacy_fits_approximation.getTransform(), pixel_frame, SkyFrame.ICRS, in_bounds=pixel_bounds

296 )

297 return Projection(

298 Transform.from_legacy(sky_wcs.getTransform(), pixel_frame, SkyFrame.ICRS, in_bounds=pixel_bounds),

299 fits_approximation=fits_approximation,

300 )

301

302

303class ProjectionAstropyView(BaseLowLevelWCS, HighLevelWCSMixin):

304 """An Astropy-interface view of a `Projection`.

305

306 Notes

307 -----

308 The constructor of this classe is considered a private implementation

309 detail; use `Projection.as_astropy` instead.

310

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

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

313 underlying `Projection` itself. It is *not* an `astropy.wcs.WCS`

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

315 only supports FITS WCS representations (see `Projection.as_fits_wcs`).

316 """

317

318 def __init__(self, ast_pixel_to_sky: Any, bbox: Box | None):

319 self._bbox = bbox

320 if bbox is not None:

321 import astshim

322

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

324 self._ast_pixel_to_sky = ast_pixel_to_sky

325

326 @property

327 def low_level_wcs(self) -> Self:

328 return self

329

330 @property

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

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

333

334 @property

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

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

337

338 @property

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

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

341

342 @property

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

344 if self._bbox is None:

345 return None

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

347

348 @property

349 def pixel_n_dim(self) -> int:

350 return 2

351

352 @property

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

354 array_shape = self.array_shape

355 return array_shape.xy if array_shape is not None else None

356

357 @property

358 def serialized_classes(self) -> bool:

359 return False

360

361 @property

362 def world_axis_names(self) -> tuple[str, str]:

363 return ("ra", "dec")

364

365 @property

366 def world_axis_object_classes(self) -> dict[str, tuple[type[SkyCoord], tuple[()], dict[str, Any]]]:

367 return {"celestial": (SkyCoord, (), {"frame": ICRS, "unit": (u.rad, u.rad)})}

368

369 @property

370 def world_axis_object_components(self) -> list[tuple[str, int, str]]:

371 return [("celestial", 0, "spherical.lon.radian"), ("celestial", 1, "spherical.lat.radian")]

372

373 @property

374 def world_axis_physical_types(self) -> tuple[str, str]:

375 return ("pos.eq.ra", "pos.eq.dec")

376

377 @property

378 def world_axis_units(self) -> tuple[str, str]:

379 return ("rad", "rad")

380

381 @property

382 def world_n_dim(self) -> int:

383 return 2

384

385 def pixel_to_world_values(self, x: np.ndarray, y: np.ndarray) -> XY[np.ndarray]:

386 return _ast_apply(self._ast_pixel_to_sky.applyForward, x=x, y=y)

387

388 def world_to_pixel_values(self, ra: np.ndarray, dec: np.ndarray) -> XY[np.ndarray]:

389 return _ast_apply(self._ast_pixel_to_sky.applyInverse, x=ra, y=dec)

390

391

392class ProjectionSerializationModel[P: pydantic.BaseModel](ArchiveTree):

393 """Serialization model for projetions."""

394

395 pixel_to_sky: TransformSerializationModel[P] = pydantic.Field(

396 description="The transform that maps pixel coordinates to the sky."

397 )

398 fits_approximation: TransformSerializationModel[P] | None = pydantic.Field(

399 default=None,

400 description=(

401 "An approximation of the pixel-to-sky transform that is exactly representable as a FITS WCS."

402 ),

403 exclude_if=is_none,

404 )

Coverage for python / lsst / images / _transforms / _projection.py: 52%

135 statements