diff --git a/docs/reference/textured-globe-glyph.md b/docs/reference/textured-globe-glyph.md index 1a7e7501..f80d3309 100644 --- a/docs/reference/textured-globe-glyph.md +++ b/docs/reference/textured-globe-glyph.md @@ -87,3 +87,33 @@ anim = globe.animate(n_frames=60, revolutions=1.0, interval=50, sun=(1.0, 0.0, 0 # from cleopatra.glyphs.base.animation import save_animation # save_animation(anim, "globe.gif") ``` + +### Aligning your own geometry with the globe (the tilt transform) + +The glyph places the sphere with a fixed transform: it spins about the polar axis, then leans that +axis `tilt_deg` from vertical about the world `x` axis — exactly the `R_tilt @ R_z` matrix +`rotation_matrix(spin)` returns. To place your own scene geometry — a marker on the surface, a ring in +the equatorial plane, an orbit plane — so it sits consistently with the rendered globe, push it through +the **same** transform with +`transform(points, spin=...)` (or grab the `(3, 3)` matrix with `rotation_matrix(spin)`). The body +frame is the unit sphere: `+z` at the north pole, so a surface point at `(lon, lat)` is +`[cos(lat)·cos(lon), cos(lat)·sin(lon), sin(lat)]` and the equatorial plane is `z = 0`. + +```python +import numpy as np +from cleopatra.basemap.reference import relief +from cleopatra.glyphs.globe.textured_globe_glyph import TexturedGlobeGlyph + +globe = TexturedGlobeGlyph(relief("low"), tilt_deg=23.44) +fig, ax = globe.draw(spin=40.0) + +# a geostationary ring in the equatorial plane, tilted+spun to match the globe +theta = np.linspace(0, 2 * np.pi, 200) +ring = np.column_stack([1.3 * np.cos(theta), 1.3 * np.sin(theta), np.zeros_like(theta)]) +ring = globe.transform(ring, spin=40.0) +ax.plot(ring[:, 0], ring[:, 1], ring[:, 2]) + +# draw() fixes the axis limits to the unit sphere; widen them so the ring is visible +for set_lim in (ax.set_xlim, ax.set_ylim, ax.set_zlim): + set_lim(-1.4, 1.4) +``` diff --git a/src/cleopatra/glyphs/globe/textured_globe_glyph.py b/src/cleopatra/glyphs/globe/textured_globe_glyph.py index c2cc6fa1..78575905 100644 --- a/src/cleopatra/glyphs/globe/textured_globe_glyph.py +++ b/src/cleopatra/glyphs/globe/textured_globe_glyph.py @@ -57,6 +57,7 @@ import matplotlib.pyplot as plt import numpy as np +import numpy.typing as npt from matplotlib.animation import FuncAnimation from matplotlib.figure import Figure from mpl_toolkits.mplot3d import Axes3D @@ -119,6 +120,8 @@ class TexturedGlobeGlyph: Methods: draw(ax=None, *, spin=0.0, sun=..., ambient=..., **kwargs): Render the globe at a given spin angle. animate(ax=None, n_frames=60, revolutions=1.0, sun=..., ...): Return a `FuncAnimation` spinning the globe. + rotation_matrix(spin=0.0): The `(3, 3)` body-to-world transform the glyph applies (tilt then spin). + transform(points, spin=0.0): Push your own `(N, 3)` scene geometry through that same transform. Notes: `TexturedGlobeGlyph` is a standalone class, not a `Glyph` subclass (like `HistogramGlyph`). The accepted option @@ -231,7 +234,6 @@ def __init__( # Filled lazily and cached by `_prepare` (sample-once contract). self._base_xyz: np.ndarray | None = None self._facecolors: np.ndarray | None = None - self._tilt_matrix: np.ndarray | None = None self._surface = None # ------------------------------------------------------------------ # @@ -292,9 +294,9 @@ def _normalize_texture(texture: np.ndarray, brightness: float) -> np.ndarray: def _prepare(self) -> None: """Sample the texture and build the base sphere mesh once, caching the results on the instance. - Computes and caches the un-spun vertex coordinates `(3, n_lat * n_lon)`, the per-face `facecolors` - `(n_lat - 1, n_lon - 1, 4)` sampled at face centres, and the fixed axial-tilt rotation matrix. Idempotent: - repeated calls (e.g. one per animation frame) return immediately. + Computes and caches the un-spun vertex coordinates `(3, n_lat * n_lon)` and the per-face `facecolors` + `(n_lat - 1, n_lon - 1, 4)` sampled at face centres. Idempotent: repeated calls (e.g. one per animation + frame) return immediately. """ if self._base_xyz is not None: return @@ -324,8 +326,6 @@ def _prepare(self) -> None: row_idx, col_idx = np.meshgrid(rows, cols, indexing="ij") self._facecolors = self._texture[row_idx, col_idx] - self._tilt_matrix = self._rotation_x(self._tilt_deg) - @staticmethod def _rotation_x(deg: float) -> np.ndarray: """Return the 3x3 matrix rotating a point cloud by `deg` degrees about the x-axis.""" @@ -353,9 +353,84 @@ def _spun_mesh(self, spin: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]: Returns: tuple: Three `(n_lat, n_lon)` arrays `(x, y, z)` for `Axes3D.plot_surface`. """ - coords = self._tilt_matrix @ (self._rotation_z(spin) @ self._base_xyz) + coords = self.rotation_matrix(spin) @ self._base_xyz return tuple(coords.reshape(3, self._n_lat, self._n_lon)) + def rotation_matrix(self, spin: float = 0.0) -> np.ndarray: + """Return the 3x3 body-to-world rotation the glyph applies at a given spin. + + This is the exact transform `draw(spin=...)` uses to place the sphere: a rotation of `spin` degrees about + the body polar axis (`z`), then the fixed axial tilt of `tilt_deg` about the world `x` axis -- + `R_tilt(x) @ R_z(spin)`. Apply it (or `transform`) to your own scene geometry so it sits consistently with + the rendered globe without reimplementing the tilt. + + Args: + spin: Rotation about the polar axis, in degrees (matching `draw`/`animate`'s `spin`). + + Returns: + numpy.ndarray: A `(3, 3)` matrix `M` such that a body-frame column vector `p` maps to world as `M @ p`. + + Examples: + - Identity at `spin=0` with no tilt: + ```python + >>> import numpy as np + >>> from cleopatra.glyphs.globe.textured_globe_glyph import TexturedGlobeGlyph + >>> globe = TexturedGlobeGlyph(np.zeros((8, 16, 3), dtype=np.uint8), tilt_deg=0.0) + >>> np.allclose(globe.rotation_matrix(0.0), np.eye(3)) + True + + ``` + + See Also: + transform: Apply this matrix to an `(N, 3)` array of points. + """ + return np.asarray(self._rotation_x(self._tilt_deg) @ self._rotation_z(spin)) + + def transform(self, points: npt.ArrayLike, spin: float = 0.0) -> np.ndarray: + """Map body-frame point(s) into world space exactly as the glyph places its mesh. + + Pushes points through `rotation_matrix(spin)` (spin about the polar axis, then the axial tilt). The body + frame is the same one the mesh is built in: a unit sphere with `+z` at the north pole, so a surface point at + `(lon, lat)` is `[cos(lat) cos(lon), cos(lat) sin(lon), sin(lat)]`, the equatorial plane is `z = 0`, and the + polar axis is `+z`. Use it to place an eclipse marker, a geostationary ring, or an orbit plane so they align + with the rendered globe. + + Args: + points: A single `(3,)` point or an `(N, 3)` array of body-frame points (any array-like). Non-finite + values (`NaN`/`inf`) are propagated, not rejected (`inf` also emits a numpy `RuntimeWarning`) -- + pass finite coordinates. + spin: Rotation about the polar axis, in degrees (matching `draw`/`animate`'s `spin`). + + Returns: + numpy.ndarray: The transformed point(s), same shape as `points` (`(3,)` or `(N, 3)`). + + Raises: + ValueError: If `points` is not `(3,)` or `(N, 3)`. + + Examples: + - The north pole maps to the tilted axis; a 90 deg tilt lays it onto `-y`: + ```python + >>> import numpy as np + >>> from cleopatra.glyphs.globe.textured_globe_glyph import TexturedGlobeGlyph + >>> globe = TexturedGlobeGlyph(np.zeros((8, 16, 3), dtype=np.uint8), tilt_deg=90.0) + >>> np.round(globe.transform([0.0, 0.0, 1.0]), 6) + array([ 0., -1., 0.]) + + ``` + + See Also: + rotation_matrix: The `(3, 3)` matrix this method applies. + """ + pts = np.asarray(points, dtype=float) + if pts.ndim not in (1, 2) or pts.shape[-1] != 3: + raise ValueError( + f"points must be a (3,) point or an (N, 3) array; got shape {pts.shape}." + ) + result = np.atleast_2d(pts) @ self.rotation_matrix(spin).T + if pts.ndim == 1: + result = result[0] + return np.asarray(result) + @staticmethod def _normalize_sun(sun: tuple[float, float, float] | None) -> np.ndarray | None: """Validate a light direction and return it as a unit vector (or `None`). diff --git a/tests/test_textured_globe_glyph.py b/tests/test_textured_globe_glyph.py index 84ea2c9a..09b4bdd8 100644 --- a/tests/test_textured_globe_glyph.py +++ b/tests/test_textured_globe_glyph.py @@ -484,6 +484,99 @@ def test_world_space_sun_honoured_under_tilt(self): ) # but the pole is not the peak -> world-space, not body-space +class TestTiltTransform: + def test_rotation_matrix_identity_without_tilt_or_spin(self, texture): + globe = TexturedGlobeGlyph(texture, tilt_deg=0.0) + assert np.allclose(globe.rotation_matrix(0.0), np.eye(3)) + + def test_rotation_matrix_is_tilt_then_spin(self, texture): + globe = TexturedGlobeGlyph(texture, tilt_deg=30.0) + expected = TexturedGlobeGlyph._rotation_x( + 30.0 + ) @ TexturedGlobeGlyph._rotation_z(47.0) + assert np.allclose(globe.rotation_matrix(47.0), expected) + + def test_transform_lands_where_the_mesh_does(self, texture): + # DoD: a point pushed through the exposed transform lands where the glyph's own mesh puts it + globe = TexturedGlobeGlyph(texture, n_lon=24, n_lat=12, tilt_deg=30.0) + globe._prepare() + spin = 47.0 + mesh = np.stack(globe._spun_mesh(spin)).reshape(3, -1).T # (N, 3) world points + out = globe.transform( + globe._base_xyz.T, spin=spin + ) # base points through the public transform + assert np.allclose(out, mesh) + + def test_transform_single_point_shape_and_value(self, texture): + globe = TexturedGlobeGlyph(texture, tilt_deg=90.0) + out = globe.transform([0.0, 0.0, 1.0]) # north pole under a 90deg x-tilt -> -y + assert out.shape == (3,) + assert np.allclose(out, [0.0, -1.0, 0.0]) + + def test_transform_array_applies_per_row(self, texture): + globe = TexturedGlobeGlyph(texture, tilt_deg=30.0) + pts = np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) + out = globe.transform(pts, spin=12.0) + assert out.shape == (3, 3) + # each row is transformed like a single point (not merely reshaped) + for row_in, row_out in zip(pts, out): + assert np.allclose(row_out, globe.transform(row_in, spin=12.0)) + + @pytest.mark.parametrize( + "bad", + [np.zeros(2), np.zeros((4, 2)), np.zeros((2, 3, 3)), 5.0, np.array(5.0)], + ) + def test_transform_bad_shape_raises(self, texture, bad): + # scalar / 0-d must raise ValueError (not IndexError from indexing shape[-1]) + globe = TexturedGlobeGlyph(texture) + with pytest.raises(ValueError): + globe.transform(bad) + + def test_default_tilt_mesh_unchanged(self, texture): + # the refactor keeps the X-axis default: the mesh equals R_x(tilt) @ R_z(spin) @ base + globe = TexturedGlobeGlyph(texture, n_lon=24, n_lat=12) + globe._prepare() + expected = TexturedGlobeGlyph._rotation_x(globe.tilt_deg) @ ( + TexturedGlobeGlyph._rotation_z(15.0) @ globe._base_xyz + ) + actual = np.stack(globe._spun_mesh(15.0)).reshape(3, -1) + assert np.allclose(actual, expected) + + def test_transform_empty_array_preserved(self, texture): + out = TexturedGlobeGlyph(texture).transform(np.zeros((0, 3))) + assert out.shape == (0, 3) + + def test_transform_1x3_not_squeezed(self, texture): + out = TexturedGlobeGlyph(texture).transform([[1.0, 0.0, 0.0]]) + assert out.shape == (1, 3) + + def test_transform_accepts_list_input(self, texture): + globe = TexturedGlobeGlyph(texture, tilt_deg=20.0) + from_list = globe.transform([0.0, 0.0, 1.0], spin=30.0) + from_array = globe.transform(np.array([0.0, 0.0, 1.0]), spin=30.0) + assert np.allclose(from_list, from_array) + + def test_rotation_matrix_orthogonal_and_fresh(self, texture): + globe = TexturedGlobeGlyph(texture, tilt_deg=23.44) + m = globe.rotation_matrix(47.0) + assert np.allclose(m @ m.T, np.eye(3)) # orthogonal + assert np.isclose(np.linalg.det(m), 1.0) # a proper rotation + m[0, 0] = 9.0 # mutating the returned matrix must not corrupt a later call + assert not np.allclose(globe.rotation_matrix(47.0), m) + + def test_transform_works_before_prepare(self, texture): + globe = TexturedGlobeGlyph(texture) + assert globe._base_xyz is None # never drawn / prepared + assert globe.transform([0.0, 0.0, 1.0], spin=10.0).shape == (3,) + + def test_transform_output_independent_of_input(self, texture): + globe = TexturedGlobeGlyph(texture, tilt_deg=0.0) + inp = np.array([1.0, 2.0, 3.0]) + out = globe.transform(inp, spin=0.0) + out[0] = 99.0 + assert inp[0] == 1.0 # the returned array does not alias the input + + def test_no_new_dependency(): """The globe uses mpl_toolkits.mplot3d, which ships with matplotlib -- no new dependency.""" import mpl_toolkits.mplot3d as m3d