Coverage for tests/test_geom.py: 100%
398 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-26 09:32 +0000
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-26 09:32 +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 pickle
15from typing import get_overloads, get_type_hints
17import astropy.units as u
18import astropy.wcs
19import numpy as np
20import pydantic
21import pytest
22from astropy.coordinates import Angle, SkyCoord
24from lsst.images import (
25 XY,
26 YX,
27 Bounds,
28 Box,
29 GeneralFrame,
30 Interval,
31 NoOverlapError,
32 SerializableXY,
33 SerializableYX,
34 SkyProjection,
35)
36from lsst.images.tests import assert_close, check_bounds_contains_broadcasting
39class IntervalModel(pydantic.BaseModel):
40 """Test Pydantic model with an Interval."""
42 interval1: Interval
43 interval2: Interval
46class BoxModel(pydantic.BaseModel):
47 """Test Pydantic model with a box."""
49 box: Box
52class PairModel(pydantic.BaseModel):
53 """Test Pydantic model with XY and YX pairs."""
55 xy: SerializableXY[float]
56 yx: SerializableYX[int]
57 bare: XY
60def _assert_type_hints_resolve(obj: object) -> None:
61 """Assert that runtime annotations and overload annotations resolve."""
62 get_type_hints(obj)
63 for overload in get_overloads(obj):
64 get_type_hints(overload)
67def test_yx() -> None:
68 """Test YX."""
69 yx = YX(5, 7)
70 assert yx == (5, 7)
71 assert yx.y == 5
72 assert yx.x == 7
74 def _plus_one(v: int) -> int:
75 return v + 1
77 new = yx.map(_plus_one)
78 assert new == (6, 8)
80 xy = yx.xy
81 assert xy == (7, 5)
84def test_xy() -> None:
85 """Test XY."""
86 xy = XY(5, 7)
87 assert xy == (5, 7)
88 assert xy.y == 7
89 assert xy.x == 5
91 def _plus_one(v: int) -> int:
92 return v + 1
94 new = xy.map(_plus_one)
95 assert new == (6, 8)
97 yx = xy.yx
98 assert yx == (7, 5)
101def test_interval_constructor() -> None:
102 """Test simple Interval construction and arithmetic."""
103 i = Interval(start=1, stop=10)
104 assert i.start == 1
105 assert i.stop == 10
106 assert i.min == 1
107 assert i.max == 9
108 assert i.size == 9
109 assert i.center == 5.0
110 assert Interval(1, 10) == i
112 assert str(i) == "1:10"
114 shifted = i + 10
115 assert shifted.start == 11
116 assert shifted.stop == 20
118 shifted = i - 10
119 assert shifted.start == -9
120 assert shifted.stop == 0
122 assert shifted == i - 10
123 assert shifted != i
125 sized = Interval.from_size(10)
126 assert sized == Interval(0, 10)
127 sized = Interval.from_size(10, 5)
128 assert sized == Interval(5, 15)
130 h = Interval.hull(2, -1, 3, 6)
131 assert h == Interval(start=-1, stop=7)
132 h2 = Interval.hull(h, 3, 40, Interval(start=-10, stop=2))
133 assert h2 == Interval(start=-10, stop=41)
136def test_interval_contains() -> None:
137 """Test containment."""
138 i = Interval(start=1, stop=10)
139 assert 5 in i
140 assert 10 not in i
142 i2 = Interval(start=2, stop=4)
143 assert i.contains(i2)
145 i3 = Interval(start=-1, stop=5)
146 assert not i.contains(i3)
148 assert i3.contains(2.5)
150 containment = i3.contains(np.array([-2, -1, 0, 1, 2, 3, 4, 5, 6]))
151 assert list(containment) == [False, True, True, True, True, True, True, False, False]
153 inter = i2.intersection(i)
154 assert inter == Interval(start=2, stop=4), f"Intersection of {i2} with {i}"
155 with pytest.raises(NoOverlapError):
156 i.intersection(Interval.factory[20:30])
157 assert i != []
160def test_interval_contains_broadcasting() -> None:
161 """Test that Interval.contains accepts array-like inputs and broadcasts."""
162 interval = Interval(start=2, stop=8)
163 # One point outside each end, one on each boundary, one in the interior.
164 xs = np.array(
165 [
166 interval.start - 1,
167 interval.start,
168 (interval.start + interval.stop) // 2,
169 interval.stop - 1,
170 interval.stop,
171 ]
172 )
173 expected_1d = np.array([interval.contains(int(x)) for x in xs])
174 # 1-D ndarray.
175 np.testing.assert_array_equal(interval.contains(xs), expected_1d)
176 # list input (array-like).
177 np.testing.assert_array_equal(interval.contains(xs.tolist()), expected_1d)
178 # 2-D array: same values reshaped to (len, 1).
179 np.testing.assert_array_equal(interval.contains(xs.reshape(-1, 1)), expected_1d.reshape(-1, 1))
182def test_runtime_type_hints_resolve() -> None:
183 """Public geometry annotations can be resolved at runtime."""
184 _assert_type_hints_resolve(Interval.contains)
185 _assert_type_hints_resolve(Box.contains)
186 _assert_type_hints_resolve(Bounds.contains)
189def test_interval_slice() -> None:
190 """Test Interval construction and subsetting with slice notation."""
191 i = Interval.factory[3:20]
192 assert i.start == 3
193 assert i.stop == 20
194 assert i.absolute[::] == i
195 assert i.local[::] == i
197 subset = i.absolute[5:10]
198 assert subset.start == 5
199 assert subset.stop == 10
201 subset = i.local[5:10]
202 assert subset.start == 8
203 assert subset.stop == 13
205 subset = i.absolute[:10]
206 assert subset.start == 3
207 assert subset.stop == 10
209 subset = i.local[:10]
210 assert subset.start == 3
211 assert subset.stop == 13
213 subset = i.absolute[10:]
214 assert subset.start == 10
215 assert subset.stop == 20
217 subset = i.local[10:]
218 assert subset.start == 13
219 assert subset.stop == 20
221 subset = i.local[3:-2]
222 assert subset.start == 6
223 assert subset.stop == 18
225 subset = i.local[-5:-2]
226 assert subset.start == 15
227 assert subset.stop == 18
229 with pytest.raises(IndexError):
230 i.absolute[:30]
232 # It might seem surprising that this does not raise, but it's exactly
233 # what what list(range(3, 20))[:30] does:
234 subset = i.local[:30]
235 assert subset.start == 3
236 assert subset.stop == 20
238 with pytest.raises(IndexError):
239 i.absolute[30:]
241 with pytest.raises(IndexError):
242 i.local[30:]
244 with pytest.raises(IndexError):
245 i.absolute[-1:10]
247 with pytest.raises(IndexError):
248 i.local[-1:10]
250 with pytest.raises(ValueError):
251 i.absolute[::2]
253 with pytest.raises(ValueError):
254 i.local[::2]
256 with pytest.raises(ValueError):
257 Interval.factory[1:2:2]
260def test_interval_usage() -> None:
261 """Test using intervals."""
262 i = Interval(start=1, stop=10)
263 d = i.dilated_by(5)
264 assert d == Interval(start=-4, stop=15)
266 s = i.slice_within(Interval(start=-1, stop=12))
267 assert s.start == 2
268 assert s.stop == 11
270 with pytest.raises(IndexError):
271 i.slice_within(Interval(start=3, stop=5))
273 val = i.linspace()
274 assert_close(val, np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]))
275 val = i.linspace(step=2.0)
276 assert_close(val, np.array([1.0, 3.0, 5.0, 7.0, 9.0]))
277 val = i.linspace(n=3)
278 assert_close(val, np.array([1.0, 5.0, 9.0]))
279 with pytest.raises(TypeError):
280 i.linspace(n=2, step=3.0)
282 assert list(i.range) == [1, 2, 3, 4, 5, 6, 7, 8, 9]
283 val = i.arange
284 assert_close(val, np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]))
287def test_interval_pydantic() -> None:
288 """Test roundtrip through pydantic serialization."""
289 i1 = Interval(start=2, stop=5)
290 i2 = Interval(start=-5, stop=10)
291 model = IntervalModel(interval1=i1, interval2=i2)
292 j_str = model.model_dump_json()
293 jmodel = IntervalModel.model_validate_json(j_str)
294 assert jmodel == model
295 assert jmodel.interval1 == i1
298def test_xy_yx_pydantic() -> None:
299 """Test that XY and YX serialize as dicts keyed by dimension name."""
300 model = PairModel(xy=XY(x=1.5, y=2.5), yx=YX(y=3, x=4), bare=XY(x=0.5, y=1))
301 data = model.model_dump()
302 assert data == {
303 "xy": {"x": 1.5, "y": 2.5},
304 "yx": {"y": 3, "x": 4},
305 "bare": {"x": 0.5, "y": 1},
306 }
308 roundtripped = PairModel.model_validate_json(model.model_dump_json())
309 assert roundtripped == model
310 assert isinstance(roundtripped.xy, XY)
311 assert isinstance(roundtripped.yx, YX)
312 assert roundtripped.yx.y == 3
313 assert roundtripped.yx.x == 4
315 json_schema = PairModel.model_json_schema()
317 def resolve(node: dict) -> dict:
318 while "$ref" in node:
319 node = json_schema["$defs"][node["$ref"].rsplit("/", 1)[-1]]
320 return node
322 xy_object, xy_array = (resolve(s) for s in json_schema["properties"]["xy"]["anyOf"])
323 assert xy_object["type"] == "object"
324 assert set(xy_object["required"]) == {"x", "y"}
325 assert xy_object["properties"]["x"]["type"] == "number"
326 assert xy_array == {"type": "array", "items": {"type": "number"}, "minItems": 2, "maxItems": 2}
327 yx_object, yx_array = (resolve(s) for s in json_schema["properties"]["yx"]["anyOf"])
328 assert yx_object["type"] == "object"
329 assert set(yx_object["required"]) == {"x", "y"}
330 assert yx_object["properties"]["y"]["type"] == "integer"
331 assert yx_array == {"type": "array", "items": {"type": "integer"}, "minItems": 2, "maxItems": 2}
334def test_xy_yx_pydantic_array_form() -> None:
335 """Test that the array serialization emitted before the object form was
336 introduced still validates, in both JSON and Python mode (already-parsed
337 JSON is validated as Python objects when reading archives).
338 """
339 from_json = PairModel.model_validate_json('{"xy": [1.5, 2.5], "yx": [3, 4], "bare": [0.5, 1.0]}')
340 assert from_json.xy == XY(x=1.5, y=2.5)
341 assert from_json.yx == YX(y=3, x=4)
342 from_python = PairModel.model_validate({"xy": [1.5, 2.5], "yx": [3, 4], "bare": [0.5, 1.0]})
343 assert from_python == from_json
345 # The array form is restricted to lists so that a transposed pair (XY
346 # and YX are both tuples) fails to validate rather than silently
347 # swapping its elements.
348 with pytest.raises(pydantic.ValidationError):
349 PairModel(xy=YX(y=2.5, x=1.5), yx=YX(y=3, x=4), bare=XY(x=0.5, y=1.0))
352def test_interval_pickle() -> None:
353 """Test pickle roundtrip."""
354 i = Interval(start=5, stop=10)
355 d = pickle.dumps(i)
356 copy = pickle.loads(d)
357 assert copy == i
360def test_box_constructor() -> None:
361 """Test Box construction and factory methods."""
362 y = Interval(-5, 5)
363 x = Interval(32, 64)
364 box = Box(y, x)
366 assert box.start == YX(-5, 32)
367 assert box.shape == YX(10, 32)
368 assert box.area == 320
369 assert box.x == x
370 assert box.y == y
371 assert box != []
373 box2 = Box.factory[-5:5, 32:64]
374 assert box2 == box
376 sbox = Box.from_shape((10, 5))
377 assert sbox == Box.factory[0:10, 0:5]
378 sbox = Box.from_shape((10, 5), start=(2, 3))
379 assert sbox == Box.factory[2:12, 3:8]
381 sbox = Box.from_shape(YX(10, 5))
382 assert sbox == Box.factory[0:10, 0:5]
383 sbox = Box.from_shape((10, 5), start=YX(2, 3))
384 assert sbox == Box.factory[2:12, 3:8]
386 sbox = Box.from_shape(XY(5, 10))
387 assert sbox == Box.factory[0:10, 0:5]
388 sbox = Box.from_shape((10, 5), start=XY(3, 2))
389 assert sbox == Box.factory[2:12, 3:8]
391 with pytest.raises(TypeError):
392 Box.from_shape(42)
393 with pytest.raises(ValueError):
394 Box.from_shape([42])
395 with pytest.raises(ValueError):
396 Box.from_shape([42, 33], start=[1, 2, 3])
398 box = Box.factory[1:2, -1:1]
399 grown = box.dilated_by(2)
400 assert grown == Box.factory[-1:4, -3:3]
403def test_box_constructor_type_checking() -> None:
404 """Verify that the Box constructor rejects non-Interval arguments."""
405 y = Interval(-5, 5)
406 x = Interval(32, 64)
407 with pytest.raises(TypeError, match="Interval"):
408 Box((-5, 5), (32, 64)) # type: ignore[arg-type]
409 with pytest.raises(TypeError, match="Interval"):
410 Box(slice(-5, 5), slice(32, 64)) # type: ignore[arg-type]
411 with pytest.raises(TypeError, match="Interval"):
412 Box(-5, 32) # type: ignore[arg-type]
413 # Each argument is checked, not just the first.
414 with pytest.raises(TypeError, match="Interval"):
415 Box(y, (32, 64)) # type: ignore[arg-type]
416 with pytest.raises(TypeError, match="Interval"):
417 Box((-5, 5), x) # type: ignore[arg-type]
418 with pytest.raises(TypeError, match="Interval"):
419 Box(y=y, x=(32, 64)) # type: ignore[arg-type]
422def test_box_contains() -> None:
423 """Test box containment against other boxes and points."""
424 box = Box.factory[0:20, 10:40]
426 assert box.contains(Box.factory[4:10, 20:25])
427 assert not box.contains(Box.factory[4:10, 35:45])
428 assert box.contains(y=4, x=15)
429 assert not box.contains(x=4, y=15)
431 contains = box.contains(
432 # Half pixel leeway.
433 x=np.array([-1, 10, 20, 30, 40, 41]),
434 y=np.array([-1, 5, 19, 20, 20, 20]),
435 )
436 assert list(contains) == [False, True, True, False, False, False]
438 with pytest.raises(TypeError):
439 box.contains(box, x=3, y=2)
440 with pytest.raises(TypeError):
441 box.contains()
444def test_box_contains_xy_yx() -> None:
445 """Verify that Box.contains accepts XY and YX positional arguments."""
446 box = Box.factory[0:20, 10:40]
447 # Scalar XY and YX — results must match the keyword form.
448 assert box.contains(XY(x=15, y=4)) == box.contains(x=15, y=4)
449 assert box.contains(YX(y=4, x=15)) == box.contains(x=15, y=4)
450 assert not box.contains(XY(x=4, y=15))
451 assert not box.contains(YX(y=15, x=4))
452 # Array XY and YX.
453 xv = np.array([-1, 10, 20, 30, 40, 41])
454 yv = np.array([-1, 5, 19, 20, 20, 20])
455 np.testing.assert_array_equal(box.contains(XY(xv, yv)), box.contains(x=xv, y=yv))
456 np.testing.assert_array_equal(box.contains(YX(yv, xv)), box.contains(x=xv, y=yv))
457 # Mixing positional point with keyword x= or y= must raise TypeError.
458 with pytest.raises(TypeError):
459 box.contains(XY(x=15, y=4), x=15)
460 with pytest.raises(TypeError):
461 box.contains(YX(y=4, x=15), y=4)
464def test_box_contains_broadcasting() -> None:
465 """Test that Box.contains broadcasts like a numpy ufunc."""
466 check_bounds_contains_broadcasting(Box.factory[0:20, 10:40])
469def test_box_intersection() -> None:
470 """Test box intersection."""
471 box1 = Box.factory[0:20, 30:50]
472 box2 = Box.factory[10:30, 40:42]
473 assert box1.intersection(box2) == Box.factory[10:20, 40:42]
474 with pytest.raises(NoOverlapError):
475 box1.intersection(Box.factory[50:70, -10:-5])
478def test_box_slicing() -> None:
479 """Test slicing."""
480 box = Box.factory[:10, :20]
481 sbox = box.absolute[4:6, :3]
482 assert sbox == Box.factory[4:6, 0:3]
483 sbox = box.local[4:6, :3]
484 assert sbox == Box.factory[4:6, 0:3]
485 sbox = box.absolute[4:, 5:]
486 assert sbox == Box.factory[4:10, 5:20]
487 sbox = box.local[4:, 5:]
488 assert sbox == Box.factory[4:10, 5:20]
489 sbox = box.absolute[XY(slice(5, None), slice(4, None))]
490 assert sbox == Box.factory[4:10, 5:20]
491 sbox = box.local[XY(slice(5, None), slice(4, None))]
492 assert sbox == Box.factory[4:10, 5:20]
494 assert Box.factory[4:10, -1:2] == Box.factory[XY(slice(-1, 2), slice(4, 10))]
496 slices = sbox.slice_within(box)
497 assert slices.x.start == 5
498 assert slices.y.start == 4
499 assert slices.x.stop == 20
500 assert slices.y.stop == 10
502 slices = Box.factory[:5, 110:119].slice_within(Box.factory[-15:12, 90:120])
503 assert slices.x.start == 20
504 assert slices.y.start == 15
505 assert slices.x.stop == 29
506 assert slices.y.stop == 20
508 with pytest.raises(IndexError):
509 box.absolute[-1:5, 3:]
510 with pytest.raises(IndexError):
511 box.local[-1:5, 3:]
512 with pytest.raises(TypeError):
513 box.absolute[3:5, :5, 4:]
514 with pytest.raises(TypeError):
515 box.local[3:5, :5, 4:]
516 with pytest.raises(TypeError):
517 Box.factory[3:5, :6, 4:]
520def test_box_mesh() -> None:
521 """Test grid creation."""
522 box = Box.factory[0:2, 0:3]
524 grid = box.meshgrid()
525 assert_close(grid.x, np.array([[0.0, 1.0, 2.0], [0.0, 1.0, 2.0]]))
526 assert_close(grid.y, np.array([[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]]))
528 grid = box.meshgrid(2)
529 assert_close(grid.x, np.array([[0.0, 2.0], [0.0, 2.0]]))
530 assert_close(grid.y, np.array([[0.0, 0.0], [1.0, 1.0]]))
532 grid = box.meshgrid([2, 1])
533 assert_close(grid.x, np.array([[0.0, 2.0]]))
534 assert_close(grid.y, np.array([[0.0, 0.0]]))
536 grid = box.meshgrid(XY(2, 1))
537 assert_close(grid.x, np.array([[0.0, 2.0]]))
538 assert_close(grid.y, np.array([[0.0, 0.0]]))
540 grid = box.meshgrid(YX(1, 2))
541 assert_close(grid.x, np.array([[0.0, 2.0]]))
542 assert_close(grid.y, np.array([[0.0, 0.0]]))
544 grid = box.meshgrid(step=3)
545 assert_close(grid.x, np.array([[0.0]]))
546 assert_close(grid.y, np.array([[0.0]]))
548 with pytest.raises(TypeError):
549 box.meshgrid(2, step=3)
551 with pytest.raises(ValueError):
552 box.meshgrid("n")
555def test_box_boundary() -> None:
556 """Test we can find the boundary."""
557 box = Box.factory[-1:9, 7:15]
558 corners = list(box.boundary())
559 assert corners[0] == (-1, 7)
560 assert corners[1] == (-1, 14)
561 assert corners[2] == (8, 14)
562 assert corners[3] == (8, 7)
565def test_box_pydantic() -> None:
566 """Test roundtrip through pydantic serialization."""
567 box = Box.factory[-1:1, 5:10]
568 model = BoxModel(box=box)
569 j_str = model.model_dump_json()
570 jmodel = BoxModel.model_validate_json(j_str)
571 assert jmodel == model
572 assert jmodel.box == box
575def test_box_pickle() -> None:
576 """Test pickle roundtrip."""
577 box = Box.factory[-1:1, 5:10]
578 d = pickle.dumps(box)
579 copy = pickle.loads(d)
580 assert copy == box
583def test_box_from_float_bounds() -> None:
584 """Test Box.from_float_bounds rounds outward to integer pixel bounds."""
585 # x in [4.6, 9.4], y in [2.6, 5.4] -> box [3:6, 5:10] in [y, x].
586 box = Box.from_float_bounds(x_min=4.6, x_max=9.4, y_min=2.6, y_max=5.4)
587 assert box == Box.factory[3:6, 5:10]
590def _make_sky_projection() -> SkyProjection:
591 """Build a gnomonic sky projection with 0.1 arcsec pixels.
593 The tangent point (12, 13) deg is at 0-based pixel (x=5, y=6).
594 """
595 wcs = astropy.wcs.WCS(naxis=2)
596 # FITS CRPIX is 1-based, so CRPIX (6, 7) is 0-based pixel (x=5, y=6).
597 wcs.wcs.crpix = [6.0, 7.0]
598 wcs.wcs.crval = [12.0, 13.0]
599 scale = 0.1 / 3600.0
600 wcs.wcs.cd = [[-scale, 0.0], [0.0, scale]]
601 wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"]
602 return SkyProjection.from_fits_wcs(wcs, GeneralFrame(unit=u.pix))
605def test_box_from_sky_circle() -> None:
606 """Test Box.from_sky_circle encloses a circular sky region in pixels."""
607 sky_projection = _make_sky_projection()
609 # A 1 arcsec circle is ~10 pixels in radius at 0.1 arcsec per pixel.
610 # Centered on the tangent point at pixel (x=5, y=6), it spans
611 # x [-5, 15] and y [-4, 16].
612 box = Box.from_sky_circle(
613 sky_projection, SkyCoord(ra=12.0 * u.deg, dec=13.0 * u.deg, frame="icrs"), Angle(1.0 * u.arcsec)
614 )
615 assert box == Box.factory[-4:17, -5:16]
617 # A center offset by +5 arcsec in dec and -5 arcsec in RA lands at
618 # (x=53.7, y=56); the RA offset scales by cos(dec).
619 box = Box.from_sky_circle(
620 sky_projection,
621 SkyCoord(ra=12.0 * u.deg - 5.0 * u.arcsec, dec=13.0 * u.deg + 5.0 * u.arcsec, frame="icrs"),
622 Angle(1.0 * u.arcsec),
623 )
624 assert box == Box.factory[46:67, 44:65]
626 # The center may be given in any frame; it is transformed to ICRS, so the
627 # same physical point expressed in galactic coordinates gives the same box.
628 galactic_center = SkyCoord(ra=12.0 * u.deg, dec=13.0 * u.deg, frame="icrs").galactic
629 box = Box.from_sky_circle(sky_projection, galactic_center, Angle(1.0 * u.arcsec))
630 assert box == Box.factory[-4:17, -5:16]
632 # Non-scalar center and radius are rejected.
633 with pytest.raises(ValueError, match="scalar SkyCoord"):
634 Box.from_sky_circle(
635 sky_projection,
636 SkyCoord(ra=[12.0, 12.1] * u.deg, dec=[13.0, 13.1] * u.deg, frame="icrs"),
637 Angle(1.0 * u.arcsec),
638 )
639 with pytest.raises(ValueError, match="scalar Angle"):
640 Box.from_sky_circle(
641 sky_projection,
642 SkyCoord(ra=12.0 * u.deg, dec=13.0 * u.deg, frame="icrs"),
643 Angle([1.0, 2.0] * u.arcsec),
644 )