# display_netcdf Guide

Plot numeric variables from NetCDF files produced by the PEM, the GCM and XIOS:
scalars, curves, vertical profiles, latitude/longitude maps, cross-sections,
vector fields, polar-stereographic views, animations, and an interactive 3D globe.

`python display_netcdf.py --help` is the authoritative option list; this page is
the practical entry point.

## Contents

| Path | Purpose |
| --- | --- |
| `display_netcdf.py` | Entry point (a thin shim over the `dispnc` package) |
| `dispnc/` | The implementation |
| `tests/` | Unit tests and the regression baseline |
| `display_netcdf.yml` | Conda environment used only by this tool |
| `MOLA_1px_per_deg.npy` | MOLA topography, 1 degree per pixel, for contour overlays and the globe |
| `molaTeam_contour_31rgb_steps.csv` | MOLA colormap steps; kept for reference, not read by the code |
| `sample_files/` | Sample files covering every supported plot type — **not versioned**, see below |

## Setup

```bash
cd PEM/toolbox/display_netcdf
conda env create -f display_netcdf.yml
conda activate disp_nc
```

Required: `netCDF4`, `numpy`, `matplotlib`, `cartopy`, `scipy`, `xarray`, `dask`.
Optional: `vedo`, needed only for `--show-3d`; `pytest`, to run the tests.

Data files are located relative to the script, so it can be run from any
directory.

## How a plot is chosen

Every dimension is given a **role** — X, Y, Z, T, unstructured, or none — and the
roles that survive slicing decide the plot. Roles come from the file's own
metadata where it exists, and from dimension names only where it does not:

| Evidence | Confidence | Typical source |
| --- | --- | --- |
| `axis = "X"/"Y"/"Z"/"T"` | 100 | XIOS output |
| CF `standard_name` | 90 | XIOS output |
| `positive = "up"/"down"` | 80 | `soildepth`, `altitude` |
| `units` (`degrees_north`, …) | 70 | anything, after normalization |
| Dimension/coordinate name | 40 | PEM and GCM files, `rlonu`/`rlatv` |
| Nothing | 10 | `nslope`, `index`, `descriptor` |

| Roles remaining | Plot |
| --- | --- |
| none | Scalar printed to the terminal |
| T | Time series |
| Z | Vertical profile — value on X, depth on Y, increasing downward |
| X, Y, or unstructured | Line plot |
| X + Y | Map, in a real projection, eligible for polar and 3D views |
| anything + Z | Cross-section, vertical coordinate on Y |
| T + X or T + Y | Hovmöller diagram |
| X + Y + Z up | 3D atmospheric shell, with `--plot-kind globe` |
| three or more | Refused, listing what is left |

`--explain` prints all of this for a given file and variable: the detected
writer, every metadata repair applied, each dimension's role and where it came
from, and the colormap decision. It is the first thing to reach for when a plot
looks wrong.

## Reading files that disagree

The three writers agree on almost nothing, so metadata is normalized on open:

- **PEM** (`diagevo.nc`, `startevo.nc`) writes only `title` and `units`, both
  free text — `"Planetary year"`, `"Degree East-West"`, `"Radian North-South"`.
- **XIOS** (`xoutdaily4pem_y*.nc`) is properly CF, but puts a sentence in
  `soildepth:standard_name`, so `standard_name` is only trusted against the CF
  list.
- **GCM/LMDZ** (`diagfi.nc`, `start.nc`, `startfi.nc`) has `long_name` on
  coordinates but no `axis`; data variables use `title`; `Time` claims
  `"days since 0000-00-0 00:00:00"`, which is not a date; `startfi.nc` stores
  latitude and longitude in radians with no units at all.

Normalization promotes `title` to `long_name`, rewrites units into CF spellings,
converts radians to degrees, repairs the unusable time units, and stamps `axis`
attributes. Radian conversion inferred from magnitude alone is always reported
under `--explain`, since it is the one repair that could mislead.

Time axes are never handed to a date library: a Mars year is not 365 days, so
decoding would be wrong rather than merely unsupported. Time stays numeric.

## Command-line mode

```bash
python display_netcdf.py FILE.nc -v VARIABLE [options]
```

| Option | Meaning |
| --- | --- |
| `-v`, `--variable` | Variable to plot |
| `-e`, `--extra-indices` | JSON object selecting dimensions (see below) |
| `-o`, `--output` | Write the figure here instead of displaying it |
| `-c`, `--cmap` | Colormap, or `auto` (default). `-c jet` restores the old look |
| `--norm` | `linear`, `log`, `symlog`, `centered` |
| `--vmin`, `--vmax` | Fix the colour limits |
| `-x`, `--x-dim` | Dimension on the X axis of a 2D plot |
| `--plot-kind` | Force a renderer instead of inferring one |
| `--anomaly` | `time` or `zonal`: plot the departure from that mean |
| `--diff` | `VAR`: plot the departure from another variable, centred on zero |
| `--reduce` | `OP:DIM[,DIM]`, repeatable. `mean`, `std`, `min`, `max`, `sum` |
| `--stats` / `--stats-only` | Print summary statistics, with or without a figure |
| `--overlay` | `VAR[:SPEC]`, repeatable: draw another variable with this one. `SPEC` is a colormap when the layer is a fill, a colour when it is a line |
| `--overlay-style` | How they share the figure: `blend` (default), `contour`, `hatch`, `panels`, `bivariate`, `glyph` — see below |
| `--overlay-threshold` | Contour levels, or the level a hatch is drawn above: values or `pNN` |
| `--overlay-scale` | Under `panels`: `own` (default) or `shared` colour scale |
| `--vector` | `U,V`: plot a vector field from two named components |
| `--vector-style` | `quiver` (default) or `stream` |
| `--vector-density` | Arrow count per axis, or streamline density |
| `--vector-background` | `magnitude` (default), `none`, or a variable name |
| `--no-remap` | Keep the raw physics grid instead of rebuilding a map |
| `--explain` | Report how the file was read and the plot chosen |
| `--list-vars` | List variables with their dimensions and exit |
| `--no-show-topo` | Drop the topography overlay |
| `--show-polar` | Also produce north and south polar views |
| `--show-3d` | Also open the 3D globe (needs `vedo`) |
| `--animate` | Step a map, its polar views and the globe through a dimension |
| `--spin` | `[TURNS]`: rotate the 3D globe camera (default: 1 turn) |
| `--fps` | Frames per second for animated output (default: 12) |
| `--figsize` | `W,H` in inches, overriding each plot kind's own default |
| `--dpi` | Resolution of saved output; the figure on screen is left alone |
| `--title` | Title the figure exactly like this |

### Animation

`--animate DIM` holds one dimension back from the plot and steps through it
instead. It applies to lat/lon maps, their polar views and the 3D globe; a
curve, a profile or a cross-section has no spare axis to spend on time and says
so rather than silently showing step zero.

Without `-o` each figure gets a slider along the bottom. With `-o` and a movie
extension (`.mp4`, `.gif`, `.avi`, `.mov`, `.mkv`, `.webm`) the frames are
written as a video: the map to that path, the polar views to
`<stem>_north<ext>` and `<stem>_south<ext>`, the globe to `<stem>_globe<ext>`.

```bash
# a Martian year of surface pressure as a movie, with both poles beside it
python display_netcdf.py sample_files/xoutdaily4pem_y1.nc -v ps_avg \
    --animate time_counter --fps 24 --show-polar -o year.mp4

# the same year on the globe, scrubbable with a slider
python display_netcdf.py sample_files/xoutdaily4pem_y1.nc -v ps_avg \
    --show-3d --animate time_counter
```

The colour scale is fixed over every frame, so the colours mean the same thing
throughout instead of sliding with each frame's own range. Frames are titled
with the coordinate value — `time_counter = 3641.67 days` — not the index.
Video needs `ffmpeg`, which `display_netcdf.yml` installs; `.gif` falls back to
Pillow when it is missing.

### Turning the globe

`--spin` rotates the camera around the planet. A fixed camera cannot tell a
feature that moves with the planet from one that moves against it, and half a
sphere is always facing away.

```bash
# one turn, written straight out as a movie
python display_netcdf.py sample_files/start.nc -v h2o_ice \
    --plot-kind globe -e '{"Time":0}' --spin -o turn.mp4

# two turns while a Martian year steps past underneath
python display_netcdf.py sample_files/xoutdaily4pem_y1.nc -v ps_avg \
    --plot-kind globe --animate time_counter --spin 2 --fps 24 -o year.mp4
```

On screen it turns by itself and stays draggable; the mouse is not held. To a
file it becomes the animation when there is nothing else to step through, and
rides on top of `--animate` when there is. One turn lasts eight seconds, so
`--fps` makes it smoother rather than faster.

### The 3D globe

`--show-3d` opens the globe after a map; `--plot-kind globe` makes it the figure,
which is what lets it keep a vertical axis and draw an atmosphere.

| Option | Meaning |
| --- | --- |
| `--globe-relief` | What shapes the globe: `topo` (MOLA, default), `none`, `data`, or a variable name |
| `--globe-surface` | What the sphere shows: `field`, `terrain` (greys), or `none`. Defaults to `terrain` under a shell |
| `--light` / `--no-light` | Shade the globe. Off by default; `--sun` turns it on |
| `--sun` | `LON,LAT`: light the globe from there and draw the terminator |
| `--sun-ls` | Sub-solar latitude from a solar longitude, in degrees |
| `--globe-export` | Write the mesh to `.vtk`, `.vtp`, `.ply`, `.obj` or `.stl` |
| `--shell-mode` | `cloud` (default), `iso`, `column`, or `layers` — see below |
| `--shell-threshold` | Isosurface level(s), comma-separated: a value or `pNN` for a percentile (default: `p99`) |
| `--shell-color` | `value` (default) colours an isosurface on the field's scale; `height` by altitude |
| `--shell-level` | Draw a single vertical level as a solid shell |
| `--shell-opacity` | Opacity of the cloud or the shell surfaces (default: 0.35) |
| `--shell-cut` | Slice the shell open with a plane you can drag. On by default for `layers` |
| `--shell-resolution` | Voxels per side of the box `cloud` is ray-cast from (default: 256) |
| `--shell-top` | Where the shell tops out when the vertical is not a length (default: 50 km) |
| `--shell-exaggeration` | Vertical exaggeration of the air (default: 30; the relief keeps its own 10) |

A corner caption names the variable, what shapes the sphere, what the shell is
and at what threshold, the altitude span and the exaggeration — the choices the
picture cannot show by itself. Clicking the globe reports longitude, latitude,
elevation and the value there, which is what `format_coord` gives the flat maps.

#### Showing an atmospheric field

Four modes, because no one of them answers every question about a volume:

| Mode | Shows | Good for |
| --- | --- | --- |
| `cloud` | The field ray-cast as a translucent volume (default) | What does the atmosphere actually look like? |
| `iso` | One or more isosurfaces | Where does the field exceed this? The shape of a cloud or a plume |
| `column` | The vertical integral on one translucent shell | Where is there any of it at all? |
| `layers` | One surface per level, cut open | How is it stacked vertically? |

```bash
# the water ice as a cloud: the default, so no shell flags at all
python display_netcdf.py sample_files/start.nc -v h2o_ice \
    --plot-kind globe -e '{"Time":0}'

# where is it, at all: the vertical integral on one shell
python display_netcdf.py sample_files/start.nc -v h2o_ice \
    --plot-kind globe -e '{"Time":0}' --shell-mode column

# its envelope, body and core at once, each in its threshold's colour
python display_netcdf.py sample_files/start.nc -v h2o_ice \
    --plot-kind globe -e '{"Time":0}' --shell-mode iso --shell-threshold p90,p99,p99.9
```

`cloud` accumulates the field along each view ray instead of drawing a surface
through it, so the atmosphere fades out at its edges rather than ending at one —
which is what an isosurface, being a decision about a threshold, cannot do. Its
transparency follows the same scale as its colour, so a field drawn on a log
scale fades on a log scale too and shows more than its peak. `--shell-opacity`
makes it denser or fainter; `--shell-resolution` makes it sharper, at the cost
of memory as the cube of N, without changing how dense it looks. It needs
`scipy` and a vedo with `Volume`, and falls back to `iso` with a note when
either is missing.

Under a shell the sphere goes neutral grey terrain by default, so the data
colours belong to the atmosphere alone; `--globe-surface field` restores the
older look where the sphere carries the field's bottom level. Nested
isosurfaces are coloured over the thresholds' own range rather than the whole
field's: `h2o_ice` spans thirty-five decades, so p90, p99 and p99.9 would
otherwise be the same yellow.

An isosurface is one flat colour over a whole body, so shape is all it has to
say — which is why `iso` is the one thing on the globe that lights itself,
without waiting for `--light`. Its surfaces are relaxed first, or the facets of
a 33x33x26 model grid would be all the lighting picked out. `--light` still
shades the globe, and the cloud with it.

`column` is a map rather than an atmosphere: the integral has no altitude of its
own, so it is drawn as a sphere — not a shell over the terrain, which gave the
map a ripple that came from the topography rather than from the field — just
high enough to clear the highest ground on the planet, about 7 km of air and
200 km once exaggerated, and interpolated onto the same 1° grid the globe is
meshed on. What it integrates is mass, `Σ q·Δp/g` from the hybrid coefficients and the
surface pressure, so a mixing ratio comes out in kg/m²; without them it falls
back on `∫q·dz` and says so. Its scale follows the plotted field's — log where
that is log — because on a linear stretch a field spanning thirty decades puts
everything but its brightest cell at zero; a log scale then gets a floor under
it, since the bottom of `h2o_ice`'s integral is a numerical zero thirty decades
below anything on the map and a scale nothing is ever near is not a scale.

`layers` is cut open, because a stack of concentric shells cannot be seen into:
the outermost covers every other, and VTK cannot even order them, since sorting
translucent actors by centroid puts every shell at the same distance. Cutting
alone would expose nothing — surfaces have no interior — so the cut carries a
cross-section of the field with it, altitude up its face. Drag the handle to
sweep it around. `--shell-cut` turns the same thing on for a cloud, which is
worth seeing from the inside.

### Several variables at once

`--overlay VAR[:SPEC]` draws another variable with the plotted one. It is
repeatable, and `--overlay-style` decides *how* they share the figure — which
matters more than it looks, because the styles answer different questions.

| Style | What it does | Reach for it when |
| --- | --- | --- |
| `blend` (default) | Translucent colour, mixing where they overlap | The overlap itself is the subject |
| `contour` | Labelled lines over the filled base | **Both fields have to stay readable** |
| `hatch` | A pattern above a threshold | The question is presence, not amount |
| `panels` | One panel per variable | They will not share a scale — works for *every* plot kind |
| `bivariate` | Both in one 2-D colour scheme, with a square key | Exactly two, and the combination is the point |
| `glyph` | The layer as sized symbols | The layer is sparse and patchy |
| `--diff VAR` | One field, the difference | The relationship is the real question |

```bash
# the classic: a shaded field under labelled contours of another
python display_netcdf.py sample_files/xoutdaily4pem_y1.nc -v tsoil_avg \
    --overlay ps_avg -e '{"time_counter":0,"soildepth":0}' --overlay-style contour

# where the CO2 frost gives way to the water ice, as translucent layers
python display_netcdf.py sample_files/start.nc -v h2o_ice \
    --overlay co2_ice -e '{"Time":0,"altitude":0}'

# where there is any of it at all, in a pattern that survives a photocopier
python display_netcdf.py sample_files/xoutdaily4pem_y1.nc -v ps_avg \
    --overlay h2o_layer1_avg -e '{"lon":0}' --overlay-style hatch

# three variables that share no scale, one panel each
python display_netcdf.py sample_files/xoutdaily4pem_y1.nc -v ps_avg \
    --overlay h2o_layer1_avg --overlay co2_layer1_avg \
    -e '{"time_counter":0}' --overlay-style panels

# and the question two layers were often standing in for
python display_netcdf.py sample_files/start.nc -v h2o_ice \
    --diff co2_ice -e '{"Time":0,"altitude":0}'
```

**Why there is a choice at all.** Blending is translucent where each field is
weak, so where two overlap the colours mix and the overlap is visible as an
overlap. That buys something real — coexistence shows at all, which drawing the
top layer opaque would hide by construction — and it costs something real: a
cell holding both CO2 and H2O comes out a colour that is in neither bar, so the
figure cannot be read back to a number exactly where it is most interesting. It
also spends the base variable's colormap, since every layer has to drop to a
single-hue ramp to stay separable.

`contour` pays neither cost. The overlay becomes *geometry* rather than paint,
so the base keeps `viridis` and both fields stay readable everywhere, including
where they cross; three sets of lines in three colours are still legible where
three washes of alpha are mud. It is how a shaded field under contoured heights
has been drawn in atmospheric science for decades, and it is the one to reach
for unless the overlap is specifically what you came to see.

`bivariate` is the principled version of mixing: both fields are cut into three
classes, every pair of classes gets one colour, and the legend is the whole
square — so every colour on the map has a cell a reader can point at. Quantized,
and exactly two variables; three have no square.

`hatch` encodes **presence, not amount**. A reader who expects a magnitude from
a hatch will misread it, which is the price of being the only style that
survives greyscale printing.

**Curves.** A 1-D plot has no colour channel to share, so `--overlay` on a line,
a time series or a profile draws several curves with a legend. Units decide the
axes: one unit is one axis, two get a twin axis opposite it — labelled, with
`[right axis]` on the legend entry — and past that there is no third side to the
figure, so the odd ones out are refused by name and `--overlay-style panels` is
suggested. Units are compared as written; a variable with no units is its own
group, because "unknown" cannot be claimed to match "K".

```bash
# two curves, one legend, a twin axis for the second unit
python display_netcdf.py sample_files/diagevo_1d.nc -v tsurf --overlay ps_avg_glob
```

**Which kinds stack.** Maps, cross-sections, the 1-D curves and the atmospheric
shell all carry several variables in one frame. A scalar, the surface globe and
`--vector` do not, and say so rather than dropping the layers in silence —
`--overlay-style panels` works for all of them.

Layers are sliced by the same `-e`, `--reduce` and `--animate` as the plotted
variable, so they cannot end up showing a different time step, and every layer
advances with it under `--animate`. One that cannot be laid over the same grid,
or that is on the same shape but different dimensions, is named and skipped, and
the figure is still drawn.

### Selecting dimensions

`--extra-indices` maps a dimension to a **0-based** index, or to `"avg"` to
average over it. Size-1 dimensions are selected automatically; anything left out
is kept whole.

Interactive prompts are **1-based** (they offer `1–<size>`) while
`--extra-indices` is 0-based. Interactive `1` and CLI `0` select the same slice.

Averaging is streamed: a mean over 669 time steps of a multi-gigabyte file uses
memory proportional to one chunk, not to the file.

### Examples

```bash
# lat/lon map of surface temperature
python display_netcdf.py sample_files/diagevo.nc -v tsurf -e '{"Time":0}' -o tsurf.png

# what did the tool make of this file?
python display_netcdf.py sample_files/startfi.nc -v tsurf --explain

# XIOS yearly output with north/south polar views
python display_netcdf.py sample_files/xoutyearly4pem_y1.nc -v co2_ice_min \
    -e '{"time_counter":0}' --show-polar -o co2.png

# time-averaged surface pressure over 669 days
python display_netcdf.py sample_files/xoutdaily4pem_y1.nc -v ps_avg -e '{"time_counter":"avg"}'

# departure of one day from the annual mean, on a diverging scale
python display_netcdf.py sample_files/xoutdaily4pem_y1.nc -v ps_avg \
    --anomaly time -e '{"time_counter":300}'

# zonal-mean soil temperature over depth and latitude
python display_netcdf.py sample_files/xoutdaily4pem_y1.nc -v tsoil_avg \
    --reduce mean:lon --reduce mean:time_counter

# statistics only, area-weighted with the file's cell_area
python display_netcdf.py sample_files/diagevo.nc -v tsurf -e '{"Time":0}' --stats-only

# soil temperature cross-section
python display_netcdf.py sample_files/xoutdaily4pem_y1.nc -v tsoil_avg -e '{"time_counter":0,"lon":0}'

# unstructured physics grid, remapped to a map automatically
python display_netcdf.py sample_files/startfi.nc -v tsurf

# vertical profile at one physics grid point
python display_netcdf.py sample_files/startfi.nc -v tsoil -e '{"physical_points":0}'

# vector field from two components on the same grid
python display_netcdf.py winds.nc --vector u,v --vector-style stream

# the four ways of showing one atmospheric field, side by side
for m in cloud iso column layers; do
  python display_netcdf.py sample_files/start.nc -v h2o_ice \
      --plot-kind globe -e '{"Time":0}' --shell-mode $m -o shell_$m.png
done

# a cloud cut open, with the cross-section on the cut face
python display_netcdf.py sample_files/start.nc -v h2o_ice \
    --plot-kind globe -e '{"Time":0}' --shell-cut

# a lit globe with the terminator, at northern summer solstice
python display_netcdf.py sample_files/start.nc -v h2o_ice \
    --plot-kind globe -e '{"Time":0}' --sun-ls 90 --light

# the same field shaped by its own values instead of by the terrain
python display_netcdf.py sample_files/start.nc -v h2o_ice \
    --plot-kind globe -e '{"Time":0}' --globe-relief data --globe-surface field

# coarser and quicker, for a first look at a big file
python display_netcdf.py sample_files/start.nc -v h2o_ice \
    --plot-kind globe -e '{"Time":0}' --shell-resolution 128

# two ices on one map, then the same pair as cloud decks
python display_netcdf.py sample_files/start.nc -v h2o_ice \
    --overlay co2_ice -e '{"Time":0,"altitude":0}' -o ices.png
python display_netcdf.py sample_files/start.nc -v h2o_ice \
    --overlay co2_ice --plot-kind globe -e '{"Time":0}' --spin -o ices.mp4
```

Exit status is `0` on success, `1` if the file, variable or plot type is
unusable, and `2` on bad arguments, so the tool composes with shell loops and
`multiple_exec.sh`.

## Interactive mode

```bash
python display_netcdf.py                     # prompts for the file too
python display_netcdf.py sample_files/diagevo.nc -c RdBu_r
```

Prompts for the file, then loops over variables, asking an index per dimension
(`a` or `average` averages, Enter or `all` keeps it whole). Every option the
command line takes is honoured inside the loop — `--norm`, `--vmin`/`--vmax`,
`--reduce`, `--anomaly`, `--stats`, `--show-polar`, `--show-3d`, `--overlay` and
the rest — with one exception: `-o` is not, because one path cannot serve a loop
over however many variables you ask for.

Two consequences worth knowing. `--reduce mean:lon` on a variable that has no
`lon` is reported and that variable is skipped, rather than silently plotted
unreduced; the loop carries on. And every figure stays on screen until you close
it, so two variables can be put side by side; answer `close` at the prompt to
clear them all at once, and Enter to quit, which closes what is left.

Tab-completion works for paths, for variable names, and — at the "which
dimension on X?" question — for exactly the dimensions that are valid answers
there. That question also accepts the coordinate name a file advertises, so
`longitude` and `lon` both select the same axis, and rejects anything else out
loud instead of quietly falling back on the default.

Yes/no questions take the whole word as well as the letter: `y`, `yes`, `oui`,
`n`, `no`, `non`, in any case. Enter takes the default shown in the prompt.

## Vector fields

`--vector U,V` names both components explicitly, because the naming convention
differs across the models this tool reads.

Both components must share a grid. `start.nc` holds `ucov(Time, altitude,
latitude, rlonu)` and `vcov(Time, altitude, rlatv, longitude)` — a real pair on a
staggered Arakawa C-grid, living on different dimensions. Drawing them together
needs de-staggering, which this tool does not do, so such a pair is refused with
a message naming both grids rather than plotted incorrectly.

## Example files

`sample_files/` is **not in SVN**: it is 62 MB of NetCDF, most of it one
file, which does not belong in the repository. Ask the toolbox maintainer for a
copy, or point the tool at your own output — nothing in `dispnc/` needs these
files. Without them the unit tests skip the cases that read a real file, and
`baseline_capture.py` says so and exits.

Any set of files works: the directory may hold more, fewer or differently named
files than the table below, and the tests adapt — see [Tests](#tests).

<!-- TODO: replace with the agreed location once one is chosen. -->

The set the table describes is the one this page was written against:

| File | Grid | Notable for |
| --- | --- | --- |
| `diagevo.nc` | 33×33 lat/lon, 1 time step | PEM diagnostics; maps and soil cross-sections |
| `diagevo_1d.nc` | single column, 275 steps | Long time series, soil profiles |
| `diagfi.nc` | 33×33 lat/lon, empty `Time` | GCM diagnostics; static fields only |
| `start.nc` | 33×33 lat/lon, 26 levels | Dynamics restart; staggered `rlonu`/`rlatv`, character variables |
| `startfi.nc` | 994 physics points | Unstructured grid; latitudes in radians with no `units` |
| `startfi_1d.nc` | 1 physics point | Single column; everything collapses to scalars and profiles |
| `startevo.nc` | 994 physics points, `nslope` | PEM restart; `units` given as `Radian North-South` |
| `xoutdaily4pem_y1.nc` | 33×32 lat/lon, 669 steps, 60 MB | XIOS daily output; largest file |
| `xoutdaily4pem_y1_1d.nc` | single column, 669 steps | XIOS single-column output |
| `xoutyearly4pem_y1.nc` | 33×32 lat/lon, 1 step | XIOS yearly averages |

Sweeping every variable of every sample plots 190 of 194; the rejections are the
three `controle_descriptor` character arrays and `diagfi.nc`'s empty `Time`.

## Behaviour notes

- Only the requested hyperslab is read, and large reductions stream through dask
  — both spellings of one, `--reduce mean:DIM` and `-e '{"DIM":"avg"}'`, so a
  reduction over a multi-gigabyte variable is accumulated a chunk at a time
  rather than read whole.
- A log colour scale has nowhere to put a zero or a negative, and a hole in a map
  reads as "the field is not here" rather than "the field is zero here". Three
  things follow. A log scale is not chosen automatically for a field that is
  mostly zero — `co2ice` is zero over 97% of a PEM step, and a log scale of it
  was 33 drawn cells out of 1089. Where one is chosen and more than 2% of the
  cells are zero or negative, `symlog` is chosen instead: log over the decades,
  linear across a window under the smallest positive value, so the zeros land at
  the bottom of the bar (this is the scale the 3D globe pre-normalizes — see the
  VTK note below). And whatever a log scale still cannot place — always the case
  under an explicit `--norm log`, which is honoured as asked and never upgraded —
  is drawn in grey, with how many said both under the figure and on the terminal.
  The grey is laid under the mesh rather than set as the colormap's "bad" colour,
  which cartopy needs transparent to wrap a map across the dateline — so missing
  data stays a hole, and only the zeros are coloured.
- `--anomaly` and `--diff` centre their colour scale on zero whatever colormap
  was named: `-c` says which colours, not what the middle colour means. Both say
  on the figure what they are, so a saved map is not mistaken for the raw field.
- `--stats` weights by area wherever the geometry allows it, including after a
  reduction: a zonal mean is still a field over latitude. Where the weights
  cannot survive a reduction the note says which one.
- An output path's directory is created if it is not there, and a figure that
  cannot be written is reported in one line with exit status 1.
- Maps are drawn on a cartopy projection, so longitudes in `[0, 360]` and
  `[-180, 180]` behave identically and the topography overlay always applies.
  A cyclic longitude axis has its seam closed.
- `--x-dim` is ignored on maps, which have a fixed orientation; use
  `--plot-kind section` for a transposed heatmap.
- Without `-o`, figures are shown; with `-o` they are written.
  `--show-polar` writes `<stem>_north<ext>` and `<stem>_south<ext>`.
- `--show-3d` with `-o` writes `<stem>_globe<ext>` beside the main figure, and
  under `--animate` that becomes a movie. Offscreen rendering needs a GL
  context; without one the globe reports it and suggests `xvfb-run`.
- The polar and 3D views reuse the colour scale of the map they follow, so
  `--norm` and `--vmin`/`--vmax` apply to them too. The globe drives a VTK
  lookup table, which knows only a linear and a log scale; under `--norm symlog`
  it colours the normalized values instead and its colour bar reads `0..1`, with
  the real range given in the bar title.
- An animation fixes its colour scale over every frame, so the colours mean the
  same thing throughout rather than sliding with each frame's own range.
- The globe is unlit by default. Shading multiplies the colormap by the angle to
  the light, so the same value reads differently depending on where it sits on
  the sphere, which is exactly what a colour scale exists to prevent. `--sun`
  asks for a lit planet with a terminator and so turns shading back on;
  `--light`/`--no-light` overrides either way.
- Where the field does not reach, the globe is grey. A globe with no data at all
  is not opened.
- Relief is exaggerated ten times; the atmosphere above it is exaggerated thirty
  times, because 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 two factors differing is the price, and the
  caption names both; `--shell-exaggeration 10` restores one shared scale.
- The camera frames the whole scene rather than the planet alone, so a stack of
  layers reaching 1.4 planet radii is in the picture instead of half out of it.
- Contours, the graticule and the terminator are placed on the *triangles* the
  globe is drawn from rather than on the smooth surface through their corners,
  so they lie in the terrain instead of standing off it. The two differ by up to
  20 km once exaggerated, which is what the lines used to have to be lifted by.
- The shell is placed at the levels' true altitudes whenever the file carries
  the hybrid coordinate (`ap`/`bp` or `aps`/`bps`) and a surface pressure. Each
  level is then a map of heights rather than a single one: it drapes the terrain
  at the bottom, where the coordinate *is* the ground, and flattens into an
  isobar at the top. The scale height is fitted to the file's own vertical axis
  and the gravity read from `controle`, so a run of the generic model on another
  planet is placed with its own numbers rather than with Mars'.
  `--shell-mode cloud` resamples into a uniform box, which has one vertical axis
  for the whole planet, so it uses the mean profile and says so.
- `--shell-mode column` integrates over mass, `Σ q·Δp/g`, which for a mixing
  ratio comes out in kg/m². Without the hybrid coefficients it falls back on the
  geometric integral `∫q·dz`, in the field's units times metres, and says which
  of the two it did — the geometric one weighs a layer by how thick it is rather
  than by how much air is in it, so most of the answer comes from the deep upper
  levels where a tracer is thin and horizontally flat.
- A vertical coordinate that is not a length — pressure, sigma, a bare level
  index — and no hybrid coefficients either has no metric altitude, so the shell
  spreads its levels evenly up to `--shell-top` and says so.
- `--plot-kind globe` over a `positive: down` vertical is refused: a soil column
  is not an atmosphere. Use `--plot-kind section` for that.
- `--vector U,V --show-3d` draws the winds as arrows tangent to the sphere.
  Staggered pairs such as `ucov`/`vcov` are still refused, as on flat maps.
- Non-numeric variables and zero-size dimensions are rejected with an explicit
  message.

## Tests

```bash
cd PEM/toolbox/display_netcdf
conda activate disp_nc
python -m pytest tests/ -q          # the whole suite, a few seconds
python -m pytest tests/ -q -rs      # ... and why each skipped case skipped
python -m pytest tests/test_coords.py -q -k rlonu    # one file, one case
```

Unit tests cover metadata normalization, axis resolution, the plot-inference
decision table, colour choice, statistics, the globe geometry and shell modes,
and the terminal prompts. They build their datasets in memory — one per writer
family, in `tests/conftest.py` — so they run in milliseconds and pass on a fresh
checkout with no sample files at all.

### Running against your own sample files

`sample_files/` may hold any set of NetCDF files. Nothing has to be renamed, and
a changed sample set never turns the suite red:

- The in-memory tests, which are nearly all of them, ignore the directory.
- `test_every_variable_of_every_sample_resolves` runs once per `.nc` file found,
  whatever it is called, and asserts what must hold of *any* file the tool can
  be pointed at: one axis per dimension, and no role claimed twice.
- A handful of cases pin down one quirk of one file — `start.nc`'s staggered
  `rlonu`, `startfi.nc`'s unstructured grid. Each of those **skips** when its
  file, or the variable it needs inside it, is not there, naming what it wanted.
  `-rs` lists them.

A test that reads a real file must go through the `tests/conftest.py` helpers,
which is what makes that skip happen rather than an error:

```python
from conftest import open_sample, sample_files

def test_soil_depth_points_down():
    # skips if startevo.nc is absent, or holds no tsoil
    with open_sample('startevo.nc', 'tsoil') as (ds, report):
        ...
```

Prefer extending the sweep, which survives the samples being replaced, to
naming a file. When a specific file really is the point, add a row to
`REAL_CASES` in `tests/test_coords.py`.

### The regression baseline

`tests/baseline_capture.py` additionally sweeps every variable of every sample
file and compares against `tests/baseline/plot_matrix.json`, which records the
exit status, plot classification, axis labels and a hash of the plotted array
for each:

```bash
python tests/baseline_capture.py --check tests/baseline/plot_matrix.json
```

The array is hashed rather than the PNG, so the baseline is immune to matplotlib
version changes and only moves when what is plotted moves.
`plot_matrix_prerefactor.json` is the same sweep taken before this refactor, kept
so the deliberate behaviour changes stay auditable.

Unlike the pytest suite, **a baseline is tied to the sample set it was taken
on**. `--check` compares the files the two have in common and reports the
difference in the sets once, at the top. After the samples change, take a fresh
baseline on the unmodified code and compare against that:

```bash
# from an unmodified checkout of the last known-good revision
python tests/baseline_capture.py --out /tmp/plot_matrix_new_samples.json
# then, from the working copy being changed
python tests/baseline_capture.py --check /tmp/plot_matrix_new_samples.json
```

Commit the new baseline over `plot_matrix.json` only once the sweep on unchanged
code is the reference — a baseline captured from modified code proves nothing.

## Related

- `../README.md`: the toolbox index and the other post-processing scripts
- `../analyse_netcdf.py`: min/max/mean and anomaly flags for the same files