Reductions¶

pypic.reduce collapses a FieldDataset along one or more axes using a chosen reduction operation. Column densities, line-of-sight integrals, slab averages, projected-max diagnostics, and peak-position maps all compose from this single verb — paired with an optional BoxSelection or SphereSelection for region restriction.

The verb is named reduce (not project) to avoid colliding with Three.js Vector3.project(camera) (camera/screen-space projection) on the webpic client side. The pypic server reduces; the viewer projects.

Reductions table¶

Reduction	xarray dispatch	Multi-axis?	Metadata behavior
`integrate`	`ds.integrate(coord=ax)` looped over axes	yes (sequential trapezoidal)	preserved
`sum`	`ds.sum(dim=axis, skipna=...)`	yes	preserved
`mean`	`ds.mean(dim=axis, skipna=...)`	yes	preserved
`median`	`ds.median(dim=axis, skipna=...)`	yes	preserved
`max` / `min`	`ds.{max,min}(dim=axis, skipna=...)`	yes	preserved
`std` / `var`	`ds.{std,var}(dim=axis, skipna=...)`	yes	preserved
`argmax` / `argmin`	`ds.{idxmax,idxmin}(dim=axis, ...)`	single axis only	overridden to `quantity_type="length"`

argmax / argmin use xarray's idxmax / idxmin to return the coordinate value of the extremum (e.g. the z-position where |J| peaks) — not a raw integer index. Multi-axis search has no single coord-value to return, so reject with ValueError.

Single-axis vs multi-axis¶

# Single axis: 3D → 2D
column = pypic.reduce(ds, "z", reduction="integrate")

# Tuple: 3D → 1D in one call (sequential trapezoidal under the hood)
line_out = pypic.reduce(ds, ("y", "z"), reduction="mean")

The single-axis and tuple forms are functionally identical for the non-integrate reductions; for integrate the multi-axis form runs a sequential trapezoidal that matches the chained two-call form to machine precision.

Unit handling after `integrate`¶

Unweighted integrate shifts the SI unit dimension by one length factor per reduced axis (number density m⁻³ → column density m⁻², energy density J/m³ → areal energy J/m²). reduce stamps the running count on attrs["reduction"]["length_axes"] so FieldDataset.in_si applies the extra length_ref**n factor at the boundary — column densities come out in the correct C/m² (or m⁻², etc.) without manual length-unit bookkeeping.

column = pypic.reduce(ds, "z", reduction="integrate")
column.in_si("rho_c")                # C/m^2 — correct SI
column.xr["rho_c"].attrs["reduction"]
# {"axis": "z", "op": "integrate", "length_axes": 1}

Chained integrates accumulate: reduce(reduce(ds, "y", "integrate"), "z", "integrate") stamps length_axes=2, equivalent to reduce(ds, ("y","z"), "integrate"). Non-integrate reductions (mean/sum/max/...) are unit-preserving and carry the count forward unchanged. Weighted integrate does not stamp length_axes because the length factor cancels in ∫ f w dx / ∫ w dx.

Post-reduction si_unit strings and openPMD 7-tuples are not yet generalized; the length_axes mechanism here is the interim fix that gets in_si() returning the right number today.

Weight semantics¶

weight=<field-name> produces yt-style density- or emission-weighted averages. Only reduction="mean" and reduction="integrate" accept a weight — other reductions raise ValueError.

Weighted mean: \(\langle f \rangle_w = \sum_i f_i w_i / \sum_i w_i\) over the reduced axes.
Weighted integrate: \(\int f w \, dx / \int w \, dx\) — the column weighted average along the line of sight (yt's weight_field convention).

# Density-weighted temperature column average
T_avg = pypic.reduce(ds, "z", reduction="integrate",
                     weight="rho_c", fields=["T_s0"])

Provenance: T_avg.xr["T_s0"].attrs["reduction"] carries {"axis": "z", "op": "integrate", "weight": "rho_c"}.

NaN handling. Under nan_policy="omit" (default) the weighted path uses joint masking: cells where the field or weight is NaN contribute zero to both numerator and denominator, skipping them consistently. Under "propagate", NaN flows through naturally. Under "raise", both field and weight are pre-checked for NaN.

Selection composition¶

selection= runs first, then the reduction operates on the selected region:

# Mean over a sub-volume
box = pypic.BoxSelection(ranges={"x": (10, 30), "y": (10, 30)})
slab = pypic.reduce(ds, "z", reduction="mean", selection=box)

# Line of sight through a sphere — NaN-mask outside, integrate inside
ball = pypic.SphereSelection(center=(0.0, 0.0, 0.0), radius=5.0,
                             keep="inside")
los = pypic.reduce(ds, "z", reduction="integrate", selection=ball)

SphereSelection pairs naturally with nan_policy="omit": outside cells become NaN and are dropped from the reduction.

Plot composition pattern¶

A reduced FieldDataset is just a smaller-dimension FieldDataset — all existing plotting helpers work unchanged:

from pypic.plotting import plot_field_slice

fig, ax = plot_field_slice(
    ds.reduce("z", reduction="integrate"),
    "rho_c",
    title="Column density",
)

reduce → plot_field_slice is the canonical pattern for column density / LOS imaging.

Time-axis reduction¶

reduce accepts any dimension on the underlying xarray dataset — not just the spatial grid.surviving_axis_names. For time-series stores written by to_zarr_timeseries, every field has time as the leading dim, so a time-mean / time-integrate works directly:

ts = pypic.from_zarr("run.zarr")     # fields shaped (nt, nx, ny, nz)
B_mean = pypic.reduce(ts, "time", reduction="mean")

Pure non-spatial reductions bypass the Cartesian-grid gate — the Jacobian only matters when integrating over a spatial axis. A mixed ("time", "r") reduction on a spherical grid still raises GeometryUnsupportedError, pending Jacobian-aware integration.

Worked examples¶

Three first-line plasma diagnostics expressible in the current API.

EDR localization via the Zenitani localizer¶

The electron-frame dissipation \(D_e\) (Zenitani EDR localizer) peaks inside the electron diffusion region. Per-column z-position of the peak gives an EDR centroid map:

ds_with_De = ds.with_derived("D_e")
# argmax returns the *coordinate value* of the maximum along z,
# not an integer index — the result is one z-position per (x, y).
z_peak = pypic.reduce(ds_with_De, "z", reduction="argmax", fields=["D_e"])
# z_peak.field_info("D_e").quantity_type == "length"

For the value at the peak (rather than the location), swap to reduction="max".

Density-weighted column temperature (synthetic LOS observable)¶

The line-of-sight temperature that an external observer would measure is the emission/density-weighted average — the yt "weight_field" convention:

T_col = pypic.reduce(ds, "z", reduction="integrate",
                     weight="rho_c", fields=["T_s0"])
# T_col["T_s0"] == ∫ T_s0 ρ_c dz / ∫ ρ_c dz

The provenance attr records the weight: T_col.xr["T_s0"].attrs["reduction"] → {"axis": "z", "op": "integrate", "weight": "rho_c"}.

Energy-budget conservation¶

For each timestep, integrate the energy densities over the full volume to get total magnetic, kinetic, and thermal energy:

totals_by_step: list[tuple[float, dict[str, float]]] = []
for step in sim.steps:
    # compute() returns one array; with_derived() attaches several.
    ds = sim.read(step).with_derived("e_B", "e_k", "e_th")
    totals = pypic.reduce(ds, ("x", "y", "z"), reduction="integrate")
    totals_by_step.append(
        (step * sim.grid.dt, {
            "E_B": float(totals["e_B"]),
            "E_k": float(totals["e_k"]),
            "E_th": float(totals["e_th"]),
        })
    )
# Plot E_B(t) + E_k(t) + E_th(t) — should be conserved in ideal MHD
# runs and slowly evolve in dissipative ones.

CLI¶

pypic reduce apply <sim-dir> --output out.zarr --axis z \
    --reduction mean --weight rho_c --fields T_s0

Supports multi-step batch writes (--step all → time-series Zarr), selection composition (--box, --plane), Icechunk versioning (--backend icechunk --tag v1), and --dry-run for plan preview.

API reference¶

`reductions` ¶

Axis reductions: collapse a FieldDataset along one or more dimensions.

The reduce function collapses a FieldDataset along one axis (axis="z") or several (axis=("y", "z")) using a chosen reduction. Supported reductions: "integrate", "sum", "mean", "median", "max", "min", "std", "var", "argmax", "argmin". Column densities, line-of-sight integrals, slab averages, projected-max diagnostics, and peak-position maps all compose from this single primitive paired with an optional BoxSelection or SphereSelection.

There is deliberately no SlabSelection: selections describe regions, not data, and a slab would fuse "pick a thick slice" with "reduce along it". Kept separate, BoxSelection stays reusable outside reductions and a viewer request is just {selection, axis, reduction}.

After an unweighted integrate the quantity_type and si_unit strings are preserved unchanged, and FieldDataset.in_si corrects the value via length_ref ** length_axes from attrs. The unit strings themselves do not shift by the reduced length factors.

`Reduction = Literal['integrate', 'sum', 'mean', 'median', 'max', 'min', 'std', 'var', 'argmax', 'argmin']` ¶

`reduce(data, axis, *, reduction='integrate', selection=None, fields=None, weight=None, nan_policy='omit')` ¶

Reduce a FieldDataset along one or more axes.

Apply selection (if given) to crop or NaN-mask, then collapse the axis dimension(s) using reduction. Returns a FieldDataset with one fewer surviving axis per name passed — the natural input to plot_field_slice, compute(...), in_si(...), and other dataset-consuming APIs.

Parameters:

Name	Type	Description	Default
`data`	`FieldDataset`	Input dataset (1-, 2-, or 3-D — must contain every name in axis).	required
`axis`	`str or tuple of str`	Dimension name(s) to reduce away (e.g. `"x"`, `("y", "z")`, `"time"`). Validated against the underlying xarray dataset's dims so non-grid dims like `time` (added by `pypic.io.to_zarr_timeseries`) are accepted alongside the spatial `grid.surviving_axis_names`. Multi-axis reduction collapses several dimensions in one call — useful for going from 3D to 1D without chaining. `"argmax"` and `"argmin"` require a single axis.	required
`reduction`	`str`	How to collapse the axis. Choices: `"integrate"` — xarray trapezoidal-rule integration over the coordinate values (looped over axes for multi-axis input). Default — matches plasma-physics convention (column densities, integrated \(\mathbf{J}\!\cdot\!\mathbf{E}\)). `"sum"` — unweighted Riemann sum (no `× dx` factor). `"mean"` / `"median"` / `"max"` / `"min"` / `"std"` / `"var"` — direct xarray equivalents. `"argmax"` / `"argmin"` — return the coordinate value of the extremum along axis (xarray's `idxmax` / `idxmin`); single-axis only. Result is a position along the reduced axis, so `quantity_type` is overridden to `"length"`.	`'integrate'`
`selection`	`BoxSelection \| SphereSelection \| None`	Optional region selector applied before reduction. `None` reduces over the full domain. Sugar for `selection.apply(data)`-then-reduce; the two forms are observationally identical.	`None`
`fields`	`Iterable[str] \| None`	Optional subset of field names (canonical or alias) to reduce. `None` reduces every data variable. Unknown names raise `KeyError` with the full list rather than being skipped.	`None`
`weight`	`str \| None`	Name of a field to weight the reduction by. Only supported for `reduction in {"mean", "integrate"}` — other reductions raise `ValueError` (`sum` is trivially achieved by pre-multiplying; `max` / `min` / `median` / `std` / `var` / `argmax` / `argmin` reject weights). `mean` returns \(\langle f \rangle_w = \sum f w / \sum w\) over the reduced axes; `integrate` returns the weighted line average \(\int f w \, dx / \int w \, dx\) (yt's emission-weighted / density-weighted column average). Resolves via `data.resolve_key` (aliases accepted); unknown names raise `KeyError`. NaN handling under `nan_policy="omit"` uses a joint mask: cells where the field or weight is NaN contribute zero to both the numerator and denominator and are skipped consistently.	`None`
`nan_policy`	`('omit', 'propagate', 'raise')`	`"omit"` (default) → `skipna=True`: NaN cells are dropped from each reduction. Pairs naturally with `SphereSelection`, which NaN-masks outside-region cells. `"propagate"` lets NaN poison the result. `"raise"` errors when any input cell is NaN. Same vocabulary and semantics as `pypic.diagnostics.l2_relative_error` — see `conventions.md § Error Norms and Divergence` for the broader rationale.	`"omit"`

Returns:

Type Description

FieldDataset

Same metadata (normalization, species, frame, transforms) with the named axis (or axes) removed from the grid and dimensions. Each surviving DataArray gains an attrs["reduction"] dict recording {axis, op} — plus result_kind: "axis_position" for argmax / argmin, weight: <canonical name> when weight is set, and length_axes: <int> after unweighted integrate (the running count of length-dimension shifts across chained reductions). The inner op key holds the reduction name ("mean", "integrate", ...) to avoid shadowing the outer reduction key.

Raises:

Type	Description
`GeometryUnsupportedError`	Non-Cartesian geometry combined with a spatial-axis reduction. Spherical / cylindrical Jacobian-aware integration is not yet implemented. Pure non-spatial reductions (e.g. along `time`) bypass this check.
`ValueError`	Unknown axis name, invalid reduction or nan_policy, multi-axis input passed with `reduction="argmax"` / `"argmin"`, or — under `nan_policy="raise"` — any NaN in a reduced field.
`KeyError`	Any name in fields fails to resolve via `FieldDataset.resolve_key`.

Notes

Unit shift after integrate. After unweighted reduction="integrate" the SI unit dimension shifts by one length factor along each reduced axis (e.g. number density m\ :sup:-3 → column density m\ :sup:-2). Per-field quantity_type and si_unit strings are preserved unchanged — but the numeric value returned by in_si() is correct: the length_axes provenance stamp here lets FieldDataset.in_si apply an extra length_ref ** length_axes factor at the boundary. Weighted integrate does not stamp length_axes because the length factor cancels between numerator and denominator. The displayed unit string and the openPMD 7-tuple do not yet shift with the reduction.

Examples:

>>> import numpy as np
>>> from pypic.coordinates import CARTESIAN
>>> from pypic.dataset import FieldDataset
>>> from pypic.grid import GridInfo
>>> from pypic.reductions import reduce
>>> from pypic.units import Normalization
>>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0),
...     origin=(0.0, 0.0, 0.0), geometry=CARTESIAN)
>>> ds = FieldDataset.from_arrays(
...     {"B_1": np.ones((4, 3, 2))}, grid, Normalization.identity())
>>> column = reduce(ds, "z", reduction="integrate")
>>> column.grid.dimensions
(4, 3)
>>> column.grid.surviving_axis_names
('x', 'y')
>>> line = reduce(ds, ("y", "z"), reduction="mean")
>>> line.grid.dimensions
(4,)

Source code in src/pypic/reductions.py

def reduce(
    data: FieldDataset,
    axis: str | tuple[str, ...],
    *,
    reduction: Reduction = "integrate",
    selection: BoxSelection | SphereSelection | None = None,
    fields: Iterable[str] | None = None,
    weight: str | None = None,
    nan_policy: NanPolicy = "omit",
) -> FieldDataset:
    r"""Reduce a FieldDataset along one or more axes.

    Apply *selection* (if given) to crop or NaN-mask, then collapse the
    *axis* dimension(s) using *reduction*.  Returns a FieldDataset with
    one fewer surviving axis per name passed — the natural input to
    ``plot_field_slice``, ``compute(...)``, ``in_si(...)``, and other
    dataset-consuming APIs.

    Parameters
    ----------
    data : FieldDataset
        Input dataset (1-, 2-, or 3-D — must contain every name in *axis*).
    axis : str or tuple of str
        Dimension name(s) to reduce away (e.g. ``"x"``, ``("y", "z")``,
        ``"time"``).  Validated against the underlying xarray dataset's
        dims so non-grid dims like ``time`` (added by
        [`pypic.io.to_zarr_timeseries`][pypic.io.to_zarr_timeseries]) are accepted
        alongside the
        spatial ``grid.surviving_axis_names``.  Multi-axis reduction
        collapses several dimensions in one call — useful for going
        from 3D to 1D without chaining.  ``"argmax"`` and ``"argmin"``
        require a single axis.
    reduction : str
        How to collapse the axis.  Choices:

        * ``"integrate"`` — xarray trapezoidal-rule integration over the
          coordinate values (looped over axes for multi-axis input).
          Default — matches plasma-physics convention (column densities,
          integrated $\mathbf{J}\!\cdot\!\mathbf{E}$).
        * ``"sum"`` — unweighted Riemann sum (no ``× dx`` factor).
        * ``"mean"`` / ``"median"`` / ``"max"`` / ``"min"`` / ``"std"`` /
          ``"var"`` — direct xarray equivalents.
        * ``"argmax"`` / ``"argmin"`` — return the *coordinate value* of
          the extremum along *axis* (xarray's ``idxmax`` / ``idxmin``);
          single-axis only.  Result is a position along the reduced
          axis, so ``quantity_type`` is overridden to ``"length"``.
    selection : BoxSelection | SphereSelection | None
        Optional region selector applied before reduction.  ``None``
        reduces over the full domain.  Sugar for
        ``selection.apply(data)``-then-reduce; the two forms are
        observationally identical.
    fields : Iterable[str] | None
        Optional subset of field names (canonical or alias) to reduce.
        ``None`` reduces every data variable.  Unknown names raise
        `KeyError` with the full list rather than being skipped.
    weight : str | None
        Name of a field to weight the reduction by.  Only supported for
        ``reduction in {"mean", "integrate"}`` — other reductions raise
        ``ValueError`` (``sum`` is trivially achieved by pre-multiplying;
        ``max`` / ``min`` / ``median`` / ``std`` / ``var`` /
        ``argmax`` / ``argmin`` reject weights).  ``mean`` returns
        $\langle f \rangle_w = \sum f w / \sum w$ over the reduced axes;
        ``integrate`` returns the weighted line average
        $\int f w \, dx / \int w \, dx$ (yt's emission-weighted /
        density-weighted column average).  Resolves via
        ``data.resolve_key`` (aliases accepted); unknown names raise
        `KeyError`.  NaN handling under ``nan_policy="omit"``
        uses a joint mask: cells where the field *or* weight is NaN
        contribute zero to both the numerator and denominator and are
        skipped consistently.
    nan_policy : {"omit", "propagate", "raise"}
        ``"omit"`` (default) → ``skipna=True``: NaN cells are dropped
        from each reduction.  Pairs naturally with
        [`SphereSelection`][pypic.selections.SphereSelection], which NaN-masks
        outside-region cells.  ``"propagate"`` lets NaN poison the
        result.  ``"raise"`` errors when any input cell is NaN.
        Same vocabulary and semantics as
        [`pypic.diagnostics.l2_relative_error`][pypic.diagnostics.l2_relative_error] —
        see
        `conventions.md § Error Norms and Divergence` for the broader
        rationale.

    Returns
    -------
    FieldDataset
        Same metadata (normalization, species, frame, transforms) with
        the named axis (or axes) removed from the grid and dimensions.
        Each surviving DataArray gains an ``attrs["reduction"]`` dict
        recording ``{axis, op}`` — plus ``result_kind: "axis_position"``
        for ``argmax`` / ``argmin``, ``weight: <canonical name>`` when
        *weight* is set, and ``length_axes: <int>`` after unweighted
        ``integrate`` (the running count of length-dimension shifts
        across chained reductions).  The inner ``op`` key holds the
        reduction name (``"mean"``, ``"integrate"``, ...) to avoid
        shadowing the outer ``reduction`` key.

    Raises
    ------
    GeometryUnsupportedError
        Non-Cartesian geometry combined with a spatial-axis reduction.
        Spherical / cylindrical Jacobian-aware integration is not yet
        implemented.  Pure non-spatial reductions (e.g. along ``time``)
        bypass this check.
    ValueError
        Unknown *axis* name, invalid *reduction* or *nan_policy*,
        multi-axis input passed with ``reduction="argmax"`` /
        ``"argmin"``, or — under ``nan_policy="raise"`` — any NaN in a
        reduced field.
    KeyError
        Any name in *fields* fails to resolve via
        ``FieldDataset.resolve_key``.

    Notes
    -----
    **Unit shift after ``integrate``.** After unweighted
    ``reduction="integrate"`` the SI unit dimension shifts by one
    length factor along each reduced axis (e.g. number density
    m\ :sup:`-3` → column density m\ :sup:`-2`).  Per-field
    ``quantity_type`` and ``si_unit`` *strings* are preserved
    unchanged — but the numeric value returned by ``in_si()`` is
    correct: the ``length_axes`` provenance stamp here lets
    ``FieldDataset.in_si`` apply an extra ``length_ref **
    length_axes`` factor at the boundary.  Weighted ``integrate``
    does not stamp ``length_axes`` because the length factor
    cancels between numerator and denominator.  The displayed
    unit string and the openPMD 7-tuple do not yet shift with the
    reduction.

    Examples
    --------
    >>> import numpy as np
    >>> from pypic.coordinates import CARTESIAN
    >>> from pypic.dataset import FieldDataset
    >>> from pypic.grid import GridInfo
    >>> from pypic.reductions import reduce
    >>> from pypic.units import Normalization
    >>> grid = GridInfo(dimensions=(4, 3, 2), spacing=(1.0, 1.0, 1.0),
    ...     origin=(0.0, 0.0, 0.0), geometry=CARTESIAN)
    >>> ds = FieldDataset.from_arrays(
    ...     {"B_1": np.ones((4, 3, 2))}, grid, Normalization.identity())
    >>> column = reduce(ds, "z", reduction="integrate")
    >>> column.grid.dimensions
    (4, 3)
    >>> column.grid.surviving_axis_names
    ('x', 'y')
    >>> line = reduce(ds, ("y", "z"), reduction="mean")
    >>> line.grid.dimensions
    (4,)
    """
    axes: tuple[str, ...] = (axis,) if isinstance(axis, str) else tuple(axis)
    _validate_reduce(axes, reduction=reduction, weight=weight, nan_policy=nan_policy)
    if selection is not None:
        data = selection.apply(data)
    _validate_axes(data, axes, reduction)

    ds_to_reduce = data.xr if fields is None else data.xr[_resolve_fields(data, fields)]
    weight_canonical = None if weight is None else _resolve_weight(data, weight)
    weight_da = None if weight_canonical is None else data.xr[weight_canonical]
    if nan_policy == "raise":
        _reject_nan(ds_to_reduce, weight_da, weight_canonical)

    reduced = _collapse(
        ds_to_reduce,
        axes,
        reduction,
        weight_da=weight_da,
        nan_policy=nan_policy,
        source=data.xr,
    )
    _stamp_provenance(reduced, data.xr, axes, reduction, weight_canonical)
    # reductions sits above dataset in the layering, so it reuses the
    # dataset's own re-wrap rather than reconstructing grid and aliases.
    return data._wrap_sliced(reduced)  # noqa: SLF001