mhd_io#

Lazy HDF readers for PSI MAS and POT3D magnetohydrodynamic model output.

This module provides a unified, unit-aware interface for reading three-dimensional field variables from Predictive Science Inc.’s MAS and POT3D solvers. Both HDF4 (.hdf) and HDF5 (.h5) files are supported through a common API. The sole public symbol exported by this module is PsiData().

Entry Point

PsiData() is a factory function that opens a PSI HDF file and returns a lazy reader object. The file extension determines the I/O backend — .h5 files are read via h5py and .hdf files via pyhdf. The model argument ('mas' or 'pot3d') selects the quantity-property and mesh-stagger metadata tables applied during metadata resolution:

from psi_io.mhd_io import PsiData

reader = PsiData('br001001.h5')                    # MAS HDF5 (default)
reader = PsiData('br001001.hdf', model='mas')      # MAS HDF4
reader = PsiData('br001.h5', model='pot3d')        # POT3D HDF5

Lazy Loading and Caching

No data are read from disk at construction time. Metadata — quantity name, sequence number, physical unit, and mesh stagger — are parsed from the filename stem and/or HDF file attributes. Data are transferred from disk only when read() or vslice() is called.

When a read or vslice call covers the entire dataset with no spatial restrictions, the result is cached on the reader object. Subsequent full-array calls return the cached copy without a second disk read. Partial reads — any call with at least one non-full-axis argument — are never cached. The cache state is exposed via the is_cached property.

Attributes

The following attributes are available on every object returned by PsiData():

quantitystr: Canonical lower-case quantity identifier (e.g. 'br', 'vr', 't'). Resolved from the filename stem, HDF file attributes, or the quantity keyword argument passed to PsiData().
sequenceint: Time-step sequence number extracted from the filename (e.g. br001001.h5 → 1001) or from HDF file attributes.
unitUnit: Conversion factor from one code unit to physical units. For MAS quantities these are the custom normalization units defined in psi_io._units (e.g. MAS_b ≈ 2.2 G for magnetic field, MAS_v ≈ 481 km s⁻¹ for velocity). For POT3D the default is dimensionless — see the PsiData() warning.
meshtuple of Mesh: Yee-grid stagger position for each spatial axis in physical (r, θ, φ) order. Each element is either HALF (offset by half a cell spacing) or MAIN (cell-centred). The integer encoding and per-quantity stagger codes are defined in psi_io._mesh; the canonical default for each quantity is stored in its Props descriptor (see psi_io._models).
propsProps: Complete property descriptor for the quantity, bundling its name, description, unit, and mesh code. Looked up from the appropriate quantity-properties mapping in psi_io._models.
descriptionstr: Human-readable description of the physical quantity (e.g. 'MAS Magnetic Field (Radial Component)'). Derived from desc.
scalesScales(r, t, p): Named tuple of coordinate scale readers. Each element wraps the one-dimensional coordinate array stored in the HDF file. The radial coordinate uses PSI_rsun (solar radii); the co-latitude and longitude coordinates use PSI_angle (radians). Each scale reader exposes the same read() interface as the main data reader.
shapetuple of int: Array dimensions in HDF storage order (Nφ, Nθ, Nr) — the radial axis is last due to Fortran column-major convention. All slicing APIs accept arguments in physical (r, θ, φ) order and reverse the indexing internally.
ndimint: Number of spatial dimensions; always 3 for MAS and POT3D field variables.
sizeint: Total element count (Nφ × Nθ × Nr).
nbytesint: Dataset size in bytes.
dtypenumpy.dtype: Element type of the stored dataset (typically float32).
attrsdict: HDF file-level attributes as a plain Python dictionary.
is_cachedbool: True once a full-array read has populated the in-memory cache; False otherwise.

Reading Data — read

odata[, r, t, p] = reader.read(*args, unit=None, mesh=None, scales=True)

Each positional argument restricts one spatial axis in physical (r, θ, φ) order:

Omitted / None — full axis.
int — single index; axis retained as a length-1 dimension.
slice — standard Python slice.
(start, stop) or (start, stop, step) — converted to a slice.
Ellipsis — expands to None for all remaining axes.

unit: Output unit. String aliases 'native' / 'code' / 'model' / 'psi' return values in MAS code units (the units defined in psi_io._units). Aliases 'real' / 'phys' / 'physical' / 'cgs' call decompose_mas_units() to express values in CGS base units. Any other string or Unit instance is forwarded to astropy.units.Quantity.to().
mesh: Target mesh stagger (MeshCodeType). Axes that are on the half mesh in the stored data but on the main mesh in mesh are averaged via remesh_array(). Up-sampling (main → half) raises ValueError.
scales: If True (default), return the corresponding coordinate slice for each axis as additional Quantity values (r, t, p) after the data array.

Note

PSI HDF files are written in Fortran column-major order so that numpy reads them with shape (Nφ, Nθ, Nr) — the radial axis is last. All positional arguments to read() and the returned coordinate scales are in physical (r, θ, φ) order regardless of on-disk layout.

Value-Space Slicing — vslice

odata[, r, t, p] = reader.vslice(*args, unit=None, mesh=None, scales=True,
                                 bounds_error=True, fill_value=None)

vslice is a superset of read() that additionally accepts physical coordinate values as positional arguments. Passing a Quantity or a bare scalar for an axis locates the two nearest grid points, loads a 2-element window, and linearly interpolates to the target value; the resulting axis has size 1. Index-space arguments (slice, int, None, Ellipsis) are handled identically to read().

bounds_error: If True (default), raise ValueError when a value is outside the coordinate range. Set to False and supply a fill_value to replace out-of-bounds points; pass fill_value=None to extrapolate silently.

When scales=True, value-interpolated axes return the target coordinate as a length-1 Quantity; index-space axes return the corresponding coordinate slice as usual.

File Lifecycle

Readers hold an open HDF file handle. Use the context-manager protocol to guarantee cleanup:

with PsiData('br001001.h5') as reader:
    data, r, t, p = reader.read(unit='Gauss')

Alternatively, call reader.close() and reader.open() to manage the handle explicitly. The handle is also released automatically when the reader is garbage-collected.

Supported Quantities

MAS (19 quantities):

Key(s)	Description	Code unit	Stagger `(r, θ, φ)`
`br`, `bt`, `bp`	Magnetic field components	`MAS_b` (≈ 2.2 G)	`HALF/MAIN/MAIN`, `MAIN/HALF/MAIN`, `MAIN/MAIN/HALF`
`vr`, `vt`, `vp`	Velocity components	`MAS_v` (≈ 481 km s⁻¹)	`MAIN/HALF/HALF`, `HALF/MAIN/HALF`, `HALF/HALF/MAIN`
`jr`, `jt`, `jp`	Current density components	`MAS_j`	`MAIN/HALF/HALF`, `HALF/MAIN/HALF`, `HALF/HALF/MAIN`
`t`, `te`, `tp`	Temperature (single / electron / proton)	`MAS_t` (≈ 28 MK)	`HALF/HALF/HALF`
`rho`	Mass density	`MAS_n` (10⁸ cm⁻³)	`HALF/HALF/HALF`
`p`	Thermal pressure	`MAS_p`	`HALF/HALF/HALF`
`ep`, `em`	Alfvén wave energy density (outward / inward)	`MAS_p`	`HALF/HALF/HALF`
`zp`, `zm`	Elsässer wave amplitudes (outward / inward)	`MAS_v`	`HALF/HALF/HALF`
`heat`	Volumetric coronal heating rate	`MAS_heat`	`HALF/HALF/HALF`

POT3D (3 quantities): br, bt, bp — magnetic field components, unit POT3D_b (dimensionless by default; see the PsiData() warning regarding unit declaration).

Examples

Full MAS radial field read with coordinate scales:

from psi_io.mhd_io import PsiData

reader = PsiData('br001001.h5')
data, r, t, p = reader.read()              # code units (MAS_b)
data, r, t, p = reader.read(unit='Gauss')  # convert to Gauss on the fly

Partial read — inner radial shell in CGS base units, coordinates suppressed:

data = reader.read(slice(0, 10), unit='physical', scales=False)

Value-space slice — extract the r = 2.5 R☉ surface:

data, r, t, p = reader.vslice(2.5, unit='Gauss')  # bare scalar → native coord unit

POT3D file — unit must be declared at construction because it is not encoded in the file:

reader = PsiData('br001.hdf', model='pot3d', unit='Gauss')
data, r, t, p = reader.read()

Inspect metadata without loading data:

reader = PsiData('rho001001.h5')
reader.quantity      # 'rho'
reader.description   # 'MAS Density'
reader.unit          # MAS_n
reader.mesh          # (Mesh.HALF, Mesh.HALF, Mesh.HALF)
reader.is_cached     # False
reader.shape         # (Nφ, Nθ, Nr) — numpy storage order

Functions

PsiData(ifile, /[, model])

Open a PSI MAS or POT3D HDF file and return the appropriate data reader.