""" The atmosphere as a cloud rather than as a surface: volume ray casting. A contour is a decision - this much and no less - and it draws a hard edge where the atmosphere has none. What reads as cloud is the field accumulated *along the view ray*, which is what a volume mapper does and what no arrangement of surfaces can imitate: concentric translucent shells all share the globe's centre, so VTK's centroid sort cannot even put them in the right order. Ray casting needs a uniform box, so this is the one place where the spherical grid is resampled into cartesian voxels. `shell.py` contours the model's own grid instead, which is the right choice there and the wrong one here. Transparency is driven by the same scale as the colour, so the two say the same thing: on a field drawn on a log scale the haze fades the way the colour bar says it does, and a field spanning decades shows more than its peak. """ import numpy as np from matplotlib.colors import Normalize from .geometry import add_bar, finite_range try: from vedo import Volume volume_available = True except ImportError: # pragma: no cover volume_available = False try: from scipy.interpolate import RegularGridInterpolator scipy_available = True except ImportError: # pragma: no cover scipy_available = False # Voxels along each side of the sampled box. 256 puts the whole globe in 67 MB # and takes about 0.7 s to fill; the shell itself is only a sixth of that box, # so most of it is empty air the sampler skips. CLOUD_RESOLUTION = 256 # The distance an opacity from the curve below is worth, as a multiple of the # atmosphere's own exaggerated depth. VTK reads a scalar opacity as "this much # over one unit distance" and corrects each sample for how far it actually # stepped, so this is a length rather than a voxel count: tied to the voxel size # it would be --shell-resolution, not the atmosphere, that decided how thick the # cloud looks. At 1.0 a ray crossing the whole atmosphere at full strength picks # up exactly --shell-opacity, and a ray grazing the limb crosses more and comes # out denser, which is what an atmosphere does. ALPHA_UNIT_DEPTHS = 1.0 # Points defining the opacity curve over the field's range. ALPHA_STOPS = 12 # Stops in a cross-section's opacity ramp. Only the first is transparent, so the # band of values it swallows is one part in this many of the whole range. SLICE_ALPHA_STOPS = 64 def available(): """ Whether a cloud can be built at all. Both parts are optional in practice: vedo is an optional dependency of the toolbox, and scipy does the resampling. """ return volume_available and scipy_available def display_scale(scale, cube): """ (transform, vmin, vmax, normalized): how the field is put into the box. VTK's volume property maps a scalar linearly onto colour and opacity and offers nothing else, so a linear scale is handed over as limits and the bar then carries the field's own numbers. Any other scale - a log one above all, which is what a field spanning decades is drawn on - is applied here instead, and the bar's title carries the real range because the bar itself runs over 0..1. This is the same split `apply_colormap` makes for a mesh. """ low, high = finite_range(cube) if scale is None or type(scale) is Normalize: vmin = low if scale is None or scale.vmin is None else float(scale.vmin) vmax = high if scale is None or scale.vmax is None else float(scale.vmax) if vmax <= vmin: vmax = vmin + 1.0 return None, vmin, vmax, False return scale, 0.0, 1.0, True def normalize(transform, cube): """ The field put on 0..1 by `transform`, or handed back as it is when there is none. The cube is flattened for the call: matplotlib's log transform takes a flat array or a column and rejects anything else, and every other caller in the toolbox happens to hand it per-vertex scalars, which are flat already. """ values = np.asarray(cube, dtype=float) if transform is None: return values scaled = np.ma.filled(np.ma.masked_invalid(transform(values.ravel())), np.nan) return np.clip(np.asarray(scaled, dtype=float), 0.0, 1.0).reshape(values.shape) def sample_volume(globe, altitudes, lats, lons, cube, transform, empty, resolution=CLOUD_RESOLUTION): """ (voxels, spacing, origin): the field resampled into a cartesian box. Everywhere there is no atmosphere - outside the level span, off the field's own latitude and longitude range, or where the data itself is missing - the box holds `empty`, the bottom of the scale, which the opacity curve pins to fully transparent. Empty air and an empty value being the same thing is what makes the cloud fade out at its edges instead of ending at one. """ values = normalize(transform, cube) if not np.any(np.isfinite(values)): return None interp = _interpolator(altitudes, lats, lons, values) low, high = float(altitudes[0]), float(altitudes[-1]) # The box holds the whole exaggerated shell, so its half-side is the top of # the atmosphere as the globe actually draws it half = globe.radius + high * globe.air_exaggeration axis = np.linspace(-half, half, int(resolution)) spacing = float(axis[1] - axis[0]) if axis.size > 1 else 2.0 * half voxels = np.full((axis.size, axis.size, axis.size), float(empty), dtype=np.float32) x, y = np.meshgrid(axis, axis, indexing='ij') flat_radius = x * x + y * y # One z-slab at a time: the coordinates of a whole 256 cube would be 400 MB # of doubles before a single value had been interpolated for k, z in enumerate(axis): r = np.sqrt(flat_radius + z * z) lat = 90.0 - np.rad2deg(np.arccos(np.clip(z / np.where(r > 0, r, 1.0), -1.0, 1.0))) lon = np.rad2deg(np.arctan2(y, x)) # The inverse of to_cartesian(..., altitude=): the relief comes off # again, so the cloud sits on the terrain the shells are lifted above height = globe.altitude_at(lat, lon, r) inside = (height >= low) & (height <= high) if not np.any(inside): continue sampled = interp(np.column_stack([height[inside], lat[inside], lon[inside]])) voxels[:, :, k][inside] = np.nan_to_num(sampled, nan=float(empty)) if not np.any(voxels > empty): return None return voxels, spacing, (-half, -half, -half) def _interpolator(altitudes, lats, lons, values): """ A (altitude, latitude, longitude) interpolator over the field. Latitudes usually run north to south, which a RegularGridInterpolator rejects, so both horizontal axes are sorted ascending and the data follows them. The longitude axis is then closed by repeating its first column a full turn later: without that the sampler falls off the end of the axis between the last longitude and the first, and the cloud comes out with a slit down one meridian. """ altitudes = np.asarray(altitudes, dtype=float) lats, lat_order = np.unique(np.asarray(lats, dtype=float), return_index=True) lons, lon_order = np.unique(np.asarray(lons, dtype=float), return_index=True) values = values[:, lat_order, :][:, :, lon_order] if lons.size > 1: lons = np.append(lons, lons[0] + 360.0) values = np.concatenate([values, values[:, :, :1]], axis=2) return RegularGridInterpolator((altitudes, lats, lons), values, bounds_error=False, fill_value=np.nan) def alpha_curve(opacity, stops=ALPHA_STOPS): """ The opacity curve over the field's range, evenly spaced from bottom to top. Zero at the bottom, so empty air stays invisible whatever the opacity, and straight from there: transparency is meant to say exactly what the colour says, and the scale has already been bent to whatever shape the field needs. `opacity` is used as it is rather than scaled up the way the surface modes scale it. A surface's alpha is applied once, where the ray crosses it; this one is applied all the way through an atmosphere, so the same number here buys a great deal more. """ return list(np.linspace(0.0, 1.0, stops) * min(1.0, float(opacity))) def build_slicer(globe, altitudes, lats, lons, cube, colormap, scale, resolution=None): """ A callable `(origin, normal) -> Mesh`: the field on any plane through it. This is what makes a cutaway say something. Clipping alone opens the scene but reveals nothing, because a stack of shells is a stack of *surfaces* - there is no interior to expose, only the far shells seen from inside. A cross-section is the thing that actually carries the vertical structure: one face, altitude up it, coloured by the field. The volume is sampled once and re-sliced per call, which is what makes the plane draggable - a slice costs about ten milliseconds against the seconds the sampling took. """ if not available(): return None altitudes = np.asarray(altitudes, dtype=float) if altitudes.size < 2 or altitudes[-1] <= altitudes[0]: return None transform, vmin, vmax, _ = display_scale(scale, cube) sampled = sample_volume(globe, altitudes, lats, lons, cube, transform, vmin, int(resolution or CLOUD_RESOLUTION)) if sampled is None: return None voxels, spacing, origin_xyz = sampled volume = Volume(voxels, spacing=(spacing,) * 3, origin=origin_xyz) # The plane spans the whole box, most of which is inside the planet or above # the atmosphere; those voxels hold the bottom of the scale, so dropping the # bottom of the scale leaves just the annulus the field actually occupies. # This has to be a ramp over the *scalar range* rather than a value per # vertex: handed a sequence, vedo spreads it across vmin..vmax, so a # per-vertex array comes out as a lookup table in point order - which is to # say, noise, and mostly invisible. fade = [0.0] + [1.0] * (SLICE_ALPHA_STOPS - 1) def make_slice(origin, normal): mesh = volume.slice_plane(origin=list(origin), normal=list(normal)) if mesh.npoints == 0: return None mesh.cmap(colormap, np.asarray(mesh.pointdata['ImageScalars'], dtype=float), vmin=vmin, vmax=vmax, alpha=fade) mesh.lighting('off') return mesh return make_slice def build_cloud(globe, altitudes, lats, lons, cube, colormap, scale, title, shell_class, opacity=0.35, lit=False, bar_pos=None, resolution=None): """ A `Shell` holding the field as one translucent cloud, or None when it cannot be built and the caller should fall back on contouring. `shell_class` is `shell.Shell`, handed in rather than imported: `shell` imports this module, and the dataclass is the only thing needed back. """ if not available(): return None # One level is a surface, not a volume: there is no depth to accumulate # along and nothing to interpolate between, and `iso` draws it properly altitudes = np.asarray(altitudes, dtype=float) if altitudes.size < 2 or altitudes[-1] <= altitudes[0]: return None resolution = int(resolution or CLOUD_RESOLUTION) transform, vmin, vmax, normalized = display_scale(scale, cube) sampled = sample_volume(globe, altitudes, lats, lons, cube, transform, vmin, resolution) if sampled is None: return None voxels, spacing, origin = sampled depth = (float(altitudes[-1]) - float(altitudes[0])) * globe.air_exaggeration volume = Volume(voxels, spacing=(spacing, spacing, spacing), origin=origin) volume.cmap(colormap, vmin=vmin, vmax=vmax) volume.alpha(alpha_curve(opacity), vmin=vmin, vmax=vmax) volume.alpha_unit(depth * ALPHA_UNIT_DEPTHS) volume.mode(0) # composite: accumulate along the ray volume.interpolation(1) # linear, which is what makes it smooth volume.jittering(True) # or the ray steps show as wood grain volume.shade(bool(lit)) low, high = finite_range(cube) bar_title = f"{title} ({low:.3g} to {high:.3g})" if normalized else title add_bar(volume, bar_title, bar_pos, use_alpha=False) detail = f"volume cloud at {resolution}^3, {low:.3g} to {high:.3g}" return shell_class(actors=[volume], note=detail, caption=[f"cloud, opacity {opacity:g}"], depth_peeling=False)