""" The 3D globe: spherical geometry, relief, and the sun. Everything here is pure numpy, so none of it needs a render window. The parts that do need vedo are skipped when it is missing, since it is an optional dependency of the toolbox. """ import numpy as np import pytest from matplotlib.colors import LogNorm, Normalize from dispnc.topography import load_topography from dispnc.render.globe import regrid, render_grid from dispnc.render.globe.geometry import (DATA_RELIEF_FRACTION, Globe, build_surface, normalize_relief) from dispnc.render.globe.scene import (OBLIQUITY_DEG, GlobeScene, subsolar_point, sun_direction, terminator_points) RADIUS = 3389500.0 def a_globe(relief=None): """ A globe on the 1 degree render grid, flat unless relief is given. """ lats, lons, lat2d, _ = render_grid() if relief is None: relief = np.zeros_like(lat2d) return Globe(relief=relief, lats=lats, lons=lons) # --- coordinates ----------------------------------------------------------- def test_cartesian_round_trip(): """ The click readout inverts to_cartesian, so the two must agree. """ globe = a_globe() rng = np.random.default_rng(0) lat = rng.uniform(-89.0, 89.0, 200) lon = rng.uniform(-179.0, 179.0, 200) back_lat, back_lon, height = globe.to_geographic(globe.to_cartesian(lat, lon)) assert np.allclose(back_lat, lat) assert np.allclose(back_lon, lon) assert np.allclose(height, 0.0, atol=1e-6) def test_relief_is_exaggerated_but_reported_unexaggerated(): """ to_cartesian multiplies elevation by the exaggeration; to_geographic must divide it back out, so a click reports metres above the datum. """ globe = a_globe() point = globe.to_cartesian(10.0, 20.0, elevation=1000.0) assert np.linalg.norm(point) == pytest.approx(RADIUS + 1000.0 * globe.exaggeration) _, _, height = globe.to_geographic(point) assert height == pytest.approx(1000.0) def test_the_air_is_stretched_harder_than_the_ground_under_it(): """ The weather is in the bottom scale height, and at the relief's own factor a 10 km cloud is a film on a 3390 km sphere - drawn, worse, inside terrain exaggerated to 182 km. The altitude therefore carries its own larger factor, while the relief keeps the one every other thing on the globe is drawn with. """ globe = a_globe(relief=np.full(render_grid()[2].shape, 1000.0)) assert globe.air_exaggeration > globe.exaggeration ground = globe.to_cartesian(10.0, 20.0) aloft = globe.to_cartesian(10.0, 20.0, altitude=10000.0) assert np.linalg.norm(ground) == pytest.approx(RADIUS + 1000.0 * globe.exaggeration) assert np.linalg.norm(aloft) == pytest.approx( np.linalg.norm(ground) + 10000.0 * globe.air_exaggeration) def test_height_above_ground_inverts_the_altitude_a_shell_was_drawn_at(): """ A point of the shell carries both exaggerations - the ground it stands on at one, its own altitude at the other - so reading it back means taking the relief off before the air factor is divided out. Read the other way, as --shell-color height used to, an isosurface over Olympus Mons came out three times as high as the same surface over the plain beside it. """ lats, lons, lat_grid, _ = render_grid() relief = np.broadcast_to(np.linspace(-8000.0, 18000.0, lats.size)[:, None], lat_grid.shape) globe = a_globe(relief=relief.copy()) for lat, lon in ((60.0, -30.0), (-45.5, 120.0), (0.5, 179.5)): point = globe.to_cartesian(lat, lon, altitude=7000.0) assert globe.height_above_ground(point) == pytest.approx(7000.0, abs=1.0) def test_poles_land_on_the_axis(): globe = a_globe() north = globe.to_cartesian(90.0, 0.0)[0] south = globe.to_cartesian(-90.0, 137.0)[0] assert north[:2] == pytest.approx([0.0, 0.0], abs=1e-6) assert north[2] == pytest.approx(RADIUS) assert south[2] == pytest.approx(-RADIUS) # --- regridding ------------------------------------------------------------ def test_regrid_drops_the_cyclic_duplicate_column(): """ add_cyclic_point on an axis already spanning -180..180 appends a copy at 191.25 degrees, which wraps straight back onto -168.75. The interpolator rejects a repeated coordinate, so regrid has to drop it. """ lons = np.append(np.linspace(-180.0, 180.0, 33), 191.25) lats = np.linspace(90.0, -90.0, 33) values = np.tile(np.arange(lons.size, dtype=float), (lats.size, 1)) _, _, lat_grid, lon_grid = render_grid() out = regrid(lats, lons, values, lat_grid, lon_grid) assert out.shape == lat_grid.shape assert np.isfinite(out).all() def test_regrid_closes_the_seam_of_a_grid_that_goes_right_round(): """ A model grid of 32 longitudes ends at 168.75, one cell short of the turn. Interpolated as it stands, the eleven degrees between that and the first longitude fall off the end of the axis and come back as missing data - an empty wedge on the surface globe, and on a shell an opaque dark red one, which is what the column mode's miscomputed slice was. """ lons = np.linspace(-180.0, 180.0, 33)[:-1] lats = np.linspace(90.0, -90.0, 33) # The field is the longitude itself, wrapped, so the seam can be read values = np.broadcast_to(lons[None, :], (lats.size, lons.size)) _, grid_lons, lat_grid, lon_grid = render_grid() out = regrid(lats, lons, values, lat_grid, lon_grid) assert np.isfinite(out).all(), "no wedge of the globe may be left uncovered" # Across the seam the field runs from 168.75 back to -180 rather than # jumping to some interpolated middle seam = out[0][(grid_lons > 168.75)] assert seam.min() >= -180.0 and seam.max() <= 168.75 def test_regrid_marks_uncovered_cells_as_gaps(): lats = np.linspace(0.0, 60.0, 20) lons = np.linspace(-60.0, 60.0, 20) values = np.ones((20, 20)) _, _, lat_grid, lon_grid = render_grid() out = regrid(lats, lons, values, lat_grid, lon_grid) assert np.isnan(out).any() assert np.nanmax(out) == pytest.approx(1.0) def test_a_regional_field_is_not_wrapped_round_the_planet(): """ Closing the seam is for a field that goes right round. A regional one has to keep its gaps, or half the planet would be filled by interpolating between its two edges. """ from dispnc.coords import close_longitude_seam, needs_seam_column assert needs_seam_column(np.linspace(-180.0, 168.75, 32)) assert not needs_seam_column(np.linspace(-60.0, 60.0, 25)) # An axis holding both -180 and +180 already closes itself; repeating its # first column would make a zero-width cell the interpolator refuses assert not needs_seam_column(np.linspace(-180.0, 180.0, 33)) # And the wrap itself: one more column, holding the first one's data lons = np.linspace(-180.0, 168.75, 32) values = np.tile(np.arange(32.0), (3, 1)) closed_lons, closed = close_longitude_seam(lons, values) assert closed_lons[-1] == pytest.approx(180.0) assert np.array_equal(closed[:, -1], values[:, 0]) assert close_longitude_seam(np.linspace(-60.0, 60.0, 25), np.ones((3, 25)))[1].shape == (3, 25) def test_regrid_handles_descending_latitudes(): """ Most of the sample files store latitude north-to-south. """ lats = np.linspace(90.0, -90.0, 19) lons = np.linspace(-180.0, 180.0, 37) values = np.broadcast_to(lats[:, None], (lats.size, lons.size)).copy() _, _, lat_grid, lon_grid = render_grid() out = regrid(lats, lons, values, lat_grid, lon_grid) # The field is latitude itself, so it must come back as the grid latitude assert np.allclose(out, lat_grid, atol=1e-6) # --- relief ---------------------------------------------------------------- def test_normalized_relief_peaks_at_the_intended_fraction(): values = np.linspace(0.0, 5.0, 100).reshape(10, 10) relief = normalize_relief(values) peak = relief.max() * Globe().exaggeration / RADIUS assert peak == pytest.approx(DATA_RELIEF_FRACTION) assert relief.min() == pytest.approx(0.0) def test_normalized_relief_follows_the_colour_scale(): """ A field spanning ten decades is flat under a linear stretch and readable under the log scale that colours it. """ values = np.geomspace(1e-10, 1e3, 64).reshape(8, 8) linear = normalize_relief(values) logged = normalize_relief(values, norm=LogNorm(vmin=1e-10, vmax=1e3)) # Linear leaves the bulk of the field on the floor; log spreads it evenly assert linear.mean() / linear.max() < 0.1 assert logged.mean() / logged.max() == pytest.approx(0.5, abs=0.1) def test_normalized_relief_survives_nan_and_a_flat_field(): assert np.all(normalize_relief(np.full((4, 4), 7.0)) == 0.0) assert np.isfinite(normalize_relief(np.array([[1.0, np.nan], [3.0, 4.0]]))).all() # --- the sun --------------------------------------------------------------- @pytest.mark.parametrize('ls, expected', [ (0.0, 0.0), (90.0, OBLIQUITY_DEG), (180.0, 0.0), (270.0, -OBLIQUITY_DEG), ]) def test_subsolar_latitude_from_solar_longitude(ls, expected): _, lat = subsolar_point(ls=ls) assert lat == pytest.approx(expected, abs=1e-6) def test_explicit_sun_wins_over_solar_longitude(): assert subsolar_point(sun=(120.0, -30.0), ls=90.0) == (120.0, -30.0) def test_no_sun_asked_for(): assert subsolar_point() is None @pytest.mark.parametrize('lon, lat', [(0.0, 0.0), (45.0, 25.0), (-120.0, -80.0), (0.0, 89.0)]) def test_terminator_is_perpendicular_to_the_sun(lon, lat): globe = a_globe() points = terminator_points(globe, lon, lat) cosines = points @ sun_direction(lon, lat) / np.linalg.norm(points, axis=1) assert np.abs(cosines).max() < 1e-12 # --- the surface mesh ------------------------------------------------------ def test_surface_is_closed_and_faces_outwards(): """ The topography grid holds cell centres, so the mesh needs a cap at each pole. Every edge must be shared by exactly two triangles, and every triangle must face away from the centre or the lighting comes out inverted. """ lats, lons, lat2d, lon2d = render_grid(nlat=30, nlon=60) elevation = np.zeros_like(lat2d) globe = Globe(relief=elevation, lats=lats, lons=lons) pts, faces, scalars = build_surface(globe, lat2d, lon2d, elevation, np.ones_like(lat2d)) assert scalars.size == pts.shape[0] == lat2d.size + 2 a, b, c = pts[faces[:, 0]], pts[faces[:, 1]], pts[faces[:, 2]] outward = np.einsum('ij,ij->i', np.cross(b - a, c - a), (a + b + c) / 3.0) assert (outward > 0).all() edges = {} for triangle in faces: for p, q in ((0, 1), (1, 2), (2, 0)): key = tuple(sorted((triangle[p], triangle[q]))) edges[key] = edges.get(key, 0) + 1 assert set(edges.values()) == {2} def test_pole_scalars_come_from_the_adjacent_ring(): lats, lons, lat2d, lon2d = render_grid(nlat=20, nlon=40) elevation = np.zeros_like(lat2d) globe = Globe(relief=elevation, lats=lats, lons=lons) values = np.tile(np.arange(20.0)[:, None], (1, 40)) values[0, ::2] = np.nan # a partly missing ring must still average _, _, scalars = build_surface(globe, lat2d, lon2d, elevation, values) assert scalars[-2] == pytest.approx(0.0) assert scalars[-1] == pytest.approx(19.0) # --- colour scale ---------------------------------------------------------- vedo = pytest.importorskip('vedo', reason="vedo is an optional dependency") @pytest.mark.parametrize('norm, probes', [ (None, [200.0, 230.0, 260.0]), (Normalize(vmin=200, vmax=260), [200.0, 215.0, 260.0]), (LogNorm(vmin=1e-3, vmax=10.0), [1e-3, 1e-2, 1e-1, 1.0, 10.0]), ]) def test_globe_colours_match_matplotlib(norm, probes): """ The globe must not re-derive a colour scale of its own: a cell has to mean the same thing as it did on the map the user just saw. """ import matplotlib from dispnc.render.globe.geometry import apply_colormap cmap = matplotlib.colormaps['viridis'] mesh = vedo.Sphere(res=12) values = np.linspace(probes[0], probes[-1], mesh.npoints) reference = norm or Normalize(vmin=values.min(), vmax=values.max()) apply_colormap(mesh, 'viridis', values, norm, 'v') lut = mesh.mapper.GetLookupTable() for probe in probes: got = [0.0, 0.0, 0.0] lut.GetColor(float(probe), got) assert got == pytest.approx(cmap(reference(probe))[:3], abs=0.01) def test_symlog_is_pre_normalized_and_says_so(): """ VTK has no symlog scale, so the values are normalized first and the bar title has to carry the real range instead. """ from matplotlib.colors import SymLogNorm from dispnc.render.globe.geometry import apply_colormap mesh = vedo.Sphere(res=8) norm = SymLogNorm(linthresh=1e-3, vmin=-5.0, vmax=5.0) title = apply_colormap(mesh, 'viridis', np.linspace(-5, 5, mesh.npoints), norm, 'v') assert title == 'v (-5 to 5)' assert mesh.mapper.GetLookupTable().GetRange() == (0.0, 1.0) def test_missing_data_is_grey_not_dark_red(): """ VTK's default NaN colour is a dark red that reads as a value. Colouring a mesh settles it in the same breath rather than leaving it to each caller to remember - the one that forgot drew a red wedge over the pole. """ from dispnc.render.globe.geometry import (NO_DATA_COLOR, NO_DATA_TRANSPARENT, apply_colormap) mesh = vedo.Sphere(res=8) values = np.linspace(0.0, 1.0, mesh.npoints) values[0] = np.nan apply_colormap(mesh, 'viridis', values, None, 'v') assert mesh.mapper.GetLookupTable().GetNanColor() == NO_DATA_COLOR # A shell has nothing to be neutral against: there, a gap is a hole apply_colormap(mesh, 'viridis', values, None, 'v', no_data=NO_DATA_TRANSPARENT) assert mesh.mapper.GetLookupTable().GetNanColor() == NO_DATA_TRANSPARENT # --- the atmospheric shell ------------------------------------------------- def test_metric_verticals_are_converted(): from dispnc.render.globe.shell import altitudes_in_metres metres, note = altitudes_in_metres(np.array([0.0, 1.0, 48.3]), 'km') assert metres == pytest.approx([0.0, 1000.0, 48300.0]) assert note is None def test_non_metric_verticals_are_spread_and_reported(): """ Pressure and sigma levels carry no altitude, so the levels are spread over a fixed depth and the user is told the vertical is indexed. """ from dispnc.render.globe.shell import altitudes_in_metres metres, note = altitudes_in_metres(np.array([1000.0, 500.0, 100.0, 10.0]), 'Pa', top_km=40.0) assert metres[0] == 0.0 and metres[-1] == pytest.approx(40000.0) assert np.all(np.diff(metres) > 0) assert 'not a length' in note def test_repeated_meridian_is_dropped(): """ Longitudes spanning -180..180 inclusive carry the same meridian twice, which makes a zero-width cell at every level. """ from dispnc.render.globe.shell import drop_duplicate_longitudes lons = np.linspace(-180.0, 180.0, 33) cube = np.arange(2 * 4 * 33, dtype=float).reshape(2, 4, 33) kept, trimmed = drop_duplicate_longitudes(lons, cube) assert kept.size == 32 and trimmed.shape == (2, 4, 32) assert kept[-1] == pytest.approx(168.75) def test_longitudes_that_do_not_wrap_are_left_alone(): from dispnc.render.globe.shell import drop_duplicate_longitudes lons = np.linspace(-180.0, 168.75, 32) cube = np.zeros((1, 1, 32)) kept, trimmed = drop_duplicate_longitudes(lons, cube) assert kept.size == 32 and trimmed.shape == (1, 1, 32) @pytest.mark.parametrize('spec, expected', [ ('1e-6', [1e-6]), (0.5, [0.5]), ('p50', [50.5]), ('p100', [100.0]), (None, [99.01]), ]) def test_threshold_accepts_values_and_percentiles(spec, expected): from dispnc.render.globe.shell import threshold_values cube = np.arange(1.0, 101.0) got = threshold_values(cube, spec) assert got == pytest.approx(expected, rel=1e-3) def test_thresholds_may_be_nested_and_come_back_sorted(): """ 'p90,p99,p99.9' is one request for three nested surfaces. They are sorted so that "inner" always means "later", whichever order the user typed. """ from dispnc.render.globe.shell import threshold_values cube = np.arange(1.0, 1001.0) assert threshold_values(cube, 'p99,p90,p99.9') == pytest.approx( [900.1, 990.01, 999.001], rel=1e-3) assert threshold_values(cube, '10, 5, 10') == [5.0, 10.0] def test_thresholds_of_an_empty_field_are_empty(): from dispnc.render.globe.shell import threshold_values assert threshold_values(np.full(10, np.nan), 'p99') == [] def test_nested_surfaces_get_rising_opacity(): """ A single surface is drawn at the opacity that was asked for - it used to be scaled up to 0.875, opaque enough to hide the planet the body sits in, back when flat colour left opacity carrying the shape. Several must still ramp, or the outer envelope hides the core it exists to reveal. """ from dispnc.render.globe.shell import _nested_alphas assert _nested_alphas(1, 0.35) == pytest.approx([0.35]) ramp = _nested_alphas(3, 0.35) assert ramp[0] == pytest.approx(0.35) and ramp[-1] == pytest.approx(0.875) assert list(ramp) == sorted(ramp) def test_layers_are_picked_evenly_over_altitude_not_index(): """ A hybrid vertical puts half its levels in the lowest few kilometres. Stepping by index would draw them all down there and describe the model's resolution instead of the atmosphere. """ from dispnc.render.globe.shell import _layer_indices # Bunched near the ground, as a sigma or hybrid coordinate is altitudes = np.concatenate([np.linspace(0, 2000, 20), np.linspace(6000, 50000, 6)]) picked = _layer_indices(altitudes, limit=5) assert len(picked) <= 5 assert picked == sorted(picked) # More than one level above the bunched bottom, which an index step misses assert sum(altitudes[i] > 2000 for i in picked) >= 3 def a_vertical(altitudes, mass=None): """ The resolved vertical a shell mode is built on, from a plain altitude axis. """ from dispnc.render.globe.shell import Vertical altitudes = np.asarray(altitudes, dtype=float) return Vertical(field=altitudes, profile=altitudes, mass=mass) def test_column_integrates_over_metres_and_keeps_gaps(): """ Without layer masses the column falls back on an integral over altitude in metres, so doubling the depth doubles the answer. Cells the field never covered stay NaN rather than integrating to a confident zero. """ from dispnc.render.globe.shell import _column lats, lons = np.array([-10.0, 10.0]), np.array([0.0, 90.0]) cube = np.ones((3, 2, 2)) cube[:, 1, 1] = np.nan shell = _column(a_globe(), a_vertical([0.0, 1000.0, 2000.0]), lats, lons, cube, 'viridis', None, 'dust', 'kg/kg', 0.35, None) # 1 everywhere over 2000 m of column assert '2e+03' in shell.note or '2000' in shell.note assert shell.actors and shell.depth_peeling assert 'kg/kg.m' in shell.note assert 'by height' in shell.note, "and it has to say it is not a mass" def test_the_column_is_weighed_by_mass_when_the_file_allows_it(): """ The physical column is the mass of the field per square metre, sum(q dp/g), which is what the hybrid coefficients and the surface pressure are for. Integrating over metres instead weighs a layer by how thick it is rather than by how much air is in it, which puts most of the answer in the deep upper levels where a tracer is thin and horizontally flat. """ from dispnc.render.globe.shell import _column lats, lons = np.array([-10.0, 10.0]), np.array([0.0, 90.0]) cube = np.full((3, 2, 2), 0.5) # 100, 200 and 300 kg/m2 of air in the three layers mass = np.stack([np.full((2, 2), m) for m in (100.0, 200.0, 300.0)]) shell = _column(a_globe(), a_vertical([0.0, 1e3, 2e3], mass=mass), lats, lons, cube, 'viridis', None, 'dust', 'kg/kg', 0.35, None) # 0.5 kg/kg over 600 kg/m2 of atmosphere assert '300' in shell.note and 'mass-weighted' in shell.note # kg/kg times kg/m2 is kg/m2; nobody writes that as 'kg/kg.kg/m2' assert 'kg/m2' in shell.note and 'kg/kg' not in shell.note def test_structured_grid_scalars_are_in_fortran_order(): """ The one that would fail silently: VTK orders structured points with the first dimension varying fastest, so a C-order ravel scrambles the cube and still produces a plausible-looking picture. """ from dispnc.render.globe.shell import build_structured_grid globe = a_globe() altitudes = np.array([0.0, 5000.0, 10000.0]) lats = np.array([-30.0, 0.0, 30.0, 60.0]) lons = np.array([-90.0, 0.0, 90.0, 180.0, -45.0]) cube = np.arange(3 * 4 * 5, dtype=float).reshape(3, 4, 5) grid = build_structured_grid(globe, altitudes, lats, lons, cube) assert grid.npoints == cube.size stored = np.asarray(grid.pointdata['Scalars']) assert np.array_equal(stored, cube.ravel(order='F')) assert not np.array_equal(stored, cube.ravel(order='C')) # And the value really does belong to the point it sits on expected = globe.to_cartesian(lats[2], lons[1], altitude=altitudes[1])[0] where = np.flatnonzero(stored == cube[1, 2, 1])[0] assert grid.coordinates[where] == pytest.approx(expected) def test_a_contoured_grid_is_closed_round_the_date_line(): """ A StructuredGrid is a box: marching cubes finds nothing outside it, so a model grid ending at 168.75 left the last eleven degrees of the planet with no surface in them at all, and every cloud crossing the date line came out as two bodies with a flat face each. """ from dispnc.render.globe.shell import build_structured_grid globe = a_globe() altitudes = np.array([0.0, 5000.0, 10000.0]) lats = np.linspace(60.0, -60.0, 9) lons = np.linspace(-180.0, 168.75, 32) # A blob straddling the date line: half at each end of the longitude axis cube = np.zeros((3, 9, 32)) cube[1, 3:6, :2] = 1.0 cube[1, 3:6, -2:] = 1.0 grid = build_structured_grid(globe, altitudes, lats, lons, cube) assert grid.npoints == 3 * 9 * 33, "one column more than the data" surface = grid.isosurface(value=0.5) lon = np.rad2deg(np.arctan2(surface.coordinates[:, 1], surface.coordinates[:, 0])) assert surface.npoints > 0 assert np.any(lon > 170.0), "the wedge before the date line has to be filled" # And a regional grid is left alone rather than swept round the planet regional = build_structured_grid(globe, altitudes, lats, np.linspace(-60.0, 60.0, 32), cube) assert regional.npoints == cube.size def test_wind_arrows_are_tangent_to_the_sphere(): """ (u, v) are eastward and northward, so both basis vectors must be unit length, perpendicular to each other, and perpendicular to the radius. """ from dispnc.render.globe.shell import wind_arrows globe = a_globe() lats = np.linspace(-80.0, 80.0, 17) lons = np.linspace(-180.0, 160.0, 18) lon2d, lat2d = np.meshgrid(lons, lats) phi = np.deg2rad(90.0 - lat2d).ravel() theta = np.deg2rad(lon2d).ravel() east = np.column_stack([-np.sin(theta), np.cos(theta), np.zeros_like(theta)]) north = np.column_stack([-np.cos(phi) * np.cos(theta), -np.cos(phi) * np.sin(theta), np.sin(phi)]) position = globe.to_cartesian(lat2d, lon2d) radial = position / np.linalg.norm(position, axis=1)[:, None] assert np.abs((east * radial).sum(axis=1)).max() < 1e-12 assert np.abs((north * radial).sum(axis=1)).max() < 1e-12 assert np.abs((east * north).sum(axis=1)).max() < 1e-12 assert np.linalg.norm(east, axis=1) == pytest.approx(1.0) assert np.linalg.norm(north, axis=1) == pytest.approx(1.0) ones = np.ones_like(lat2d) assert wind_arrows(globe, lats, lons, ones, ones) is not None assert wind_arrows(globe, lats, lons, ones * 0, ones * 0) is None # --- animation ------------------------------------------------------------- def an_axis(role, dim, size): from dispnc.coords import Axis return Axis(role=role, dim=dim, size=size, coord=dim, values=np.arange(size, dtype=float), source='axis-attr', confidence=100) def test_animated_dimension_is_held_back_from_the_plan(): """ The animated axis must leave the axis list, so what is inferred is the shape of one frame rather than a three-dimensional field. """ from dispnc.pipeline import _split_frames axes = [an_axis('T', 'Time', 5), an_axis('Y', 'lat', 3), an_axis('X', 'lon', 4)] data = np.arange(5 * 3 * 4, dtype=float).reshape(5, 3, 4) first, kept, frames, axis, note = _split_frames(data, axes, 'Time') assert [a.dim for a in kept] == ['lat', 'lon'] assert frames.shape == (5, 3, 4) assert np.array_equal(first, data[0]) assert axis.dim == 'Time' assert '5 frames' in note def test_animated_axis_is_moved_to_the_front(): """ The animated dimension is not always first in the file. """ from dispnc.pipeline import _split_frames axes = [an_axis('Y', 'lat', 3), an_axis('T', 'Time', 5), an_axis('X', 'lon', 4)] data = np.arange(3 * 5 * 4, dtype=float).reshape(3, 5, 4) first, kept, frames, axis, _ = _split_frames(data, axes, 'Time') assert frames.shape == (5, 3, 4) assert [a.dim for a in kept] == ['lat', 'lon'] assert np.array_equal(first, data[:, 0, :]) assert axis.dim == 'Time' @pytest.mark.parametrize('dim, size, expected', [ ('nope', 5, 'not among'), ('Time', 1, 'single step'), ]) def test_animation_refusals_leave_the_data_alone(dim, size, expected): from dispnc.pipeline import _split_frames axes = [an_axis('T', 'Time', size), an_axis('Y', 'lat', 3), an_axis('X', 'lon', 4)] data = np.zeros((size, 3, 4)) out, kept, frames, axis, note = _split_frames(data, axes, dim) assert frames is None and axis is None assert out is data and kept is axes assert expected in note def test_a_map_is_animated_without_asking_for_the_globe(): """ Frames used to be refused unless --show-3d or --plot-kind globe was given. A lat/lon map animates on its own now, so the split must not care. """ from dispnc.pipeline import ANIMATABLE, _split_frames axes = [an_axis('T', 'Time', 5), an_axis('Y', 'lat', 3), an_axis('X', 'lon', 4)] _, _, frames, _, _ = _split_frames(np.zeros((5, 3, 4)), axes, 'Time') assert frames is not None assert 'geomap' in ANIMATABLE and 'globe' in ANIMATABLE assert 'section' not in ANIMATABLE def test_animation_pins_the_colour_scale_across_frames(): """ The one that flickers if it regresses. `choose_style` leaves the norm unset for an ordinary field, and an unset norm makes vedo rescale to whatever array it is handed. Handing it one frame at a time would slide the colours with that frame's own minimum and maximum, so the animation has to pin the limits over every frame first. """ from dispnc.colors import choose_style from dispnc.render.globe import GlobeOptions, _attach_animation, render_grid frames = np.stack([np.full((4, 4), 10.0), np.full((4, 4), 200.0)]) _, norm, _ = choose_style(frames, 'tsurf') assert norm is None, "this test is only meaningful while the norm starts unset" lats = np.linspace(80.0, -80.0, 4) lons = np.linspace(-180.0, 135.0, 4) lon2d, lat2d = np.meshgrid(lons, lats) grid_lats, grid_lons, lat_grid, lon_grid = render_grid(nlat=12, nlon=24) # The mesh has to match the render grid, or cmap is handed the wrong number # of scalars and quietly falls back to a 0..1 scale globe = Globe(relief=np.zeros_like(lat_grid), lats=grid_lats, lons=grid_lons) pts, faces, _ = build_surface(globe, lat_grid, lon_grid, np.zeros_like(lat_grid), np.zeros_like(lat_grid)) mesh = vedo.Mesh([pts, faces]) scene = GlobeScene(globe) _attach_animation(scene, mesh, GlobeOptions(frames=frames, frame_dim='Time'), lat2d, lon2d, lat_grid, lon_grid, 'viridis', None, 'tsurf', 'K') ranges = [] for index in range(len(frames)): scene._animation['on_frame'](index) ranges.append(mesh.mapper.GetLookupTable().GetRange()) assert ranges[0] == ranges[1] == (10.0, 200.0) def test_video_paths_are_recognised(): from dispnc.render.extras import globe_output_path assert globe_output_path('out.mp4', animated=True) == 'out_globe.mp4' assert globe_output_path('out.png', animated=True) == 'out_globe.mp4' assert globe_output_path('out.pdf', animated=False) == 'out_globe.png' assert globe_output_path(None) is None def test_matplotlib_never_asked_to_write_a_movie(): """ '--animate ... -o movie.mp4' names the movie; the flat figure beside it has to fall back to a format matplotlib supports. """ from dispnc.figure import still_path assert still_path('movie.mp4') == 'movie.png' assert still_path('fig.pdf') == 'fig.pdf' assert still_path(None) is None def test_layers_are_thinned_and_transparent_where_empty(): """ Two things keep the stack from becoming a solid ball that hides the globe: the levels are thinned to MAX_LAYERS, and alpha follows the value so empty air disappears instead of contributing another opaque lid. """ from dispnc.render.globe.shell import MAX_LAYERS, build_shell class Axis: def __init__(self, values, units=None): self.values, self.units = values, units class Plan: pass levels = np.linspace(0.1, 40.0, 26) lats = np.linspace(80.0, -80.0, 9) lons = np.linspace(-180.0, 160.0, 12) plan = Plan() plan.z = Axis(levels, 'km') plan.y = Axis(lats) plan.x = Axis(lons) # Empty everywhere except one blob halfway up cube = np.zeros((26, 9, 12)) cube[13, 4:6, 5:7] = 1.0 shell = build_shell(a_globe(), plan, cube, 'viridis', 'dust', mode='layers', opacity=0.5) actors = shell.actors assert 0 < len(actors) <= MAX_LAYERS assert shell.depth_peeling, "a translucent stack needs order-independent blending" alphas = np.concatenate([np.asarray(a.pointdata['Scalars_alpha']) if a.pointdata['Scalars_alpha'] is not None else np.zeros(a.npoints) for a in actors]) assert alphas.min() == pytest.approx(0.0), "empty air must be fully transparent" def a_cloud_plan(levels=None): """ A minimal plan over a small atmospheric cube, for the shell modes. """ class Axis: def __init__(self, values, units=None): self.values, self.units = values, units class Plan: pass plan = Plan() plan.z = Axis(np.linspace(0.1, 40.0, 26) if levels is None else levels, 'km') plan.y = Axis(np.linspace(80.0, -80.0, 9)) plan.x = Axis(np.linspace(-180.0, 160.0, 12)) return plan def test_nested_isosurfaces_are_drawn_from_inside_out(): """ Three thresholds must give three bodies, each less transparent than the one around it, all on the field's own colour scale rather than a second one. """ from dispnc.render.globe.shell import build_shell plan = a_cloud_plan() cube = np.zeros((26, 9, 12)) # A blob whose core is denser than its edge, so every level has a surface cube[10:16, 3:7, 4:8] = 1.0 cube[12:14, 4:6, 5:7] = 5.0 shell = build_shell(a_globe(), plan, cube, 'viridis', 'dust', mode='iso', threshold='0.5,2,4', opacity=0.3) assert len(shell.actors) == 3 alphas = [a.alpha() for a in shell.actors] assert alphas == sorted(alphas), "the core must be the least transparent" assert 'iso at' in shell.caption[0] def test_a_shell_says_what_it_is_on_screen(): """ The caption is the only record of the threshold, the vertical span and the exaggeration once the terminal has scrolled away. """ from dispnc.render.globe.shell import build_shell globe = a_globe() shell = build_shell(globe, a_cloud_plan(), np.ones((26, 9, 12)), 'viridis', 'dust', mode='column') text = ' '.join(shell.caption) assert 'column' in text assert 'km' in text # Both factors, because they are not the same number assert f"relief x{globe.exaggeration:g}" in text assert f"altitude x{globe.air_exaggeration:g}" in text def test_the_levels_of_a_shell_follow_the_terrain_when_the_file_says_where(): """ Handed the hybrid altitudes, a level is a map of heights rather than one height: it drapes the ground at the bottom and flattens into an isobar at the top. Handed nothing, the shell keeps the vertical axis' own values and every level is a sphere, which is what it always used to be. """ from dispnc.render.globe.shell import build_shell from dispnc.vertical import Altitudes plan = a_cloud_plan(levels=np.linspace(0.1, 40.0, 26)) nlat, nlon = plan.y.values.size, plan.x.values.size cube = np.ones((26, nlat, nlon)) # A ridge down one meridian: the levels over it start lower above the ground ridge = np.zeros((nlat, nlon)) ridge[:, 5] = 3000.0 above_ground = np.maximum(np.linspace(10.0, 40e3, 26)[:, None, None] - ridge, 1.0) air = Altitudes(above_ground=above_ground, interfaces=None, gravity=3.72, note='hybrid') shaped = build_shell(a_globe(), plan, cube, 'viridis', 'dust', mode='level', level=20, altitudes=air) plain = build_shell(a_globe(), plan, cube, 'viridis', 'dust', mode='level', level=20) radii = np.linalg.norm(shaped.actors[0].coordinates, axis=1) assert np.ptp(radii) > 0.0, "a hybrid level is not a sphere" assert 'hybrid' in shaped.note flat = np.linalg.norm(plain.actors[0].coordinates, axis=1) assert np.ptp(flat) == pytest.approx(0.0, abs=1e-6), \ "without the coefficients every level stays one height" # --- the cloud ------------------------------------------------------------- def a_blob(shape=(26, 9, 12)): """ An empty atmosphere with one blob halfway up, which is what a cloud is. """ cube = np.zeros(shape) cube[12:15, 3:6, 4:7] = 1.0 return cube def test_a_cloud_is_one_volume_and_needs_no_depth_peeling(): """ The whole point of ray casting the field is that the mapper composites it itself: one actor, in the right order, without the peeling a stack of translucent surfaces has to ask for. """ from vedo import Volume from dispnc.render.globe.shell import build_shell shell = build_shell(a_globe(), a_cloud_plan(), a_blob(), 'viridis', 'dust', mode='cloud', resolution=48) assert len(shell.actors) == 1 volume = shell.actors[0] assert isinstance(volume, Volume) assert list(volume.dimensions()) == [48, 48, 48] assert not shell.depth_peeling assert 'cloud' in ' '.join(shell.caption) def test_empty_air_stays_at_the_bottom_of_the_scale(): """ A voxel the atmosphere never reaches has to read as the bottom of the scale, which the opacity curve pins to fully transparent. Anything else and the cloud ends at the box instead of fading out. """ from dispnc.render.globe.cloud import sample_volume voxels, spacing, origin = sample_volume( a_globe(), np.linspace(0.0, 40e3, 26), np.linspace(80, -80, 9), np.linspace(-180, 160, 12), a_blob(), None, 0.0, resolution=48) assert voxels.min() == 0.0, "empty air must sit at the bottom of the scale" assert voxels.max() > 0.0, "and the blob must survive the resampling" # The corner of the box is outside the planet altogether assert voxels[0, 0, 0] == 0.0 # And so is its centre, which is inside the planet assert voxels[24, 24, 24] == 0.0 assert spacing > 0 and len(origin) == 3 def test_the_cloud_is_faded_on_the_same_scale_it_is_coloured_on(): """ The user's choice: transparency follows the plot's own norm. h2o_ice spans thirty-five decades on a log scale, so a linear alpha would leave everything but the peak invisible. """ from dispnc.render.globe.cloud import display_scale, normalize cube = np.array([[[1e-30, 1e-15, 1e-5]]]) transform, vmin, vmax, normalized = display_scale(LogNorm(vmin=1e-30, vmax=1e-5), cube) assert normalized and (vmin, vmax) == (0.0, 1.0) # A (level, lat, lon) cube, which is the shape matplotlib's log transform # refuses unless it is flattened first scaled = normalize(transform, cube) assert scaled.shape == cube.shape assert 0.4 < scaled[0, 0, 1] < 0.7, "a mid-decade value belongs in the middle" assert scaled[0, 0, 0] == 0.0 and scaled[0, 0, 2] == pytest.approx(1.0) # A linear scale needs no transform: VTK takes the limits, so the bar keeps # the field's own numbers transform, vmin, vmax, normalized = display_scale(Normalize(vmin=0.0, vmax=4.0), cube) assert transform is None and not normalized assert (vmin, vmax) == (0.0, 4.0) def test_the_cloud_has_no_slit_down_the_date_line(): """ The longitude axis stops at 160 while the field wraps all the way round, so without closing the circle the sampler falls off the end and the cloud comes out with a meridian-wide gap. """ from dispnc.render.globe.cloud import _interpolator values = np.ones((26, 9, 12)) interp = _interpolator(np.linspace(0.0, 40e3, 26), np.linspace(80, -80, 9), np.linspace(-180, 160, 12), values) seam = interp([[20e3, 0.0, lon] for lon in (160.0, 175.0, 179.9, 180.0)]) assert np.all(np.isfinite(seam)), "the seam must carry values, not gaps" assert np.allclose(seam, 1.0) def test_the_opacity_unit_is_the_depth_of_the_atmosphere(): """ VTK reads an opacity as "this much over one unit distance" and corrects each sample for the distance it stepped, so the unit has to be a length. Tied to the voxel size instead, --shell-resolution rather than the atmosphere would decide how thick the cloud looks. """ from dispnc.render.globe.shell import build_shell coarse, fine = [build_shell(a_globe(), a_cloud_plan(), a_blob(), 'viridis', 'dust', mode='cloud', resolution=n).actors[0] for n in (32, 64)] assert coarse.alpha_unit() == pytest.approx(fine.alpha_unit()), \ "resolution must not change how dense the cloud is" # a_cloud_plan spans 0.1 to 40 km, stretched by the air's own exaggeration assert coarse.alpha_unit() == pytest.approx( (40.0 - 0.1) * 1000.0 * a_globe().air_exaggeration) def test_the_column_follows_a_log_scale_instead_of_flattening_it(): """ The integral is not the field, so its range is its own - but the *shape* of the scale carries over. Drawn linearly, a field spanning decades put all but its brightest cell at the floor and the globe came back empty. """ from dispnc.render.globe.shell import _column_scale amount = np.geomspace(1e-4, 1e-2, 100) logged = _column_scale(LogNorm(vmin=1e-30, vmax=1e-5), amount, 1e-4, 1e-2) assert isinstance(logged, LogNorm) assert (logged.vmin, logged.vmax) == (1e-4, 1e-2), "the integral's own range" plain = _column_scale(Normalize(vmin=0, vmax=1), amount, 2.0, 8.0) assert type(plain) is Normalize and (plain.vmin, plain.vmax) == (2.0, 8.0) # A field with no range at all cannot be scaled by it flat = _column_scale(None, amount, 5.0, 5.0) assert flat.vmin == 5.0 and flat.vmax > 5.0 def test_a_log_column_gets_a_floor_under_its_numerical_zero(): """ `h2o_ice` integrates to 2e-35 in the driest cell and 4e-2 in the wettest. A log scale over all thirty-four decades puts every value that matters within a per cent of the top, and the globe comes back one flat yellow - which is what it did. The floor lifts the bottom to where the values are. """ from dispnc.render.globe.shell import MAX_COLUMN_DECADES, _column_scale # A hundred cells between 1e-4 and 4e-2, and one numerical zero amount = np.append(np.geomspace(1e-4, 4e-2, 100), 2e-35) scale = _column_scale(LogNorm(vmin=1e-40, vmax=1e-5), amount, float(amount.min()), float(amount.max())) assert isinstance(scale, LogNorm) assert scale.vmax == pytest.approx(4e-2) assert scale.vmin > 1e-30, "the numerical zero must not set the floor" assert scale.vmax / scale.vmin <= 10.0 ** MAX_COLUMN_DECADES # The bulk of the map still has room to spread out over the scale assert scale.vmin <= np.percentile(amount[amount > 0], 25) def test_a_column_that_spans_little_keeps_its_own_floor(): """ The floor is for a scale that has run away, not for every log scale: a field spanning four decades is drawn over exactly those four. """ from dispnc.render.globe.shell import _column_scale amount = np.geomspace(1e-6, 1e-2, 100) scale = _column_scale(LogNorm(vmin=1e-8, vmax=1.0), amount, 1e-6, 1e-2) assert (scale.vmin, scale.vmax) == (1e-6, 1e-2) def test_the_column_clears_the_highest_ground_on_the_planet(): """ The column shell was once laid against the terrain, to read as the map it is. Ten kilometres above the datum sounds close enough to the ground until you remember the relief on this globe is exaggerated tenfold and reaches a hundred and eighty: better than a quarter of Mars grew straight through the shell and the map came out in pieces. It no longer has to go all the way to the model top to be safe - the air is stretched three times harder than the relief, so seven kilometres of it already clear every summit on the planet. And it is a sphere, not a shell over the terrain. A column integral is a map of the atmosphere; drawing it on the ground gave it a ripple that came from the topography rather than from the field. """ from dispnc.render.globe.geometry import RELIEF_EXAGGERATION from dispnc.render.globe.shell import _column topo = load_topography() if topo is None: # pragma: no cover pytest.skip("topography is needed to have terrain to clear") # A globe with real relief: a flat one cannot reproduce any of this globe = a_globe(relief=regrid(topo.lats, topo.lons, topo.values, *render_grid()[2:])) shell = _column(globe, a_vertical([0.0, 48e3]), np.array([-10.0, 10.0]), np.array([0.0, 90.0]), np.ones((2, 2, 2)), 'viridis', None, 'dust', 'kg/kg', 0.35, None) radii = np.linalg.norm(shell.actors[0].coordinates, axis=1) highest_ground = globe.radius + topo.values.max() * RELIEF_EXAGGERATION # A metre of slack: vedo keeps its points in float32, so a 3390 km radius # is quantized to about a quarter of a metre assert np.ptp(radii) == pytest.approx(0.0, abs=1.0), "the map is a sphere" assert radii.min() > highest_ground, "no summit may reach through the shell" assert radii.max() < globe.radius + 48e3 * globe.air_exaggeration, \ "and it no longer has to sit at the model top to manage it" assert 'km up' in shell.note, "and it has to say where it put it" def test_a_layer_still_rides_over_the_ground_the_column_no_longer_does(): """ The distinction the column depends on: a model level belongs over the ground it sits on, a map of a column integral does not. """ from dispnc.render.globe.shell import _level_mesh lats, lons, lat_grid, _ = render_grid() relief = np.broadcast_to(np.linspace(-8000.0, 18000.0, lats.size)[:, None], lat_grid.shape) globe = a_globe(relief=relief.copy()) over_ground = _level_mesh(globe, 5000.0, lats[::10], lons[::10]) on_datum = _level_mesh(globe, 5000.0, lats[::10], lons[::10], on_datum=True) assert np.ptp(np.linalg.norm(over_ground.coordinates, axis=1)) > 100e3 assert np.ptp(np.linalg.norm(on_datum.coordinates, axis=1)) == pytest.approx(0.0, abs=1.0) def test_isosurfaces_are_smoothed_and_lit_without_being_asked(): """ An isosurface is one flat colour over a whole body, so shape is all it has to say - and a facetted, unlit surface says none of it. This is the one thing on the globe that does not take the globe's lighting answer. """ from dispnc.render.globe.shell import build_shell plan = a_cloud_plan() cube = np.zeros((26, 9, 12)) cube[10:16, 3:7, 4:8] = 1.0 shell = build_shell(a_globe(), plan, cube, 'viridis', 'dust', mode='iso', threshold='0.5', lit=False) surface = shell.actors[0] assert surface.pointdata['Normals'] is not None, "shading needs normals" assert surface.properties.GetInterpolation() > 0, "lit despite lit=False" def test_the_shell_and_the_surface_do_not_share_a_colour_bar_slot(): """ They were the same rectangle, so a terrain globe under a stack of layers drew the elevation bar and the field bar on top of each other and neither could be read. """ from dispnc.render.globe import SURFACE_BAR_POS from dispnc.render.globe.shell import SHELL_BAR_POS assert SURFACE_BAR_POS != SHELL_BAR_POS def test_layers_ask_to_be_cut_open_and_bring_a_cross_section(): """ Concentric shells cannot be seen into: the outer one covers every other, and VTK cannot even order them, since sorting translucent actors by centroid puts every shell at the same distance. The cut is what makes the mode readable, and the cross-section is what the cut reveals - clipping alone exposes nothing, because surfaces have no interior. """ from dispnc.render.globe.shell import build_shell shell = build_shell(a_globe(), a_cloud_plan(), a_blob(), 'viridis', 'dust', mode='layers', resolution=32) assert shell.cut, "layers cannot be read any other way" assert shell.slicer is not None section = shell.slicer(np.zeros(3), np.array([-1.0, 0.0, 0.0])) assert section is not None and section.npoints > 0 def test_the_cut_is_off_for_the_modes_that_do_not_need_it_and_forceable(): """ A cloud is legible whole, so it is not cut unless asked; asking must work all the same, since a cloud is worth looking inside. """ from dispnc.render.globe.shell import build_shell plain = build_shell(a_globe(), a_cloud_plan(), a_blob(), 'viridis', 'dust', mode='cloud', resolution=32) assert not plain.cut and plain.slicer is None forced = build_shell(a_globe(), a_cloud_plan(), a_blob(), 'viridis', 'dust', mode='cloud', resolution=32, cut=True) assert forced.cut and forced.slicer is not None quiet = build_shell(a_globe(), a_cloud_plan(), a_blob(), 'viridis', 'dust', mode='layers', resolution=32, cut=False) assert not quiet.cut def test_a_cutter_clips_every_actor_it_is_given(): """ Cutting the shell and leaving the planet whole swaps one lid for another: the globe fills the opening. Everything with a mapper has to end at the same plane. """ import vedo from dispnc.render.globe import scene as scene_mod globe = a_globe() scene = scene_mod.GlobeScene(globe) ball, ring = vedo.Sphere(res=8), vedo.Sphere(res=8) scene.add(ball, ring, vedo.Text2D("not a mesh")) scene.add_cutter() for mesh in (ball, ring): assert mesh.mapper.GetNumberOfClippingPlanes() == 1 def test_the_peel_count_travels_with_the_peeling_flag(monkeypatch): """ Turning peeling on and leaving VTK at four passes resolves the first two shells of a stack and guesses the rest. Both are process-global, so both have to be put back. """ from vedo import settings from dispnc.render.globe import scene as scene_mod seen = [] class FakePlotter: def __init__(self, **kwargs): seen.append(settings.max_number_of_peels) def show(self, *actors, **kwargs): return self def close(self): pass monkeypatch.setattr(scene_mod, 'Plotter', FakePlotter) monkeypatch.setattr(settings, 'use_depth_peeling', False) monkeypatch.setattr(settings, 'max_number_of_peels', 4) scene_mod.GlobeScene(a_globe()).want_depth_peeling(True).show() assert seen == [scene_mod.DEPTH_PEELS] assert settings.max_number_of_peels == 4, "and must give it back" def test_a_single_level_says_so_instead_of_taking_the_process_down(): """ One level has no depth to accumulate along, so the cloud stands aside - and the contouring it stands aside for has no cells to march either. VTK does not merely fail at that, it segfaults, so the fallback has to stop first and name a mode that does work. """ from dispnc.render.globe.shell import build_shell plan = a_cloud_plan(levels=np.array([12.0])) shell = build_shell(a_globe(), plan, np.ones((1, 9, 12)), 'viridis', 'dust', mode='cloud', threshold='0.5', resolution=32) assert not shell.actors assert 'single level' in shell.note assert '--shell-level' in shell.note or '--shell-mode' in shell.note def test_a_cloud_falls_back_on_contouring_when_it_cannot_be_built(): """ vedo is optional and scipy does the resampling, so the default mode has to degrade to one that needs neither rather than draw nothing at all. """ from dispnc.render.globe import cloud as cloud_mod from dispnc.render.globe.shell import build_shell with pytest.MonkeyPatch.context() as patch: patch.setattr(cloud_mod, 'volume_available', False) shell = build_shell(a_globe(), a_cloud_plan(), a_blob(), 'viridis', 'dust', mode='cloud', threshold='0.5') assert shell.actors, "the fallback has to draw something" assert 'iso at' in shell.caption[0] assert 'contouring instead' in shell.note def test_the_alpha_curve_starts_transparent_and_follows_the_scale(): """ Zero at the bottom so empty air disappears, and straight from there, so the transparency says what the colour says. --shell-opacity is the top of it unscaled: a volume applies its alpha along the whole ray, not once. """ from dispnc.render.globe.cloud import alpha_curve curve = alpha_curve(0.35, stops=5) assert curve[0] == 0.0 assert curve[-1] == pytest.approx(0.35) assert curve == sorted(curve) assert curve[2] == pytest.approx(0.175), "straight, not bent a second time" def test_unlit_is_the_default_until_a_sun_is_asked_for(): """ Shading multiplies the colormap by the angle to the light, so the same value reads differently across the sphere. --sun is the request for the opposite. """ from dispnc.render.globe import GlobeOptions assert GlobeOptions().lit(None) is False assert GlobeOptions().lit((0.0, 20.0)) is True assert GlobeOptions(lighting=True).lit(None) is True assert GlobeOptions(lighting=False).lit((0.0, 20.0)) is False def test_the_globe_window_is_closed_when_it_is_dismissed(monkeypatch): """ The reported bug: the interactor returns when the user shuts the window, but the VTK render window was never torn down, so it lingered and the next globe of an interactive session opened behind it. """ from dispnc.render.globe import scene as scene_mod closed = [] class FakePlotter: def __init__(self, **kwargs): self.kwargs = kwargs def show(self, *actors, **kwargs): return self def close(self): closed.append(self) def add_callback(self, *a, **k): pass monkeypatch.setattr(scene_mod, 'Plotter', FakePlotter) scene = scene_mod.GlobeScene(a_globe()) assert scene.show() == 0 assert len(closed) == 1 assert scene._plotter is None, "a stale plotter keeps the click callback alive" def test_the_window_is_closed_even_when_showing_raises(monkeypatch): from dispnc.render.globe import scene as scene_mod closed = [] class ExplodingPlotter: def __init__(self, **kwargs): pass def show(self, *actors, **kwargs): raise RuntimeError("no GL context") def close(self): closed.append(self) monkeypatch.setattr(scene_mod, 'Plotter', ExplodingPlotter) with pytest.raises(RuntimeError): scene_mod.GlobeScene(a_globe()).show() assert len(closed) == 1 def test_depth_peeling_does_not_leak_into_later_views(monkeypatch): """ vedo reads settings.use_depth_peeling when a Plotter is built, and it is process-global. Setting it where the translucent actor is made leaked it into every later view, flat ones included. """ from vedo import settings from dispnc.render.globe import scene as scene_mod class FakePlotter: def __init__(self, **kwargs): # What vedo does: read the flag at construction time seen.append(settings.use_depth_peeling) def show(self, *actors, **kwargs): return self def close(self): pass seen = [] monkeypatch.setattr(scene_mod, 'Plotter', FakePlotter) monkeypatch.setattr(settings, 'use_depth_peeling', False) scene_mod.GlobeScene(a_globe()).want_depth_peeling(True).show() assert seen == [True], "the translucent scene must get depth peeling" assert settings.use_depth_peeling is False, "and must give it back" # --- the spin and the encoder ---------------------------------------------- def test_spin_becomes_its_own_animation_and_rides_on_an_existing_one(): """ A still globe with --spin has nothing to step through, so the turn is the animation. With --animate already set it must ride along instead of replacing it, or asking for both silently loses the data. """ from dispnc.render.globe import scene as scene_mod alone = scene_mod.GlobeScene(a_globe()).spin(1.0, fps=10) assert alone._animation is not None assert alone._animation['count'] == int(10 * scene_mod.SPIN_SECONDS) stepped = [] both = scene_mod.GlobeScene(a_globe()) both.animate(5, stepped.append, label='time') both.spin(1.0, fps=10) assert both._animation['count'] == 5, "the data's frames, not the turn's" both._animation['on_frame'](3) assert stepped == [3], "the original callback still has to run" def test_a_whole_turn_ends_where_it_began(): """ The first frame is drawn where the camera already is, so the steps have to fall in the gaps between frames rather than before the first one. """ from dispnc.render.globe.scene import GlobeScene scene = GlobeScene(a_globe()).spin(1.0, fps=12) spin = scene._spin assert spin['step'] * spin['frames'] == pytest.approx(360.0) def test_the_encoder_refuses_clearly_without_ffmpeg(monkeypatch, capsys): """ A raw pipe needs the binary. Saying so beats writing a file nobody can play. """ from dispnc.render import movie monkeypatch.setattr(movie.shutil, 'which', lambda name: None) assert movie.open_encoder('out.mp4', 12, (100, 200, 3)) is None assert 'ffmpeg' in capsys.readouterr().out def test_the_camera_frames_the_atmosphere_and_not_just_the_planet(): """ A fixed 4.5 radii frames 1.21 R of half-height, and the air is drawn at thirty times its own depth: a 48 km model top reaches 1.43 R, so the outer half of every stack of layers and the whole rim of every cloud fell outside the picture. The distance follows the scene. """ from dispnc.render.globe.scene import CAMERA_DISTANCE, GlobeScene globe = a_globe() bare = GlobeScene(globe) assert np.linalg.norm(bare.camera()['pos']) == pytest.approx( globe.radius * CAMERA_DISTANCE), "a surface globe is framed as it was" # A shell out at 1.4 R, the size --shell-mode layers reaches reach = globe.radius * 1.4 with_shell = GlobeScene(globe).add(vedo.Sphere(r=reach)) distance = np.linalg.norm(with_shell.camera()['pos']) assert distance == pytest.approx(reach * CAMERA_DISTANCE, rel=0.01) assert distance > globe.radius * CAMERA_DISTANCE # An explicit distance still wins, and is still in planet radii assert np.linalg.norm(with_shell.camera(distance=3.0)['pos']) == pytest.approx( globe.radius * 3.0) # --- lines that lie on the ground ------------------------------------------ def a_real_globe(): """ A globe carrying MOLA, since a flat one cannot reproduce anything here. """ topo = load_topography() if topo is None: # pragma: no cover pytest.skip("topography is needed to have terrain to clear") lats, lons, _, _ = render_grid() # Exactly what render_globe hands the Globe for --globe-relief topo return Globe(relief=topo.values, lats=lats, lons=lons), topo.values def _mesh_elevation(relief, lats, lons): """ The height of the surface mesh itself at arbitrary lat/lon, in metres. Not an interpolation of the relief - the actual triangle. `build_surface` cuts each cell along the p1-p2 diagonal, so a point on one side of that diagonal sits on the plane through (p0, p1, p2) and on the other side on the plane through (p1, p2, p3). Anything drawn on the globe has to clear this, and nothing else: the cell's highest corner is a bound, not the answer, and demanding a line clear that would ask for a lift the mesh never needs. """ grid_lats, grid_lons, _, _ = render_grid() # build_surface closes the seam with (j + 1) % nlon, so the cell spanning # the date line is a cell like any other and the test has to see it that way grid_lons = np.append(grid_lons, grid_lons[0] + 360.0) relief = np.concatenate([relief, relief[:, :1]], axis=1) nlat, nlon = relief.shape lons = grid_lons[0] + np.mod(lons - grid_lons[0], 360.0) i = np.clip(np.searchsorted(-grid_lats, -lats) - 1, 0, nlat - 2) j = np.clip(np.searchsorted(grid_lons, lons) - 1, 0, nlon - 2) # Where in the cell, as fractions along longitude and latitude u = np.clip((lons - grid_lons[j]) / (grid_lons[j + 1] - grid_lons[j]), 0.0, 1.0) v = np.clip((lats - grid_lats[i]) / (grid_lats[i + 1] - grid_lats[i]), 0.0, 1.0) z0, z1 = relief[i, j], relief[i, j + 1] z2, z3 = relief[i + 1, j], relief[i + 1, j + 1] lower = z0 + u * (z1 - z0) + v * (z2 - z0) upper = z3 + (1.0 - u) * (z2 - z3) + (1.0 - v) * (z1 - z3) return np.where(u + v <= 1.0, lower, upper) def _sunk(points, relief, globe): """ How far each point sits below the mesh under it, in metres. Negative is clear of it. """ lats, lons, _ = globe.to_geographic(points) ground = _mesh_elevation(relief, np.atleast_1d(lats), np.atleast_1d(lons)) radius = np.linalg.norm(points, axis=1) return (globe.radius + ground * globe.exaggeration) - radius def test_the_graticule_never_sinks_into_the_terrain(): """ The regression this was written for: meridians drawn by interpolating the relief at points *between* the mesh nodes dipped below the triangles built from those nodes, so every line vanished into every escarpment and came out the far side. Twenty-one kilometres under, at worst. """ from dispnc.render.globe import graticule globe, relief = a_real_globe() grid_lats, grid_lons, _, _ = render_grid() for line in graticule(globe, grid_lats, grid_lons): sunk = _sunk(line.coordinates, relief, globe) assert sunk.max() <= 0.0, f"buried by {sunk.max() / 1000:.1f} km" def test_the_topography_contours_never_sink_either(): """ Contours land wherever the field puts them, so they cannot be snapped to the grid the way the graticule is - the dilated relief is what carries them. """ from dispnc.render.globe import topography_contours globe, relief = a_real_globe() _, _, lat_grid, lon_grid = render_grid() lines = topography_contours(globe, lat_grid, lon_grid, relief, levels=6) assert lines, "there should be contours to check" worst = max(_sunk(line.coordinates, relief, globe).max() for line in lines) assert worst <= 0.0, f"buried by {worst / 1000:.1f} km" def test_lines_hug_the_ground_rather_than_floating_over_it(): """ Clearing the terrain is half of it; the other half is not standing off it. A single global offset big enough for the worst cliff on Mars would put every line sixty kilometres up, which is not "on the topography" either. Since `on_surface` reads the triangles rather than the smooth surface through their corners, the standoff is the lift *exactly*, everywhere - not the lift plus however much the two surfaces happened to differ by, which was up to twenty kilometres and was what made the lines read as hovering. """ from dispnc.render.globe import graticule, topography_contours from dispnc.render.globe.scene import CONTOUR_LIFT, GRATICULE_LIFT globe, relief = a_real_globe() grid_lats, grid_lons, lat_grid, lon_grid = render_grid() for lines, lift in ( (graticule(globe, grid_lats, grid_lons), GRATICULE_LIFT), (topography_contours(globe, lat_grid, lon_grid, relief, levels=6), CONTOUR_LIFT)): gaps = np.concatenate([-_sunk(line.coordinates, relief, globe) for line in lines]) ceiling = lift * globe.exaggeration assert gaps.max() == pytest.approx(ceiling, rel=1e-3) assert gaps.min() == pytest.approx(ceiling, rel=1e-3) assert ceiling < globe.radius * 0.001, "and the lift itself stays small" def test_the_mesh_elevation_is_the_drawn_surface_not_the_smooth_one(): """ The two differ by the cell's twist, which over MOLA reaches 1977 m - twenty kilometres once exaggerated, and the reason everything laid on the globe used to be lifted so far clear of it. Checked against the test's own reading of `build_surface`'s triangles, written independently of the one in `geometry`. """ globe, relief = a_real_globe() rng = np.random.default_rng(4) lats = rng.uniform(-89.4, 89.4, 2000) lons = rng.uniform(-180.0, 180.0, 2000) on_mesh = globe.mesh_elevation_at(lats, lons) assert on_mesh == pytest.approx(_mesh_elevation(relief, lats, lons), abs=1e-6) # And it is genuinely a different answer from the smooth surface assert np.abs(on_mesh - globe.elevation_at(lats, lons)).max() > 100.0 def test_the_mesh_and_the_smooth_surface_agree_at_the_nodes(): """ Every node of the grid is a corner of four triangles and a knot of the bilinear surface, so the two readings have to meet there. """ globe, _ = a_real_globe() grid_lats, grid_lons, _, _ = render_grid() lat2d, lon2d = np.meshgrid(grid_lats[::7], grid_lons[::11], indexing='ij') assert globe.mesh_elevation_at(lat2d, lon2d) == pytest.approx( globe.elevation_at(lat2d, lon2d), abs=1e-6) # --- colours that do not collide ------------------------------------------- def test_truncated_uses_only_the_range_it_was_given(): """ Narrowing the colours a scale is drawn in is not the same as rescaling the data, so the bar still reads in real units. """ import matplotlib.pyplot as plt from dispnc.colors import truncated base = plt.get_cmap('viridis') part = truncated('viridis', 0.45, 1.0) assert part(0.0) == pytest.approx(base(0.45), abs=0.01) assert part(1.0) == pytest.approx(base(1.0), abs=0.01) # and it takes a Colormap as readily as a name, since --cmap passes either assert truncated(part, 0.0, 1.0)(0.5) == pytest.approx(part(0.5), abs=0.01) def test_the_terrain_under_a_shell_carries_no_hue_to_compete_with_it(): """ It used to be `bone`, which runs black to white through blue - exactly where viridis's dark end lives - so an isosurface and the globe beneath it came out the same colour and the picture was one purple-grey mush. """ from dispnc.render.globe import TERRAIN_COLORMAP for position in (0.0, 0.25, 0.5, 0.75, 1.0): red, green, blue = TERRAIN_COLORMAP(position)[:3] assert red == pytest.approx(green, abs=0.01), "a grey has no hue" assert green == pytest.approx(blue, abs=0.01) # And it stops short of both ends, leaving black and white to everything else assert TERRAIN_COLORMAP(0.0)[0] > 0.05 assert TERRAIN_COLORMAP(1.0)[0] < 0.95 def test_the_outermost_isosurface_is_not_the_darkest_colour_on_the_scale(): """ With several thresholds the scale is stretched across them, which puts the *outer* envelope at the bottom of the range - and the outer envelope is the one in view. On the full colormap that made it a near-black purple. """ import matplotlib.pyplot as plt from dispnc.render.globe.shell import build_shell plan = a_cloud_plan() cube = np.zeros((26, 9, 12)) cube[10:16, 3:7, 4:8] = 1.0 cube[12:14, 4:6, 5:7] = 5.0 shell = build_shell(a_globe(), plan, cube, 'viridis', 'dust', mode='iso', threshold='0.5,2,4') outer = shell.actors[0] # The surface carries one constant scalar; its colour is what the lookup # table makes of that value value = float(np.asarray(outer.pointdata['Scalars'])[0]) colour = np.zeros(3) outer.mapper.GetLookupTable().GetColor(value, colour) darkest = np.asarray(plt.get_cmap('viridis')(0.0)[:3]) assert colour.sum() > darkest.sum() * 1.5, "the envelope has to be visible" # And it has a hue, which is what actually tells it from the ground: the # terrain is grey at every point of its scale, so any colour at all reads as # "this is the shell". Brightness would not do - the pale end of the terrain # is brighter than a mid-scale isosurface and still unmistakably not one. from dispnc.render.globe import TERRAIN_COLORMAP assert colour.max() - colour.min() > 0.2, "the surface must not be grey" for position in (0.0, 0.5, 1.0): ground = np.asarray(TERRAIN_COLORMAP(position)[:3]) assert ground.max() - ground.min() < 0.02