Coverage for tests/test_transforms.py: 83%
529 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-14 07:48 +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)
53EXTERNAL_DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None)
56@pytest.fixture(scope="session")
57def legacy_camera() -> Any:
58 """Return a legacy Camera loaded from camera.fits.
60 Skips if TESTDATA_IMAGES_DIR is unset or lsst.afw.cameraGeom is
61 unavailable.
62 """
63 if EXTERNAL_DATA_DIR is None: 63 ↛ 65line 63 didn't jump to line 65 because the condition on line 63 was always true
64 pytest.skip("TESTDATA_IMAGES_DIR is not in the environment.")
65 try:
66 from lsst.afw.cameraGeom import Camera
67 except ImportError:
68 pytest.skip("'lsst.afw.cameraGeom' could not be imported.")
69 filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "camera.fits")
70 return Camera.readFits(filename)
73@pytest.fixture(scope="session")
74def legacy_detector_wcs() -> dict[str, Any]:
75 """Return WCS-related objects read from visit_image.fits.
77 Skips if TESTDATA_IMAGES_DIR is unset or lsst.afw.image is unavailable.
78 """
79 if EXTERNAL_DATA_DIR is None: 79 ↛ 81line 79 didn't jump to line 81 because the condition on line 79 was always true
80 pytest.skip("TESTDATA_IMAGES_DIR is not in the environment.")
81 try:
82 from lsst.afw.image import ExposureFitsReader
83 except ImportError:
84 pytest.skip("'lsst.afw.image' could not be imported.")
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):
165 SkyProjection(Transform(ICRS, ICRS, astshim.UnitMap(2)))
167 with pytest.raises(ValueError):
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_detector_wcs(legacy_detector_wcs: dict[str, Any]) -> None:
364 """Test the Transform/SkyProjection representation of a detector WCS."""
365 legacy_wcs = legacy_detector_wcs["legacy_wcs"]
366 wcs_bbox = legacy_detector_wcs["wcs_bbox"]
367 subimage_bbox = legacy_detector_wcs["subimage_bbox"]
368 detector_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=wcs_bbox)
369 sky_projection = SkyProjection.from_legacy(legacy_wcs, detector_frame)
370 assert sky_projection.fits_approximation is not None
371 compare_sky_projection_to_legacy_wcs(sky_projection, legacy_wcs, detector_frame, subimage_bbox)
372 # When we convert from a legacy SkyWcs, the internal AST Mapping needs
373 # to really be an AST FrameSet in order to be able to convert back.
374 assert "Begin FrameSet" in sky_projection.show()
375 compare_sky_projection_to_legacy_wcs(
376 sky_projection, sky_projection.to_legacy(), detector_frame, subimage_bbox
377 )
378 assert "Begin FrameSet" in sky_projection.fits_approximation.show()
379 compare_sky_projection_to_legacy_wcs(
380 sky_projection.fits_approximation,
381 sky_projection.fits_approximation.to_legacy(),
382 detector_frame,
383 subimage_bbox,
384 is_fits=True,
385 )
386 with RoundtripJson(sky_projection, "SkyProjection") as roundtrip:
387 pass
388 compare_sky_projection_to_legacy_wcs(roundtrip.result, legacy_wcs, detector_frame, subimage_bbox)
389 # The AST FrameSet-ness needs to propagate through serialization.
390 assert "Begin FrameSet" in roundtrip.result.show()
391 compare_sky_projection_to_legacy_wcs(
392 sky_projection, roundtrip.result.to_legacy(), detector_frame, subimage_bbox
393 )
394 with RoundtripJson(sky_projection.fits_approximation, "SkyProjection") as roundtrip:
395 pass
396 compare_sky_projection_to_legacy_wcs(
397 roundtrip.result,
398 legacy_wcs.getFitsApproximation(),
399 detector_frame,
400 subimage_bbox,
401 is_fits=True,
402 )
403 assert "Begin FrameSet" in roundtrip.result.show()
404 compare_sky_projection_to_legacy_wcs(
405 sky_projection.fits_approximation,
406 roundtrip.result.to_legacy(),
407 detector_frame,
408 subimage_bbox,
409 is_fits=True,
410 )
413@dataclasses.dataclass
414class FrameSetTestHolder:
415 """A top-level object that holds a CameraFrameSet and a transform
416 extracted from it, for testing archive pointers and frame set references.
417 """
419 frames: CameraFrameSet
420 pixels_to_fp: Transform[DetectorFrame, FocalPlaneFrame]
422 def serialize[P: pydantic.BaseModel](self, archive: OutputArchive[P]) -> FrameSetTestHolderModel[P]:
423 frames_model = archive.serialize_frame_set(
424 "frames", self.frames, self.frames.serialize, key=id(self.frames)
425 )
426 pixels_to_fp_model = archive.serialize_direct(
427 "pixels_to_fp", functools.partial(self.pixels_to_fp.serialize, use_frame_sets=True)
428 )
429 return FrameSetTestHolderModel[P](frames=frames_model, pixels_to_fp=pixels_to_fp_model)
431 @staticmethod
432 def _get_archive_tree_type[P: pydantic.BaseModel](
433 pointer_type: type[P],
434 ) -> type[FrameSetTestHolderModel[P]]:
435 return FrameSetTestHolderModel[pointer_type] # type: ignore
438class FrameSetTestHolderModel[P: pydantic.BaseModel](ArchiveTree):
439 """The serialization model for FrameSetTestHolder."""
441 SCHEMA_NAME: ClassVar[str] = "_test_frame_set_holder"
442 SCHEMA_VERSION: ClassVar[str] = "1.0.0"
443 MIN_READ_VERSION: ClassVar[int] = 1
444 PUBLIC_TYPE: ClassVar[type] = FrameSetTestHolder
446 frames: CameraFrameSetSerializationModel | P
447 pixels_to_fp: TransformSerializationModel[P]
449 def deserialize(self, archive: InputArchive[Any]) -> FrameSetTestHolder:
450 assert not isinstance(self.frames, CameraFrameSetSerializationModel), "Archive pointer expected."
451 frames = archive.deserialize_pointer(
452 self.frames, CameraFrameSetSerializationModel, CameraFrameSetSerializationModel.deserialize
453 )
454 pixels_to_fp = self.pixels_to_fp.deserialize(archive)
455 return FrameSetTestHolder(frames, pixels_to_fp)
458@dataclasses.dataclass
459class _BroadcastTestData:
460 """Shared inputs for broadcasting/scalar tests on Transform and
461 SkyProjection.
462 """
464 in_frame: DetectorFrame
465 out_frame: GeneralFrame
466 matrix: np.ndarray
467 scalar_x: float
468 scalar_y: float
469 xv: list[int]
470 yv: list[int]
471 sky_proj: SkyProjection[DetectorFrame]
474@pytest.fixture
475def broadcast_test_data() -> _BroadcastTestData:
476 """Return shared inputs for broadcasting/scalar tests."""
477 in_frame = DetectorFrame(instrument="Inst", visit=1, detector=1, bbox=Box.factory[0:20, 0:20])
478 return _BroadcastTestData(
479 in_frame=in_frame,
480 out_frame=GeneralFrame(unit=u.pix),
481 matrix=np.array([[2.0, 0.5], [-0.3, 1.5]]),
482 scalar_x=3.0,
483 scalar_y=7.0,
484 xv=[1, 2, 3],
485 yv=[4, 5, 6],
486 sky_proj=make_random_sky_projection(np.random.default_rng(42), in_frame, in_frame.bbox),
487 )
490def test_apply_forward_scalar(broadcast_test_data: _BroadcastTestData) -> None:
491 """Verify that apply_forward and apply_inverse accept scalar x/y and return
492 scalar floats, and that the _q variants accept scalar Quantity inputs.
493 """
494 t = broadcast_test_data.sky_proj.pixel_to_sky_transform
495 # apply_forward with Python float scalars should return XY of floats.
496 result_fwd = t.apply_forward(x=broadcast_test_data.scalar_x, y=broadcast_test_data.scalar_y)
497 assert type(result_fwd.x) is float
498 assert type(result_fwd.y) is float
499 # Values must match the corresponding single-element array call.
500 ref_fwd = t.apply_forward(
501 x=np.array([broadcast_test_data.scalar_x]), y=np.array([broadcast_test_data.scalar_y])
502 )
503 assert result_fwd.x == ref_fwd.x[0]
504 assert result_fwd.y == ref_fwd.y[0]
505 # apply_inverse round-trips back to the original scalars.
506 result_inv = t.apply_inverse(x=result_fwd.x, y=result_fwd.y)
507 assert type(result_inv.x) is float
508 assert type(result_inv.y) is float
509 np.testing.assert_allclose(result_inv.x, broadcast_test_data.scalar_x, atol=1e-12)
510 np.testing.assert_allclose(result_inv.y, broadcast_test_data.scalar_y, atol=1e-12)
511 # apply_forward_q / apply_inverse_q with scalar Quantity inputs.
512 x_q = broadcast_test_data.scalar_x * t.in_frame.unit
513 y_q = broadcast_test_data.scalar_y * t.in_frame.unit
514 result_fwd_q = t.apply_forward_q(x=x_q, y=y_q)
515 assert result_fwd_q.x.shape == ()
516 assert result_fwd_q.y.shape == ()
517 np.testing.assert_allclose(result_fwd_q.x.to_value(t.out_frame.unit), result_fwd.x, atol=1e-12)
518 result_inv_q = t.apply_inverse_q(x=result_fwd_q.x, y=result_fwd_q.y)
519 assert result_inv_q.x.shape == ()
520 np.testing.assert_allclose(
521 result_inv_q.x.to_value(t.in_frame.unit), broadcast_test_data.scalar_x, atol=1e-12
522 )
525def test_apply_array_like_and_integer_input(broadcast_test_data: _BroadcastTestData) -> None:
526 """Verify that apply_forward accepts Python lists and integer-dtype
527 arrays, returning float64 ndarray results consistent with float64
528 array input.
529 """
530 t = broadcast_test_data.sky_proj.pixel_to_sky_transform
531 # Python list input should return an ndarray.
532 result_list = t.apply_forward(x=broadcast_test_data.xv, y=broadcast_test_data.yv)
533 assert isinstance(result_list.x, np.ndarray)
534 assert isinstance(result_list.y, np.ndarray)
535 ref = t.apply_forward(x=np.array(broadcast_test_data.xv), y=np.array(broadcast_test_data.yv))
536 np.testing.assert_array_equal(result_list.x, ref.x)
537 np.testing.assert_array_equal(result_list.y, ref.y)
538 # Integer dtype arrays should not raise and should return float64.
539 xi = np.array(broadcast_test_data.xv, dtype=np.int32)
540 yi = np.array(broadcast_test_data.yv, dtype=np.int32)
541 result_int = t.apply_forward(x=xi, y=yi)
542 assert result_int.x.dtype == np.float64
543 assert result_int.y.dtype == np.float64
544 np.testing.assert_array_equal(result_int.x, ref.x)
545 np.testing.assert_array_equal(result_int.y, ref.y)
548def test_apply_broadcast(broadcast_test_data: _BroadcastTestData) -> None:
549 """Verify that apply_forward and apply_inverse broadcast x and y like
550 a NumPy ufunc, in both 1-D and 2-D cases.
551 """
552 t = broadcast_test_data.sky_proj.pixel_to_sky_transform
553 xv = np.array(broadcast_test_data.xv)
554 yv = np.array(broadcast_test_data.yv + [7])
555 # 1-D broadcast: array x, scalar y.
556 result_1d = t.apply_forward(x=xv, y=broadcast_test_data.scalar_y)
557 assert isinstance(result_1d.x, np.ndarray)
558 assert result_1d.x.shape == xv.shape
559 ref_1d = t.apply_forward(x=xv, y=np.full_like(xv, broadcast_test_data.scalar_y))
560 np.testing.assert_array_equal(result_1d.x, ref_1d.x)
561 np.testing.assert_array_equal(result_1d.y, ref_1d.y)
562 # 2-D broadcast: column x (M,1) × row y (1,N) -> (M,N).
563 x2d = xv[:, np.newaxis] # shape (3, 1)
564 y2d = yv[np.newaxis, :] # shape (1, 4)
565 result_2d = t.apply_forward(x=x2d, y=y2d)
566 assert result_2d.x.shape == (3, 4)
567 assert result_2d.y.shape == (3, 4)
568 # Values must match the fully expanded meshgrid call.
569 xmesh, ymesh = np.meshgrid(xv, yv, indexing="ij")
570 ref_2d = t.apply_forward(x=xmesh, y=ymesh)
571 np.testing.assert_array_equal(result_2d.x, ref_2d.x)
572 np.testing.assert_array_equal(result_2d.y, ref_2d.y)
573 # apply_inverse also broadcasts.
574 result_inv_2d = t.apply_inverse(x=result_2d.x, y=result_2d.y)
575 assert result_inv_2d.x.shape == (3, 4)
576 np.testing.assert_allclose(result_inv_2d.x, xmesh, atol=1e-12)
577 np.testing.assert_allclose(result_inv_2d.y, ymesh, atol=1e-12)
580def test_sky_projection_broadcast(broadcast_test_data: _BroadcastTestData) -> None:
581 """Verify that SkyProjection.pixel_to_sky, sky_to_pixel, and the
582 Astropy view broadcast x and y like a NumPy ufunc.
583 """
584 p = broadcast_test_data
585 xv = np.array(p.xv)
586 yv = np.array(p.yv + [7])
587 # 1-D broadcast: array x, scalar y.
588 sc_1d = p.sky_proj.pixel_to_sky(x=xv, y=p.scalar_y)
589 assert sc_1d.shape == xv.shape
590 ref_1d = p.sky_proj.pixel_to_sky(x=xv, y=np.full_like(xv, p.scalar_y))
591 np.testing.assert_allclose(sc_1d.ra.rad, ref_1d.ra.rad, atol=1e-12)
592 np.testing.assert_allclose(sc_1d.dec.rad, ref_1d.dec.rad, atol=1e-12)
593 # 2-D broadcast: column x (M,1) × row y (1,N) -> (M,N).
594 x2d = xv[:, np.newaxis] # shape (3, 1)
595 y2d = yv[np.newaxis, :] # shape (1, 4)
596 sc_2d = p.sky_proj.pixel_to_sky(x=x2d, y=y2d)
597 assert sc_2d.shape == (3, 4)
598 xmesh, ymesh = np.meshgrid(xv, yv, indexing="ij")
599 ref_2d = p.sky_proj.pixel_to_sky(x=xmesh, y=ymesh)
600 np.testing.assert_allclose(sc_2d.ra.rad, ref_2d.ra.rad, atol=1e-12)
601 np.testing.assert_allclose(sc_2d.dec.rad, ref_2d.dec.rad, atol=1e-12)
602 # sky_to_pixel round-trips back to the original grid.
603 pix_2d = p.sky_proj.sky_to_pixel(sc_2d)
604 assert pix_2d.x.shape == (3, 4)
605 np.testing.assert_allclose(pix_2d.x, xmesh, atol=1e-9)
606 np.testing.assert_allclose(pix_2d.y, ymesh, atol=1e-9)
607 # SkyProjectionAstropyView.pixel_to_world_values also broadcasts.
608 view = p.sky_proj.as_astropy()
609 world_2d = view.pixel_to_world_values(x2d, y2d)
610 assert world_2d[0].shape == (3, 4)
611 assert world_2d[1].shape == (3, 4)
612 np.testing.assert_allclose(world_2d[0], ref_2d.ra.rad, atol=1e-12)
613 np.testing.assert_allclose(world_2d[1], ref_2d.dec.rad, atol=1e-12)
614 # SkyProjectionAstropyView.world_to_pixel_values also broadcasts.
615 ra_2d = ref_2d.ra.rad[:, np.newaxis, :] # (3, 1, 4) — over-broadcast to check
616 dec_2d = ref_2d.dec.rad[np.newaxis, :, :] # (1, 3, 4)
617 pix_world = view.world_to_pixel_values(ra_2d, dec_2d)
618 assert pix_world[0].shape == (3, 3, 4)
621def test_apply_xy_yx(broadcast_test_data: _BroadcastTestData) -> None:
622 """Verify that apply_forward, apply_inverse, and the _q variants accept
623 XY and YX positional arguments, producing results identical to the
624 equivalent x=/y= keyword calls.
625 """
626 p = broadcast_test_data
627 t = p.sky_proj.pixel_to_sky_transform
628 sx, sy = p.scalar_x, p.scalar_y
629 xv, yv = np.array(p.xv, dtype=float), np.array(p.yv, dtype=float)
631 # --- apply_forward: scalar ---
632 ref_fwd = t.apply_forward(x=sx, y=sy)
633 assert t.apply_forward(XY(sx, sy)) == ref_fwd
634 assert t.apply_forward(YX(sy, sx)) == ref_fwd
636 # --- apply_forward: array ---
637 ref_fwd_arr = t.apply_forward(x=xv, y=yv)
638 np.testing.assert_array_equal(t.apply_forward(XY(xv, yv)).x, ref_fwd_arr.x)
639 np.testing.assert_array_equal(t.apply_forward(YX(yv, xv)).x, ref_fwd_arr.x)
641 # --- apply_inverse: scalar ---
642 ref_inv = t.apply_inverse(x=ref_fwd.x, y=ref_fwd.y)
643 assert t.apply_inverse(XY(ref_fwd.x, ref_fwd.y)) == ref_inv
644 assert t.apply_inverse(YX(ref_fwd.y, ref_fwd.x)) == ref_inv
646 # --- apply_forward_q: scalar ---
647 x_q = sx * t.in_frame.unit
648 y_q = sy * t.in_frame.unit
649 ref_fwd_q = t.apply_forward_q(x=x_q, y=y_q)
650 result_q = t.apply_forward_q(XY(x_q, y_q))
651 np.testing.assert_allclose(result_q.x.value, ref_fwd_q.x.value, atol=1e-12)
652 result_q_yx = t.apply_forward_q(YX(y_q, x_q))
653 np.testing.assert_allclose(result_q_yx.x.value, ref_fwd_q.x.value, atol=1e-12)
655 # --- apply_inverse_q: scalar ---
656 ref_inv_q = t.apply_inverse_q(x=ref_fwd_q.x, y=ref_fwd_q.y)
657 result_inv_q = t.apply_inverse_q(XY(ref_fwd_q.x, ref_fwd_q.y))
658 np.testing.assert_allclose(result_inv_q.x.value, ref_inv_q.x.value, atol=1e-12)
660 # --- TypeError on bad combinations ---
661 with pytest.raises(TypeError):
662 t.apply_forward(XY(sx, sy), x=sx)
663 with pytest.raises(TypeError):
664 t.apply_forward(YX(sy, sx), y=sy)
665 with pytest.raises(TypeError):
666 t.apply_forward()
669def test_pixel_to_sky_xy_yx(broadcast_test_data: _BroadcastTestData) -> None:
670 """Verify that SkyProjection.pixel_to_sky accepts XY and YX positional
671 arguments, producing results identical to the x=/y= keyword form.
672 """
673 p = broadcast_test_data
674 sx, sy = p.scalar_x, p.scalar_y
675 xv, yv = np.array(p.xv, dtype=float), np.array(p.yv, dtype=float)
677 # Scalar XY and YX.
678 ref_scalar = p.sky_proj.pixel_to_sky(x=sx, y=sy)
679 result_xy = p.sky_proj.pixel_to_sky(XY(sx, sy))
680 result_yx = p.sky_proj.pixel_to_sky(YX(sy, sx))
681 np.testing.assert_allclose(result_xy.ra.rad, ref_scalar.ra.rad, atol=1e-12)
682 np.testing.assert_allclose(result_yx.ra.rad, ref_scalar.ra.rad, atol=1e-12)
684 # Array XY and YX.
685 ref_array = p.sky_proj.pixel_to_sky(x=xv, y=yv)
686 result_xy_arr = p.sky_proj.pixel_to_sky(XY(xv, yv))
687 result_yx_arr = p.sky_proj.pixel_to_sky(YX(yv, xv))
688 np.testing.assert_allclose(result_xy_arr.ra.rad, ref_array.ra.rad, atol=1e-12)
689 np.testing.assert_allclose(result_yx_arr.ra.rad, ref_array.ra.rad, atol=1e-12)
691 # TypeError on bad combinations.
692 with pytest.raises(TypeError):
693 p.sky_proj.pixel_to_sky(XY(sx, sy), x=sx)
694 with pytest.raises(TypeError):
695 p.sky_proj.pixel_to_sky(YX(sy, sx), y=sy)
696 with pytest.raises(TypeError):
697 p.sky_proj.pixel_to_sky()
700def _rotated_tan(rot_deg: float, *, crval2: float = 30.0, scale_y: float = 0.2) -> SkyProjection:
701 """Return a rotated TAN projection with given pixel scales."""
702 cx = (0.2 * u.arcsec).to_value(u.deg)
703 cy = (scale_y * u.arcsec).to_value(u.deg)
704 t = np.deg2rad(rot_deg)
705 header = {
706 "CTYPE1": "RA---TAN",
707 "CTYPE2": "DEC--TAN",
708 "CRPIX1": 50,
709 "CRPIX2": 100,
710 "CRVAL1": 45.0,
711 "CRVAL2": crval2,
712 "CD1_1": -cx * np.cos(t),
713 "CD1_2": cy * np.sin(t),
714 "CD2_1": -cx * np.sin(t),
715 "CD2_2": -cy * np.cos(t),
716 }
717 return SkyProjection.from_fits_wcs(astropy.wcs.WCS(header), GeneralFrame(unit=u.pix))
720def test_sky_projection_nominal_pixel_scale() -> None:
721 """_nominal_pixel_scale reports per sky axis [longitude, latitude].
723 Faithful KPG1_PXSCL port: the scale attaches to the sky axis, so a 90 deg
724 rotation swaps the returned [RA, Dec] scales. Great-circle distances keep
725 it correct near the poles.
726 """
727 bbox = Box.factory[0:200, 0:100]
729 # Unrotated anisotropic WCS: RA scale 0.2, Dec scale 0.3.
730 np.testing.assert_allclose(
731 _rotated_tan(0.0, scale_y=0.3)._nominal_pixel_scale(bbox), [0.2, 0.3], rtol=1e-3
732 )
733 # Rotated 90 deg: the sky-axis scales swap.
734 np.testing.assert_allclose(
735 _rotated_tan(90.0, scale_y=0.3)._nominal_pixel_scale(bbox), [0.3, 0.2], rtol=1e-3
736 )
737 # Reference pixel ~2 arcsec from the north pole: great-circle scale holds.
738 np.testing.assert_allclose(
739 _rotated_tan(30.0, crval2=89.9995)._nominal_pixel_scale(bbox), [0.2, 0.2], rtol=1e-3
740 )
743def test_sky_projection_describe() -> None:
744 """SkyProjection._describe reports pixels, pixel scale, and a corners
745 table.
746 """
747 rng = np.random.default_rng(43)
748 bbox = Box.factory[0:200, 0:100]
749 pixel_frame = GeneralFrame(unit=u.pix)
750 sky_projection = make_random_sky_projection(rng, pixel_frame, bbox)
752 # Without a bbox this projection still has pixel_bounds, so the report is
753 # characterized over those rather than falling back to the pixel origin.
754 assert sky_projection.pixel_bounds is not None
755 report = sky_projection.describe()
756 assert isinstance(report, Report)
757 assert report.type_name == "SkyProjection"
758 assert report.title == "ICRS coordinates"
759 # The title carries the sky frame, so it is not repeated as a field.
760 assert not any(f.label == "domain" for f in report.fields)
761 # The uninformative pixel_to_sky field is gone, so repr has no ARG fields.
762 assert not any(f.role is FieldRole.ARG for f in report.fields)
763 assert not any(f.label == "origin pixel" for f in report.fields)
764 assert any(f.label == "center pixel" for f in report.fields)
765 assert any(t.title == "Corners" for t in report.tables)
766 # Exactly one nominal-pixel-scale field is present (single or dual form).
767 scale_fields = [f for f in report.fields if f.label.startswith("Nominal pixel scale")]
768 assert len(scale_fields) == 1
770 # With a bbox: a Corners table and a center-pixel field for the box
771 # center, and still no origin-pixel fallback.
772 report = sky_projection.describe(bbox=bbox)
773 corners = next(t for t in report.tables if t.title == "Corners")
774 assert corners.columns == ["x", "y", "RA", "Dec", "RA (°)", "Dec (°)"]
775 assert len(corners.rows) == 4
776 # The first two columns carry the pixel coordinates of each corner. These
777 # are the box corners expanded by half a pixel, so that they bound the full
778 # area the image covers rather than the centers of the outermost pixels.
779 assert [(row[0], row[1]) for row in corners.rows] == [
780 ("-0.5", "-0.5"),
781 ("-0.5", "199.5"),
782 ("99.5", "199.5"),
783 ("99.5", "-0.5"),
784 ]
785 assert not any(f.label == "origin pixel" for f in report.fields)
786 center = next(f for f in report.fields if f.label == "center pixel")
787 assert center.role is FieldRole.DERIVED
788 assert center.value.startswith("(x=")
790 # FITS-WCS availability is reported (projection is FITS-representable).
791 fits_field = next(f for f in report.fields if f.label == "fits_wcs")
792 assert fits_field.value == "available"
794 # A projection cannot be rebuilt from a string, so repr is descriptive
795 # rather than an eval-ish call. It does not depend on a bbox and does not
796 # evaluate the mapping.
797 assert repr(sky_projection) == "<SkyProjection: GeneralFrame → ICRS>"
800def test_sky_projection_describe_pixel_scale_forms() -> None:
801 """The pixel-scale field collapses to one value when the axes agree."""
802 bbox = Box.factory[0:200, 0:100]
804 # Isotropic scale: a single "Nominal pixel scale" field to 0.01 arcsec.
805 report = _rotated_tan(0.0).describe(bbox=bbox)
806 scale = next(f for f in report.fields if f.label == "Nominal pixel scale")
807 assert scale.value == "0.20 arcsec"
808 assert not any(f.label == "Nominal pixel scales" for f in report.fields)
810 # Anisotropic scale: axes are named with the AST SkyFrame sky-axis labels.
811 report = _rotated_tan(0.0, scale_y=0.3).describe(bbox=bbox)
812 scales = next(f for f in report.fields if f.label == "Nominal pixel scales")
813 assert scales.value == "Right ascension 0.20 arcsec; Declination 0.30 arcsec"
814 assert not any(f.label == "Nominal pixel scale" for f in report.fields)
817def test_sky_projection_describe_fits_wcs_availability() -> None:
818 """fits_wcs is probed against a box, never blindly reported available."""
819 # No supplied bbox and no pixel bounds: representability cannot be checked.
820 projection = _rotated_tan(0.0)
821 assert projection.pixel_bounds is None
822 report = projection.describe()
823 fits_field = next(f for f in report.fields if f.label == "fits_wcs")
824 assert fits_field.value == "unknown"
826 # A projection with pixel bounds probes those bounds when no bbox is given.
827 bounded = SkyProjection.from_fits_wcs(
828 _rotated_tan(0.0).pixel_to_sky_transform.as_fits_wcs(Box.factory[0:32, 0:32]),
829 GeneralFrame(unit=u.pix),
830 pixel_bounds=Box.factory[0:32, 0:32],
831 )
832 report = bounded.describe()
833 fits_field = next(f for f in report.fields if f.label == "fits_wcs")
834 assert fits_field.value in ("available", "none")
837def test_sky_projection_astropy_view_repr_str() -> None:
838 """The Astropy view has informative str and repr, not the default object
839 form.
840 """
841 bbox = Box.factory[0:200, 0:100]
842 sky_projection = _rotated_tan(0.0)
844 bounded = sky_projection.as_astropy(bbox)
845 assert str(bounded) == "SkyProjectionAstropyView([y=0:200, x=0:100] → ICRS)"
846 bounded_repr = repr(bounded)
847 assert bounded_repr.startswith("SkyProjectionAstropyView\n")
848 assert "ICRS (ra, dec)" in bounded_repr
849 assert str(bbox.shape) in bounded_repr
850 # The reported pixel (0, 0) sky position matches the view's own transform.
851 ra, dec = bounded.pixel_to_world_values(0.0, 0.0)
852 ra_hms = Longitude(float(ra) * u.rad).to_string(unit=u.hour, sep="hms", pad=True, precision=1)
853 assert ra_hms in bounded_repr
855 unbounded = sky_projection.as_astropy()
856 assert str(unbounded) == "SkyProjectionAstropyView(unbounded → ICRS)"
857 # An unbounded view omits the array-shape line but keeps the reference.
858 assert "array shape" not in repr(unbounded)
859 assert "pixel (0, 0)" in repr(unbounded)
862def test_sky_projection_describe_origin_pixel_is_a_last_resort() -> None:
863 """The origin pixel appears only when no bounding box can be had.
865 Pixel (0, 0) is not a reference point of a projection, so it is reported
866 only when neither the caller nor the projection itself supplies a box to
867 characterize over.
868 """
869 rng = np.random.default_rng(43)
870 bbox = Box.factory[0:200, 0:100]
871 bounded = make_random_sky_projection(rng, GeneralFrame(unit=u.pix), bbox)
872 assert bounded.pixel_bounds is not None
873 # A box from either source displaces the origin fallback.
874 for report in (bounded.describe(), bounded.describe(bbox=bbox)):
875 assert not any(f.label == "origin pixel" for f in report.fields)
877 # With neither, the origin is all that is left, and the corners table and
878 # FITS-WCS probe drop out with it.
879 unbounded = _rotated_tan(0.0)
880 assert unbounded.pixel_bounds is None
881 report = unbounded.describe()
882 assert not any(f.label == "center pixel" for f in report.fields)
883 assert not any(t.title == "Corners" for t in report.tables)
884 assert next(f for f in report.fields if f.label == "fits_wcs").value == "unknown"
885 origin = next(f for f in report.fields if f.label == "origin pixel")
886 assert origin.role is FieldRole.DERIVED
887 assert origin.value.startswith("(x=0, y=0) →")
890def test_sky_projection_describe_origin_pixel_matches_transform() -> None:
891 """The origin-pixel field reports pixel (0, 0)'s actual sky position."""
892 sky_projection = _rotated_tan(0.0)
893 assert sky_projection.pixel_bounds is None
895 report = sky_projection.describe()
896 ref = next(f for f in report.fields if f.label == "origin pixel")
897 sky00 = sky_projection.pixel_to_sky(x=0, y=0)
898 # Sexagesimal (hms/dms) and labeled decimal degrees both appear.
899 assert sky00.ra.to_string(unit=u.hour, sep="hms", pad=True, precision=1) in ref.value
900 assert sky00.dec.to_string(sep="dms", pad=True, alwayssign=True, precision=0) in ref.value
901 assert f"RA {sky00.ra.deg:.6f}°" in ref.value
902 assert f"Dec {sky00.dec.deg:+.6f}°" in ref.value
905def test_frame_describe_preserves_pydantic_repr() -> None:
906 """Frames expose describe() while retaining pydantic's repr."""
907 frame = GeneralFrame(unit=u.pix)
908 report = frame.describe()
909 assert report.type_name == "GeneralFrame"
910 assert {f.label for f in report.fields} >= {"unit"}
911 # The mixin must not shadow pydantic's repr.
912 assert repr(frame).startswith("GeneralFrame(")
915def test_transform_describe() -> None:
916 """Transform._describe reports its frames and bounds."""
917 pixel_frame = DetectorFrame(**DP2_VISIT_DETECTOR_DATA_ID, bbox=Box.factory[:5, :4])
918 transform = Transform(pixel_frame, ICRS, astshim.UnitMap(2))
919 report = transform.describe()
920 assert isinstance(report, Report)
921 assert report.type_name == "Transform"
922 labels = {f.label for f in report.fields}
923 assert {"in_frame", "out_frame"} <= labels
924 assert any(f.label == "mapping" for f in report.fields)
927def test_sky_projection_describe_extent() -> None:
928 """The extent gives the box's angular size and its orientation.
930 The size is the great-circle distance across the area the box covers, and
931 the orientation is the position angle of the pixel y axis.
932 """
933 projection = _rotated_tan(40.0)
934 # 275 pixels at 0.2 arcsec/pixel spans 55 arcsec.
935 report = projection.describe(bbox=Box.factory[0:275, 0:275])
936 extent = next(f for f in report.fields if f.label == "Image extent")
937 assert extent.role is FieldRole.DERIVED
938 assert extent.value == "55 x 55 arcsec @ 140 deg E of N"
940 # The unit follows the scale of the box.
941 def size(pixels: int) -> str:
942 report = projection.describe(bbox=Box.factory[0:pixels, 0:pixels])
943 return next(f.value for f in report.fields if f.label == "Image extent")
945 assert "arcsec" in size(275)
946 assert "arcmin" in size(4000)
947 assert "deg" in size(60000)
949 # A long, thin box has no unit that suits both sides, so each gets its
950 # own; a square one names the shared unit once.
951 def extent_of(x_pixels: int, y_pixels: int) -> str:
952 report = projection.describe(bbox=Box.factory[0:y_pixels, 0:x_pixels])
953 return next(f.value for f in report.fields if f.label == "Image extent")
955 assert extent_of(10, 1000).startswith("2 arcsec x 3.33 arcmin ")
956 assert extent_of(1000, 10).startswith("3.33 arcmin x 2 arcsec ")
957 assert extent_of(275, 275).startswith("55 x 55 arcsec ")
959 # No box, so no extent to report.
960 assert projection.pixel_bounds is None
961 assert not any(f.label == "Image extent" for f in projection.describe().fields)
964def test_sky_projection_describe_extent_position_angle() -> None:
965 """The reported angle is the position angle of the pixel y axis.
967 ``_rotated_tan`` builds a CD matrix whose y axis lands at ``180 - rot``
968 degrees East of North, which pins both the axis and the direction of
969 increasing angle.
970 """
971 bbox = Box.factory[0:200, 0:100]
972 for rot in (0.0, 30.0, 90.0, 270.0):
973 report = _rotated_tan(rot).describe(bbox=bbox)
974 value = next(f.value for f in report.fields if f.label == "Image extent")
975 assert value.endswith(f"@ {(180 - rot) % 360:g} deg E of N"), (rot, value)
978def test_sky_projection_describe_extent_is_measured_at_the_box() -> None:
979 """The orientation is sampled at the box, not at the pixel origin.
981 Meridians converge near the pole, so the position angle at a pixel far
982 outside the box says nothing about the box: here the origin differs from
983 the box by nearly 200 degrees.
984 """
985 projection = _rotated_tan(40.0, crval2=89.9995)
986 bbox = Box.factory[0:200, 0:100]
988 def position_angle_at(x: float, y: float) -> float:
989 here = projection.pixel_to_sky(x=x, y=y)
990 up = projection.pixel_to_sky(x=x, y=y + 1.0)
991 return here.position_angle(up).to_value(u.deg) % 360.0
993 at_origin = position_angle_at(0.0, 0.0)
994 at_center = position_angle_at(bbox.x.center, bbox.y.center)
995 assert abs(at_origin - at_center) > 100.0
997 value = next(f.value for f in projection.describe(bbox=bbox).fields if f.label == "Image extent")
998 assert value.endswith(f"@ {round(at_center) % 360} deg E of N")