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