clawhdf5-netcdf4: variables' dimensions come from the file

Variables got the first unused dimension of equal size, so a variable on
an unlimited dimension with fewer records got an anonymous dim_<n>, and
dimensions of one size could be swapped. Resolve them as netCDF-C does
(libhdf5/hdf5open.c): _Netcdf4Coordinates ids, else the scales
DIMENSION_LIST references (the last one attached to an axis), searched in
the variable's group and its parents; a coordinate variable is on its own
scale. Size matching remains only for axes the file names nothing for.

variables()/variable_names() leave out dimension scales that are only
dimensions, and _nc4_non_coord_<name> is the variable <name>.
Variable::shape is the netCDF shape (an unlimited dimension's length) and
the reads pad unwritten records with the fill value (_FillValue, else
NC_FILL_*; NaN from read_f64); Variable::stored_shape is the HDF5 extent.
New NetCDF4File::variable_names.

Tests compare with netCDF4-python variable by variable: the known-issues
reproducer, equal sizes, (p, p), scalars, inherited dimensions, unwritten
records, h5py dimension scales, h5netcdf and xarray files. CI installs
h5netcdf. known-issues entry moved to Fixed (history); stale open-table
row for the unlimited-size fix removed.

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>

crates/clawhdf5-netcdf4/src/dimension.rs

@@ -22,73 +22,123 @@ pub struct Dimension {
     pub is_unlimited: bool,
 }
 
-/// Extract dimensions from an HDF5 group (root or subgroup).
+/// A dimension scale of one group: the dataset that defines a dimension.
+#[derive(Debug, Clone)]
+pub(crate) struct Scale {
+    /// Object header address of the scale's dataset (what a variable's
+    /// `DIMENSION_LIST` references).
+    pub address: u64,
+    /// Its `_Netcdf4Dimid` (what a variable's `_Netcdf4Coordinates` lists).
+    pub dimid: Option<i64>,
+    /// Index of its dimension in [`GroupDims::dims`].
+    pub dim: usize,
+}
+
+/// The dimensions a group defines, with the scales that define them.
+#[derive(Debug, Clone, Default)]
+pub(crate) struct GroupDims {
+    /// The group's dimensions, in `_Netcdf4Dimid` order (then discovery order).
+    pub dims: Vec<Dimension>,
+    /// The dimension scales behind `dims`; empty when the group has no
+    /// dimension scale and `dims` were inferred from 1-D datasets.
+    pub scales: Vec<Scale>,
+}
+
+impl GroupDims {
+    /// The dimension defined by the scale at `address`.
+    pub fn by_address(&self, address: u64) -> Option<&Dimension> {
+        self.scales
+            .iter()
+            .find(|s| s.address == address)
+            .map(|s| &self.dims[s.dim])
+    }
+
+    /// The dimension whose scale has `_Netcdf4Dimid` `id`.
+    pub fn by_dimid(&self, id: i64) -> Option<&Dimension> {
+        self.scales
+            .iter()
+            .find(|s| s.dimid == Some(id))
+            .map(|s| &self.dims[s.dim])
+    }
+}
+
+/// The dimensions of an HDF5 group (root or subgroup).
 ///
 /// NetCDF-4 stores dimensions as datasets with `CLASS=DIMENSION_SCALE`. A fixed
 /// dimension's size is the dataset's first (and typically only) shape extent.
 /// Unlimited dimensions have `max_dimensions[0] == u64::MAX` in the HDF5 dataspace;
-/// their size is computed by `unlimited_len`.
-pub(crate) fn extract_dimensions(
+/// their size is computed by `unlimited_len`. A group with no dimension
+/// scale at all (not written by a netCDF library) gets one dimension per
+/// 1-D dataset instead.
+pub(crate) fn group_dims(
     file: &clawhdf5::File,
     group: &clawhdf5::Group<'_>,
-) -> Result<Vec<Dimension>, Error> {
+) -> Result<GroupDims, Error> {
+    let addresses: HashMap<String, u64> = group.entries()?.into_iter().collect();
     let dataset_names = group.datasets()?;
-    let mut dims = Vec::new();
-    let mut seen_dimids: HashMap<i64, usize> = HashMap::new();
+    // (dimid, dimension, scale address), in discovery order.
+    let mut found: Vec<(Option<i64>, Dimension, u64)> = Vec::new();
 
     for ds_name in &dataset_names {
         let ds = group.dataset(ds_name)?;
         let attrs = ds.attrs()?;
-
-        // Check if this is a dimension scale
         if !is_dimension_scale(&attrs) {
             continue;
         }
-
+        let Some(&address) = addresses.get(ds_name) else {
+            continue;
+        };
         let shape = ds.shape()?;
-        let is_unlimited = check_unlimited(file, group, ds_name);
+        let is_unlimited = is_unlimited(&ds);
         let size = if is_unlimited {
             unlimited_len(file, &attrs, &shape)
         } else {
             shape.first().copied().unwrap_or(0)
         };
-
-        let dimid = get_dimid(&attrs);
-
         let dim = Dimension {
             name: ds_name.clone(),
             size,
             is_unlimited,
         };
-
-        if let Some(id) = dimid {
-            seen_dimids.insert(id, dims.len());
-        }
-        dims.push(dim);
+        found.push((get_dimid(&attrs), dim, address));
     }
 
-    // Sort by dimid if available, otherwise keep discovery order
-    if !seen_dimids.is_empty() {
-        let mut pairs: Vec<(i64, Dimension)> = Vec::new();
-        let mut unordered = Vec::new();
-
-        for (i, dim) in dims.into_iter().enumerate() {
-            let id = seen_dimids
-                .iter()
-                .find(|(_, idx)| **idx == i)
-                .map(|(k, _)| *k);
-            if let Some(id) = id {
-                pairs.push((id, dim));
-            } else {
-                unordered.push(dim);
+    if found.is_empty() {
+        // Fallback: infer dimensions from dataset shapes and names.
+        // In NetCDF-4, coordinate variables are datasets whose name matches
+        // a dimension name. If there are no explicit DIMENSION_SCALE attributes,
+        // we look for 1-D datasets that might be coordinate variables.
+        let mut dims = Vec::new();
+        for ds_name in &dataset_names {
+            let ds = group.dataset(ds_name)?;
+            let shape = ds.shape()?;
+            if shape.len() == 1 {
+                dims.push(Dimension {
+                    name: ds_name.clone(),
+                    size: shape[0],
+                    is_unlimited: is_unlimited(&ds),
+                });
             }
         }
-        pairs.sort_by_key(|(id, _)| *id);
-        dims = pairs.into_iter().map(|(_, d)| d).collect();
-        dims.extend(unordered);
+        return Ok(GroupDims {
+            dims,
+            scales: Vec::new(),
+        });
     }
 
-    Ok(dims)
+    // By dimid; scales without one keep their discovery order after those
+    // with one (the sort is stable).
+    found.sort_by_key(|(id, ..)| (id.is_none(), id.unwrap_or(0)));
+    let mut out = GroupDims::default();
+    for (i, (dimid, dim, address)) in found.into_iter().enumerate() {
+        out.dims.push(dim);
+        out.scales.push(Scale {
+            address,
+            dimid,
+            dim: i,
+        });
+    }
+    Ok(out)
 }
 
 /// The start of the `NAME` attribute netCDF-C gives a dimension scale that
@@ -107,10 +157,7 @@ const PURE_DIMENSION_NAME: &str = "This is a netCDF dimension but not a netCDF v
 /// scale's own extent, as before.
 fn unlimited_len(file: &clawhdf5::File, attrs: &HashMap<String, AttrValue>, shape: &[u64]) -> u64 {
     let own = shape.first().copied().unwrap_or(0);
-    let is_variable = !matches!(
-        attrs.get("NAME"),
-        Some(AttrValue::String(n)) if n.starts_with(PURE_DIMENSION_NAME)
-    );
+    let is_variable = !is_pure_dimension(attrs);
     let Some(refs) = reference_list(file, attrs) else {
         return own;
     };
@@ -171,8 +218,59 @@ fn reference_list(
     Some(addresses.into_iter().map(|r| r.address).zip(axes).collect())
 }
 
+/// The dimension scale attached to each axis of a variable, from its
+/// `DIMENSION_LIST` attribute (HDF5 dimension scales: one variable-length
+/// sequence of object references per axis) — the address of the scale
+/// netCDF-C takes for the axis, or `None` for an axis with none. netCDF-C's
+/// `dimscale_visitor` lets `H5DSiterate_scales` visit every scale attached
+/// to the axis and keeps the last, so with several (h5py's `attach_scale`
+/// twice) the last one is the axis's dimension. `None` overall when the
+/// attribute is missing or not in that form.
+pub(crate) fn dimension_list(
+    file: &clawhdf5::File,
+    attrs: &HashMap<String, AttrValue>,
+) -> Option<Vec<Option<u64>>> {
+    use clawhdf5_format::data_read::read_object_references;
+    use clawhdf5_format::datatype::Datatype;
+    use clawhdf5_format::vl_data::VlResolver;
+    let Some(AttrValue::Raw { datatype, data, .. }) = attrs.get("DIMENSION_LIST") else {
+        return None;
+    };
+    let Datatype::VariableLength {
+        is_string: false,
+        base_type,
+        ..
+    } = datatype
+    else {
+        return None;
+    };
+    let sb = file.superblock();
+    let base_size = usize::try_from(base_type.type_size()).ok()?;
+    let sequences = VlResolver::new_in(file.storage(), sb.offset_size, sb.length_size)
+        .sequences(data, base_size)
+        .ok()?;
+    sequences
+        .iter()
+        .map(|refs| {
+            let refs = read_object_references(refs, base_type, sb.offset_size).ok()?;
+            Some(refs.iter().rev().find(|r| !r.is_null()).map(|r| r.address))
+        })
+        .collect()
+}
+
+/// Whether a dataset is a dimension scale that is only a dimension, not a
+/// netCDF variable: netCDF-C and h5netcdf give it this `NAME`, and netCDF-C
+/// does not list it among the variables.
+pub(crate) fn is_pure_dimension(attrs: &HashMap<String, AttrValue>) -> bool {
+    is_dimension_scale(attrs)
+        && matches!(
+            attrs.get("NAME"),
+            Some(AttrValue::String(n)) if n.starts_with(PURE_DIMENSION_NAME)
+        )
+}
+
 /// Check if a dataset's attributes mark it as a dimension scale.
-fn is_dimension_scale(attrs: &HashMap<String, AttrValue>) -> bool {
+pub(crate) fn is_dimension_scale(attrs: &HashMap<String, AttrValue>) -> bool {
     if let Some(AttrValue::String(class)) = attrs.get("CLASS") {
         return class == "DIMENSION_SCALE";
     }
@@ -180,7 +278,7 @@ fn is_dimension_scale(attrs: &HashMap<String, AttrValue>) -> bool {
 }
 
 /// Get the _Netcdf4Dimid attribute value if present.
-fn get_dimid(attrs: &HashMap<String, AttrValue>) -> Option<i64> {
+pub(crate) fn get_dimid(attrs: &HashMap<String, AttrValue>) -> Option<i64> {
     match attrs.get("_Netcdf4Dimid") {
         Some(AttrValue::I64(id)) => Some(*id),
         Some(AttrValue::U64(id)) => Some(*id as i64),
@@ -188,56 +286,8 @@ fn get_dimid(attrs: &HashMap<String, AttrValue>) -> Option<i64> {
     }
 }
 
-/// Check if a dimension is unlimited by inspecting the HDF5 dataspace max_dimensions.
-///
-/// A dimension is unlimited when `max_dimensions[0] == u64::MAX` in the HDF5 dataspace.
-fn check_unlimited(_file: &clawhdf5::File, group: &clawhdf5::Group<'_>, ds_name: &str) -> bool {
-    let ds = match group.dataset(ds_name) {
-        Ok(ds) => ds,
-        Err(_) => return false,
-    };
-
-    match ds.max_dimensions() {
-        Ok(Some(max_dims)) => max_dims.first().copied() == Some(u64::MAX),
-        _ => false,
-    }
-}
-
-/// Extract dimensions from an HDF5 group using both dimension scale attributes
-/// and variable DIMENSION_LIST references.
-///
-/// This is a more robust approach that also discovers dimensions from variables
-/// that reference them, even when dimension scales aren't explicitly set.
-pub(crate) fn extract_dimensions_from_datasets(
-    group: &clawhdf5::Group<'_>,
-    file: &clawhdf5::File,
-) -> Result<Vec<Dimension>, Error> {
-    // First try the standard approach with DIMENSION_SCALE
-    let mut dims = extract_dimensions(file, group)?;
-
-    // If we found dimensions, return them
-    if !dims.is_empty() {
-        return Ok(dims);
-    }
-
-    // Fallback: infer dimensions from dataset shapes and names.
-    // In NetCDF-4, coordinate variables are datasets whose name matches
-    // a dimension name. If there are no explicit DIMENSION_SCALE attributes,
-    // we look for 1-D datasets that might be coordinate variables.
-    let dataset_names = group.datasets()?;
-    for ds_name in &dataset_names {
-        let ds = group.dataset(ds_name)?;
-        let shape = ds.shape()?;
-        if shape.len() == 1 {
-            // This 1-D dataset could be a coordinate variable / dimension
-            let is_unlimited = check_unlimited(file, group, ds_name);
-            dims.push(Dimension {
-                name: ds_name.clone(),
-                size: shape[0],
-                is_unlimited,
-            });
-        }
-    }
-
-    Ok(dims)
+/// Whether a dataset's first axis is unlimited (`max_dimensions[0] ==
+/// u64::MAX` in its dataspace).
+fn is_unlimited(ds: &clawhdf5::Dataset<'_>) -> bool {
+    matches!(ds.max_dimensions(), Ok(Some(max_dims)) if max_dims.first() == Some(&u64::MAX))
 }

crates/clawhdf5-netcdf4/src/group.rs

@@ -9,12 +9,15 @@ use clawhdf5::AttrValue;
 
 use crate::dimension::{self, Dimension};
 use crate::error::Error;
-use crate::variable::{self, Variable};
+use crate::scope::{self, Scope};
+use crate::variable::Variable;
 
 /// A NetCDF-4 group corresponding to an HDF5 group.
 pub struct NetCDF4Group<'f> {
     /// Group name.
     name: String,
+    /// Path of the group from the root (`/`-separated).
+    path: String,
     /// Underlying HDF5 file.
     file: &'f clawhdf5::File,
     /// Underlying HDF5 group.
@@ -25,11 +28,13 @@ impl<'f> NetCDF4Group<'f> {
     /// Create a new NetCDF4Group from an HDF5 group.
     pub(crate) fn new(
         name: String,
+        path: String,
         file: &'f clawhdf5::File,
         hdf5_group: clawhdf5::Group<'f>,
     ) -> Self {
         Self {
             name,
+            path,
             file,
             hdf5_group,
         }
@@ -40,27 +45,22 @@ impl<'f> NetCDF4Group<'f> {
         &self.name
     }
 
-    /// List dimensions defined in this group.
+    /// List dimensions defined in this group (not those of its parent
+    /// groups, which its variables can also use).
     pub fn dimensions(&self) -> Result<Vec<Dimension>, Error> {
-        dimension::extract_dimensions_from_datasets(&self.hdf5_group, self.file)
+        Ok(dimension::group_dims(self.file, &self.hdf5_group)?.dims)
     }
 
-    /// List variables in this group.
+    /// List variables in this group: its datasets, except the dimension
+    /// scales that are only dimensions. Their dimensions can be defined in
+    /// this group or a parent group.
     pub fn variables(&self) -> Result<Vec<Variable<'f>>, Error> {
-        let dims = self.dimensions()?;
-        variable::build_variables(&self.hdf5_group, &dims)
+        Scope::new(self.file, &self.path)?.variables()
     }
 
     /// Get a specific variable by name.
     pub fn variable(&self, name: &str) -> Result<Variable<'f>, Error> {
-        let dims = self.dimensions()?;
-        let ds = self
-            .hdf5_group
-            .dataset(name)
-            .map_err(|_| Error::VariableNotFound(name.to_string()))?;
-        let shape = ds.shape()?;
-        let var_dims = crate::variable::match_dimensions_to_variable(&shape, &dims);
-        Ok(Variable::new(name.to_string(), ds, var_dims))
+        scope::variable_at(self.file, &self.path, name)
     }
 
     /// Read all attributes of this group.
@@ -79,12 +79,18 @@ impl<'f> NetCDF4Group<'f> {
             .hdf5_group
             .group(name)
             .map_err(|_| Error::GroupNotFound(name.to_string()))?;
-        Ok(NetCDF4Group::new(name.to_string(), self.file, hdf5_group))
+        Ok(NetCDF4Group::new(
+            name.to_string(),
+            format!("{}/{name}", self.path),
+            self.file,
+            hdf5_group,
+        ))
     }
 
-    /// List dataset (variable) names in this group.
+    /// The names of this group's variables (see
+    /// [`variables`](Self::variables)).
     pub fn variable_names(&self) -> Result<Vec<String>, Error> {
-        Ok(self.hdf5_group.datasets()?)
+        Scope::new(self.file, &self.path)?.variable_names()
     }
 }
 

crates/clawhdf5-netcdf4/src/lib.rs

@@ -27,6 +27,7 @@ pub mod cf;
 pub mod dimension;
 pub mod error;
 pub mod group;
+mod scope;
 pub mod types;
 pub mod variable;
 
@@ -75,25 +76,25 @@ impl NetCDF4File {
 
     /// List dimensions defined in the root group.
     pub fn dimensions(&self) -> Result<Vec<Dimension>, Error> {
-        dimension::extract_dimensions_from_datasets(&self.hdf5.root(), &self.hdf5)
+        Ok(dimension::group_dims(&self.hdf5, &self.hdf5.root())?.dims)
     }
 
-    /// List all variables in the root group.
+    /// List all variables in the root group: its datasets, except the
+    /// dimension scales that are only dimensions (netCDF-C does not list
+    /// them either).
     pub fn variables(&self) -> Result<Vec<Variable<'_>>, Error> {
-        let dims = self.dimensions()?;
-        variable::build_variables(&self.hdf5.root(), &dims)
+        scope::Scope::new(&self.hdf5, "/")?.variables()
+    }
+
+    /// The names of the root group's variables (see
+    /// [`variables`](Self::variables)).
+    pub fn variable_names(&self) -> Result<Vec<String>, Error> {
+        scope::Scope::new(&self.hdf5, "/")?.variable_names()
     }
 
     /// Get a specific variable by name from the root group.
     pub fn variable(&self, name: &str) -> Result<Variable<'_>, Error> {
-        let dims = self.dimensions()?;
-        let ds = self
-            .hdf5
-            .dataset(name)
-            .map_err(|_| Error::VariableNotFound(name.to_string()))?;
-        let shape = ds.shape()?;
-        let var_dims = variable::match_dimensions_to_variable(&shape, &dims);
-        Ok(Variable::new(name.to_string(), ds, var_dims))
+        scope::variable_at(&self.hdf5, "", name)
     }
 
     /// Read all global (root group) attributes.
@@ -112,7 +113,12 @@ impl NetCDF4File {
             .hdf5
             .group(name)
             .map_err(|_| Error::GroupNotFound(name.to_string()))?;
-        Ok(NetCDF4Group::new(name.to_string(), &self.hdf5, hdf5_group))
+        Ok(NetCDF4Group::new(
+            name.to_string(),
+            name.to_string(),
+            &self.hdf5,
+            hdf5_group,
+        ))
     }
 
     /// Access the underlying HDF5 file for advanced operations.

crates/clawhdf5-netcdf4/src/scope.rs

@@ -0,0 +1,231 @@
+//! A group's variables and the dimensions they are defined on.
+//!
+//! netCDF-C (`libhdf5/hdf5open.c`) gives a variable its dimensions from the
+//! file, never by size: the dimension ids in its `_Netcdf4Coordinates`
+//! attribute (each dimension scale's `_Netcdf4Dimid`), else the dimension
+//! scales its `DIMENSION_LIST` attribute references, looked up in the
+//! variable's group and then each parent group up to the root. Only an axis
+//! with neither (a file not written by a netCDF library) gets a dimension
+//! by size. Dimension scales that are only dimensions are not variables, and
+//! a variable stored as `_nc4_non_coord_<name>` (a variable sharing a
+//! dimension's name without being its coordinate variable) is `<name>`.
+
+use std::collections::{HashMap, HashSet};
+
+use clawhdf5::AttrValue;
+
+use crate::dimension::{self, Dimension, GroupDims};
+use crate::error::Error;
+use crate::variable::Variable;
+
+/// The prefix netCDF-C gives the dataset of a variable that has a
+/// dimension's name but is not that dimension's coordinate variable (the
+/// dimension's scale holds the name).
+const NON_COORD_PREFIX: &str = "_nc4_non_coord_";
+
+/// A group, with the dimensions visible from it.
+pub(crate) struct Scope<'f> {
+    file: &'f clawhdf5::File,
+    group: clawhdf5::Group<'f>,
+    /// This group's dimensions, then its parent's, and so on to the root's.
+    levels: Vec<GroupDims>,
+}
+
+impl<'f> Scope<'f> {
+    /// The group at `path` (`/`-separated from the root; `""` or `"/"` is
+    /// the root).
+    pub fn new(file: &'f clawhdf5::File, path: &str) -> Result<Self, Error> {
+        let parts: Vec<&str> = path.split('/').filter(|p| !p.is_empty()).collect();
+        let mut levels = Vec::with_capacity(parts.len() + 1);
+        for n in (0..=parts.len()).rev() {
+            let group = file.group(&parts[..n].join("/"))?;
+            levels.push(dimension::group_dims(file, &group)?);
+        }
+        let group = file.group(&parts.join("/"))?;
+        Ok(Self {
+            file,
+            group,
+            levels,
+        })
+    }
+
+    /// The group's datasets, as `(dataset name, object header address)` in
+    /// listing order.
+    fn datasets(&self) -> Result<Vec<(String, u64)>, Error> {
+        let datasets: HashSet<String> = self.group.datasets()?.into_iter().collect();
+        Ok(self
+            .group
+            .entries()?
+            .into_iter()
+            .filter(|(name, _)| datasets.contains(name))
+            .collect())
+    }
+
+    /// The group's variables: every dataset but the dimension scales that
+    /// are only dimensions.
+    pub fn variables(&self) -> Result<Vec<Variable<'f>>, Error> {
+        let mut variables = Vec::new();
+        for (ds_name, address) in self.datasets()? {
+            let ds = self.file.dataset_at(address)?;
+            let attrs = ds.attrs()?;
+            if dimension::is_pure_dimension(&attrs) {
+                continue;
+            }
+            variables.push(self.variable_from(nc_name(&ds_name), address, ds, attrs)?);
+        }
+        Ok(variables)
+    }
+
+    /// The names of the group's variables.
+    pub fn variable_names(&self) -> Result<Vec<String>, Error> {
+        let mut names = Vec::new();
+        for (ds_name, address) in self.datasets()? {
+            let attrs = self.file.dataset_at(address)?.attrs()?;
+            if !dimension::is_pure_dimension(&attrs) {
+                names.push(nc_name(&ds_name));
+            }
+        }
+        Ok(names)
+    }
+
+    /// The variable called `name`: the dataset `_nc4_non_coord_<name>` if
+    /// there is one, else the dataset `<name>` unless it is only a
+    /// dimension.
+    pub fn variable(&self, name: &str) -> Result<Variable<'f>, Error> {
+        let not_found = || Error::VariableNotFound(name.to_string());
+        let datasets = self.datasets()?;
+        let prefixed = format!("{NON_COORD_PREFIX}{name}");
+        let address = datasets
+            .iter()
+            .find(|(n, _)| *n == prefixed)
+            .or_else(|| datasets.iter().find(|(n, _)| n == name))
+            .map(|&(_, address)| address)
+            .ok_or_else(not_found)?;
+        let ds = self.file.dataset_at(address)?;
+        let attrs = ds.attrs()?;
+        if dimension::is_pure_dimension(&attrs) {
+            return Err(not_found());
+        }
+        self.variable_from(nc_name(name), address, ds, attrs)
+    }
+
+    fn variable_from(
+        &self,
+        name: String,
+        address: u64,
+        ds: clawhdf5::Dataset<'f>,
+        attrs: HashMap<String, AttrValue>,
+    ) -> Result<Variable<'f>, Error> {
+        let shape = ds.shape()?;
+        let dims = self.variable_dims(address, &attrs, &shape);
+        Ok(Variable::new(name, ds, dims, attrs))
+    }
+
+    /// The first dimension, searching this group and then its ancestors,
+    /// that `find` picks.
+    fn find<'a>(
+        &'a self,
+        find: impl Fn(&'a GroupDims) -> Option<&'a Dimension>,
+    ) -> Option<Dimension> {
+        self.levels.iter().find_map(find).cloned()
+    }
+
+    /// The dimensions of the dataset at `address`, one per axis of `shape`,
+    /// as netCDF-C resolves them (see the module docs).
+    fn variable_dims(
+        &self,
+        address: u64,
+        attrs: &HashMap<String, AttrValue>,
+        shape: &[u64],
+    ) -> Vec<Dimension> {
+        let rank = shape.len();
+        let mut dims: Vec<Option<Dimension>> = vec![None; rank];
+        if rank == 0 {
+            return Vec::new();
+        }
+        // A coordinate variable is the scale of its (first) dimension.
+        dims[0] = self.levels[0].by_address(address).cloned();
+
+        if let Some(ids) = coordinates(attrs).filter(|ids| ids.len() == rank) {
+            for (slot, id) in dims.iter_mut().zip(ids) {
+                if slot.is_none() {
+                    *slot = self.find(|level| level.by_dimid(id));
+                }
+            }
+        }
+        if dims.iter().any(Option::is_none)
+            && let Some(scales) =
+                dimension::dimension_list(self.file, attrs).filter(|s| s.len() == rank)
+        {
+            for (slot, scale) in dims.iter_mut().zip(scales) {
+                if slot.is_none()
+                    && let Some(scale) = scale
+                {
+                    *slot = self.find(|level| level.by_address(scale));
+                }
+            }
+        }
+
+        // Neither: the first dimension of this group of the same size not
+        // already taken by another such axis, else an anonymous one.
+        let own = &self.levels[0].dims;
+        let mut used = vec![false; own.len()];
+        dims.into_iter()
+            .zip(shape)
+            .map(|(dim, &size)| {
+                dim.unwrap_or_else(|| {
+                    match own
+                        .iter()
+                        .enumerate()
+                        .find(|&(i, d)| !used[i] && d.size == size)
+                    {
+                        Some((i, d)) => {
+                            used[i] = true;
+                            d.clone()
+                        }
+                        None => Dimension {
+                            name: format!("dim_{size}"),
+                            size,
+                            is_unlimited: false,
+                        },
+                    }
+                })
+            })
+            .collect()
+    }
+}
+
+/// The variable `name` of the group at `group_path`; `name` may itself be
+/// a path (`"sub/var"`), relative to that group.
+pub(crate) fn variable_at<'f>(
+    file: &'f clawhdf5::File,
+    group_path: &str,
+    name: &str,
+) -> Result<Variable<'f>, Error> {
+    match name.trim_start_matches('/').rsplit_once('/') {
+        Some((dir, leaf)) => Scope::new(file, &format!("{group_path}/{dir}"))
+            .map_err(|_| Error::VariableNotFound(name.to_string()))?
+            .variable(leaf),
+        None => Scope::new(file, group_path)?.variable(name.trim_start_matches('/')),
+    }
+}
+
+/// The netCDF name of the dataset `ds_name`.
+fn nc_name(ds_name: &str) -> String {
+    ds_name
+        .strip_prefix(NON_COORD_PREFIX)
+        .unwrap_or(ds_name)
+        .to_string()
+}
+
+/// A variable's `_Netcdf4Coordinates`: the `_Netcdf4Dimid` of the dimension
+/// of each axis.
+fn coordinates(attrs: &HashMap<String, AttrValue>) -> Option<Vec<i64>> {
+    match attrs.get("_Netcdf4Coordinates")? {
+        AttrValue::I64Array(ids) => Some(ids.clone()),
+        AttrValue::I64(id) => Some(vec![*id]),
+        AttrValue::U64Array(ids) => ids.iter().map(|&id| i64::try_from(id).ok()).collect(),
+        AttrValue::U64(id) => Some(vec![i64::try_from(*id).ok()?]),
+        _ => None,
+    }
+}

crates/clawhdf5-netcdf4/src/variable.rs

@@ -2,12 +2,18 @@
 //!
 //! Variables in NetCDF-4 are HDF5 datasets. This module wraps them with
 //! dimension associations and CF attribute support.
+//!
+//! A variable along an unlimited dimension has that dimension's length in
+//! netCDF, even when fewer records of it have been written (its HDF5 dataset
+//! is shorter): [`Variable::shape`] is the netCDF shape and the reads return
+//! that many values, the unwritten ones as the fill value, as netCDF-C does.
+//! [`Variable::stored_shape`] is the dataset's extent.
 
 use std::collections::HashMap;
 
 use clawhdf5::AttrValue;
 
-use crate::cf::{self, CfAttributes};
+use crate::cf::{self, CfAttributes, FillValue};
 use crate::dimension::Dimension;
 use crate::error::Error;
 use crate::types::{NcType, dtype_to_nctype};
@@ -20,18 +26,23 @@ pub struct Variable<'f> {
     dataset: clawhdf5::Dataset<'f>,
     /// Dimensions associated with this variable.
     dims: Vec<Dimension>,
-    /// Cached attributes.
-    attrs_cache: Option<HashMap<String, AttrValue>>,
+    /// The dataset's attributes.
+    attrs: HashMap<String, AttrValue>,
 }
 
 impl<'f> Variable<'f> {
     /// Create a new Variable wrapping an HDF5 dataset.
-    pub(crate) fn new(name: String, dataset: clawhdf5::Dataset<'f>, dims: Vec<Dimension>) -> Self {
+    pub(crate) fn new(
+        name: String,
+        dataset: clawhdf5::Dataset<'f>,
+        dims: Vec<Dimension>,
+        attrs: HashMap<String, AttrValue>,
+    ) -> Self {
         Self {
             name,
             dataset,
             dims,
-            attrs_cache: None,
+            attrs,
         }
     }
 
@@ -40,13 +51,27 @@ impl<'f> Variable<'f> {
         &self.name
     }
 
-    /// The dimensions of this variable.
+    /// The dimensions of this variable, one per axis: the ones the file
+    /// gives it (`_Netcdf4Coordinates`, else `DIMENSION_LIST`), found in its
+    /// group or a parent group. An axis the file gives no dimension (a file
+    /// not written by a netCDF library) gets the first dimension of the
+    /// variable's group of the same size, else an anonymous `dim_<size>`.
     pub fn dimensions(&self) -> &[Dimension] {
         &self.dims
     }
 
-    /// The shape of this variable (dimension sizes).
+    /// The shape of this variable as netCDF reports it: along an unlimited
+    /// dimension, the dimension's current length (the longest variable on
+    /// it), even if fewer records of this variable have been written;
+    /// otherwise the dataset's extent. The reads return this many values.
     pub fn shape(&self) -> Result<Vec<u64>, Error> {
+        Ok(nc_shape(&self.dataset.shape()?, &self.dims))
+    }
+
+    /// The extent of the HDF5 dataset: what has been written. It differs
+    /// from [`shape`](Self::shape) only along an unlimited dimension that
+    /// another variable has more records of.
+    pub fn stored_shape(&self) -> Result<Vec<u64>, Error> {
         Ok(self.dataset.shape()?)
     }
 
@@ -58,59 +83,88 @@ impl<'f> Variable<'f> {
 
     /// Read all attributes as a HashMap.
     pub fn attrs(&mut self) -> Result<&HashMap<String, AttrValue>, Error> {
-        if self.attrs_cache.is_none() {
-            self.attrs_cache = Some(self.dataset.attrs()?);
-        }
-        Ok(self
-            .attrs_cache
-            .as_ref()
-            .expect("invariant: attrs_cache is Some after initialization"))
+        Ok(&self.attrs)
     }
 
     /// Extract CF convention attributes.
     pub fn cf_attributes(&mut self) -> Result<CfAttributes, Error> {
-        let attrs = self.attrs()?;
-        Ok(cf::extract_cf_attributes(attrs))
+        Ok(cf::extract_cf_attributes(&self.attrs))
     }
 
     /// Read data as f64 with scale_factor/add_offset applied.
     ///
     /// Missing values (matching `_FillValue` or `missing_value`) become NaN.
     /// If no scale_factor or add_offset attributes exist, returns the raw f64 data.
+    /// Records along an unlimited dimension that this variable has not
+    /// written (see [`shape`](Self::shape)) are NaN.
     pub fn read_f64(&mut self) -> Result<Vec<f64>, Error> {
         let raw = self.dataset.read_f64()?;
-        let cf = self.cf_attributes()?;
-        Ok(cf::apply_scale_offset(&raw, &cf))
+        let cf = cf::extract_cf_attributes(&self.attrs);
+        self.padded(cf::apply_scale_offset(&raw, &cf), || Ok(f64::NAN))
     }
 
     /// Read raw data as f64 without any scale/offset transformation.
+    ///
+    /// Unwritten records along an unlimited dimension read as the fill
+    /// value (`_FillValue`, else netCDF's default for the type), as in the
+    /// other `read_raw_*` methods and [`read_string`](Self::read_string).
     pub fn read_raw_f64(&self) -> Result<Vec<f64>, Error> {
-        Ok(self.dataset.read_f64()?)
+        self.padded_read(self.dataset.read_f64()?)
     }
 
     /// Read raw data as f32 without any scale/offset transformation.
     pub fn read_raw_f32(&self) -> Result<Vec<f32>, Error> {
-        Ok(self.dataset.read_f32()?)
+        self.padded_read(self.dataset.read_f32()?)
     }
 
     /// Read raw data as i32 without any scale/offset transformation.
     pub fn read_raw_i32(&self) -> Result<Vec<i32>, Error> {
-        Ok(self.dataset.read_i32()?)
+        self.padded_read(self.dataset.read_i32()?)
     }
 
     /// Read raw data as i64 without any scale/offset transformation.
     pub fn read_raw_i64(&self) -> Result<Vec<i64>, Error> {
-        Ok(self.dataset.read_i64()?)
+        self.padded_read(self.dataset.read_i64()?)
     }
 
     /// Read raw data as u64 without any scale/offset transformation.
     pub fn read_raw_u64(&self) -> Result<Vec<u64>, Error> {
-        Ok(self.dataset.read_u64()?)
+        self.padded_read(self.dataset.read_u64()?)
     }
 
     /// Read raw data as strings.
     pub fn read_string(&self) -> Result<Vec<String>, Error> {
-        Ok(self.dataset.read_string()?)
+        self.padded_read(self.dataset.read_string()?)
+    }
+
+    /// The fill value netCDF-C gives the variable's unwritten values: its
+    /// `_FillValue`, else the default fill value of its type (`NC_FILL_*`).
+    fn fill_value(&self) -> Result<FillValue, Error> {
+        if let Some(fill) = cf::extract_cf_attributes(&self.attrs).fill_value {
+            return Ok(fill);
+        }
+        Ok(default_fill(self.nc_type()?))
+    }
+
+    /// `data`, read in the dataset's extent, laid out in the variable's
+    /// netCDF shape with the fill value in the positions not written.
+    fn padded_read<T: Clone + FromFill>(&self, data: Vec<T>) -> Result<Vec<T>, Error> {
+        self.padded(data, || Ok(T::from_fill(&self.fill_value()?)))
+    }
+
+    /// Like [`padded_read`](Self::padded_read), padding with what `fill`
+    /// returns (called only when there is something to pad).
+    fn padded<T: Clone>(
+        &self,
+        data: Vec<T>,
+        fill: impl FnOnce() -> Result<T, Error>,
+    ) -> Result<Vec<T>, Error> {
+        let extent = self.dataset.shape()?;
+        let shape = nc_shape(&extent, &self.dims);
+        if shape == extent {
+            return Ok(data);
+        }
+        pad(data, &extent, &shape, fill()?)
     }
 
     /// Read raw bytes without any type conversion.
@@ -122,23 +176,23 @@ impl<'f> Variable<'f> {
         let dtype = self.dataset.dtype()?;
         match dtype {
             clawhdf5::DType::F64 => {
-                let vals = self.dataset.read_f64()?;
+                let vals = self.read_raw_f64()?;
                 Ok(vals.iter().flat_map(|v| v.to_le_bytes()).collect())
             }
             clawhdf5::DType::F32 => {
-                let vals = self.dataset.read_f32()?;
+                let vals = self.read_raw_f32()?;
                 Ok(vals.iter().flat_map(|v| v.to_le_bytes()).collect())
             }
             clawhdf5::DType::I32 => {
-                let vals = self.dataset.read_i32()?;
+                let vals = self.read_raw_i32()?;
                 Ok(vals.iter().flat_map(|v| v.to_le_bytes()).collect())
             }
             clawhdf5::DType::I64 => {
-                let vals = self.dataset.read_i64()?;
+                let vals = self.read_raw_i64()?;
                 Ok(vals.iter().flat_map(|v| v.to_le_bytes()).collect())
             }
             clawhdf5::DType::U64 => {
-                let vals = self.dataset.read_u64()?;
+                let vals = self.read_raw_u64()?;
                 Ok(vals.iter().flat_map(|v| v.to_le_bytes()).collect())
             }
             other => {
@@ -169,64 +223,137 @@ impl std::fmt::Debug for Variable<'_> {
     }
 }
 
-/// Build variables from a group's datasets and associated dimensions.
-pub(crate) fn build_variables<'f>(
-    group: &clawhdf5::Group<'f>,
-    available_dims: &[Dimension],
-) -> Result<Vec<Variable<'f>>, Error> {
-    let dataset_names = group.datasets()?;
-    let mut variables = Vec::new();
-
-    for ds_name in &dataset_names {
-        let ds = group.dataset(ds_name)?;
-        let shape = ds.shape()?;
-
-        // Associate dimensions with this variable.
-        // First try DIMENSION_LIST attribute, then fall back to shape matching.
-        let var_dims = match_dimensions_to_variable(&shape, available_dims);
-
-        variables.push(Variable::new(ds_name.clone(), ds, var_dims));
+/// The netCDF shape of a variable whose dataset has `extent`: along an
+/// unlimited dimension the dimension's length, which is at least the extent.
+fn nc_shape(extent: &[u64], dims: &[Dimension]) -> Vec<u64> {
+    if dims.len() != extent.len() {
+        return extent.to_vec();
     }
-
-    Ok(variables)
+    extent
+        .iter()
+        .zip(dims)
+        .map(|(&e, d)| if d.is_unlimited { e.max(d.size) } else { e })
+        .collect()
 }
 
-/// Match dimensions to a variable based on shape.
-///
-/// For each axis of the variable, find a dimension with matching size.
-/// If multiple dimensions have the same size, prefer exact name matching
-/// from the convention order.
-pub(crate) fn match_dimensions_to_variable(
-    shape: &[u64],
-    available_dims: &[Dimension],
-) -> Vec<Dimension> {
-    let mut result = Vec::with_capacity(shape.len());
-
-    // Track which dimensions have been used to avoid duplicates
-    let mut used = vec![false; available_dims.len()];
-
-    for &dim_size in shape {
-        let mut matched = false;
-
-        // Find a dimension with matching size that hasn't been used yet
-        for (i, dim) in available_dims.iter().enumerate() {
-            if !used[i] && dim.size == dim_size {
-                result.push(dim.clone());
-                used[i] = true;
-                matched = true;
+/// `data`, row-major in `extent`, placed in a row-major array of `shape`
+/// (as many axes, each at least as long) filled with `fill`.
+fn pad<T: Clone>(data: Vec<T>, extent: &[u64], shape: &[u64], fill: T) -> Result<Vec<T>, Error> {
+    let too_big = || Error::TypeError(format!("variable of shape {shape:?} is too large"));
+    let to_usize = |dims: &[u64]| -> Result<Vec<usize>, Error> {
+        dims.iter()
+            .map(|&d| usize::try_from(d).map_err(|_| too_big()))
+            .collect()
+    };
+    let (extent, shape) = (to_usize(extent)?, to_usize(shape)?);
+    let total = shape
+        .iter()
+        .try_fold(1usize, |n, &d| n.checked_mul(d))
+        .ok_or_else(too_big)?;
+    if extent.len() != shape.len()
+        || extent.iter().zip(&shape).any(|(e, s)| e > s)
+        || extent.iter().product::<usize>() != data.len()
+    {
+        return Err(Error::TypeError(format!(
+            "{} values of extent {extent:?} do not fit shape {shape:?}",
+            data.len()
+        )));
+    }
+    let (Some((&row, outer)), Some(&row_stride)) = (extent.split_last(), shape.last()) else {
+        return Ok(data);
+    };
+    let mut out = vec![fill; total];
+    if row == 0 {
+        return Ok(out);
+    }
+    // The position of the current row along each outer axis.
+    let mut index = vec![0usize; outer.len()];
+    for chunk in data.chunks_exact(row) {
+        let offset = index.iter().zip(&shape).fold(0, |o, (&i, &n)| o * n + i);
+        out[offset * row_stride..][..row].clone_from_slice(chunk);
+        for (i, &n) in index.iter_mut().zip(outer).rev() {
+            *i += 1;
+            if *i < n {
                 break;
             }
-        }
-
-        if !matched {
-            // Create an anonymous dimension for unmatched sizes
-            result.push(Dimension {
-                name: format!("dim_{dim_size}"),
-                size: dim_size,
-                is_unlimited: false,
-            });
+            *i = 0;
         }
     }
+    Ok(out)
+}
 
-    result
+/// netCDF's default fill value for a type (`NC_FILL_*` in `netcdf.h`).
+fn default_fill(nc_type: NcType) -> FillValue {
+    match nc_type {
+        NcType::Byte => FillValue::Int(-127),
+        NcType::UByte => FillValue::UInt(255),
+        NcType::Short => FillValue::Int(-32767),
+        NcType::UShort => FillValue::UInt(65535),
+        NcType::Int => FillValue::Int(-2_147_483_647),
+        NcType::UInt => FillValue::UInt(4_294_967_295),
+        NcType::Int64 => FillValue::Int(-9_223_372_036_854_775_806),
+        NcType::UInt64 => FillValue::UInt(18_446_744_073_709_551_614),
+        NcType::Float => FillValue::Float(f64::from(9.969_21e36_f32)),
+        NcType::Double => FillValue::Float(9.969_209_968_386_869e36),
+        NcType::String => FillValue::String(String::new()),
+        NcType::Char => FillValue::Int(0),
+    }
+}
+
+/// A fill value converted to the element type of a read, as the read
+/// converts the stored values.
+trait FromFill {
+    fn from_fill(fill: &FillValue) -> Self;
+}
+
+macro_rules! numeric_from_fill {
+    ($($t:ty),*) => {$(
+        impl FromFill for $t {
+            fn from_fill(fill: &FillValue) -> Self {
+                match fill {
+                    FillValue::Float(v) => *v as $t,
+                    FillValue::Int(v) => *v as $t,
+                    FillValue::UInt(v) => *v as $t,
+                    FillValue::String(_) => <$t>::default(),
+                }
+            }
+        }
+    )*};
+}
+numeric_from_fill!(f64, f32, i32, i64, u64);
+
+impl FromFill for String {
+    fn from_fill(fill: &FillValue) -> Self {
+        match fill {
+            FillValue::String(s) => s.clone(),
+            _ => String::new(),
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn pad_places_rows() {
+        // (2, 1) written of (2, 4): each row padded, not the tail.
+        let out = pad(vec![1, 2], &[2, 1], &[2, 4], 0).unwrap();
+        assert_eq!(out, vec![1, 0, 0, 0, 2, 0, 0, 0]);
+        // Leading axis short.
+        let out = pad(vec![1, 2, 3, 4], &[2, 2], &[3, 2], -1).unwrap();
+        assert_eq!(out, vec![1, 2, 3, 4, -1, -1]);
+        // Nothing written.
+        let out = pad(Vec::<i32>::new(), &[0, 3], &[2, 3], 7).unwrap();
+        assert_eq!(out, vec![7; 6]);
+        // 3-D, middle axis short.
+        let out = pad(vec![1, 2, 3, 4], &[2, 1, 2], &[2, 2, 2], 0).unwrap();
+        assert_eq!(out, vec![1, 2, 0, 0, 3, 4, 0, 0]);
+    }
+
+    #[test]
+    fn pad_rejects_wrong_length() {
+        assert!(pad(vec![1, 2, 3], &[2, 2], &[3, 2], 0).is_err());
+        assert!(pad(vec![1, 2, 3, 4], &[2, 2], &[1, 4], 0).is_err());
+    }
 }