clawhdf5/crates/clawhdf5-netcdf4/tests/interop_tests.rs

//! NetCDF-4 interop tests: Python creates NetCDF-4 files, clawhdf5-netcdf4 reads them.
//!
//! Tests are skipped if python3 or netCDF4/xarray Python packages are not available.

use std::process::Command;

use clawhdf5_netcdf4::{AttrValue, NcType, NetCDF4File};

// ---------------------------------------------------------------------------
// Helpers
// ---------------------------------------------------------------------------
/// The Python interpreter to drive interop checks with.
///
/// `CLAWHDF5_PYTHON` lets these run against a virtualenv holding h5py, which
/// on a PEP 668 "externally managed" system is the only place it can be
/// installed. Without it the suite silently skips, and a silent skip here is
/// how a datatype bug once reached a release.
fn python() -> String {
    std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".to_string())
}

/// When `CLAWHDF5_REQUIRE_INTEROP=1` (set in CI), a missing Python dependency
/// is a test failure instead of a silent skip.
fn interop_required() -> bool {
    std::env::var("CLAWHDF5_REQUIRE_INTEROP").is_ok_and(|v| v == "1")
}

fn netcdf4_python_available() -> bool {
    Command::new(python())
        .args(["-c", "import netCDF4; print(netCDF4.__version__)"])
        .output()
        .map(|o| o.status.success())
        .unwrap_or(false)
}

fn xarray_available() -> bool {
    Command::new(python())
        .args(["-c", "import xarray; print(xarray.__version__)"])
        .output()
        .map(|o| o.status.success())
        .unwrap_or(false)
}

macro_rules! skip_if_no_netcdf4 {
    () => {
        if !netcdf4_python_available() {
            assert!(
                !interop_required(),
                "CLAWHDF5_REQUIRE_INTEROP=1 but python3 with netCDF4 is not available"
            );
            eprintln!("SKIP: python3 with netCDF4 not available");
            return;
        }
    };
}

macro_rules! skip_if_no_xarray {
    () => {
        if !xarray_available() {
            assert!(
                !interop_required(),
                "CLAWHDF5_REQUIRE_INTEROP=1 but python3 with xarray is not available"
            );
            eprintln!("SKIP: python3 with xarray not available");
            return;
        }
    };
}

fn run_python(script: &str) {
    let output = Command::new(python())
        .args(["-c", script])
        .output()
        .expect("failed to run python3");
    if !output.status.success() {
        let stderr = String::from_utf8_lossy(&output.stderr);
        let stdout = String::from_utf8_lossy(&output.stdout);
        panic!("Python script failed:\nSTDOUT: {stdout}\nSTDERR: {stderr}");
    }
}

// ===========================================================================
// 1. Python netCDF4 creates file with dims, vars, CF attrs -> read with clawhdf5
// ===========================================================================

#[test]
fn netcdf4_python_creates_cf_file_clawhdf5_reads() {
    skip_if_no_netcdf4!();
    let dir = tempfile::tempdir().unwrap();
    let path = dir.path().join("cf_test.nc");
    let path_str = path.display().to_string();

    let script = format!(
        r#"
import netCDF4 as nc
import numpy as np

ds = nc.Dataset("{path_str}", "w", format="NETCDF4")
ds.Conventions = "CF-1.8"
ds.title = "Test CF Dataset"

# Dimensions
lat_dim = ds.createDimension("lat", 3)
lon_dim = ds.createDimension("lon", 4)
time_dim = ds.createDimension("time", None)  # unlimited

# Coordinate variables
lat = ds.createVariable("lat", "f4", ("lat",))
lat.units = "degrees_north"
lat.standard_name = "latitude"
lat[:] = [10.0, 20.0, 30.0]

lon = ds.createVariable("lon", "f4", ("lon",))
lon.units = "degrees_east"
lon.standard_name = "longitude"
lon[:] = [-120.0, -110.0, -100.0, -90.0]

time = ds.createVariable("time", "f8", ("time",))
time.units = "hours since 2000-01-01"
time.calendar = "standard"
time[:] = [0.0, 6.0, 12.0]

# Data variable
temp = ds.createVariable("temperature", "f4", ("time", "lat", "lon"),
                         fill_value=-9999.0)
temp.units = "K"
temp.long_name = "Air Temperature"
temp.standard_name = "air_temperature"
data = np.arange(36, dtype=np.float32).reshape(3, 3, 4) + 270.0
temp[:] = data

ds.close()
"#
    );
    run_python(&script);

    let file = NetCDF4File::open(&path).unwrap();

    // Check dimensions
    let dims = file.dimensions().unwrap();
    let dim_names: Vec<&str> = dims.iter().map(|d| d.name.as_str()).collect();
    assert!(dim_names.contains(&"lat"));
    assert!(dim_names.contains(&"lon"));
    assert!(dim_names.contains(&"time"));

    let lat_dim = dims.iter().find(|d| d.name == "lat").unwrap();
    assert_eq!(lat_dim.size, 3);

    let lon_dim = dims.iter().find(|d| d.name == "lon").unwrap();
    assert_eq!(lon_dim.size, 4);

    let time_dim = dims.iter().find(|d| d.name == "time").unwrap();
    assert_eq!(time_dim.size, 3);

    // Check variables
    let variables = file.variables().unwrap();
    let var_names: Vec<String> = variables.iter().map(|v| v.name().to_string()).collect();
    assert!(var_names.contains(&"lat".to_string()));
    assert!(var_names.contains(&"lon".to_string()));
    assert!(var_names.contains(&"time".to_string()));
    assert!(var_names.contains(&"temperature".to_string()));

    // Read lat values
    let lat_var = file.variable("lat").unwrap();
    let lat_vals = lat_var.read_raw_f32().unwrap();
    assert_eq!(lat_vals, vec![10.0f32, 20.0, 30.0]);

    // Check CF attributes on temperature
    let mut temp_var = file.variable("temperature").unwrap();
    let cf = temp_var.cf_attributes().unwrap();
    assert_eq!(cf.units.as_deref(), Some("K"));
    assert_eq!(cf.long_name.as_deref(), Some("Air Temperature"));
    assert_eq!(cf.standard_name.as_deref(), Some("air_temperature"));

    // Read temperature data
    let temp_vals = temp_var.read_raw_f32().unwrap();
    assert_eq!(temp_vals.len(), 36); // 3 * 3 * 4
    assert!((temp_vals[0] - 270.0).abs() < 0.01);
    assert!((temp_vals[35] - 305.0).abs() < 0.01);

    // Check global attributes
    let global_attrs = file.global_attrs().unwrap();
    assert!(matches!(global_attrs.get("Conventions"), Some(AttrValue::String(s)) if s == "CF-1.8"));
    assert!(
        matches!(global_attrs.get("title"), Some(AttrValue::String(s)) if s == "Test CF Dataset")
    );
}

// ===========================================================================
// 2. Python xarray creates file -> clawhdf5 reads dimensions and variables
// ===========================================================================

#[test]
fn xarray_creates_file_clawhdf5_reads() {
    skip_if_no_xarray!();
    let dir = tempfile::tempdir().unwrap();
    let path = dir.path().join("xarray_test.nc");
    let path_str = path.display().to_string();

    let script = format!(
        r#"
import xarray as xr
import numpy as np
import pandas as pd

# Create xarray Dataset
times = pd.date_range("2020-01-01", periods=5, freq="D")
lats = [10.0, 20.0, 30.0]
lons = [-120.0, -110.0]

temp = np.random.RandomState(42).randn(5, 3, 2).astype(np.float64) * 10 + 280
precip = np.random.RandomState(123).rand(5, 3, 2).astype(np.float64) * 50

ds = xr.Dataset(
    {{
        "temperature": (["time", "lat", "lon"], temp, {{"units": "K", "long_name": "Temperature"}}),
        "precipitation": (["time", "lat", "lon"], precip, {{"units": "mm/day", "long_name": "Precipitation"}}),
    }},
    coords={{
        "time": times,
        "lat": lats,
        "lon": lons,
    }},
    attrs={{"Conventions": "CF-1.8", "source": "xarray test"}},
)
ds.to_netcdf("{path_str}", engine="netcdf4")
"#
    );
    run_python(&script);

    let file = NetCDF4File::open(&path).unwrap();

    // Check dimensions
    let dims = file.dimensions().unwrap();
    let dim_names: Vec<&str> = dims.iter().map(|d| d.name.as_str()).collect();
    assert!(dim_names.contains(&"lat"));
    assert!(dim_names.contains(&"lon"));
    assert!(dim_names.contains(&"time"));

    let lat_dim = dims.iter().find(|d| d.name == "lat").unwrap();
    assert_eq!(lat_dim.size, 3);

    let lon_dim = dims.iter().find(|d| d.name == "lon").unwrap();
    assert_eq!(lon_dim.size, 2);

    // Check variables exist
    let variables = file.variables().unwrap();
    let var_names: Vec<String> = variables.iter().map(|v| v.name().to_string()).collect();
    assert!(var_names.contains(&"temperature".to_string()));
    assert!(var_names.contains(&"precipitation".to_string()));

    // Read temperature variable
    let mut temp_var = file.variable("temperature").unwrap();
    let shape = temp_var.shape().unwrap();
    assert_eq!(shape, vec![5, 3, 2]); // time=5, lat=3, lon=2
    let temp_vals = temp_var.read_raw_f64().unwrap();
    assert_eq!(temp_vals.len(), 30); // 5*3*2

    // Read precipitation variable
    let precip_var = file.variable("precipitation").unwrap();
    let precip_vals = precip_var.read_raw_f64().unwrap();
    assert_eq!(precip_vals.len(), 30);

    // Check CF attributes
    let cf = temp_var.cf_attributes().unwrap();
    assert_eq!(cf.units.as_deref(), Some("K"));
    assert_eq!(cf.long_name.as_deref(), Some("Temperature"));

    // Check global attributes
    let global_attrs = file.global_attrs().unwrap();
    assert!(matches!(global_attrs.get("Conventions"), Some(AttrValue::String(s)) if s == "CF-1.8"));
}

// ===========================================================================
// 3. Python netCDF4 creates file with groups -> clawhdf5 reads groups
// ===========================================================================

#[test]
fn netcdf4_python_creates_grouped_file_clawhdf5_reads() {
    skip_if_no_netcdf4!();
    let dir = tempfile::tempdir().unwrap();
    let path = dir.path().join("grouped_test.nc");
    let path_str = path.display().to_string();

    let script = format!(
        r#"
import netCDF4 as nc
import numpy as np

ds = nc.Dataset("{path_str}", "w", format="NETCDF4")
ds.title = "Grouped NetCDF4 file"

# Root-level dimension and variable
ds.createDimension("x", 5)
x_var = ds.createVariable("x", "f8", ("x",))
x_var[:] = [1.0, 2.0, 3.0, 4.0, 5.0]

# Group: surface
surface = ds.createGroup("surface")
surface.description = "Surface observations"
surface.createDimension("station", 3)
temp = surface.createVariable("temperature", "f4", ("station",))
temp.units = "K"
temp[:] = [288.0, 290.0, 285.0]

# Group: upper_air
upper = ds.createGroup("upper_air")
upper.description = "Upper air soundings"
upper.createDimension("level", 4)
press = upper.createVariable("pressure", "f4", ("level",))
press.units = "hPa"
press[:] = [1000.0, 850.0, 500.0, 200.0]

ds.close()
"#
    );
    run_python(&script);

    let file = NetCDF4File::open(&path).unwrap();

    // Check root variable
    let x_var = file.variable("x").unwrap();
    let x_vals = x_var.read_raw_f64().unwrap();
    assert_eq!(x_vals, vec![1.0, 2.0, 3.0, 4.0, 5.0]);

    // Check group names
    let group_names = file.group_names().unwrap();
    assert!(group_names.contains(&"surface".to_string()));
    assert!(group_names.contains(&"upper_air".to_string()));

    // Check surface group
    let surface = file.group("surface").unwrap();
    let surface_attrs = surface.attrs().unwrap();
    assert!(matches!(
        surface_attrs.get("description"),
        Some(AttrValue::String(s)) if s == "Surface observations"
    ));

    let surf_vars = surface.variables().unwrap();
    let surf_var_names: Vec<String> = surf_vars.iter().map(|v| v.name().to_string()).collect();
    assert!(surf_var_names.contains(&"temperature".to_string()));

    let temp_var = surface.variable("temperature").unwrap();
    let temp_vals = temp_var.read_raw_f32().unwrap();
    assert_eq!(temp_vals, vec![288.0f32, 290.0, 285.0]);

    // Check upper_air group
    let upper = file.group("upper_air").unwrap();
    let press_var = upper.variable("pressure").unwrap();
    let press_vals = press_var.read_raw_f32().unwrap();
    assert_eq!(press_vals, vec![1000.0f32, 850.0, 500.0, 200.0]);
}

#[test]
fn netcdf4_python_string_variable_clawhdf5_reads() {
    // NC_STRING variables are HDF5 variable-length strings, which
    // `read_string` refused ("expected String, got VariableLength") until
    // 2026-09-26.
    skip_if_no_netcdf4!();
    let dir = tempfile::tempdir().unwrap();
    let path = dir.path().join("strings.nc");
    let path_str = path.display().to_string();
    let script = format!(
        r#"
import netCDF4 as nc
import numpy as np
ds = nc.Dataset("{path_str}", "w", format="NETCDF4")
ds.createDimension("station", 4)
v = ds.createVariable("name", str, ("station",))
v[:] = np.array(["Oslo", "", "São Paulo", "x"], dtype=object)
ds.close()
"#
    );
    run_python(&script);

    let file = NetCDF4File::open(&path).unwrap();
    let names = file.variable("name").unwrap().read_string().unwrap();
    assert_eq!(names, vec!["Oslo", "", "São Paulo", "x"]);
}

// ===========================================================================
// Unlimited dimensions: the length netCDF-C reports
// ===========================================================================

/// An unlimited dimension's length is the largest extent of the variables
/// using it, in any group (netCDF-C's `nc4_find_dim_len`), not its dimension
/// scale's extent (which netCDF-C leaves at 0): variables of different
/// lengths, one in a subgroup, a dimension no variable has written, a
/// coordinate variable, a subgroup's own unlimited dimension. Compared with
/// what netCDF4-python reports for the same file.
#[test]
fn unlimited_dimension_lengths_match_netcdf4_python() {
    skip_if_no_netcdf4!();
    let dir = tempfile::tempdir().unwrap();
    let path = dir.path().join("unlimited.nc");
    let script = format!(
        r#"
import netCDF4 as nc
import numpy as np
with nc.Dataset({path:?}, "w", format="NETCDF4") as f:
    f.createDimension("time", None)
    f.createDimension("empty", None)
    f.createDimension("rec", None)
    f.createDimension("x", 3)
    f.createVariable("t", "f8", ("time", "x"))[0:2, :] = np.ones((2, 3))
    f.createVariable("a", "i4", ("time",))[0:4] = np.arange(4)
    f.createVariable("e", "i4", ("empty",))
    f.createVariable("rec", "f4", ("rec",))[0:3] = [1, 2, 3]
    f.createVariable("r", "f4", ("x", "rec"))[:, 0:5] = np.ones((3, 5))
    g = f.createGroup("sub")
    g.createVariable("c", "i4", ("time",))[0:6] = np.arange(6)
    g.createDimension("srec", None)
    g.createVariable("s", "i4", ("srec", "x"))[0:1, :] = np.ones((1, 3))
with nc.Dataset({path:?}) as f:
    for grp in (f, f.groups["sub"]):
        for name, d in grp.dimensions.items():
            print(grp.path, name, len(d), d.isunlimited())
"#,
        path = path.display().to_string()
    );
    let out = Command::new(python())
        .args(["-c", &script])
        .output()
        .expect("failed to run python3");
    assert!(
        out.status.success(),
        "{}",
        String::from_utf8_lossy(&out.stderr)
    );
    let expected: Vec<String> = String::from_utf8(out.stdout)
        .unwrap()
        .lines()
        .map(str::to_string)
        .collect();
    // time: t has 2 records, a 4 and sub/c 6; rec: the coordinate variable
    // has 3, r 5; empty: nothing written.
    for want in [
        "/ time 6 True",
        "/ empty 0 True",
        "/ rec 5 True",
        "/ x 3 False",
        "/sub srec 1 True",
    ] {
        assert!(
            expected.iter().any(|l| l == want),
            "netCDF4 reports {expected:?}"
        );
    }

    let file = NetCDF4File::open(&path).unwrap();
    let sub = file.group("sub").unwrap();
    let mut got = Vec::new();
    for (path, dims) in [
        ("/", file.dimensions().unwrap()),
        ("/sub", sub.dimensions().unwrap()),
    ] {
        for d in dims {
            let unlimited = if d.is_unlimited { "True" } else { "False" };
            got.push(format!("{path} {} {} {unlimited}", d.name, d.size));
        }
    }
    assert_eq!(got, expected);
}

// ===========================================================================
// Variables' dimensions, shapes and values as netCDF4-python reports them
// ===========================================================================

/// Whether python can import `module`.
fn python_has(module: &str) -> bool {
    Command::new(python())
        .args(["-c", &format!("import {module}")])
        .output()
        .map(|o| o.status.success())
        .unwrap_or(false)
}

/// h5netcdf is not in every interop environment (CI installs it; a local
/// `.venv` may not have it), so its tests skip without it even under
/// `CLAWHDF5_REQUIRE_INTEROP=1`.
macro_rules! skip_if_no_h5netcdf {
    () => {
        if !python_has("h5netcdf") {
            eprintln!("SKIP: python3 with h5netcdf not available");
            return;
        }
    };
}

/// Every variable of the file at `path`, in every group, as netCDF4-python
/// reports it: `"<group path> <name> (<dims>) (<shape>)"` and its values
/// (numeric variables; element by element with masking off, so unwritten
/// records are the fill value), sorted by the description.
///
/// Values are read one element at a time because netCDF-C 4.9.3 lays out a
/// whole-variable read of a variable shorter than an unlimited dimension
/// that is not its first wrongly (the written values first, then the fill);
/// element reads, and reads of one index of the leading axis, are right.
fn netcdf4_view(path: &std::path::Path) -> Vec<(String, Vec<f64>)> {
    let script = r#"
import sys
import numpy as np
import netCDF4 as nc
def walk(g):
    for name, v in g.variables.items():
        v.set_auto_mask(False)
        head = "%s %s (%s) (%s)" % (g.path, name, ",".join(v.dimensions), ",".join(map(str, v.shape)))
        vals = []
        if v.dtype != str and v.dtype.kind in "iuf":
            vals = [repr(float(v[i])) for i in np.ndindex(v.shape)]
        print(head + "|" + " ".join(vals))
    for sub in g.groups.values():
        walk(sub)
with nc.Dataset(sys.argv[1]) as f:
    walk(f)
"#;
    let out = Command::new(python())
        .args(["-c", script, &path.display().to_string()])
        .output()
        .expect("failed to run python3");
    assert!(
        out.status.success(),
        "{}",
        String::from_utf8_lossy(&out.stderr)
    );
    let mut view: Vec<(String, Vec<f64>)> = String::from_utf8(out.stdout)
        .unwrap()
        .lines()
        .map(|line| {
            let (head, vals) = line.split_once('|').unwrap();
            let vals = vals
                .split_whitespace()
                .map(|v| v.parse().unwrap())
                .collect();
            (head.to_string(), vals)
        })
        .collect();
    view.sort_by(|a, b| a.0.cmp(&b.0));
    view
}

/// The same view of the file through clawhdf5-netcdf4.
fn clawhdf5_view(path: &std::path::Path) -> Vec<(String, Vec<f64>)> {
    fn describe(
        group_path: &str,
        vars: Vec<clawhdf5_netcdf4::Variable<'_>>,
    ) -> Vec<(String, Vec<f64>)> {
        vars.into_iter()
            .map(|v| {
                let dims: Vec<&str> = v.dimensions().iter().map(|d| d.name.as_str()).collect();
                let shape: Vec<String> = v.shape().unwrap().iter().map(u64::to_string).collect();
                let head = format!(
                    "{group_path} {} ({}) ({})",
                    v.name(),
                    dims.join(","),
                    shape.join(",")
                );
                let vals = match v.nc_type().unwrap() {
                    NcType::String | NcType::Char => Vec::new(),
                    _ => v.read_raw_f64().unwrap(),
                };
                (head, vals)
            })
            .collect()
    }
    fn walk(
        group_path: &str,
        group: &clawhdf5_netcdf4::NetCDF4Group<'_>,
        out: &mut Vec<(String, Vec<f64>)>,
    ) {
        out.extend(describe(group_path, group.variables().unwrap()));
        for name in group.group_names().unwrap() {
            walk(
                &format!("{group_path}/{name}"),
                &group.group(&name).unwrap(),
                out,
            );
        }
    }
    let file = NetCDF4File::open(path).unwrap();
    let mut view = describe("/", file.variables().unwrap());
    for name in file.group_names().unwrap() {
        walk(&format!("/{name}"), &file.group(&name).unwrap(), &mut view);
    }
    view.sort_by(|a, b| a.0.cmp(&b.0));
    view
}

/// clawhdf5-netcdf4 reports the same variables, dimensions, shapes and
/// values (bit for bit, NaN equal to NaN) as netCDF4-python.
fn assert_same_view(path: &std::path::Path) {
    let want = netcdf4_view(path);
    let got = clawhdf5_view(path);
    let heads = |v: &[(String, Vec<f64>)]| v.iter().map(|(h, _)| h.clone()).collect::<Vec<_>>();
    assert_eq!(heads(&got), heads(&want), "variables differ from netCDF4's");
    for ((head, got), (_, want)) in got.iter().zip(&want) {
        let same = got.len() == want.len()
            && got
                .iter()
                .zip(want)
                .all(|(a, b)| a.to_bits() == b.to_bits() || (a.is_nan() && b.is_nan()));
        assert!(same, "{head}: got {got:?}, netCDF4 reads {want:?}");
    }
}

/// The reproducer of the known-issues entry: `a` is on the unlimited `time`
/// (5 long through `b`) with 2 records, not on an anonymous `dim_2`; the
/// pure dimension scales `time` and `empty` are not variables; `a` has
/// shape (5,) and reads its 3 unwritten records as the fill value.
#[test]
fn variable_dimensions_come_from_the_file() {
    skip_if_no_netcdf4!();
    let dir = tempfile::tempdir().unwrap();
    let path = dir.path().join("repro.nc");
    run_python(&format!(
        r#"
import netCDF4 as nc
import numpy as np
with nc.Dataset({path:?}, "w") as f:
    f.createDimension("time", None)
    f.createDimension("empty", None)
    f.createDimension("x", 3)
    f.createVariable("a", "i4", ("time",))[0:2] = [1, 2]
    f.createVariable("b", "f4", ("time", "x"))[0:5, :] = np.arange(15).reshape(5, 3)
    f.createVariable("e", "i4", ("empty",))
    f.createVariable("c", "i4", ("x",))[:] = [7, 8, 9]
"#,
        path = path.display().to_string()
    ));
    assert_same_view(&path);

    let file = NetCDF4File::open(&path).unwrap();
    let mut names = file.variable_names().unwrap();
    names.sort();
    assert_eq!(names, ["a", "b", "c", "e"]);
    assert!(matches!(
        file.variable("time"),
        Err(clawhdf5_netcdf4::Error::VariableNotFound(_))
    ));
    let a = file.variable("a").unwrap();
    assert_eq!(a.dimensions()[0].name, "time");
    assert_eq!(a.shape().unwrap(), [5]);
    assert_eq!(a.stored_shape().unwrap(), [2]);
    assert_eq!(
        a.read_raw_i32().unwrap(),
        [1, 2, -2_147_483_647, -2_147_483_647, -2_147_483_647]
    );
}

/// Dimensions of one size are told apart by the file, not by order: `p`
/// and `q` are both 2 long, and `v(q, p)`, `same(p, p)` (one dimension
/// twice), a scalar, `q`'s coordinate variable, a variable called `p` that
/// is not `p`'s coordinate variable (stored as `_nc4_non_coord_p`), and
/// variables in a subgroup and a sub-subgroup on dimensions of their
/// ancestors.
#[test]
fn equal_size_and_inherited_dimensions_match_netcdf4_python() {
    skip_if_no_netcdf4!();
    let dir = tempfile::tempdir().unwrap();
    let path = dir.path().join("dims.nc");
    run_python(&format!(
        r#"
import netCDF4 as nc
import numpy as np
with nc.Dataset({path:?}, "w") as f:
    f.createDimension("p", 2)
    f.createDimension("q", 2)
    f.createVariable("v", "i4", ("q", "p"))[:] = np.array([[1, 2], [3, 4]])
    f.createVariable("same", "i4", ("p", "p"))[:] = np.array([[5, 6], [7, 8]])
    f.createVariable("s", "f8", ())[...] = 3.5
    f.createVariable("q", "f4", ("q",))[:] = [0, 1]
    f.createVariable("p", "f4", ("q", "p"))[:] = np.array([[0, 1], [2, 3]])
    g = f.createGroup("g")
    g.createDimension("r", 2)
    g.createVariable("w", "i4", ("r", "q", "p"))[:] = np.arange(8).reshape(2, 2, 2)
    h = g.createGroup("h")
    h.createVariable("z", "i4", ("p", "r"))[:] = np.array([[1, 2], [3, 4]])
"#,
        path = path.display().to_string()
    ));
    assert_same_view(&path);

    let file = NetCDF4File::open(&path).unwrap();
    let v = file.variable("v").unwrap();
    let dims: Vec<&str> = v.dimensions().iter().map(|d| d.name.as_str()).collect();
    assert_eq!(dims, ["q", "p"]);
    let p = file.variable("p").unwrap();
    assert!(!p.is_coordinate());
    assert!(file.variable("q").unwrap().is_coordinate());
    let s = file.variable("s").unwrap();
    assert!(s.dimensions().is_empty());
    assert_eq!(s.shape().unwrap(), Vec::<u64>::new());
    let z = file
        .group("g")
        .unwrap()
        .group("h")
        .unwrap()
        .variable("z")
        .unwrap();
    let dims: Vec<&str> = z.dimensions().iter().map(|d| d.name.as_str()).collect();
    assert_eq!(dims, ["p", "r"]);
}

/// Variables shorter than their unlimited dimension have its length and
/// read the fill value (`_FillValue`, else netCDF's default for the type)
/// where nothing was written — also when the unlimited dimension is not
/// the first; `read_f64` gives NaN there.
#[test]
fn unwritten_records_read_as_fill_like_netcdf4_python() {
    skip_if_no_netcdf4!();
    let dir = tempfile::tempdir().unwrap();
    let path = dir.path().join("pad.nc");
    run_python(&format!(
        r#"
import netCDF4 as nc
import numpy as np
with nc.Dataset({path:?}, "w") as f:
    f.createDimension("t", None)
    f.createDimension("x", 2)
    f.createVariable("a", "i4", ("t",))[0:2] = [1, 2]
    f.createVariable("f", "f4", ("x", "t"), fill_value=-5.0)[:, 0:1] = np.array([[1], [2]])
    f.createVariable("d", "f8", ("t",))[0:4] = [1, 2, 3, 4]
    f.createVariable("u", "u8", ("t",))[0:1] = [1]
    f.createVariable("b", "i1", ("t", "x"))[0:3, :] = np.ones((3, 2))
    f.createVariable("st", str, ("t",))[0] = "hi"
    g = f.createGroup("g")
    g.createVariable("k", "f4", ("t",))[0:1] = [9]
"#,
        path = path.display().to_string()
    ));
    assert_same_view(&path);

    let file = NetCDF4File::open(&path).unwrap();
    let mut f = file.variable("f").unwrap();
    assert_eq!(f.shape().unwrap(), [2, 4]);
    assert_eq!(f.stored_shape().unwrap(), [2, 1]);
    assert_eq!(
        f.read_raw_f32().unwrap(),
        [1.0, -5.0, -5.0, -5.0, 2.0, -5.0, -5.0, -5.0]
    );
    let read = f.read_f64().unwrap();
    assert_eq!(read[0], 1.0);
    assert!(read[1].is_nan() && read[7].is_nan());
    let st = file.variable("st").unwrap();
    assert_eq!(st.read_string().unwrap(), ["hi", "", "", ""]);
    assert_eq!(st.shape().unwrap(), [4]);
}

/// A file with HDF5 dimension scales but none of netCDF's own attributes
/// (h5py's `dims` API): the dimensions come from `DIMENSION_LIST`, so
/// `v(q, p)` is not `v(p, q)` although both are 2 long; with two scales
/// attached to one axis (`w`), netCDF-C takes the last.
#[test]
fn h5py_dimension_scales_match_netcdf4_python() {
    skip_if_no_netcdf4!();
    let dir = tempfile::tempdir().unwrap();
    let path = dir.path().join("scales.h5");
    run_python(&format!(
        r#"
import h5py
import numpy as np
with h5py.File({path:?}, "w") as f:
    f["p"] = np.arange(2.0)
    f["q"] = np.arange(2.0) + 10
    f["p"].make_scale("p")
    f["q"].make_scale("q")
    f["v"] = np.arange(4).reshape(2, 2)
    f["v"].dims[0].attach_scale(f["q"])
    f["v"].dims[1].attach_scale(f["p"])
    f["w"] = np.arange(2)
    f["w"].dims[0].attach_scale(f["p"])
    f["w"].dims[0].attach_scale(f["q"])
"#,
        path = path.display().to_string()
    ));
    assert_same_view(&path);
}

/// Files h5netcdf writes (its own implementation of the netCDF-4
/// conventions over h5py): an unlimited dimension, equal sizes, a subgroup
/// on inherited dimensions, a scalar.
#[test]
fn h5netcdf_file_matches_netcdf4_python() {
    skip_if_no_netcdf4!();
    skip_if_no_h5netcdf!();
    let dir = tempfile::tempdir().unwrap();
    let path = dir.path().join("h5netcdf.nc");
    run_python(&format!(
        r#"
import h5netcdf
import numpy as np
with h5netcdf.File({path:?}, "w") as f:
    f.dimensions = {{"p": 2, "q": 2, "t": None}}
    f.create_variable("v", ("q", "p"), "i4")[...] = np.array([[1, 2], [3, 4]])
    f.create_variable("q", ("q",), "f4")[...] = [0, 1]
    f.create_variable("same", ("p", "p"), "i4")[...] = np.array([[5, 6], [7, 8]])
    a = f.create_variable("a", ("t", "p"), "f8")
    f.resize_dimension("t", 3)
    a[...] = np.ones((3, 2))
    f.create_variable("short", ("t",), "i4")
    g = f.create_group("g")
    g.dimensions = {{"r": 2}}
    g.create_variable("w", ("r", "q", "p"), "i4")[...] = np.arange(8).reshape(2, 2, 2)
    g.create_variable("s", (), "f8")[...] = 2.5
"#,
        path = path.display().to_string()
    ));
    assert_same_view(&path);
}

/// Files xarray writes, through netCDF4 and (when installed) h5netcdf:
/// coordinates, two dimensions of one size, an unlimited dimension.
#[test]
fn xarray_files_match_netcdf4_python() {
    skip_if_no_netcdf4!();
    skip_if_no_xarray!();
    let dir = tempfile::tempdir().unwrap();
    let mut engines = vec!["netcdf4"];
    if python_has("h5netcdf") {
        engines.push("h5netcdf");
    } else {
        eprintln!("SKIP: xarray with engine h5netcdf (h5netcdf not available)");
    }
    for engine in engines {
        let path = dir.path().join(format!("xarray_{engine}.nc"));
        run_python(&format!(
            r#"
import numpy as np
import xarray as xr
ds = xr.Dataset(
    {{
        "temp": (("time", "lat", "lon"), np.arange(12.0).reshape(3, 2, 2)),
        "grid": (("lon", "lat"), np.array([[1, 2], [3, 4]], dtype="i4")),
        "scalar": ((), 1.5),
    }},
    coords={{"time": [0.0, 6.0, 12.0], "lat": [10.0, 20.0], "lon": [5.0, 6.0]}},
)
ds.to_netcdf({path:?}, engine={engine:?}, unlimited_dims=["time"])
"#,
            path = path.display().to_string()
        ));
        assert_same_view(&path);
    }
}