Coverage for tests/test_transforms.py: 83%
538 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-26 09:47 +0000
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.
12from __future__ import annotations
14import dataclasses
15import functools
16import os
17from typing import Any, ClassVar
19import astropy.units as u
20import astropy.wcs
21import numpy as np
22import pydantic
23import pytest
24from astropy.coordinates import Longitude
26from lsst.images import (
27 ICRS,
28 XY,
29 YX,
30 Box,
31 DetectorFrame,
32 FocalPlaneFrame,
33 GeneralFrame,
34 SkyProjection,
35 Transform,
36 TransformSerializationModel,
37)
38from lsst.images._transforms import _ast as astshim
39from lsst.images.cameras import CameraFrameSet, CameraFrameSetSerializationModel
40from lsst.images.describe import FieldRole, Report
41from lsst.images.fits import PointerModel
42from lsst.images.serialization import ArchiveTree, InputArchive, JsonRef, OutputArchive
43from lsst.images.tests import (
44 DP2_VISIT_DETECTOR_DATA_ID,
45 RoundtripFits,
46 RoundtripJson,
47 check_transform,
48 compare_sky_projection_to_legacy_wcs,
49 legacy_points_to_xy_array,
50 make_random_sky_projection,
51 reset_afw_mask_planes, # noqa: F401
52)
54EXTERNAL_DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None)
57@pytest.fixture(scope="session")
58def legacy_camera() -> Any:
59 """Return a legacy Camera loaded from camera.fits.
61 Skips if TESTDATA_IMAGES_DIR is unset or lsst.afw.cameraGeom is
62 unavailable.
63 """
64 if EXTERNAL_DATA_DIR is None: 64 ↛ 66line 64 didn't jump to line 66 because the condition on line 64 was always true
65 pytest.skip("TESTDATA_IMAGES_DIR is not in the environment.")
66 try:
67 from lsst.afw.cameraGeom import Camera
68 except ImportError:
69 pytest.skip("'lsst.afw.cameraGeom' could not be imported.")
70 filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "camera.fits")
71 return Camera.readFits(filename)
74@pytest.fixture
75def legacy_detector_wcs(reset_afw_mask_planes: None) -> dict[str, Any]: # noqa: F811
76 """Return WCS-related objects read from visit_image.fits.
78 Skips if TESTDATA_IMAGES_DIR is unset or lsst.afw.image is unavailable.
79 """
80 # reset_afw_mask_planes will have already skipped if afw is not available.
81 from lsst.afw.image import ExposureFitsReader
83 if EXTERNAL_DATA_DIR is None: 83 ↛ 85line 83 didn't jump to line 85 because the condition on line 83 was always true
84 pytest.skip("TESTDATA_IMAGES_DIR is not in the environment.")
85 filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "visit_image.fits")
86 reader = ExposureFitsReader(filename)
87 return {
88 "legacy_wcs": reader.readWcs(),
89 "wcs_bbox": Box.from_legacy(reader.readDetector().getBBox()),
90 "subimage_bbox": Box.from_legacy(reader.readBBox()),
91 }
94def test_identity() -> None:
95 """Test an identity transform."""
96 frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.factory[:5, :4])
97 xy = frame.bbox.meshgrid().map(np.ravel)
98 identity = Transform.identity(frame)
99 check_transform(identity, xy, xy, frame, frame)
100 assert identity.decompose() == []
101 with RoundtripJson(identity) as roundtrip:
102 pass
103 check_transform(roundtrip.result, xy, xy, frame, frame)
106def test_transform_equality() -> None:
107 """Test Transform.__eq__ across all of its comparison branches."""
108 pixel_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.factory[:5, :4])
109 focal_plane = FocalPlaneFrame(instrument="LSSTCam", visit=1, unit=u.mm)
110 # A distinct frame for the in-frame and out-frame branches.
111 alt_frame = DetectorFrame(instrument="LSSTCam", visit=1, detector=12, bbox=Box.factory[:5, :4])
112 in_bounds = Box.factory[:5, :4]
113 out_bounds = Box.factory[:10, :8]
115 def make(
116 *,
117 in_frame: Any = pixel_frame,
118 out_frame: Any = focal_plane,
119 ast_mapping: astshim.Mapping | None = None,
120 in_bounds_: Box | None = in_bounds,
121 out_bounds_: Box | None = out_bounds,
122 components: Any = (),
123 ) -> Transform[Any, Any]:
124 return Transform(
125 in_frame,
126 out_frame,
127 ast_mapping if ast_mapping is not None else astshim.UnitMap(2),
128 in_bounds=in_bounds_,
129 out_bounds=out_bounds_,
130 components=components,
131 )
133 base = make()
135 # Identity short-circuit: an object is always equal to itself.
136 assert base == base
138 # Two independently constructed but equivalent transforms are equal,
139 # and equality is symmetric.
140 assert base == make()
141 assert make() == base
143 # Comparison against a non-Transform yields NotImplemented, so Python
144 # falls back to identity: the objects are unequal and != is True.
145 assert not (base == "not a transform")
146 assert base != "not a transform"
147 assert base != None # noqa: E711
148 assert base != 42
150 # Each remaining branch differs from base in exactly one attribute.
151 assert base != make(ast_mapping=astshim.ShiftMap([1.0, 2.0]))
152 assert base != make(in_bounds_=out_bounds)
153 assert base != make(out_bounds_=in_bounds)
154 assert base != make(in_frame=alt_frame)
155 assert base != make(out_frame=alt_frame)
156 assert base != make(components=[Transform.identity(alt_frame)])
159def test_sky_projection_equality() -> None:
160 """Test SkyProjection.__eq__ across all of its comparison branches."""
161 pixel_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.factory[:5, :4])
163 # Check the two failure modes.
164 with pytest.raises(ValueError, match="not a mapping from pixel coordinates"):
165 SkyProjection(Transform(ICRS, ICRS, astshim.UnitMap(2)))
167 with pytest.raises(ValueError, match="not a mapping to ICRS"):
168 SkyProjection(Transform(pixel_frame, pixel_frame, astshim.UnitMap(2)))
170 def make_pixel_to_sky(ast_mapping: astshim.Mapping | None = None) -> Transform[Any, Any]:
171 mapping = ast_mapping if ast_mapping is not None else astshim.UnitMap(2)
172 return Transform(pixel_frame, ICRS, mapping)
174 base = SkyProjection(make_pixel_to_sky())
176 # Identity short-circuit: an object is always equal to itself.
177 assert base == base
179 # Two independently constructed but equivalent projections are equal.
180 assert base == SkyProjection(make_pixel_to_sky())
182 # Comparison against a non-SkyProjection yields NotImplemented.
183 assert not (base == "not a projection")
184 assert base != "not a projection"
185 assert base != None # noqa: E711
187 # Differ only in the pixel-to-sky transform.
188 assert base != SkyProjection(make_pixel_to_sky(astshim.ShiftMap([1.0, 2.0])))
190 # The fits_approximation branch: absent on base but present here.
191 with_approx = SkyProjection(
192 make_pixel_to_sky(), fits_approximation=make_pixel_to_sky(astshim.ShiftMap([0.1, 0.2]))
193 )
194 assert base != with_approx
196 # Equal pixel-to-sky and equal fits_approximations are equal.
197 with_approx_again = SkyProjection(
198 make_pixel_to_sky(), fits_approximation=make_pixel_to_sky(astshim.ShiftMap([0.1, 0.2]))
199 )
200 assert with_approx == with_approx_again
202 # Same pixel-to-sky transform but a different fits_approximation.
203 other_approx = SkyProjection(
204 make_pixel_to_sky(), fits_approximation=make_pixel_to_sky(astshim.ShiftMap([0.3, 0.4]))
205 )
206 assert with_approx != other_approx
209def test_affine_2x2() -> None:
210 """Test an affine transform constructed from a 2x2 matrix."""
211 in_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.factory[:5, :4])
212 out_frame = GeneralFrame(unit=u.pix)
213 transform_matrix = np.array([[2.0, 0.25], [-0.75, 0.8]])
214 in_xy = in_frame.bbox.meshgrid().map(np.ravel)
215 in_matrix = np.array([in_xy.x, in_xy.y])
216 out_matrix = np.dot(transform_matrix, in_matrix)
217 check_transform(
218 Transform.affine(in_frame, out_frame, transform_matrix),
219 in_xy,
220 XY(x=out_matrix[0, :], y=out_matrix[1, :]),
221 in_frame,
222 out_frame,
223 in_atol=1e-15 * u.pix,
224 out_atol=1e-15 * u.pix,
225 )
228def test_affine_3x3() -> None:
229 """Test an affine transform constructed from a 3x3 matrix."""
230 in_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.factory[:5, :4])
231 out_frame = GeneralFrame(unit=u.pix)
232 transform_matrix = np.array([[2.0, 0.25, -0.5], [-0.75, 0.8, 0.4], [0.0, 0.0, 1.0]])
233 in_xy = in_frame.bbox.meshgrid().map(np.ravel)
234 in_matrix = np.array([in_xy.x, in_xy.y, np.ones(in_xy.x.shape)])
235 out_matrix = np.dot(transform_matrix, in_matrix)
236 check_transform(
237 Transform.affine(in_frame, out_frame, transform_matrix),
238 in_xy,
239 XY(x=out_matrix[0, :], y=out_matrix[1, :]),
240 in_frame,
241 out_frame,
242 in_atol=1e-15 * u.pix,
243 out_atol=1e-15 * u.pix,
244 )
247def compare_to_legacy_camera(legacy_camera: Any, frame_set: CameraFrameSet) -> None:
248 """Assert that transforms extracted from a CameraFrameSet match the
249 legacy afw implementations.
250 """
251 from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, PIXELS
252 from lsst.geom import Point2D
254 legacy_detector = legacy_camera[16]
255 pixel_legacy_points = [Point2D(50.0, 60.0), Point2D(801.2, 322.8), Point2D(33.5, 22.1)]
256 fp_legacy_points = [legacy_detector.transform(p, PIXELS, FOCAL_PLANE) for p in pixel_legacy_points]
257 fa_legacy_points = [legacy_detector.transform(p, PIXELS, FIELD_ANGLE) for p in pixel_legacy_points]
258 pixel_xy_array = legacy_points_to_xy_array(pixel_legacy_points)
259 fp_xy_array = legacy_points_to_xy_array(fp_legacy_points)
260 fa_xy_array = legacy_points_to_xy_array(fa_legacy_points)
261 # Test transforms extracted directly from the frame set.
262 pixel_to_fp = frame_set[frame_set.detector(16), frame_set.focal_plane()]
263 check_transform(pixel_to_fp, pixel_xy_array, fp_xy_array, frame_set.detector(16), frame_set.focal_plane())
264 pixel_to_fa = frame_set[frame_set.detector(16), frame_set.field_angle()]
265 check_transform(pixel_to_fa, pixel_xy_array, fa_xy_array, frame_set.detector(16), frame_set.field_angle())
266 fp_to_fa = frame_set[frame_set.focal_plane(), frame_set.field_angle()]
267 check_transform(fp_to_fa, fp_xy_array, fa_xy_array, frame_set.focal_plane(), frame_set.field_angle())
268 # Test a composition.
269 pixel_to_fa_indirect = pixel_to_fp.then(fp_to_fa)
270 check_transform(
271 pixel_to_fa_indirect,
272 pixel_xy_array,
273 fa_xy_array,
274 frame_set.detector(16),
275 frame_set.field_angle(),
276 )
277 pixel_to_fp_d, fp_to_fa_d = pixel_to_fa_indirect.decompose()
278 check_transform(
279 pixel_to_fp_d, pixel_xy_array, fp_xy_array, frame_set.detector(16), frame_set.focal_plane()
280 )
281 check_transform(fp_to_fa_d, fp_xy_array, fa_xy_array, frame_set.focal_plane(), frame_set.field_angle())
282 fa_to_fp_d, fp_to_pixel_d = pixel_to_fa_indirect.inverted().decompose()
283 check_transform(fa_to_fp_d, fa_xy_array, fp_xy_array, frame_set.field_angle(), frame_set.focal_plane())
284 check_transform(
285 fp_to_pixel_d, fp_xy_array, pixel_xy_array, frame_set.focal_plane(), frame_set.detector(16)
286 )
289def test_camera(legacy_camera: Any) -> None:
290 """Test CameraFrameSet construction, transforms, and FITS/JSON
291 serialization round-trips.
293 Also verifies the archive system's pointer and frame-set reference
294 machinery.
295 """
296 legacy_camera = legacy_camera
297 frame_set = CameraFrameSet.from_legacy(legacy_camera)
298 detector_id: int = DP2_VISIT_DETECTOR_DATA_ID["detector"]
299 compare_to_legacy_camera(legacy_camera, frame_set)
300 test_holder = FrameSetTestHolder(
301 frames=frame_set,
302 pixels_to_fp=frame_set[frame_set.detector(detector_id), frame_set.focal_plane()],
303 )
304 with RoundtripFits(test_holder) as roundtrip1:
305 assert len(roundtrip1.serialized.pixels_to_fp.frames) == 2
306 assert len(roundtrip1.serialized.pixels_to_fp.bounds) == 2
307 assert len(roundtrip1.serialized.pixels_to_fp.mappings) == 1
308 # Instead of storing the AST mapping directly, we should have
309 # stored a reference to the frame set:
310 assert isinstance(roundtrip1.serialized.pixels_to_fp.mappings[0], PointerModel)
311 compare_to_legacy_camera(legacy_camera, roundtrip1.result.frames)
312 assert roundtrip1.result.pixels_to_fp.in_frame == frame_set.detector(detector_id)
313 assert roundtrip1.result.pixels_to_fp.out_frame == frame_set.focal_plane()
314 assert (
315 roundtrip1.result.pixels_to_fp._ast_mapping.simplified().show()
316 == test_holder.pixels_to_fp._ast_mapping.simplified().show()
317 )
318 with RoundtripJson(test_holder) as roundtrip2:
319 assert len(roundtrip2.serialized.pixels_to_fp.frames) == 2
320 assert len(roundtrip2.serialized.pixels_to_fp.bounds) == 2
321 assert len(roundtrip2.serialized.pixels_to_fp.mappings) == 1
322 # Instead of storing the AST mapping directly, we should have
323 # stored a reference to the frame set:
324 assert isinstance(roundtrip2.serialized.pixels_to_fp.mappings[0], JsonRef)
325 raw_data = roundtrip2.inspect()
326 assert len(raw_data["indirect"]) == 1
327 assert raw_data["frames"] == {"$ref": "#/indirect/0"}
328 compare_to_legacy_camera(legacy_camera, roundtrip2.result.frames)
329 assert roundtrip2.result.pixels_to_fp.in_frame == frame_set.detector(detector_id)
330 assert roundtrip2.result.pixels_to_fp.out_frame == frame_set.focal_plane()
331 assert (
332 roundtrip2.result.pixels_to_fp._ast_mapping.simplified().show()
333 == test_holder.pixels_to_fp._ast_mapping.simplified().show()
334 )
337def test_fits_wcs_projection_to_legacy() -> None:
338 """Verify that a projection created by from_fits_wcs can be converted
339 to a legacy SkyWcs.
341 The AST pixel frame uses the domain PIXEL, while lsst.afw.geom.SkyWcs
342 requires PIXELS, so to_legacy has to rename it.
343 """
344 pytest.importorskip("lsst.afw.geom")
345 rng = np.random.default_rng(43)
346 bbox = Box.factory[75:275, 25:225]
347 pixel_frame = GeneralFrame(unit=u.pix)
348 sky_projection = make_random_sky_projection(rng, pixel_frame, bbox)
349 legacy_wcs = sky_projection.to_legacy()
350 compare_sky_projection_to_legacy_wcs(sky_projection, legacy_wcs, pixel_frame, bbox, is_fits=True)
351 # The conversion must not modify the projection in place: its own AST
352 # mapping keeps the PIXEL domain, and the conversion is repeatable.
353 frame_set = sky_projection.pixel_to_sky_transform._ast_mapping
354 assert isinstance(frame_set, astshim.FrameSet)
355 domains = {frame_set.getFrame(i, copy=False).domain for i in range(1, frame_set.nFrame + 1)}
356 assert "PIXEL" in domains
357 assert "PIXELS" not in domains
358 compare_sky_projection_to_legacy_wcs(
359 sky_projection, sky_projection.to_legacy(), pixel_frame, bbox, is_fits=True
360 )
363def test_as_fits_wcs_unrepresentable_returns_none_not_keyerror() -> None:
364 """Test that a transform that is not exactly FITS representable returns
365 None from ``as_fits_wcs`` rather than crashing on the partial header AST
366 writes.
368 A pixel->sky mapping that AST can only partially encode (here a polynomial
369 that produces no valid primary celestial WCS) writes a nonzero number of
370 FITS cards whose header lacks a primary ``CTYPE``; constructing
371 ``astropy.wcs.WCS`` from it previously raised ``KeyError``.
372 """
373 pixel_frame = GeneralFrame(unit=u.pix)
374 coeff_f = np.array(
375 [
376 [1.0, 1, 0, 0], # out1 += 1.0
377 [1.0, 2, 0, 0], # out2 += 1.0
378 [1.0, 1, 1, 0], # out1 += x
379 [1.0, 2, 0, 1], # out2 += y
380 [1e-6, 1, 2, 0], # out1 += 1e-6 x^2
381 [1e-6, 2, 0, 2], # out2 += 1e-6 y^2
382 ]
383 )
384 sky_projection = SkyProjection(Transform(pixel_frame, ICRS, astshim.PolyMap(coeff_f, 2, "")))
385 bbox = Box.factory[0:32, 0:32]
386 assert sky_projection.as_fits_wcs(bbox, allow_approximation=False) is None
389def test_detector_wcs(legacy_detector_wcs: dict[str, Any]) -> None:
390 """Test the Transform/SkyProjection representation of a detector WCS."""
391 legacy_wcs = legacy_detector_wcs["legacy_wcs"]
392 wcs_bbox = legacy_detector_wcs["wcs_bbox"]
393 subimage_bbox = legacy_detector_wcs["subimage_bbox"]
394 detector_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=wcs_bbox)
395 sky_projection = SkyProjection.from_legacy(legacy_wcs, detector_frame)
396 assert sky_projection.fits_approximation is not None
397 compare_sky_projection_to_legacy_wcs(sky_projection, legacy_wcs, detector_frame, subimage_bbox)
398 # When we convert from a legacy SkyWcs, the internal AST Mapping needs
399 # to really be an AST FrameSet in order to be able to convert back.
400 assert "Begin FrameSet" in sky_projection.show()
401 compare_sky_projection_to_legacy_wcs(
402 sky_projection, sky_projection.to_legacy(), detector_frame, subimage_bbox
403 )
404 assert "Begin FrameSet" in sky_projection.fits_approximation.show()
405 compare_sky_projection_to_legacy_wcs(
406 sky_projection.fits_approximation,
407 sky_projection.fits_approximation.to_legacy(),
408 detector_frame,
409 subimage_bbox,
410 is_fits=True,
411 )
412 with RoundtripJson(sky_projection, "SkyProjection") as roundtrip:
413 pass
414 compare_sky_projection_to_legacy_wcs(roundtrip.result, legacy_wcs, detector_frame, subimage_bbox)
415 # The AST FrameSet-ness needs to propagate through serialization.
416 assert "Begin FrameSet" in roundtrip.result.show()
417 compare_sky_projection_to_legacy_wcs(
418 sky_projection, roundtrip.result.to_legacy(), detector_frame, subimage_bbox
419 )
420 with RoundtripJson(sky_projection.fits_approximation, "SkyProjection") as roundtrip:
421 pass
422 compare_sky_projection_to_legacy_wcs(
423 roundtrip.result,
424 legacy_wcs.getFitsApproximation(),
425 detector_frame,
426 subimage_bbox,
427 is_fits=True,
428 )
429 assert "Begin FrameSet" in roundtrip.result.show()
430 compare_sky_projection_to_legacy_wcs(
431 sky_projection.fits_approximation,
432 roundtrip.result.to_legacy(),
433 detector_frame,
434 subimage_bbox,
435 is_fits=True,
436 )
439def test_detector_wcs_full_bbox(legacy_detector_wcs: dict[str, Any]) -> None:
440 """Test that a from_legacy sky projection matchesthe legacy WCS over the
441 whole detector bbox, not just a subregion away from the origin.
443 We use a smaller box in most other tests for performance reasons, but this
444 one is more senstive to the precision floor close to the origin in
445 particular.
446 """
447 legacy_wcs = legacy_detector_wcs["legacy_wcs"]
448 wcs_bbox = legacy_detector_wcs["wcs_bbox"]
449 detector_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=wcs_bbox)
450 sky_projection = SkyProjection.from_legacy(legacy_wcs, detector_frame)
451 compare_sky_projection_to_legacy_wcs(sky_projection, legacy_wcs, detector_frame, wcs_bbox)
454@dataclasses.dataclass
455class FrameSetTestHolder:
456 """A top-level object that holds a CameraFrameSet and a transform
457 extracted from it, for testing archive pointers and frame set references.
458 """
460 frames: CameraFrameSet
461 pixels_to_fp: Transform[DetectorFrame, FocalPlaneFrame]
463 def serialize[P: pydantic.BaseModel](self, archive: OutputArchive[P]) -> FrameSetTestHolderModel[P]:
464 frames_model = archive.serialize_frame_set(
465 "frames", self.frames, self.frames.serialize, key=id(self.frames)
466 )
467 pixels_to_fp_model = archive.serialize_direct(
468 "pixels_to_fp", functools.partial(self.pixels_to_fp.serialize, use_frame_sets=True)
469 )
470 return FrameSetTestHolderModel[P](frames=frames_model, pixels_to_fp=pixels_to_fp_model)
472 @staticmethod
473 def _get_archive_tree_type[P: pydantic.BaseModel](
474 pointer_type: type[P],
475 ) -> type[FrameSetTestHolderModel[P]]:
476 return FrameSetTestHolderModel[pointer_type] # type: ignore
479class FrameSetTestHolderModel[P: pydantic.BaseModel](ArchiveTree):
480 """The serialization model for FrameSetTestHolder."""
482 SCHEMA_NAME: ClassVar[str] = "_test_frame_set_holder"
483 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
484 MIN_READ_VERSION: ClassVar[int] = 1
485 PUBLIC_TYPE: ClassVar[type] = FrameSetTestHolder
487 frames: CameraFrameSetSerializationModel | P
488 pixels_to_fp: TransformSerializationModel[P]
490 def deserialize(self, archive: InputArchive[Any]) -> FrameSetTestHolder:
491 assert not isinstance(self.frames, CameraFrameSetSerializationModel), "Archive pointer expected."
492 frames = archive.deserialize_pointer(
493 self.frames, CameraFrameSetSerializationModel, CameraFrameSetSerializationModel.deserialize
494 )
495 pixels_to_fp = self.pixels_to_fp.deserialize(archive)
496 return FrameSetTestHolder(frames, pixels_to_fp)
499@dataclasses.dataclass
500class _BroadcastTestData:
501 """Shared inputs for broadcasting/scalar tests on Transform and
502 SkyProjection.
503 """
505 in_frame: DetectorFrame
506 out_frame: GeneralFrame
507 matrix: np.ndarray
508 scalar_x: float
509 scalar_y: float
510 xv: list[int]
511 yv: list[int]
512 sky_proj: SkyProjection[DetectorFrame]
515@pytest.fixture
516def broadcast_test_data() -> _BroadcastTestData:
517 """Return shared inputs for broadcasting/scalar tests."""
518 in_frame = DetectorFrame(instrument="Inst", visit=1, detector=1, bbox=Box.factory[0:20, 0:20])
519 return _BroadcastTestData(
520 in_frame=in_frame,
521 out_frame=GeneralFrame(unit=u.pix),
522 matrix=np.array([[2.0, 0.5], [-0.3, 1.5]]),
523 scalar_x=3.0,
524 scalar_y=7.0,
525 xv=[1, 2, 3],
526 yv=[4, 5, 6],
527 sky_proj=make_random_sky_projection(np.random.default_rng(42), in_frame, in_frame.bbox),
528 )
531def test_apply_forward_scalar(broadcast_test_data: _BroadcastTestData) -> None:
532 """Verify that apply_forward and apply_inverse accept scalar x/y and return
533 scalar floats, and that the _q variants accept scalar Quantity inputs.
534 """
535 t = broadcast_test_data.sky_proj.pixel_to_sky_transform
536 # apply_forward with Python float scalars should return XY of floats.
537 result_fwd = t.apply_forward(x=broadcast_test_data.scalar_x, y=broadcast_test_data.scalar_y)
538 assert type(result_fwd.x) is float
539 assert type(result_fwd.y) is float
540 # Values must match the corresponding single-element array call.
541 ref_fwd = t.apply_forward(
542 x=np.array([broadcast_test_data.scalar_x]), y=np.array([broadcast_test_data.scalar_y])
543 )
544 assert result_fwd.x == ref_fwd.x[0]
545 assert result_fwd.y == ref_fwd.y[0]
546 # apply_inverse round-trips back to the original scalars.
547 result_inv = t.apply_inverse(x=result_fwd.x, y=result_fwd.y)
548 assert type(result_inv.x) is float
549 assert type(result_inv.y) is float
550 np.testing.assert_allclose(result_inv.x, broadcast_test_data.scalar_x, atol=1e-12)
551 np.testing.assert_allclose(result_inv.y, broadcast_test_data.scalar_y, atol=1e-12)
552 # apply_forward_q / apply_inverse_q with scalar Quantity inputs.
553 x_q = broadcast_test_data.scalar_x * t.in_frame.unit
554 y_q = broadcast_test_data.scalar_y * t.in_frame.unit
555 result_fwd_q = t.apply_forward_q(x=x_q, y=y_q)
556 assert result_fwd_q.x.shape == ()
557 assert result_fwd_q.y.shape == ()
558 np.testing.assert_allclose(result_fwd_q.x.to_value(t.out_frame.unit), result_fwd.x, atol=1e-12)
559 result_inv_q = t.apply_inverse_q(x=result_fwd_q.x, y=result_fwd_q.y)
560 assert result_inv_q.x.shape == ()
561 np.testing.assert_allclose(
562 result_inv_q.x.to_value(t.in_frame.unit), broadcast_test_data.scalar_x, atol=1e-12
563 )
566def test_apply_array_like_and_integer_input(broadcast_test_data: _BroadcastTestData) -> None:
567 """Verify that apply_forward accepts Python lists and integer-dtype
568 arrays, returning float64 ndarray results consistent with float64
569 array input.
570 """
571 t = broadcast_test_data.sky_proj.pixel_to_sky_transform
572 # Python list input should return an ndarray.
573 result_list = t.apply_forward(x=broadcast_test_data.xv, y=broadcast_test_data.yv)
574 assert isinstance(result_list.x, np.ndarray)
575 assert isinstance(result_list.y, np.ndarray)
576 ref = t.apply_forward(x=np.array(broadcast_test_data.xv), y=np.array(broadcast_test_data.yv))
577 np.testing.assert_array_equal(result_list.x, ref.x)
578 np.testing.assert_array_equal(result_list.y, ref.y)
579 # Integer dtype arrays should not raise and should return float64.
580 xi = np.array(broadcast_test_data.xv, dtype=np.int32)
581 yi = np.array(broadcast_test_data.yv, dtype=np.int32)
582 result_int = t.apply_forward(x=xi, y=yi)
583 assert result_int.x.dtype == np.float64
584 assert result_int.y.dtype == np.float64
585 np.testing.assert_array_equal(result_int.x, ref.x)
586 np.testing.assert_array_equal(result_int.y, ref.y)
589def test_apply_broadcast(broadcast_test_data: _BroadcastTestData) -> None:
590 """Verify that apply_forward and apply_inverse broadcast x and y like
591 a NumPy ufunc, in both 1-D and 2-D cases.
592 """
593 t = broadcast_test_data.sky_proj.pixel_to_sky_transform
594 xv = np.array(broadcast_test_data.xv)
595 yv = np.array(broadcast_test_data.yv + [7])
596 # 1-D broadcast: array x, scalar y.
597 result_1d = t.apply_forward(x=xv, y=broadcast_test_data.scalar_y)
598 assert isinstance(result_1d.x, np.ndarray)
599 assert result_1d.x.shape == xv.shape
600 ref_1d = t.apply_forward(x=xv, y=np.full_like(xv, broadcast_test_data.scalar_y))
601 np.testing.assert_array_equal(result_1d.x, ref_1d.x)
602 np.testing.assert_array_equal(result_1d.y, ref_1d.y)
603 # 2-D broadcast: column x (M,1) × row y (1,N) -> (M,N).
604 x2d = xv[:, np.newaxis] # shape (3, 1)
605 y2d = yv[np.newaxis, :] # shape (1, 4)
606 result_2d = t.apply_forward(x=x2d, y=y2d)
607 assert result_2d.x.shape == (3, 4)
608 assert result_2d.y.shape == (3, 4)
609 # Values must match the fully expanded meshgrid call.
610 xmesh, ymesh = np.meshgrid(xv, yv, indexing="ij")
611 ref_2d = t.apply_forward(x=xmesh, y=ymesh)
612 np.testing.assert_array_equal(result_2d.x, ref_2d.x)
613 np.testing.assert_array_equal(result_2d.y, ref_2d.y)
614 # apply_inverse also broadcasts.
615 result_inv_2d = t.apply_inverse(x=result_2d.x, y=result_2d.y)
616 assert result_inv_2d.x.shape == (3, 4)
617 np.testing.assert_allclose(result_inv_2d.x, xmesh, atol=1e-12)
618 np.testing.assert_allclose(result_inv_2d.y, ymesh, atol=1e-12)
621def test_sky_projection_broadcast(broadcast_test_data: _BroadcastTestData) -> None:
622 """Verify that SkyProjection.pixel_to_sky, sky_to_pixel, and the
623 Astropy view broadcast x and y like a NumPy ufunc.
624 """
625 p = broadcast_test_data
626 xv = np.array(p.xv)
627 yv = np.array(p.yv + [7])
628 # 1-D broadcast: array x, scalar y.
629 sc_1d = p.sky_proj.pixel_to_sky(x=xv, y=p.scalar_y)
630 assert sc_1d.shape == xv.shape
631 ref_1d = p.sky_proj.pixel_to_sky(x=xv, y=np.full_like(xv, p.scalar_y))
632 np.testing.assert_allclose(sc_1d.ra.rad, ref_1d.ra.rad, atol=1e-12)
633 np.testing.assert_allclose(sc_1d.dec.rad, ref_1d.dec.rad, atol=1e-12)
634 # 2-D broadcast: column x (M,1) × row y (1,N) -> (M,N).
635 x2d = xv[:, np.newaxis] # shape (3, 1)
636 y2d = yv[np.newaxis, :] # shape (1, 4)
637 sc_2d = p.sky_proj.pixel_to_sky(x=x2d, y=y2d)
638 assert sc_2d.shape == (3, 4)
639 xmesh, ymesh = np.meshgrid(xv, yv, indexing="ij")
640 ref_2d = p.sky_proj.pixel_to_sky(x=xmesh, y=ymesh)
641 np.testing.assert_allclose(sc_2d.ra.rad, ref_2d.ra.rad, atol=1e-12)
642 np.testing.assert_allclose(sc_2d.dec.rad, ref_2d.dec.rad, atol=1e-12)
643 # sky_to_pixel round-trips back to the original grid.
644 pix_2d = p.sky_proj.sky_to_pixel(sc_2d)
645 assert pix_2d.x.shape == (3, 4)
646 np.testing.assert_allclose(pix_2d.x, xmesh, atol=1e-9)
647 np.testing.assert_allclose(pix_2d.y, ymesh, atol=1e-9)
648 # SkyProjectionAstropyView.pixel_to_world_values also broadcasts.
649 view = p.sky_proj.as_astropy()
650 world_2d = view.pixel_to_world_values(x2d, y2d)
651 assert world_2d[0].shape == (3, 4)
652 assert world_2d[1].shape == (3, 4)
653 np.testing.assert_allclose(world_2d[0], ref_2d.ra.rad, atol=1e-12)
654 np.testing.assert_allclose(world_2d[1], ref_2d.dec.rad, atol=1e-12)
655 # SkyProjectionAstropyView.world_to_pixel_values also broadcasts.
656 ra_2d = ref_2d.ra.rad[:, np.newaxis, :] # (3, 1, 4) — over-broadcast to check
657 dec_2d = ref_2d.dec.rad[np.newaxis, :, :] # (1, 3, 4)
658 pix_world = view.world_to_pixel_values(ra_2d, dec_2d)
659 assert pix_world[0].shape == (3, 3, 4)
662def test_apply_xy_yx(broadcast_test_data: _BroadcastTestData) -> None:
663 """Verify that apply_forward, apply_inverse, and the _q variants accept
664 XY and YX positional arguments, producing results identical to the
665 equivalent x=/y= keyword calls.
666 """
667 p = broadcast_test_data
668 t = p.sky_proj.pixel_to_sky_transform
669 sx, sy = p.scalar_x, p.scalar_y
670 xv, yv = np.array(p.xv, dtype=float), np.array(p.yv, dtype=float)
672 # --- apply_forward: scalar ---
673 ref_fwd = t.apply_forward(x=sx, y=sy)
674 assert t.apply_forward(XY(sx, sy)) == ref_fwd
675 assert t.apply_forward(YX(sy, sx)) == ref_fwd
677 # --- apply_forward: array ---
678 ref_fwd_arr = t.apply_forward(x=xv, y=yv)
679 np.testing.assert_array_equal(t.apply_forward(XY(xv, yv)).x, ref_fwd_arr.x)
680 np.testing.assert_array_equal(t.apply_forward(YX(yv, xv)).x, ref_fwd_arr.x)
682 # --- apply_inverse: scalar ---
683 ref_inv = t.apply_inverse(x=ref_fwd.x, y=ref_fwd.y)
684 assert t.apply_inverse(XY(ref_fwd.x, ref_fwd.y)) == ref_inv
685 assert t.apply_inverse(YX(ref_fwd.y, ref_fwd.x)) == ref_inv
687 # --- apply_forward_q: scalar ---
688 x_q = sx * t.in_frame.unit
689 y_q = sy * t.in_frame.unit
690 ref_fwd_q = t.apply_forward_q(x=x_q, y=y_q)
691 result_q = t.apply_forward_q(XY(x_q, y_q))
692 np.testing.assert_allclose(result_q.x.value, ref_fwd_q.x.value, atol=1e-12)
693 result_q_yx = t.apply_forward_q(YX(y_q, x_q))
694 np.testing.assert_allclose(result_q_yx.x.value, ref_fwd_q.x.value, atol=1e-12)
696 # --- apply_inverse_q: scalar ---
697 ref_inv_q = t.apply_inverse_q(x=ref_fwd_q.x, y=ref_fwd_q.y)
698 result_inv_q = t.apply_inverse_q(XY(ref_fwd_q.x, ref_fwd_q.y))
699 np.testing.assert_allclose(result_inv_q.x.value, ref_inv_q.x.value, atol=1e-12)
701 # --- TypeError on bad combinations ---
702 with pytest.raises(TypeError):
703 t.apply_forward(XY(sx, sy), x=sx)
704 with pytest.raises(TypeError):
705 t.apply_forward(YX(sy, sx), y=sy)
706 with pytest.raises(TypeError):
707 t.apply_forward()
710def test_pixel_to_sky_xy_yx(broadcast_test_data: _BroadcastTestData) -> None:
711 """Verify that SkyProjection.pixel_to_sky accepts XY and YX positional
712 arguments, producing results identical to the x=/y= keyword form.
713 """
714 p = broadcast_test_data
715 sx, sy = p.scalar_x, p.scalar_y
716 xv, yv = np.array(p.xv, dtype=float), np.array(p.yv, dtype=float)
718 # Scalar XY and YX.
719 ref_scalar = p.sky_proj.pixel_to_sky(x=sx, y=sy)
720 result_xy = p.sky_proj.pixel_to_sky(XY(sx, sy))
721 result_yx = p.sky_proj.pixel_to_sky(YX(sy, sx))
722 np.testing.assert_allclose(result_xy.ra.rad, ref_scalar.ra.rad, atol=1e-12)
723 np.testing.assert_allclose(result_yx.ra.rad, ref_scalar.ra.rad, atol=1e-12)
725 # Array XY and YX.
726 ref_array = p.sky_proj.pixel_to_sky(x=xv, y=yv)
727 result_xy_arr = p.sky_proj.pixel_to_sky(XY(xv, yv))
728 result_yx_arr = p.sky_proj.pixel_to_sky(YX(yv, xv))
729 np.testing.assert_allclose(result_xy_arr.ra.rad, ref_array.ra.rad, atol=1e-12)
730 np.testing.assert_allclose(result_yx_arr.ra.rad, ref_array.ra.rad, atol=1e-12)
732 # TypeError on bad combinations.
733 with pytest.raises(TypeError):
734 p.sky_proj.pixel_to_sky(XY(sx, sy), x=sx)
735 with pytest.raises(TypeError):
736 p.sky_proj.pixel_to_sky(YX(sy, sx), y=sy)
737 with pytest.raises(TypeError):
738 p.sky_proj.pixel_to_sky()
741def _rotated_tan(rot_deg: float, *, crval2: float = 30.0, scale_y: float = 0.2) -> SkyProjection:
742 """Return a rotated TAN projection with given pixel scales."""
743 cx = (0.2 * u.arcsec).to_value(u.deg)
744 cy = (scale_y * u.arcsec).to_value(u.deg)
745 t = np.deg2rad(rot_deg)
746 header = {
747 "CTYPE1": "RA---TAN",
748 "CTYPE2": "DEC--TAN",
749 "CRPIX1": 50,
750 "CRPIX2": 100,
751 "CRVAL1": 45.0,
752 "CRVAL2": crval2,
753 "CD1_1": -cx * np.cos(t),
754 "CD1_2": cy * np.sin(t),
755 "CD2_1": -cx * np.sin(t),
756 "CD2_2": -cy * np.cos(t),
757 }
758 return SkyProjection.from_fits_wcs(astropy.wcs.WCS(header), GeneralFrame(unit=u.pix))
761def test_sky_projection_nominal_pixel_scale() -> None:
762 """_nominal_pixel_scale reports per sky axis [longitude, latitude].
764 Faithful KPG1_PXSCL port: the scale attaches to the sky axis, so a 90 deg
765 rotation swaps the returned [RA, Dec] scales. Great-circle distances keep
766 it correct near the poles.
767 """
768 bbox = Box.factory[0:200, 0:100]
770 # Unrotated anisotropic WCS: RA scale 0.2, Dec scale 0.3.
771 np.testing.assert_allclose(
772 _rotated_tan(0.0, scale_y=0.3)._nominal_pixel_scale(bbox), [0.2, 0.3], rtol=1e-3
773 )
774 # Rotated 90 deg: the sky-axis scales swap.
775 np.testing.assert_allclose(
776 _rotated_tan(90.0, scale_y=0.3)._nominal_pixel_scale(bbox), [0.3, 0.2], rtol=1e-3
777 )
778 # Reference pixel ~2 arcsec from the north pole: great-circle scale holds.
779 np.testing.assert_allclose(
780 _rotated_tan(30.0, crval2=89.9995)._nominal_pixel_scale(bbox), [0.2, 0.2], rtol=1e-3
781 )
784def test_sky_projection_describe() -> None:
785 """SkyProjection._describe reports pixels, pixel scale, and a corners
786 table.
787 """
788 rng = np.random.default_rng(43)
789 bbox = Box.factory[0:200, 0:100]
790 pixel_frame = GeneralFrame(unit=u.pix)
791 sky_projection = make_random_sky_projection(rng, pixel_frame, bbox)
793 # Without a bbox this projection still has pixel_bounds, so the report is
794 # characterized over those rather than falling back to the pixel origin.
795 assert sky_projection.pixel_bounds is not None
796 report = sky_projection.describe()
797 assert isinstance(report, Report)
798 assert report.type_name == "SkyProjection"
799 assert report.title == "ICRS coordinates"
800 # The title carries the sky frame, so it is not repeated as a field.
801 assert not any(f.label == "domain" for f in report.fields)
802 # The uninformative pixel_to_sky field is gone, so repr has no ARG fields.
803 assert not any(f.role is FieldRole.ARG for f in report.fields)
804 assert not any(f.label == "origin pixel" for f in report.fields)
805 assert any(f.label == "center pixel" for f in report.fields)
806 assert any(t.title == "Corners" for t in report.tables)
807 # Exactly one nominal-pixel-scale field is present (single or dual form).
808 scale_fields = [f for f in report.fields if f.label.startswith("Nominal pixel scale")]
809 assert len(scale_fields) == 1
811 # With a bbox: a Corners table and a center-pixel field for the box
812 # center, and still no origin-pixel fallback.
813 report = sky_projection.describe(bbox=bbox)
814 corners = next(t for t in report.tables if t.title == "Corners")
815 assert corners.columns == ["x", "y", "RA", "Dec", "RA (°)", "Dec (°)"]
816 assert len(corners.rows) == 4
817 # The first two columns carry the pixel coordinates of each corner. These
818 # are the box corners expanded by half a pixel, so that they bound the full
819 # area the image covers rather than the centers of the outermost pixels.
820 assert [(row[0], row[1]) for row in corners.rows] == [
821 ("-0.5", "-0.5"),
822 ("-0.5", "199.5"),
823 ("99.5", "199.5"),
824 ("99.5", "-0.5"),
825 ]
826 assert not any(f.label == "origin pixel" for f in report.fields)
827 center = next(f for f in report.fields if f.label == "center pixel")
828 assert center.role is FieldRole.DERIVED
829 assert center.value.startswith("(x=")
831 # FITS-WCS availability is reported (projection is FITS-representable).
832 fits_field = next(f for f in report.fields if f.label == "fits_wcs")
833 assert fits_field.value == "available"
835 # A projection cannot be rebuilt from a string, so repr is descriptive
836 # rather than an eval-ish call. It does not depend on a bbox and does not
837 # evaluate the mapping.
838 assert repr(sky_projection) == "<SkyProjection: GeneralFrame → ICRS>"
841def test_sky_projection_describe_pixel_scale_forms() -> None:
842 """The pixel-scale field collapses to one value when the axes agree."""
843 bbox = Box.factory[0:200, 0:100]
845 # Isotropic scale: a single "Nominal pixel scale" field to 0.01 arcsec.
846 report = _rotated_tan(0.0).describe(bbox=bbox)
847 scale = next(f for f in report.fields if f.label == "Nominal pixel scale")
848 assert scale.value == "0.20 arcsec"
849 assert not any(f.label == "Nominal pixel scales" for f in report.fields)
851 # Anisotropic scale: axes are named with the AST SkyFrame sky-axis labels.
852 report = _rotated_tan(0.0, scale_y=0.3).describe(bbox=bbox)
853 scales = next(f for f in report.fields if f.label == "Nominal pixel scales")
854 assert scales.value == "Right ascension 0.20 arcsec; Declination 0.30 arcsec"
855 assert not any(f.label == "Nominal pixel scale" for f in report.fields)
858def test_sky_projection_describe_fits_wcs_availability() -> None:
859 """fits_wcs is probed against a box, never blindly reported available."""
860 # No supplied bbox and no pixel bounds: representability cannot be checked.
861 projection = _rotated_tan(0.0)
862 assert projection.pixel_bounds is None
863 report = projection.describe()
864 fits_field = next(f for f in report.fields if f.label == "fits_wcs")
865 assert fits_field.value == "unknown"
867 # A projection with pixel bounds probes those bounds when no bbox is given.
868 bounded = SkyProjection.from_fits_wcs(
869 _rotated_tan(0.0).pixel_to_sky_transform.as_fits_wcs(Box.factory[0:32, 0:32]),
870 GeneralFrame(unit=u.pix),
871 pixel_bounds=Box.factory[0:32, 0:32],
872 )
873 report = bounded.describe()
874 fits_field = next(f for f in report.fields if f.label == "fits_wcs")
875 assert fits_field.value in ("available", "none")
878def test_sky_projection_astropy_view_repr_str() -> None:
879 """The Astropy view has informative str and repr, not the default object
880 form.
881 """
882 bbox = Box.factory[0:200, 0:100]
883 sky_projection = _rotated_tan(0.0)
885 bounded = sky_projection.as_astropy(bbox)
886 assert str(bounded) == "SkyProjectionAstropyView([y=0:200, x=0:100] → ICRS)"
887 bounded_repr = repr(bounded)
888 assert bounded_repr.startswith("SkyProjectionAstropyView\n")
889 assert "ICRS (ra, dec)" in bounded_repr
890 assert str(bbox.shape) in bounded_repr
891 # The reported pixel (0, 0) sky position matches the view's own transform.
892 ra, dec = bounded.pixel_to_world_values(0.0, 0.0)
893 ra_hms = Longitude(float(ra) * u.rad).to_string(unit=u.hour, sep="hms", pad=True, precision=1)
894 assert ra_hms in bounded_repr
896 unbounded = sky_projection.as_astropy()
897 assert str(unbounded) == "SkyProjectionAstropyView(unbounded → ICRS)"
898 # An unbounded view omits the array-shape line but keeps the reference.
899 assert "array shape" not in repr(unbounded)
900 assert "pixel (0, 0)" in repr(unbounded)
903def test_sky_projection_describe_origin_pixel_is_a_last_resort() -> None:
904 """The origin pixel appears only when no bounding box can be had.
906 Pixel (0, 0) is not a reference point of a projection, so it is reported
907 only when neither the caller nor the projection itself supplies a box to
908 characterize over.
909 """
910 rng = np.random.default_rng(43)
911 bbox = Box.factory[0:200, 0:100]
912 bounded = make_random_sky_projection(rng, GeneralFrame(unit=u.pix), bbox)
913 assert bounded.pixel_bounds is not None
914 # A box from either source displaces the origin fallback.
915 for report in (bounded.describe(), bounded.describe(bbox=bbox)):
916 assert not any(f.label == "origin pixel" for f in report.fields)
918 # With neither, the origin is all that is left, and the corners table and
919 # FITS-WCS probe drop out with it.
920 unbounded = _rotated_tan(0.0)
921 assert unbounded.pixel_bounds is None
922 report = unbounded.describe()
923 assert not any(f.label == "center pixel" for f in report.fields)
924 assert not any(t.title == "Corners" for t in report.tables)
925 assert next(f for f in report.fields if f.label == "fits_wcs").value == "unknown"
926 origin = next(f for f in report.fields if f.label == "origin pixel")
927 assert origin.role is FieldRole.DERIVED
928 assert origin.value.startswith("(x=0, y=0) →")
931def test_sky_projection_describe_origin_pixel_matches_transform() -> None:
932 """The origin-pixel field reports pixel (0, 0)'s actual sky position."""
933 sky_projection = _rotated_tan(0.0)
934 assert sky_projection.pixel_bounds is None
936 report = sky_projection.describe()
937 ref = next(f for f in report.fields if f.label == "origin pixel")
938 sky00 = sky_projection.pixel_to_sky(x=0, y=0)
939 # Sexagesimal (hms/dms) and labeled decimal degrees both appear.
940 assert sky00.ra.to_string(unit=u.hour, sep="hms", pad=True, precision=1) in ref.value
941 assert sky00.dec.to_string(sep="dms", pad=True, alwayssign=True, precision=0) in ref.value
942 assert f"RA {sky00.ra.deg:.6f}°" in ref.value
943 assert f"Dec {sky00.dec.deg:+.6f}°" in ref.value
946def test_frame_describe_preserves_pydantic_repr() -> None:
947 """Frames expose describe() while retaining pydantic's repr."""
948 frame = GeneralFrame(unit=u.pix)
949 report = frame.describe()
950 assert report.type_name == "GeneralFrame"
951 assert {f.label for f in report.fields} >= {"unit"}
952 # The mixin must not shadow pydantic's repr.
953 assert repr(frame).startswith("GeneralFrame(")
956def test_transform_describe() -> None:
957 """Transform._describe reports its frames and bounds."""
958 pixel_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.factory[:5, :4])
959 transform = Transform(pixel_frame, ICRS, astshim.UnitMap(2))
960 report = transform.describe()
961 assert isinstance(report, Report)
962 assert report.type_name == "Transform"
963 labels = {f.label for f in report.fields}
964 assert {"in_frame", "out_frame"} <= labels
965 assert any(f.label == "mapping" for f in report.fields)
968def test_sky_projection_describe_extent() -> None:
969 """The extent gives the box's angular size and its orientation.
971 The size is the great-circle distance across the area the box covers, and
972 the orientation is the position angle of the pixel y axis.
973 """
974 projection = _rotated_tan(40.0)
975 # 275 pixels at 0.2 arcsec/pixel spans 55 arcsec.
976 report = projection.describe(bbox=Box.factory[0:275, 0:275])
977 extent = next(f for f in report.fields if f.label == "Image extent")
978 assert extent.role is FieldRole.DERIVED
979 assert extent.value == "55 x 55 arcsec @ 140 deg E of N"
981 # The unit follows the scale of the box.
982 def size(pixels: int) -> str:
983 report = projection.describe(bbox=Box.factory[0:pixels, 0:pixels])
984 return next(f.value for f in report.fields if f.label == "Image extent")
986 assert "arcsec" in size(275)
987 assert "arcmin" in size(4000)
988 assert "deg" in size(60000)
990 # A long, thin box has no unit that suits both sides, so each gets its
991 # own; a square one names the shared unit once.
992 def extent_of(x_pixels: int, y_pixels: int) -> str:
993 report = projection.describe(bbox=Box.factory[0:y_pixels, 0:x_pixels])
994 return next(f.value for f in report.fields if f.label == "Image extent")
996 assert extent_of(10, 1000).startswith("2 arcsec x 3.33 arcmin ")
997 assert extent_of(1000, 10).startswith("3.33 arcmin x 2 arcsec ")
998 assert extent_of(275, 275).startswith("55 x 55 arcsec ")
1000 # No box, so no extent to report.
1001 assert projection.pixel_bounds is None
1002 assert not any(f.label == "Image extent" for f in projection.describe().fields)
1005def test_sky_projection_describe_extent_position_angle() -> None:
1006 """The reported angle is the position angle of the pixel y axis.
1008 ``_rotated_tan`` builds a CD matrix whose y axis lands at ``180 - rot``
1009 degrees East of North, which pins both the axis and the direction of
1010 increasing angle.
1011 """
1012 bbox = Box.factory[0:200, 0:100]
1013 for rot in (0.0, 30.0, 90.0, 270.0):
1014 report = _rotated_tan(rot).describe(bbox=bbox)
1015 value = next(f.value for f in report.fields if f.label == "Image extent")
1016 assert value.endswith(f"@ {(180 - rot) % 360:g} deg E of N"), (rot, value)
1019def test_sky_projection_describe_extent_is_measured_at_the_box() -> None:
1020 """The orientation is sampled at the box, not at the pixel origin.
1022 Meridians converge near the pole, so the position angle at a pixel far
1023 outside the box says nothing about the box: here the origin differs from
1024 the box by nearly 200 degrees.
1025 """
1026 projection = _rotated_tan(40.0, crval2=89.9995)
1027 bbox = Box.factory[0:200, 0:100]
1029 def position_angle_at(x: float, y: float) -> float:
1030 here = projection.pixel_to_sky(x=x, y=y)
1031 up = projection.pixel_to_sky(x=x, y=y + 1.0)
1032 return here.position_angle(up).to_value(u.deg) % 360.0
1034 at_origin = position_angle_at(0.0, 0.0)
1035 at_center = position_angle_at(bbox.x.center, bbox.y.center)
1036 assert abs(at_origin - at_center) > 100.0
1038 value = next(f.value for f in projection.describe(bbox=bbox).fields if f.label == "Image extent")
1039 assert value.endswith(f"@ {round(at_center) % 360} deg E of N")