From 86bddaa0e4b714a1ad457c82799a10fb0645c118 Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Thu, 9 Apr 2026 13:14:21 +0200 Subject: [PATCH 1/5] partial review --- dascore/data_registry.txt | 1 + dascore/io/core.py | 20 +- dascore/io/netcdf/__init__.py | 5 + dascore/io/netcdf/core.py | 239 ++++ dascore/io/netcdf/utils.py | 627 ++++++++++ dascore/utils/io.py | 7 +- pyproject.toml | 2 + tests/test_io/test_common_io.py | 10 +- tests/test_io/test_io_core.py | 38 + tests/test_io/test_netcdf/test_netcdf.py | 1440 ++++++++++++++++++++++ 10 files changed, 2378 insertions(+), 11 deletions(-) create mode 100644 dascore/io/netcdf/__init__.py create mode 100644 dascore/io/netcdf/core.py create mode 100644 dascore/io/netcdf/utils.py create mode 100644 tests/test_io/test_netcdf/test_netcdf.py diff --git a/dascore/data_registry.txt b/dascore/data_registry.txt index 879d1509d..caf9c2492 100644 --- a/dascore/data_registry.txt +++ b/dascore/data_registry.txt @@ -39,3 +39,4 @@ das_vader_1.jld2 cbeace6bfcc17711796886ba072ae68bb792903c38a052d592ed968885933d8 febus_2.h5 c118960a94e37fbff0eb5c33856d34cdfe81609902c4feaedab9949498d31c23 https://github.com/dasdae/test_data/raw/master/das/febus_2.h5 febg1_C1_2023-05-10T12.25.03+0000.bsl e1a8ff72f3ec1805129267df916f41419bf7fa3a4993602e2b85e721cad922ae https://github.com/dasdae/test_data/raw/master/dss/febg1_C1_2023-05-10T12.25.03+0000.bsl febg1_C1_2023-05-10T12.27.33+0000.bsl 233df0c184796944442ae19beddcf962aba1c9fab337fd1c74971f0c4d513a36 https://github.com/dasdae/test_data/raw/master/dss/febg1_C1_2023-05-10T12.27.33+0000.bsl +xdas_netcdf.nc 9e53fa1ce8395fedbb195048b3eb2832b87cb6883867cdff30be93078e8027f7 https://github.com/dasdae/test_data/raw/master/das/xdas_netcdf.nc diff --git a/dascore/io/core.py b/dascore/io/core.py index 0bca95863..621366700 100644 --- a/dascore/io/core.py +++ b/dascore/io/core.py @@ -376,6 +376,8 @@ def _get_prioritized_list(self, input_type="file"): second_class_fiber_ios = [] for format_name in self.known_formats: unsorted = self._format_version[format_name] + if not unsorted: + continue keys = sorted(unsorted, reverse=True) fiber_ios = [unsorted[key] for key in keys] priority_fiber_ios.append(fiber_ios[0]) @@ -399,10 +401,20 @@ def load_plugins(self, format: str | None = None): # Load one, or all, formats for form in formats: entries = [name for name in self._eps.index if name.startswith(form)] - for eps in self._eps.loc[entries]: - self.register_fiberio(eps()()) - # The selected format(s) should now be loaded - assert set(formats).isdisjoint(self.unloaded_formats) + for name, loader in self._eps.loc[entries].items(): + fiberio = self._load_entry_point(name, loader) + if fiberio is not None: + self.register_fiberio(fiberio) + return + + def _load_entry_point(self, name: str, loader) -> FiberIO | None: + """Load one FiberIO entry point, skipping broken registrations.""" + try: + return loader()() + except (ImportError, MissingOptionalDependencyError) as exc: + msg = f"Failed to load FiberIO plugin {name!r}: {exc}" + warnings.warn(msg, UserWarning, stacklevel=2) + return None def register_fiberio(self, fiberio: FiberIO): """Register a new fiber IO to manage.""" diff --git a/dascore/io/netcdf/__init__.py b/dascore/io/netcdf/__init__.py new file mode 100644 index 000000000..f5c3a7a48 --- /dev/null +++ b/dascore/io/netcdf/__init__.py @@ -0,0 +1,5 @@ +"""NetCDF IO support for DASCore using CF (Climate and Forecast) conventions.""" + +from __future__ import annotations + +from dascore.io.netcdf.core import NetCDFCFV18 diff --git a/dascore/io/netcdf/core.py b/dascore/io/netcdf/core.py new file mode 100644 index 000000000..7da3f4e42 --- /dev/null +++ b/dascore/io/netcdf/core.py @@ -0,0 +1,239 @@ +"""Core NetCDF IO implementation with CF conventions.""" + +from __future__ import annotations + +from pathlib import Path + +import h5py + +import dascore as dc +from dascore.constants import SpoolType +from dascore.io import FiberIO +from dascore.io.core import ScanPayload, _make_scan_payload +from dascore.utils.hdf5 import H5Reader +from dascore.utils.io import patch_to_xarray +from dascore.utils.misc import optional_import + +from .utils import ( + extract_patch_attrs_from_netcdf, + find_main_data_variable, + get_xarray_data_var_name, + get_xarray_engine, + get_cf_data_attrs, + get_cf_global_attrs, + get_cf_version, + is_netcdf4_file, + iter_written_aux_coords, + parse_cf_version, + read_netcdf_coordinates, + XDAS_PAYLOAD_VARIABLE, + coord_attrs, +) + + +class NetCDFCFV18(FiberIO): + """NetCDF-4 IO with CF-1.8 conventions, using xarray/netcdf4 for IO.""" + + name = "NETCDF_CF" + version = "1.8" + preferred_extensions = ("nc", "nc4", "netcdf") + + def get_format(self, resource: H5Reader, **kwargs) -> tuple[str, str] | bool: + """Return format tuple if file is a CF-convention NetCDF-4, else False.""" + if not is_netcdf4_file(resource): + return False + cf_version = get_cf_version(resource) + if not cf_version: + return False + try: + if parse_cf_version(cf_version) >= (1, 6): + return self.name, self.version + except (TypeError, ValueError): + pass + return False + + def read(self, resource: Path, **kwargs) -> SpoolType: + """Read a NetCDF-4 file into a Spool.""" + xr = optional_import("xarray") + engine = get_xarray_engine() + # Read structural metadata first so xarray only has to resolve the final + # payload variable and dense coordinate values. + data_var_name, attrs_dict, coords = self._read_metadata(resource) + data_var_name, data_array, data = self._read_data_array( + resource, xr, engine, data_var_name + ) + patch = self._build_patch( + data=data, + data_array=data_array, + coords=coords, + attrs_dict=attrs_dict, + data_var_name=data_var_name, + ) + patch = self._apply_coordinate_filtering(patch, kwargs) + if not patch.data.size: + return dc.spool([]) + return dc.spool([patch]) + + def _read_metadata(self, resource: Path): + """Read NetCDF metadata needed before opening through xarray.""" + with h5py.File(resource, "r") as h5file: + data_var_name = self._get_data_variable_name(h5file) + attrs_dict = self._get_patch_attrs(h5file) + coords = read_netcdf_coordinates(h5file, data_var_name) + return data_var_name, attrs_dict, coords + + def _read_data_array(self, resource: Path, xr, engine: str, data_var_name: str): + """Load the selected xarray data variable from disk.""" + # TODO: consider a lazy read path for large NetCDF arrays. + with xr.open_dataset(resource, engine=engine) as dataset: + data_array = dataset.get(data_var_name) + if data_array is None: + data_var_name = get_xarray_data_var_name(dataset) + data_array = dataset[data_var_name] + data = data_array.load().data + return data_var_name, data_array, data + + def _build_patch(self, *, data, data_array, coords, attrs_dict, data_var_name: str): + """Merge coordinate metadata and construct the output patch.""" + source_patch_id = ( + XDAS_PAYLOAD_VARIABLE if data_var_name is None else data_var_name + ) + # Start from the HDF-derived coordinates, then add any extra xarray-only + # coordinates that were materialized during decode. + coords_dict = { + name: ( + coord.values + if name in coords.dims + else (coords.dim_map[name], coord.values) + ) + for name, coord in coords.coord_map.items() + } + for name, coord in data_array.coords.items(): + if name not in coords_dict: + coords_dict[name] = (coord.dims, coord.values) + coords = dc.get_coord_manager(coords=coords_dict, dims=coords.dims) + return dc.Patch( + data=data, + coords=coords, + dims=coords.dims, + attrs=attrs_dict | {"_source_patch_id": source_patch_id}, + ) + + def _build_data_array(self, patch: dc.Patch): + """Convert a patch to an xarray DataArray with coordinate attrs.""" + data_array = patch_to_xarray(patch).rename("data") + for name, coord in patch.coords.coord_map.items(): + if coord._partial: + continue + data_array.coords[name].attrs.update(coord_attrs(name, coord)) + return data_array + + def _build_dataset(self, patch: dc.Patch, data_array): + """Create the xarray Dataset and attach CF metadata.""" + global_attrs = get_cf_global_attrs(patch.attrs, self.version) + if patch.attrs.data_type: + global_attrs["source_data_type"] = patch.attrs.data_type + dataset = data_array.to_dataset() + # NetCDF attrs cannot safely preserve DASCore's null-ish metadata, so + # filter those out before handing the dataset to xarray. + dataset.attrs.update( + { + key: value + for key, value in global_attrs.items() + if value not in (None, "") + } + ) + dataset["data"].attrs.update( + get_cf_data_attrs(patch.attrs.data_type or "acoustic_signal") + ) + aux_coord_names = tuple(iter_written_aux_coords(patch)) + if aux_coord_names: + dataset["data"].attrs["coordinates"] = " ".join(aux_coord_names) + return dataset + + def _get_write_encoding(self, **kwargs): + """Translate explicit write options into xarray encoding hints.""" + compression = kwargs.get("compression") + if compression not in ("gzip", None, False): + msg = "xarray netcdf4 writing currently supports only gzip compression." + raise ValueError(msg) + chunks = kwargs.get("chunks") + encoding: dict[str, object] = {} + if chunks not in (None, False, True): + encoding["chunksizes"] = tuple(chunks) + if compression == "gzip": + encoding["zlib"] = True + encoding["complevel"] = kwargs.get("compression_opts", 4) + encoding["shuffle"] = True + return encoding + + def write(self, spool: SpoolType, resource: Path, **kwargs) -> None: + """ + Write a Spool to NetCDF-4 with CF-1.8 conventions. + + Parameters + ---------- + kwargs + compression: 'gzip', None, or False + compression_opts: gzip level 1-9 (default 4) + chunks: True to defer chunking to xarray/backend defaults, or an + explicit tuple of chunk sizes + """ + patch = self._validate_and_extract_patch(spool) + optional_import("xarray") + engine = get_xarray_engine() + # Build the xarray object first, then attach CF metadata and optional + # storage hints in separate steps to keep write policy isolated. + data_array = self._build_data_array(patch) + dataset = self._build_dataset(patch, data_array) + encoding = self._get_write_encoding(**kwargs) + dataset.to_netcdf( + resource, + engine=engine, + encoding={"data": encoding} if encoding else None, + ) + + def scan(self, resource: H5Reader, **kwargs) -> list[ScanPayload]: + """Scan NetCDF file to extract metadata without loading data.""" + data_var_name = self._get_data_variable_name(resource) + data_var = resource[data_var_name] + coords = read_netcdf_coordinates(resource, data_var_name) + attrs_dict = self._get_patch_attrs(resource) + return [ + _make_scan_payload( + attrs=attrs_dict, + coords=coords, + dims=coords.dims, + shape=data_var.shape, + dtype=str(data_var.dtype), + source_patch_id=data_var_name, + ) + ] + + def _get_data_variable_name(self, resource: H5Reader) -> str: + """Return the main NetCDF data variable name or raise.""" + data_var_name = find_main_data_variable(resource) + if data_var_name is None: + msg = "No suitable data variable found in NetCDF file" + raise ValueError(msg) + return data_var_name + + def _get_patch_attrs(self, resource: H5Reader) -> dict: + """Extract patch attrs from a NetCDF resource.""" + return extract_patch_attrs_from_netcdf(resource) + + def _apply_coordinate_filtering(self, patch: dc.Patch, kwargs: dict) -> dc.Patch: + """Apply coordinate selection kwargs to a loaded patch.""" + coord_kwargs = {k: v for k, v in kwargs.items() if k in patch.coords.coord_map} + return patch.select(**coord_kwargs) if coord_kwargs else patch + + def _validate_and_extract_patch(self, spool: SpoolType) -> dc.Patch: + """Validate write input and return the single supported patch.""" + patches = [spool] if isinstance(spool, dc.Patch) else list(spool) + if len(patches) == 0: + msg = "Cannot write empty spool" + raise ValueError(msg) + if len(patches) > 1: + msg = "Multi-patch spools not yet supported for NetCDF output" + raise NotImplementedError(msg) + return patches[0] diff --git a/dascore/io/netcdf/utils.py b/dascore/io/netcdf/utils.py new file mode 100644 index 000000000..6e1d5680d --- /dev/null +++ b/dascore/io/netcdf/utils.py @@ -0,0 +1,627 @@ +"""Utilities for NetCDF IO with CF conventions support. + +See https://cfconventions.org/ for the Climate and Forecast metadata standard. +""" + +from __future__ import annotations + +import datetime +from functools import cache +from typing import TYPE_CHECKING + +import h5py +import numpy as np + +import dascore as dc +from dascore.exceptions import MissingOptionalDependencyError +from dascore.utils.misc import optional_import +from dascore.utils.time import to_datetime64, to_float + +if TYPE_CHECKING: + from dascore.core.attrs import PatchAttrs + from dascore.core.coordmanager import CoordManager + +# CF-compliant time reference (Unix epoch is standard) +CF_TIME_REFERENCE = "seconds since 1970-01-01 00:00:00" +CF_CALENDAR = "proleptic_gregorian" +XDAS_PAYLOAD_VARIABLE = "__values__" + +# CF standard names for DAS data types +CF_STANDARD_NAMES = { + "strain": "strain", + "strain_rate": "strain_rate", + "velocity": "velocity", + "acceleration": "acceleration", + "temperature": "air_temperature", + "pressure": "air_pressure", + "acoustic": "acoustic_signal", +} + + +def coord_attrs(name: str, coord) -> dict[str, str]: + """Return CF-ish coordinate attrs for xarray-backed NetCDF output.""" + # Keep time explicit because CF consumers treat it as a special semantic + # axis, not just another coordinate with datetime-like values. + if name == "time": + return { + "standard_name": "time", + "long_name": "Time", + "axis": "T", + } + return { + "long_name": name.replace("_", " ").title(), + "standard_name": name.lower(), + "units": str(coord.units.units) if coord.units else ("m" if "depth" in name else "1"), + } + + +def get_xarray_data_var_name(dataset) -> str: + """Return the main xarray data variable name.""" + if "data" in dataset.data_vars: + return "data" + # XDAS-style files can surface the primary payload under a None key while + # exposing coord helper arrays (for example *_indices/*_values) as data vars. + if None in dataset.data_vars: + return None + if len(dataset.data_vars) == 1: + return next(iter(dataset.data_vars)) + msg = "No suitable data variable found in NetCDF file" + raise ValueError(msg) + + +def parse_cf_version(cf_version: str) -> tuple[int, int]: + """Parse a CF version string into comparable major/minor integers.""" + parts = cf_version.split(".") + major = int(parts[0]) + minor = int(parts[1]) if len(parts) > 1 else 0 + return major, minor + + +@cache +def get_xarray_engine(on_missing: str = "raise") -> str | None: + """Return the preferred xarray engine for NetCDF-4 files.""" + # Map importable module names onto the engine strings xarray expects. + module_to_engine = { + "netCDF4": "netcdf4", + "h5netcdf": "h5netcdf", + } + for module_name, engine_name in module_to_engine.items(): + mod = optional_import(module_name, on_missing="ignore") + if mod is not None: + return engine_name + # If no backend available ignore or raise. + if on_missing == "ignore": + return None + msg = ( + "Either netCDF4 or h5netcdf is required for NetCDF-4 read/write " + "functionality." + ) + raise MissingOptionalDependencyError(msg) + + +def iter_written_aux_coords(patch: dc.Patch): + """Yield names of auxiliary coordinates serialized to NetCDF.""" + for name, coord in patch.coords.coord_map.items(): + if coord._partial or name in patch.dims: + continue + yield name + + +def is_netcdf4_file(h5file: h5py.File) -> bool: + """ + Check if an HDF5 file follows NetCDF-4 conventions. + + Parameters + ---------- + h5file + Open h5py.File object + + Returns + ------- + bool + True if file appears to be NetCDF-4 format + + Notes + ----- + NetCDF-4 files are identified by: + - _NCProperties attribute (NetCDF-4 specific) + - Conventions attribute starting with "CF-" + """ + try: + # Check for NetCDF-specific markers + if "_NCProperties" in h5file.attrs: + return True + + # Check for Conventions attribute indicating CF compliance + conventions = h5file.attrs.get("Conventions", "") + if isinstance(conventions, bytes): + conventions = conventions.decode("utf-8", errors="ignore") + if conventions and "CF" in conventions: + return True + + return False + + except (AttributeError, KeyError): + return False + + +def get_cf_version(h5file: h5py.File) -> str | None: + """ + Extract CF convention version from NetCDF file. + + Parameters + ---------- + h5file + Open h5py.File object + + Returns + ------- + str | None + CF version string (e.g., "1.8") or None if not found + """ + conventions = h5file.attrs.get("Conventions", "") + if isinstance(conventions, bytes): + conventions = conventions.decode("utf-8", errors="ignore") + + # Handle various CF convention formats + if "CF-" in conventions: + # Format: "CF-1.8" or "CF-1.8, ACDD-1.3" + parts = conventions.split("CF-", 1)[1].split()[0].rstrip(",;") + return parts + elif conventions.startswith("CF "): + # Format: "CF 1.8" + return conventions.split()[1].rstrip(",;") + + return None + + +def datetime64_to_cf_time(dt_array: np.ndarray) -> np.ndarray: + """Convert datetime64 array to CF time (seconds since 1970-01-01).""" + return to_float(dt_array) + + +def cf_time_to_datetime64(time_array: np.ndarray, units: str) -> np.ndarray: + """Convert CF time array with units string to datetime64[ns]. + + Parameters + ---------- + time_array + Numeric time values. + units + CF units string, e.g. "seconds since 1970-01-01 00:00:00". + """ + if " since " not in units: + msg = f"Invalid CF time units format: {units}" + raise ValueError(msg) + time_unit, ref_str = units.split(" since ", 1) + time_unit = time_unit.strip().lower() + unit_to_seconds = { + "days": 86400.0, + "hours": 3600.0, + "minutes": 60.0, + "seconds": 1.0, + "milliseconds": 1e-3, + "microseconds": 1e-6, + } + if time_unit not in unit_to_seconds: + msg = f"Unsupported time unit: {time_unit}" + raise ValueError(msg) + # Convert to seconds since 1970 then to datetime64 + ref_epoch = to_float(to_datetime64(ref_str.strip())) + seconds_from_epoch = ( + np.asarray(time_array, dtype=np.float64) * unit_to_seconds[time_unit] + + ref_epoch + ) + return to_datetime64(seconds_from_epoch) + + +def get_cf_data_attrs(data_type: str = "acoustic_signal") -> dict[str, str | float]: + """ + Get CF-compliant attributes for DAS data variable. + + Parameters + ---------- + data_type + Type of DAS data + + Returns + ------- + dict + CF-compliant data attributes + """ + attrs = { + "long_name": "Distributed Acoustic Sensing data", + "_FillValue": np.nan, + } + + # Map data type to CF standard name + data_type_lower = data_type.lower() if data_type else "acoustic" + + # Check for exact match first + if data_type_lower in CF_STANDARD_NAMES: + attrs["standard_name"] = CF_STANDARD_NAMES[data_type_lower] + else: + # Check for partial matches + for key, std_name in CF_STANDARD_NAMES.items(): + if key in data_type_lower: + attrs["standard_name"] = std_name + break + else: + # Default to generic acoustic signal + attrs["standard_name"] = "acoustic_signal" + + # Add units based on data type + if "strain_rate" in data_type_lower: + attrs["units"] = "1/s" + attrs["long_name"] = "Strain rate" + elif "strain" in data_type_lower: + attrs["units"] = "1" + attrs["long_name"] = "Strain" + elif "velocity" in data_type_lower: + attrs["units"] = "m/s" + attrs["long_name"] = "Velocity" + elif "temperature" in data_type_lower: + attrs["units"] = "K" + attrs["long_name"] = "Temperature" + elif "pressure" in data_type_lower: + attrs["units"] = "Pa" + attrs["long_name"] = "Pressure" + else: + # Generic units for acoustic data + attrs["units"] = "1" + + return attrs + + +def get_cf_global_attrs( + patch_attrs: PatchAttrs, cf_version: str = "1.8" +) -> dict[str, str]: + """ + Get CF-compliant global attributes from PatchAttrs. + + Parameters + ---------- + patch_attrs + DASCore patch attributes + cf_version + CF convention version + + Returns + ------- + dict + CF-compliant global attributes + """ + now = datetime.datetime.now(datetime.timezone.utc) + attrs = { + "Conventions": f"CF-{cf_version}", + "title": "DAS data from DASCore", + "source": f"DASCore v{dc.__version__}", + "history": f"{now.isoformat()}: Created by DASCore", + "references": "https://dascore.org", + "comment": "Distributed Acoustic Sensing data", + "date_created": now.isoformat(), + } + + # Add optional attributes from patch + if patch_attrs.station: + attrs["station"] = patch_attrs.station + if patch_attrs.network: + attrs["network"] = patch_attrs.network + if patch_attrs.instrument_id: + attrs["instrument"] = patch_attrs.instrument_id + if patch_attrs.acquisition_id: + attrs["acquisition"] = patch_attrs.acquisition_id + if patch_attrs.tag: + attrs["tag"] = patch_attrs.tag + + # Add data category/type info + if patch_attrs.data_category: + attrs["data_category"] = patch_attrs.data_category + if hasattr(patch_attrs, "category") and patch_attrs.category: + attrs["category"] = patch_attrs.category + if patch_attrs.data_type: + attrs["data_type"] = patch_attrs.data_type + + # Add processing history if available + if hasattr(patch_attrs, "history") and patch_attrs.history: + history_str = " | ".join(str(h) for h in patch_attrs.history) + attrs["processing_history"] = history_str + + return attrs + + +def extract_patch_attrs_from_netcdf(h5file: h5py.File) -> dict: + """ + Extract patch attributes from NetCDF global attributes. + + Parameters + ---------- + h5file + Open h5py.File object + + Returns + ------- + dict + Dictionary of patch attributes + """ + attrs = {} + + # Map CF global attributes to patch attributes + attr_mapping = { + "station": "station", + "network": "network", + "instrument": "instrument_id", + "acquisition": "acquisition_id", + "tag": "tag", + "data_category": "data_category", + "category": "category", + } + + for cf_name, patch_name in attr_mapping.items(): + if cf_name in h5file.attrs: + value = h5file.attrs[cf_name] + if isinstance(value, bytes): + value = value.decode("utf-8", errors="ignore") + if value: # Only add non-empty values + attrs[patch_name] = value + for cf_name in ("data_type", "source_data_type"): + if cf_name not in h5file.attrs: + continue + value = h5file.attrs[cf_name] + if isinstance(value, bytes): + value = value.decode("utf-8", errors="ignore") + if value: + attrs.setdefault("data_type", value) + + return attrs + + +def _handle_time_interpolation( + h5file: h5py.File, coord_name: str, coord_data: np.ndarray +) -> np.ndarray | None: + """ + Handle coordinate interpolation for time coordinates. + + Parameters + ---------- + h5file + Open h5py.File object + coord_name + Name of the coordinate (e.g., "time") + coord_data + Raw coordinate data array + + Returns + ------- + np.ndarray | None + Converted datetime64 array if interpolation data found, None otherwise + """ + # Compatibility path for xdas-produced NetCDF files that store time as CF + # tie points in `_values` plus optional `_indices`. This is + # not generic CF behavior and does not require an xdas dependency. + # Look for time_values with proper CF units + time_values_name = f"{coord_name}_values" + if time_values_name not in h5file: + return None + + time_values_var = h5file[time_values_name] + units = time_values_var.attrs.get("units", "") + if isinstance(units, bytes): + units = units.decode("utf-8", errors="ignore") + + if "since" not in units: + return None + + # Get the time values and interpolate to full coordinate array + time_values = time_values_var[:] + time_indices_name = f"{coord_name}_indices" + + if time_indices_name in h5file: + # Use indices for interpolation + time_indices = h5file[time_indices_name][:] + if len(time_values) >= 2 and len(time_indices) >= 2: + # Linear interpolation from tie points to full coordinate + interpolated_values = np.interp(coord_data, time_indices, time_values) + else: + # Fall back to using time_values directly + interpolated_values = time_values + else: + # Use time_values directly + interpolated_values = time_values + + # Convert to datetime64 + try: + return cf_time_to_datetime64(interpolated_values, units) + except (ValueError, KeyError): + return None + + +def read_netcdf_coordinates( + h5file: h5py.File, data_var_name: str | None = None +) -> CoordManager: + """Read coordinate information from a NetCDF file into a CoordManager.""" + coords = {} + coord_order = [] + + # Use provided data variable name, or discover it + main_data_var = data_var_name or find_main_data_variable(h5file) + expected_shape = None + if main_data_var and main_data_var in h5file: + data_var = h5file[main_data_var] + expected_shape = data_var.shape + + # Try to get dimension order from DIMENSION_LIST + if "DIMENSION_LIST" in data_var.attrs: + dim_list = data_var.attrs["DIMENSION_LIST"] + for ref_array in dim_list: + try: + ref = ref_array[0] + dim_scale = h5file[ref] + dim_name = dim_scale.name.strip("/") + coord_order.append(dim_name) + except (IndexError, KeyError, TypeError): + # If we can't resolve reference, we'll fall back to discovery + pass + + # Collect dimension-scale coordinates + dim_coords = {} + for name, dataset in h5file.items(): + if not isinstance(dataset, h5py.Dataset): + continue + if not (dataset.is_scale or "NAME" in dataset.attrs): + continue + # Skip auxiliaries that don't match any data dimension size + if expected_shape and dataset.shape[0] not in expected_shape: + continue + if any( + skip in name.lower() + for skip in ("_points", "_indices", "_values", "_interpolation") + ): + continue + data = dataset[:] + if name == "time" or dataset.attrs.get("axis") == "T": + units = dataset.attrs.get("units", "") + if isinstance(units, bytes): + units = units.decode("utf-8", errors="ignore") + if "since" in units: + try: + data = cf_time_to_datetime64(data, units) + dim_coords["time"] = data + continue + except (ValueError, KeyError): + pass + time_coord = _handle_time_interpolation(h5file, name, data) + if time_coord is not None: + dim_coords["time"] = time_coord + continue + dim_coords[name] = data + + # Build coords dict preserving dimension order from DIMENSION_LIST + coords = {} + if coord_order: + for dim_name in coord_order: + if dim_name in dim_coords: + coords[dim_name] = dim_coords[dim_name] + for name, data in dim_coords.items(): + if name not in coords: + coords[name] = data + else: + coords = dim_coords + + # Collect non-dimension coordinates stored with _DASCORE_DIMS attribute + if main_data_var and main_data_var in h5file: + coord_attr = h5file[main_data_var].attrs.get("coordinates", "") + if isinstance(coord_attr, bytes): + coord_attr = coord_attr.decode("utf-8", errors="ignore") + for name in coord_attr.split(): + if name in coords or name not in h5file: + continue + dataset = h5file[name] + if not isinstance(dataset, h5py.Dataset): + continue + dims_attr = dataset.attrs.get("_DASCORE_DIMS", "") + if isinstance(dims_attr, bytes): + dims_attr = dims_attr.decode("utf-8", errors="ignore") + if dims_attr: + associated_dims = tuple(dims_attr.split(",")) + coords[name] = (associated_dims, dataset[:]) + + return dc.core.coordmanager.get_coord_manager(coords) + + +def validate_cf_compliance(h5file: h5py.File) -> list[str]: + """ + Validate CF compliance and return list of issues. + + Parameters + ---------- + h5file + Open h5py.File object + + Returns + ------- + list[str] + List of CF compliance issues found + """ + issues = [] + + # Check for required global attributes + if "Conventions" not in h5file.attrs: + issues.append("Missing required 'Conventions' global attribute") + + # Check coordinate variables + for name, dataset in h5file.items(): + if isinstance(dataset, h5py.Dataset) and dataset.is_scale: + # Check for required coordinate attributes + if "units" not in dataset.attrs: + issues.append(f"Coordinate '{name}' missing 'units' attribute") + + # Check time coordinate + if name == "time" or dataset.attrs.get("standard_name") == "time": + units = dataset.attrs.get("units", "") + if not units or "since" not in str(units): + issues.append(f"Time coordinate '{name}' has invalid units") + + # Check data variables + for name, dataset in h5file.items(): + if isinstance(dataset, h5py.Dataset) and not dataset.is_scale: + # Check for required data variable attributes + if "units" not in dataset.attrs: + issues.append(f"Data variable '{name}' missing 'units' attribute") + if "long_name" not in dataset.attrs: + issues.append(f"Data variable '{name}' missing 'long_name' attribute") + + return issues + + +def find_main_data_variable(h5file: h5py.File) -> str | None: + """ + Find the main data variable in the NetCDF file. + + Looks for 2D+ datasets that are not dimension scales. + Prioritizes variables with standard names suggesting DAS data. + + Parameters + ---------- + h5file + Open h5py.File object + + Returns + ------- + str | None + Name of the main data variable, or None if not found + """ + priority_names = [ + "data", + "acoustic_data", + "das_data", + "strain", + "strain_rate", + "velocity", + "amplitude", + ] + + candidates = [] + + for name, item in h5file.items(): + if not _is_data_variable_candidate(item): + continue + + # Check for priority matches first + if name in priority_names or _has_priority_standard_name(item, priority_names): + return name + + candidates.append(name) + + return candidates[0] if candidates else None + + +def _is_data_variable_candidate(item) -> bool: + """Check if item is a candidate for main data variable.""" + return isinstance(item, h5py.Dataset) and not item.is_scale and item.ndim >= 2 + + +def _has_priority_standard_name(item, priority_names: list[str]) -> bool: + """Check if dataset has a standard_name matching priority names.""" + std_name = str(item.attrs.get("standard_name", "")).lower() + return any(name in std_name for name in priority_names) diff --git a/dascore/utils/io.py b/dascore/utils/io.py index db10cf8a0..8eb1abe31 100644 --- a/dascore/utils/io.py +++ b/dascore/utils/io.py @@ -255,7 +255,12 @@ def __del__(self): def patch_to_xarray(patch: PatchType): """Return a data array with patch contents.""" xr = optional_import("xarray") - attrs = dict(patch.attrs) + # Omit None-valued attrs because xarray backends may reject them during + # NetCDF serialization, while a missing attr round-trips cleanly. + attrs = { + key: value for key, value in dict(patch.attrs).items() + if value is not None + } patch_dims = patch.dims coords = {} for name, coord in patch.coords.coord_map.items(): diff --git a/pyproject.toml b/pyproject.toml index 006f4b7d5..400ad97d9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -64,6 +64,7 @@ dependencies = [ extras = [ "xarray", + "netCDF4", "findiff", "obspy", "numba", @@ -148,6 +149,7 @@ RSF__V1 = "dascore.io.rsf.core:RSFV1" WAV = "dascore.io.wav.core:WavIO" XMLBINARY__V1 = "dascore.io.xml_binary.core:XMLBinaryV1" GDR_DAS__V1 = "dascore.io.gdr.core:GDR_V1" +NETCDF_CF__V1_8 = "dascore.io.netcdf.core:NetCDFCFV18" [project.entry-points."dascore.patch_namespace"] io = "dascore.io:PatchIO" diff --git a/tests/test_io/test_common_io.py b/tests/test_io/test_common_io.py index 2dfcaef15..31ef65f92 100644 --- a/tests/test_io/test_common_io.py +++ b/tests/test_io/test_common_io.py @@ -9,7 +9,6 @@ """ from __future__ import annotations - from contextlib import suppress from functools import cache from io import BytesIO, UnsupportedOperation @@ -29,6 +28,7 @@ from dascore.io.febus import Febus1, Febus2 from dascore.io.gdr import GDR_V1 from dascore.io.h5simple import H5Simple +from dascore.io.netcdf import NetCDFCFV18 from dascore.io.neubrex import NeubrexDASV1, NeubrexRFSV1 from dascore.io.optodas import OptoDASV8 from dascore.io.pickle import PickleIO @@ -45,11 +45,8 @@ ) from dascore.utils.downloader import fetch, get_registry_df from dascore.utils.misc import all_close, iterate -from tests.test_io._common_io_test_utils import ( - get_flat_io_test, - skip_missing, - skip_timeout, -) +from tests.test_io._common_io_test_utils import get_flat_io_test, skip_missing, skip_timeout + # --- Fixtures @@ -87,6 +84,7 @@ ), Terra15FormatterV5(): ("terra15_v5_test_file.hdf5",), Terra15FormatterV6(): ("terra15_v6_test_file.hdf5",), + NetCDFCFV18(): ("xdas_netcdf.nc",), } # This tuple is for fiber io which support a write method and can write diff --git a/tests/test_io/test_io_core.py b/tests/test_io/test_io_core.py index 09a7f1b02..d68009e0b 100644 --- a/tests/test_io/test_io_core.py +++ b/tests/test_io/test_io_core.py @@ -372,6 +372,44 @@ def test_load_plugins_empty_entry_points(self, format_manager): format_manager._eps = pd.Series(dtype=object) format_manager.load_plugins() + def test_load_entry_point_warns_on_broken_plugin(self, format_manager): + """Broken plugin loaders should warn and return None.""" + + def loader(): + raise ImportError("boom") + + with pytest.warns(UserWarning, match="Failed to load FiberIO plugin 'broken'"): + out = format_manager._load_entry_point("broken", loader) + + assert out is None + + def test_load_entry_point_raises_on_runtime_plugin_error(self, format_manager): + """Runtime plugin construction failures should not be silently skipped.""" + + def loader(): + raise RuntimeError("boom") + + with pytest.raises(RuntimeError, match="boom"): + format_manager._load_entry_point("broken", loader) + + def test_prioritized_list_skips_unloaded_formats(self, format_manager): + """Formats with no registered versions should not break prioritization.""" + manager = type(format_manager)("dascore.fiber_io") + manager.__dict__.pop("known_formats", None) + manager._eps = pd.Series( + { + "BROKEN__V1": lambda: (_ for _ in ()).throw(ImportError("boom")), + "GOOD__V1": lambda: _ReadOnlySummaryFormatter, + } + ) + + prioritized = manager._get_prioritized_list() + + registered = manager._format_version[_ReadOnlySummaryFormatter.name.upper()] + assert "1" in registered + assert isinstance(registered["1"], _ReadOnlySummaryFormatter) + assert isinstance(prioritized, tuple) + class TestFormatter: """Tests for adding file supports through Formatter.""" diff --git a/tests/test_io/test_netcdf/test_netcdf.py b/tests/test_io/test_netcdf/test_netcdf.py new file mode 100644 index 000000000..b13fb369d --- /dev/null +++ b/tests/test_io/test_netcdf/test_netcdf.py @@ -0,0 +1,1440 @@ +"""Tests for NetCDF IO with CF conventions.""" + +from __future__ import annotations + +import importlib.util + +import h5py +import numpy as np +import pytest + +import dascore as dc +from dascore.exceptions import MissingOptionalDependencyError +from dascore.io.netcdf import core as netcdf_core +from dascore.io.netcdf import utils as netcdf_utils +from dascore.io.netcdf.utils import ( + cf_time_to_datetime64, + datetime64_to_cf_time, + extract_patch_attrs_from_netcdf, + find_main_data_variable, + get_cf_data_attrs, + get_cf_global_attrs, + get_cf_version, + is_netcdf4_file, + read_netcdf_coordinates, + validate_cf_compliance, +) +from dascore.utils.downloader import fetch + + +def _get_xarray_netcdf_engine() -> str | None: + """Return an xarray engine that can open NetCDF-4 files, if available.""" + if importlib.util.find_spec("netCDF4"): + return "netcdf4" + if importlib.util.find_spec("h5netcdf"): + return "h5netcdf" + return None + + +def _require_xarray_netcdf_engine() -> str: + """Return an xarray NetCDF backend or skip the test.""" + if not importlib.util.find_spec("xarray"): + pytest.skip("xarray not installed") + engine = _get_xarray_netcdf_engine() + if engine is None: + pytest.skip("xarray NetCDF-4 backend not installed") + return engine + + +def _assert_patch_round_trip_equal(expected: dc.Patch, observed: dc.Patch) -> None: + """Assert two patches are equivalent for round-trip testing.""" + assert expected.equals(observed) + assert expected.dims == observed.dims + assert set(expected.coords.coord_map) == set(observed.coords.coord_map) + expected_attrs = expected.attrs.model_dump() + observed_attrs = observed.attrs.model_dump() + expected_attrs.pop("_source_patch_id", None) + observed_attrs.pop("_source_patch_id", None) + assert expected_attrs == observed_attrs + + +def _assert_patch_compatible_with_xarray_output( + expected: dc.Patch, observed: dc.Patch +) -> None: + """Assert xarray-read DASCore output preserves patch content plus CF attrs.""" + assert expected.equals(observed) + assert expected.dims == observed.dims + assert set(expected.coords.coord_map) == set(observed.coords.coord_map) + observed_attrs = observed.attrs.model_dump() + for key, value in expected.attrs.model_dump().items(): + assert observed_attrs.get(key) == value + + +@pytest.fixture +def example_patch(): + """Create an example patch for NetCDF round-trip tests.""" + return dc.get_example_patch("random_das") + + +@pytest.fixture +def patch_with_attrs(example_patch): + """Create a patch with representative string attrs preserved in NetCDF.""" + return example_patch.update( + attrs={ + "station": "TEST_STATION", + "network": "TEST_NET", + "instrument_id": "TEST_INST", + "acquisition_id": "ACQ_01", + "tag": "netcdf_roundtrip", + "data_type": "strain_rate", + "data_category": "DAS", + } + ) + + +@pytest.fixture +def patch_with_non_dim_coords(example_patch): + """Create a patch with 1D and 2D non-dimensional coordinates.""" + shape = example_patch.shape + latitude = np.linspace(40.0, 41.0, shape[0]) + quality = np.broadcast_to(np.linspace(0.0, 1.0, shape[1]), shape) + coords = example_patch.coords.update( + latitude=("distance", latitude), + quality=(("distance", "time"), quality), + ) + return example_patch.new(coords=coords) + + +@pytest.fixture +def minimal_cf_netcdf_path(tmp_path): + """Create a minimal CF-compliant NetCDF file using only h5py.""" + path = tmp_path / "minimal_cf.nc" + with h5py.File(path, "w") as h5file: + h5file.attrs["Conventions"] = "CF-1.8" + h5file.attrs["_NCProperties"] = "version=2,netcdf=4.9.0" + time_ds = h5file.create_dataset("time", data=np.arange(10)) + time_ds.make_scale("time") + distance_ds = h5file.create_dataset("distance", data=np.arange(5)) + distance_ds.make_scale("distance") + data_ds = h5file.create_dataset("data", data=np.zeros((5, 10))) + data_ds.dims[0].attach_scale(distance_ds) + data_ds.dims[1].attach_scale(time_ds) + return path + + +@pytest.fixture +def rng(): + """Return a seeded generator for synthetic NetCDF test data.""" + return np.random.default_rng(0) + + +class TestNetCDFUtils: + """Tests for NetCDF utility functions.""" + + @pytest.fixture + def test_datetime_array(self): + """Create test datetime array for CF time conversion tests.""" + return np.array( + [ + "2023-01-01T00:00:00", + "2023-01-01T01:00:00", + "2023-01-01T02:00:00", + ], + dtype="datetime64[ns]", + ) + + @pytest.fixture + def expected_cf_times(self): + """Expected CF time values for test datetime array.""" + expected_base = 1672531200.0 # 2023-01-01T00:00:00 in epoch seconds + return np.array([expected_base, expected_base + 3600, expected_base + 7200]) + + @pytest.fixture + def test_patch_attrs(self): + """Create test patch attributes for global attrs tests.""" + return dc.PatchAttrs( + station="TEST_STATION", + network="TEST_NET", + instrument_id="TEST_INST", + ) + + def test_datetime64_to_cf_time(self, test_datetime_array, expected_cf_times): + """Test conversion of datetime64 to CF time format.""" + cf_times = datetime64_to_cf_time(test_datetime_array) + np.testing.assert_array_almost_equal(cf_times, expected_cf_times) + + def test_get_cf_data_attrs(self): + """Test CF data attributes generation.""" + attrs = get_cf_data_attrs("strain_rate") + assert attrs["standard_name"] == "strain_rate" + assert attrs["units"] == "1/s" + assert "_FillValue" in attrs + + def test_get_cf_global_attrs(self, test_patch_attrs): + """Test CF global attributes generation.""" + attrs = get_cf_global_attrs(test_patch_attrs, "1.8") + assert attrs["Conventions"] == "CF-1.8" + assert attrs["station"] == "TEST_STATION" + assert attrs["network"] == "TEST_NET" + assert attrs["instrument"] == "TEST_INST" + assert "institution" not in attrs + + def test_cf_time_to_datetime64(self): + """Test CF time to datetime64 conversion.""" + # Test with seconds since epoch + cf_times = np.array([0, 3600, 7200]) # 0, 1, 2 hours since epoch + units = "seconds since 1970-01-01 00:00:00" + + result = cf_time_to_datetime64(cf_times, units) + + expected = np.array( + ["1970-01-01T00:00:00", "1970-01-01T01:00:00", "1970-01-01T02:00:00"], + dtype="datetime64[ns]", + ) + + np.testing.assert_array_equal(result, expected) + + def test_cf_time_to_datetime64_different_units(self): + """Test CF time conversion with different time units.""" + # Test with days since epoch + cf_times = np.array([0, 1, 2]) + units = "days since 2023-01-01 00:00:00" + + result = cf_time_to_datetime64(cf_times, units) + + expected = np.array( + ["2023-01-01T00:00:00", "2023-01-02T00:00:00", "2023-01-03T00:00:00"], + dtype="datetime64[ns]", + ) + + np.testing.assert_array_equal(result, expected) + + def test_cf_time_to_datetime64_invalid_units(self): + """Test CF time conversion with invalid units.""" + cf_times = np.array([0, 1, 2]) + invalid_units = "invalid format" + + with pytest.raises(ValueError, match="Invalid CF time units format"): + cf_time_to_datetime64(cf_times, invalid_units) + + def test_find_main_data_variable(self, tmp_path, rng): + """Test finding main data variable in NetCDF file.""" + path = tmp_path / "test_data_var.nc" + + with h5py.File(path, "w") as h5file: + # Create dimension scales + time_ds = h5file.create_dataset("time", data=np.arange(100)) + time_ds.make_scale("time") + + # Create various datasets. + h5file.create_dataset("metadata", data=np.array([1, 2, 3])) + h5file.create_dataset("other_data", data=rng.random((50, 50))) + + # Create priority data variable + strain_data = h5file.create_dataset( + "strain_data", data=rng.random((100, 50)) + ) + strain_data.attrs["standard_name"] = b"strain" + + with h5py.File(path, "r") as h5file: + main_var = find_main_data_variable(h5file) + assert main_var == "strain_data" + + def test_find_main_data_variable_no_priority(self, tmp_path, rng): + """Test finding main data variable when no priority match.""" + path = tmp_path / "test_no_priority.nc" + + with h5py.File(path, "w") as h5file: + h5file.create_dataset("first_candidate", data=rng.random((50, 50))) + h5file.create_dataset("second_candidate", data=rng.random((60, 60))) + + with h5py.File(path, "r") as h5file: + main_var = find_main_data_variable(h5file) + assert main_var == "first_candidate" # Returns first candidate + + def test_find_main_data_variable_none_found(self, tmp_path): + """Test finding main data variable when none found.""" + path = tmp_path / "test_none_found.nc" + + with h5py.File(path, "w") as h5file: + h5file.create_dataset("only_1d", data=np.array([1, 2, 3])) # Only 1D data + + with h5py.File(path, "r") as h5file: + main_var = find_main_data_variable(h5file) + assert main_var is None + + def test_get_cf_version(self, tmp_path): + """Test extracting CF version from NetCDF file.""" + path = tmp_path / "test_cf_version.nc" + + # Test CF-1.8 format + with h5py.File(path, "w") as h5file: + h5file.attrs["Conventions"] = np.bytes_("CF-1.8") + + with h5py.File(path, "r") as h5file: + version = get_cf_version(h5file) + assert version == "1.8" + + def test_get_cf_version_space_format(self, tmp_path): + """Test extracting CF version with space format.""" + path = tmp_path / "test_cf_space.nc" + + with h5py.File(path, "w") as h5file: + h5file.attrs["Conventions"] = b"CF 1.7" + + with h5py.File(path, "r") as h5file: + version = get_cf_version(h5file) + assert version == "1.7" + + def test_get_cf_version_with_comma_separated_conventions(self, tmp_path): + """Version parsing should tolerate additional comma-separated conventions.""" + path = tmp_path / "test_cf_comma.nc" + + with h5py.File(path, "w") as h5file: + h5file.attrs["Conventions"] = "CF-1.8, ACDD-1.3" + + with h5py.File(path, "r") as h5file: + version = get_cf_version(h5file) + assert version == "1.8" + + def test_get_cf_version_none(self, tmp_path): + """Test extracting CF version when none present.""" + path = tmp_path / "test_no_cf.nc" + + with h5py.File(path, "w") as h5file: + h5file.attrs["other_attr"] = b"value" + + with h5py.File(path, "r") as h5file: + version = get_cf_version(h5file) + assert version is None + + def test_is_netcdf4_file_from_conventions_bytes(self, tmp_path): + """CF conventions alone should identify a NetCDF-like file.""" + path = tmp_path / "test_conventions.nc" + with h5py.File(path, "w") as h5file: + h5file.attrs["Conventions"] = b"CF-1.8" + + with h5py.File(path, "r") as h5file: + assert is_netcdf4_file(h5file) + + def test_is_netcdf4_file_from_dimension_scales(self, tmp_path): + """Dimension scales alone should not identify a NetCDF-like file.""" + path = tmp_path / "test_dimension_scale.nc" + with h5py.File(path, "w") as h5file: + time_ds = h5file.create_dataset("time", data=np.arange(5)) + time_ds.make_scale("time") + + with h5py.File(path, "r") as h5file: + assert not is_netcdf4_file(h5file) + + def test_get_cf_data_attrs_partial_match(self): + """Partial data-type matches should map to known CF names.""" + attrs = get_cf_data_attrs("my_acceleration_trace") + assert attrs["standard_name"] == "acceleration" + + def test_get_cf_global_attrs_optional_fields(self): + """Optional patch attrs should be propagated to global metadata.""" + attrs = dc.PatchAttrs( + station="TEST_STATION", + network="TEST_NET", + instrument_id="TEST_INST", + acquisition_id="ACQ_01", + tag="taggy", + data_category="DAS", + data_type="strain_rate", + history=("step-1", "step-2"), + category="processed", + ) + + out = get_cf_global_attrs(attrs, "1.8") + + assert out["acquisition"] == "ACQ_01" + assert out["tag"] == "taggy" + assert out["data_category"] == "DAS" + assert out["data_type"] == "strain_rate" + assert out["category"] == "processed" + assert out["processing_history"] == "step-1 | step-2" + + def test_extract_patch_attrs_source_data_type_fallback(self, tmp_path): + """Source data type should populate data_type when primary attr is absent.""" + path = tmp_path / "source_data_type.nc" + with h5py.File(path, "w") as h5file: + h5file.attrs["network"] = np.bytes_("TEST_NET") + h5file.attrs["source_data_type"] = np.bytes_("strain_rate") + + with h5py.File(path, "r") as h5file: + attrs = extract_patch_attrs_from_netcdf(h5file) + + assert attrs["network"] == "TEST_NET" + assert attrs["data_type"] == "strain_rate" + + def test_handle_time_interpolation_without_indices(self, tmp_path): + """Time interpolation should fall back to raw time_values when needed.""" + path = tmp_path / "time_interp.nc" + with h5py.File(path, "w") as h5file: + h5file.create_dataset("time", data=np.arange(3)) + time_values = h5file.create_dataset("time_values", data=np.arange(3)) + time_values.attrs["units"] = "seconds since 1970-01-01 00:00:00" + + with h5py.File(path, "r") as h5file: + out = netcdf_utils._handle_time_interpolation(h5file, "time", np.arange(3)) + + expected = np.array( + ["1970-01-01T00:00:00", "1970-01-01T00:00:01", "1970-01-01T00:00:02"], + dtype="datetime64[ns]", + ) + np.testing.assert_array_equal(out, expected) + + def test_handle_time_interpolation_invalid_units_returns_none(self, tmp_path): + """Interpolation should fail quietly when auxiliary units are invalid.""" + path = tmp_path / "time_interp_invalid.nc" + with h5py.File(path, "w") as h5file: + h5file.create_dataset("time", data=np.arange(3)) + time_values = h5file.create_dataset("time_values", data=np.arange(3)) + time_values.attrs["units"] = "invalid" + + with h5py.File(path, "r") as h5file: + out = netcdf_utils._handle_time_interpolation(h5file, "time", np.arange(3)) + + assert out is None + + def test_read_netcdf_coordinates_with_non_dim_coord_attrs(self, tmp_path): + """Non-dimensional coordinates should honor _DASCORE_DIMS metadata.""" + path = tmp_path / "non_dim_coords.nc" + with h5py.File(path, "w") as h5file: + time_ds = h5file.create_dataset("time", data=np.arange(4)) + time_ds.make_scale("time") + dist_ds = h5file.create_dataset("distance", data=np.arange(3)) + dist_ds.make_scale("distance") + data_ds = h5file.create_dataset("data", data=np.ones((3, 4))) + data_ds.dims[0].attach_scale(dist_ds) + data_ds.dims[1].attach_scale(time_ds) + data_ds.attrs["coordinates"] = np.bytes_("latitude") + lat_ds = h5file.create_dataset("latitude", data=np.linspace(1.0, 2.0, 3)) + lat_ds.attrs["_DASCORE_DIMS"] = np.bytes_("distance") + + with h5py.File(path, "r") as h5file: + coords = read_netcdf_coordinates(h5file) + + assert coords.dims == ("distance", "time") + assert "latitude" in coords.coord_map + assert coords.dim_map["latitude"] == ("distance",) + + def test_validate_cf_compliance_invalid_time_units(self, tmp_path): + """Invalid time units should be reported as a CF issue.""" + path = tmp_path / "invalid_time_units.nc" + with h5py.File(path, "w") as h5file: + h5file.attrs["Conventions"] = "CF-1.8" + time_ds = h5file.create_dataset("time", data=np.arange(4)) + time_ds.make_scale("time") + time_ds.attrs["units"] = "seconds" + + with h5py.File(path, "r") as h5file: + issues = validate_cf_compliance(h5file) + + assert any("invalid units" in issue for issue in issues) + + def test_is_netcdf4_file_handles_attribute_error(self): + """Attribute errors during detection should return False.""" + + class BrokenFile: + @property + def attrs(self): + raise AttributeError + + assert not is_netcdf4_file(BrokenFile()) + + def test_is_netcdf4_file_decodes_byte_conventions_from_mapping(self): + """Conventions bytes from a mapping-like attrs object should decode.""" + + class FakeFile: + def __init__(self): + self.attrs = {"Conventions": b"CF-1.8"} + + def values(self): + return [] + + assert is_netcdf4_file(FakeFile()) + + def test_is_netcdf4_file_from_dimension_list_attr(self, tmp_path): + """DIMENSION_LIST attrs alone should not mark a file as NetCDF-like.""" + path = tmp_path / "test_dimension_list.nc" + with h5py.File(path, "w") as h5file: + data_ds = h5file.create_dataset("data", data=np.ones((2, 2))) + data_ds.attrs["DIMENSION_LIST"] = "present" + + with h5py.File(path, "r") as h5file: + assert not is_netcdf4_file(h5file) + + def test_get_cf_version_decodes_bytes(self, tmp_path): + """CF version extraction should decode byte attrs.""" + path = tmp_path / "cf_version_bytes.nc" + with h5py.File(path, "w") as h5file: + h5file.attrs["Conventions"] = np.bytes_("CF-1.9") + + with h5py.File(path, "r") as h5file: + assert get_cf_version(h5file) == "1.9" + + def test_handle_time_interpolation_decodes_units_and_handles_bad_units( + self, tmp_path + ): + """Interpolation should decode bytes and handle unsupported units.""" + decoded_path = tmp_path / "time_interp_bytes.nc" + with h5py.File(decoded_path, "w") as h5file: + h5file.create_dataset("time", data=np.arange(3)) + time_values = h5file.create_dataset("time_values", data=np.arange(3)) + time_values.attrs["units"] = np.bytes_("seconds since 1970-01-01 00:00:00") + time_indices = h5file.create_dataset("time_indices", data=np.array([0])) + time_indices[...] = np.array([0]) + + with h5py.File(decoded_path, "r") as h5file: + out = netcdf_utils._handle_time_interpolation(h5file, "time", np.arange(3)) + + assert out.dtype == "datetime64[ns]" + + invalid_path = tmp_path / "time_interp_unsupported.nc" + with h5py.File(invalid_path, "w") as h5file: + h5file.create_dataset("time", data=np.arange(3)) + time_values = h5file.create_dataset("time_values", data=np.arange(3)) + time_values.attrs["units"] = "fortnights since 2023-01-01" + + with h5py.File(invalid_path, "r") as h5file: + out = netcdf_utils._handle_time_interpolation(h5file, "time", np.arange(3)) + + assert out is None + + def test_read_netcdf_coordinates_handles_bad_refs_and_auxiliaries(self, tmp_path): + """Coordinate reader should tolerate bad refs and skip auxiliaries.""" + path = tmp_path / "coords_bad_refs.nc" + with h5py.File(path, "w") as h5file: + h5file.create_group("aux_group") + h5file.create_dataset("time_values", data=np.arange(3)) + data_ds = h5file.create_dataset("data", data=np.ones((2, 3))) + data_ds.attrs["DIMENSION_LIST"] = np.array([[b"/missing"]], dtype="S8") + data_ds.attrs["coordinates"] = np.bytes_("missing aux_group latitude") + lat_ds = h5file.create_dataset("latitude", data=np.arange(2)) + lat_ds.attrs["_DASCORE_DIMS"] = np.bytes_("distance") + dist_ds = h5file.create_dataset("distance", data=np.arange(2)) + dist_ds.make_scale("distance") + time_ds = h5file.create_dataset("time", data=np.arange(3)) + time_ds.make_scale("time") + time_ds.attrs["units"] = np.bytes_("fortnights since 2023-01-01") + + with h5py.File(path, "r") as h5file: + coords = read_netcdf_coordinates(h5file) + + assert coords.dims == ("distance", "time") + assert "latitude" in coords.coord_map + + def test_read_netcdf_coordinates_without_dimension_list(self, tmp_path): + """Coordinate reader should fall back to discovered order when needed.""" + path = tmp_path / "coords_no_dimension_list.nc" + with h5py.File(path, "w") as h5file: + dist_ds = h5file.create_dataset("distance", data=np.arange(2)) + dist_ds.make_scale("distance") + time_ds = h5file.create_dataset("time", data=np.arange(3)) + time_ds.make_scale("time") + time_ds.attrs["units"] = "seconds since 1970-01-01 00:00:00" + h5file.create_dataset("data", data=np.ones((2, 3))) + + with h5py.File(path, "r") as h5file: + coords = read_netcdf_coordinates(h5file) + + assert coords.dims == ("distance", "time") + np.testing.assert_array_equal(coords.coord_map["distance"].values, np.arange(2)) + + def test_read_netcdf_coordinates_adds_unordered_dim_coords(self, tmp_path): + """Coords discovered outside DIMENSION_LIST should still be retained.""" + path = tmp_path / "coords_extra_scale.nc" + with h5py.File(path, "w") as h5file: + dist_ds = h5file.create_dataset("distance", data=np.arange(2)) + dist_ds.make_scale("distance") + time_ds = h5file.create_dataset("time", data=np.arange(3)) + time_ds.make_scale("time") + time_ds.attrs["units"] = np.bytes_("seconds since 1970-01-01 00:00:00") + extra_ds = h5file.create_dataset("channel", data=np.arange(2)) + extra_ds.make_scale("channel") + data_ds = h5file.create_dataset("data", data=np.ones((2, 3))) + data_ds.dims[0].attach_scale(dist_ds) + data_ds.dims[1].attach_scale(time_ds) + + with h5py.File(path, "r") as h5file: + coords = read_netcdf_coordinates(h5file) + + assert "channel" in coords.coord_map + + def test_read_netcdf_coordinates_skips_auxiliary_named_scales(self, tmp_path): + """Auxiliary scale names should be skipped before coord collection.""" + path = tmp_path / "coords_aux_skip.nc" + with h5py.File(path, "w") as h5file: + dist_ds = h5file.create_dataset("distance", data=np.arange(2)) + dist_ds.make_scale("distance") + time_ds = h5file.create_dataset("time", data=np.arange(3)) + time_ds.make_scale("time") + time_ds.attrs["units"] = "seconds since 1970-01-01 00:00:00" + aux_ds = h5file.create_dataset("time_values", data=np.arange(3)) + aux_ds.make_scale("time_values") + h5file.create_dataset("data", data=np.ones((2, 3))) + + with h5py.File(path, "r") as h5file: + coords = read_netcdf_coordinates(h5file) + + assert "time_values" not in coords.coord_map + + +class TestNetCDFCoreHelpers: + """Direct tests for lightweight NetCDF core helpers.""" + + def test_coord_attrs_cover_distance_and_depth(self, example_patch): + """Coordinate helper should populate distance and depth CF attrs.""" + distance_attrs = netcdf_utils.coord_attrs( + "distance", example_patch.coords.coord_map["distance"] + ) + depth_coord = dc.get_coord(data=np.arange(3), units="ft") + depth_attrs = netcdf_utils.coord_attrs("sensor_depth", depth_coord) + + assert distance_attrs["standard_name"] == "distance" + assert distance_attrs["units"] == "m" + assert depth_attrs["units"] == "ft" + + def test_get_xarray_data_var_name(self): + """Dataset helper should find expected data variables or raise.""" + xr = pytest.importorskip("xarray") + ds_with_data = xr.Dataset({"data": (("x",), [1, 2])}) + ds_single = xr.Dataset({"signal": (("x",), [1, 2])}) + ds_multi = xr.Dataset({"signal": (("x",), [1, 2]), "other": (("x",), [3, 4])}) + + class _DatasetWithNone: + data_vars = {None: object(), "distance_indices": object()} + + assert netcdf_utils.get_xarray_data_var_name(ds_with_data) == "data" + assert netcdf_utils.get_xarray_data_var_name(_DatasetWithNone()) is None + assert netcdf_utils.get_xarray_data_var_name(ds_single) == "signal" + with pytest.raises(ValueError, match="No suitable data variable found"): + netcdf_utils.get_xarray_data_var_name(ds_multi) + + def test_get_format_false_cases(self, minimal_cf_netcdf_path, tmp_path): + """Format detection should reject missing and invalid CF versions.""" + formatter = netcdf_core.NetCDFCFV18() + + no_cf_path = tmp_path / "no_cf_version.nc" + with h5py.File(no_cf_path, "w") as h5file: + h5file.attrs["_NCProperties"] = "version=2" + + invalid_cf_path = tmp_path / "invalid_cf_version.nc" + with h5py.File(invalid_cf_path, "w") as h5file: + h5file.attrs["Conventions"] = "CF-not-a-number" + + assert dc.get_format(minimal_cf_netcdf_path) == ("NETCDF_CF", "1.8") + with h5py.File(no_cf_path, "r") as h5file: + assert formatter.get_format(h5file) is False + with h5py.File(invalid_cf_path, "r") as h5file: + assert formatter.get_format(h5file) is False + + def test_get_format_accepts_comma_separated_conventions(self, tmp_path): + """Format detection should accept CF versions followed by extra conventions.""" + formatter = netcdf_core.NetCDFCFV18() + path = tmp_path / "comma_conventions.nc" + + with h5py.File(path, "w") as h5file: + h5file.attrs["Conventions"] = "CF-1.8, ACDD-1.3" + h5file.attrs["_NCProperties"] = "version=2" + time_ds = h5file.create_dataset("time", data=np.arange(3)) + time_ds.make_scale("time") + dist_ds = h5file.create_dataset("distance", data=np.arange(2)) + dist_ds.make_scale("distance") + data_ds = h5file.create_dataset("data", data=np.ones((2, 3))) + data_ds.dims[0].attach_scale(dist_ds) + data_ds.dims[1].attach_scale(time_ds) + + with h5py.File(path, "r") as h5file: + assert formatter.get_format(h5file) == ("NETCDF_CF", "1.8") + + def test_get_data_variable_name_raises_for_missing_data(self, tmp_path): + """Formatter should raise when no main data variable exists.""" + formatter = netcdf_core.NetCDFCFV18() + path = tmp_path / "missing_data.nc" + with h5py.File(path, "w") as h5file: + h5file.create_dataset("time", data=np.arange(3)) + + with h5py.File(path, "r") as h5file: + with pytest.raises(ValueError, match="No suitable data variable found"): + formatter._get_data_variable_name(h5file) + + def test_apply_coordinate_filtering_handles_selected_and_unrelated_kwargs( + self, example_patch + ): + """Filtering should only use known coordinate kwargs.""" + formatter = netcdf_core.NetCDFCFV18() + + unfiltered = formatter._apply_coordinate_filtering(example_patch, {"tag": "x"}) + filtered = formatter._apply_coordinate_filtering( + example_patch, {"distance": (10, 20)} + ) + + assert unfiltered.equals(example_patch) + assert filtered.data.shape[0] < example_patch.data.shape[0] + + def test_read_uses_xarray_dataset_and_merges_missing_coords( + self, minimal_cf_netcdf_path, monkeypatch + ): + """Read should accept xarray-backed data and merge extra coords.""" + + class FakeCoord: + def __init__(self, dims, values): + self.dims = dims + self.values = values + + class FakeDataArray: + def __init__(self, data): + self.data = data + self.coords = { + "distance": FakeCoord(("distance",), np.arange(data.shape[0])), + "time": FakeCoord(("time",), np.arange(data.shape[1])), + "latitude": FakeCoord( + ("distance",), np.linspace(1.0, 2.0, data.shape[0]) + ), + } + + def load(self): + return self + + class FakeDataset: + def __init__(self, data_array): + self.data_vars = {None: data_array} + self._data_array = data_array + + def __enter__(self): + return self + + def __exit__(self, *args): + return False + + def get(self, name): + return None + + def __getitem__(self, item): + assert item is None + return self._data_array + + fake_coords = dc.get_coord_manager( + coords={"distance": np.arange(2), "time": np.arange(3)}, + dims=("distance", "time"), + ) + fake_data_array = FakeDataArray(np.arange(6).reshape(2, 3)) + fake_dataset = FakeDataset(fake_data_array) + fake_xarray = type( + "FakeXarray", + (), + {"open_dataset": staticmethod(lambda *args, **kwargs: fake_dataset)}, + ) + formatter = netcdf_core.NetCDFCFV18() + + def _optional_import(name, on_missing="raise"): + if name == "xarray": + return fake_xarray + if name == "netCDF4": + return object() + return None + + monkeypatch.setattr(netcdf_core, "optional_import", _optional_import) + monkeypatch.setattr( + netcdf_core, + "read_netcdf_coordinates", + lambda *args: fake_coords, + ) + monkeypatch.setattr( + formatter, "_get_data_variable_name", lambda resource: "__values__" + ) + monkeypatch.setattr( + formatter, "_get_patch_attrs", lambda resource: {"tag": "fake"} + ) + + spool = formatter.read(minimal_cf_netcdf_path) + patch = spool[0] + + np.testing.assert_array_equal(patch.data, fake_data_array.data) + assert patch.attrs.tag == "fake" + assert "latitude" in patch.coords.coord_map + + def test_read_falls_back_to_single_xarray_data_var( + self, minimal_cf_netcdf_path, monkeypatch + ): + """Read should use the only xarray data var when HDF-derived name misses.""" + + class FakeDataArray: + def __init__(self): + self.data = np.ones((2, 2)) + self.coords = {} + + def load(self): + return self + + class FakeDataset: + def __init__(self): + self.data_vars = {"signal": FakeDataArray()} + + def __enter__(self): + return self + + def __exit__(self, *args): + return False + + def get(self, name): + return None + + def __getitem__(self, item): + return self.data_vars[item] + + fake_coords = dc.get_coord_manager( + coords={"distance": np.arange(2), "time": np.arange(2)}, + dims=("distance", "time"), + ) + fake_xarray = type( + "FakeXarray", + (), + {"open_dataset": staticmethod(lambda *args, **kwargs: FakeDataset())}, + ) + formatter = netcdf_core.NetCDFCFV18() + + def _optional_import(name, on_missing="raise"): + if name == "xarray": + return fake_xarray + if name == "netCDF4": + return object() + return None + + monkeypatch.setattr(netcdf_core, "optional_import", _optional_import) + monkeypatch.setattr( + netcdf_core, "read_netcdf_coordinates", lambda *args: fake_coords + ) + monkeypatch.setattr( + formatter, "_get_data_variable_name", lambda resource: "missing" + ) + monkeypatch.setattr( + formatter, "_get_patch_attrs", lambda resource: {"tag": "fallback"} + ) + + spool = formatter.read(minimal_cf_netcdf_path) + assert spool[0].attrs.tag == "fallback" + + def test_read_returns_empty_spool_for_empty_filtered_patch( + self, minimal_cf_netcdf_path, monkeypatch + ): + """Read should return an empty spool after filtering removes all data.""" + + class FakeDataArray: + def __init__(self): + self.data = np.ones((1, 1)) + self.coords = {} + + def load(self): + return self + + class FakeDataset: + def __init__(self): + self.data_vars = {"data": FakeDataArray()} + + def __enter__(self): + return self + + def __exit__(self, *args): + return False + + def get(self, name): + return self.data_vars["data"] + + formatter = netcdf_core.NetCDFCFV18() + fake_xarray = type( + "FakeXarray", + (), + {"open_dataset": staticmethod(lambda *args, **kwargs: FakeDataset())}, + ) + fake_coords = dc.get_coord_manager( + coords={"distance": np.arange(1), "time": np.arange(1)}, + dims=("distance", "time"), + ) + empty_patch = dc.Patch( + data=np.empty((0, 0)), + coords=dc.get_coord_manager( + coords={"distance": np.array([]), "time": np.array([])}, + dims=("distance", "time"), + ), + dims=("distance", "time"), + attrs={"tag": "empty"}, + ) + + def _optional_import(name, on_missing="raise"): + if name == "xarray": + return fake_xarray + if name == "netCDF4": + return object() + return None + + monkeypatch.setattr(netcdf_core, "optional_import", _optional_import) + monkeypatch.setattr( + netcdf_core, "read_netcdf_coordinates", lambda *args: fake_coords + ) + monkeypatch.setattr( + formatter, "_get_data_variable_name", lambda resource: "data" + ) + monkeypatch.setattr( + formatter, "_get_patch_attrs", lambda resource: {"tag": "empty"} + ) + monkeypatch.setattr( + formatter, "_apply_coordinate_filtering", lambda patch, kwargs: empty_patch + ) + + spool = formatter.read(minimal_cf_netcdf_path) + assert len(spool) == 0 + + def test_write_uses_xarray_dataset_path(self, example_patch, tmp_path, monkeypatch): + """Write should pass CF attrs and encoding through the xarray path.""" + + class FakeCoord: + def __init__(self): + self.attrs = {} + + class FakeDataVar: + def __init__(self): + self.attrs = {} + + class FakeDataset: + def __init__(self): + self.attrs = {} + self.data = FakeDataVar() + self.to_netcdf_calls = [] + + def __getitem__(self, item): + assert item == "data" + return self.data + + def to_netcdf(self, *args, **kwargs): + self.to_netcdf_calls.append((args, kwargs)) + + class FakeDataArray: + def __init__(self): + self.coords = { + "distance": FakeCoord(), + "time": FakeCoord(), + "latitude": FakeCoord(), + } + self.dataset = FakeDataset() + + def rename(self, name): + assert name == "data" + return self + + def to_dataset(self): + return self.dataset + + fake_data_array = FakeDataArray() + formatter = netcdf_core.NetCDFCFV18() + + def _optional_import(name, on_missing="raise"): + if name == "xarray": + return object() + if name == "netCDF4": + return object() + return None + + monkeypatch.setattr( + netcdf_core, + "optional_import", + _optional_import, + ) + monkeypatch.setattr( + netcdf_core, + "patch_to_xarray", + lambda patch: fake_data_array, + ) + + patch_with_partial = example_patch.new( + coords=example_patch.coords.update( + latitude=("distance", np.linspace(0.0, 1.0, example_patch.shape[0])), + quality=("distance", dc.get_coord(shape=(example_patch.shape[0],))), + ) + ) + out_path = tmp_path / "write_stub.nc" + formatter.write( + patch_with_partial.update(attrs={"data_type": "strain_rate"}), out_path + ) + + dataset = fake_data_array.dataset + assert dataset.attrs["Conventions"] == "CF-1.8" + assert dataset.attrs["source_data_type"] == "strain_rate" + assert dataset.data.attrs["standard_name"] == "strain_rate" + assert dataset.data.attrs["coordinates"] == "latitude" + assert dataset.to_netcdf_calls + _args, kwargs = dataset.to_netcdf_calls[0] + assert kwargs["encoding"] is None + + +class TestNetCDFIO: + """Tests for NetCDF IO functionality.""" + + @pytest.fixture + def example_patch(self): + """Create an example patch for testing.""" + return dc.get_example_patch("random_das") + + @pytest.fixture + def netcdf_path(self, example_patch, tmp_path): + """Create a test NetCDF file.""" + _require_xarray_netcdf_engine() + path = tmp_path / "test.nc" + # Write patch to NetCDF format + dc.write(example_patch, path, file_format="netcdf_cf") + return path + + def test_write_netcdf(self, example_patch, tmp_path): + """Test writing a patch to NetCDF format.""" + _require_xarray_netcdf_engine() + path = tmp_path / "test_write.nc" + + # Write patch + dc.write(example_patch, path, file_format="netcdf_cf") + + # Check file exists and is valid HDF5/NetCDF + assert path.exists() + + # Open with h5py to verify structure + with h5py.File(path, "r") as h5file: + # Check it's detected as NetCDF + assert is_netcdf4_file(h5file) + + # Check global attributes + assert "Conventions" in h5file.attrs + conventions = h5file.attrs["Conventions"] + if isinstance(conventions, bytes): + conventions = conventions.decode() + assert "CF-" in conventions + + # Check coordinates exist + assert "time" in h5file + assert "distance" in h5file + + # Check data variable exists + assert "data" in h5file + + # Check dimension scales + assert h5file["time"].is_scale + assert h5file["distance"].is_scale + + def test_read_netcdf(self, netcdf_path): + """Test reading a NetCDF file.""" + spool = dc.read(netcdf_path, file_format="netcdf_cf") + patch = spool[0] + assert isinstance(patch, dc.Patch) + assert "time" in patch.coords + assert "distance" in patch.coords + assert patch.data.ndim == 2 + assert patch.attrs["_source_patch_id"] == "data" + + def test_scan_netcdf(self, netcdf_path): + """Test scanning a NetCDF file for metadata.""" + summary_list = dc.scan(netcdf_path, file_format="netcdf_cf") + assert len(summary_list) == 1 + summary = summary_list[0] + assert summary.source_format == "NETCDF_CF" + assert "time" in summary.coords + assert "distance" in summary.coords + assert summary.source_patch_id == "data" + + def test_get_format(self, minimal_cf_netcdf_path): + """Test format detection.""" + assert dc.get_format(minimal_cf_netcdf_path) == ("NETCDF_CF", "1.8") + + def test_get_format_without_xarray_import( + self, minimal_cf_netcdf_path, monkeypatch + ): + """Format detection should not depend on xarray being importable.""" + import importlib + + original_import_module = importlib.import_module + + def _import_module(name, package=None): + if name == "xarray": + raise ImportError("xarray disabled for test") + return original_import_module(name, package) + + monkeypatch.setattr(importlib, "import_module", _import_module) + + assert dc.get_format(minimal_cf_netcdf_path) == ("NETCDF_CF", "1.8") + + def test_round_trip(self, example_patch, tmp_path): + """Test round-trip: patch -> NetCDF -> patch.""" + _require_xarray_netcdf_engine() + path = tmp_path / "roundtrip.nc" + + # Write and read back + dc.write(example_patch, path, file_format="netcdf_cf") + spool = dc.read(path, file_format="netcdf_cf") + recovered_patch = spool[0] + + # Check data preservation + np.testing.assert_array_almost_equal( + example_patch.data, recovered_patch.data, decimal=6 + ) + + # Check coordinate preservation + for coord_name in example_patch.coords.coord_map: + orig_coord = example_patch.coords.get_array(coord_name) + recovered_coord = recovered_patch.coords.get_array(coord_name) + + if coord_name == "time": + # Time coordinates might have slight precision differences + # due to CF time conversion (float64 seconds -> datetime64[ns]) + time_diff = np.abs(orig_coord - recovered_coord) + assert np.all( + time_diff < np.timedelta64(200, "us") + ) # 200 microsecond tolerance + else: + np.testing.assert_array_almost_equal( + orig_coord, recovered_coord, decimal=6 + ) + + +class TestNetCDFXarrayCompatibility: + """Tests for compatibility between DASCore NetCDF output and xarray.""" + + @pytest.fixture( + params=("example_patch", "patch_with_attrs", "patch_with_non_dim_coords") + ) + def patch_variant(self, request): + """Parametrize representative patch variants for compatibility checks.""" + return request.getfixturevalue(request.param) + + def test_patch_to_xarray_can_be_serialized_with_xarray( + self, patch_variant, tmp_path + ): + """A patch converted to xarray should round-trip through xarray IO.""" + xr = pytest.importorskip("xarray") + path = tmp_path / "xarray_roundtrip.nc" + + data_array = dc.io.patch_to_xarray(patch_variant) + data_array.to_netcdf(path, engine="scipy") + + reopened = xr.open_dataarray(path, engine="scipy") + round_tripped = dc.io.xarray_to_patch(reopened) + reopened.close() + + _assert_patch_round_trip_equal(patch_variant, round_tripped) + + def test_current_netcdf_output_is_readable_by_xarray(self, patch_variant, tmp_path): + """DASCore NetCDF output should be consumable by xarray backends.""" + xr = pytest.importorskip("xarray") + engine = _require_xarray_netcdf_engine() + + path = tmp_path / "dascore_netcdf.nc" + dc.write(patch_variant, path, file_format="netcdf_cf") + + data_array = xr.open_dataarray(path, engine=engine) + round_tripped = dc.io.xarray_to_patch(data_array) + data_array.close() + + _assert_patch_compatible_with_xarray_output(patch_variant, round_tripped) + + def test_written_file_exposes_expected_metadata_to_xarray( + self, patch_with_attrs, tmp_path + ): + """Xarray should see the expected CF metadata on DASCore output.""" + xr = pytest.importorskip("xarray") + engine = _require_xarray_netcdf_engine() + + path = tmp_path / "metadata_visible.nc" + dc.write(patch_with_attrs, path, file_format="netcdf_cf") + + with xr.open_dataset(path, engine=engine) as dataset: + assert dataset.attrs["Conventions"] == "CF-1.8" + assert dataset.attrs["source_data_type"] == "strain_rate" + assert dataset.attrs["station"] == "TEST_STATION" + assert dataset.attrs["network"] == "TEST_NET" + + assert list(dataset.data_vars) == ["data"] + assert list(dataset.coords) == ["distance", "time"] + assert dataset["data"].attrs["standard_name"] == "strain_rate" + assert dataset["data"].attrs["units"] == "1/s" + assert dataset["time"].attrs["axis"] == "T" + assert dataset["time"].attrs["standard_name"] == "time" + assert dataset["distance"].attrs["standard_name"] == "distance" + assert dataset["distance"].attrs["units"] == "m" + assert "featureType" not in dataset.attrs + assert "coordinates" not in dataset["data"].attrs + + def test_dascore_and_xarray_roundtrip_agree_for_non_dim_coords( + self, patch_with_non_dim_coords, tmp_path + ): + """The same patch should round-trip through both IO paths identically.""" + xr = pytest.importorskip("xarray") + _require_xarray_netcdf_engine() + path = tmp_path / "xarray_non_dim.nc" + + dascore_path = tmp_path / "dascore_non_dim.nc" + dc.write(patch_with_non_dim_coords, dascore_path, file_format="netcdf_cf") + dascore_patch = dc.read(dascore_path, file_format="netcdf_cf")[0] + + data_array = dc.io.patch_to_xarray(patch_with_non_dim_coords) + data_array.to_netcdf(path, engine="scipy") + reopened = xr.open_dataarray(path, engine="scipy") + xarray_patch = dc.io.xarray_to_patch(reopened) + reopened.close() + + _assert_patch_round_trip_equal(patch_with_non_dim_coords, dascore_patch) + _assert_patch_round_trip_equal(patch_with_non_dim_coords, xarray_patch) + _assert_patch_round_trip_equal(dascore_patch, xarray_patch) + + +class TestNetCDFEdgeCases: + """Test edge cases and error conditions.""" + + @pytest.fixture + def multi_patch_spool(self): + """Create a spool with multiple patches for testing.""" + patch1 = dc.get_example_patch("random_das") + patch2 = dc.get_example_patch("random_das") + return dc.spool([patch1, patch2]) + + @pytest.fixture + def invalid_hdf5_file(self, tmp_path): + """Create an invalid HDF5 file (not NetCDF).""" + path = tmp_path / "invalid.nc" + with h5py.File(path, "w") as h5file: + rng = np.random.default_rng() + h5file.create_dataset("random_data", data=rng.standard_normal((100, 50))) + return path + + @pytest.fixture + def compressed_netcdf_file(self, tmp_path): + """Create a compressed NetCDF file for testing.""" + _require_xarray_netcdf_engine() + patch = dc.get_example_patch("random_das") + path = tmp_path / "compressed.nc" + dc.write( + patch, + path, + file_format="netcdf_cf", + compression="gzip", + compression_opts=9, + ) + return path, patch + + def test_empty_spool_write_error(self, tmp_path): + """Test that writing empty spool raises error.""" + path = tmp_path / "empty.nc" + empty_spool = dc.spool([]) + + with pytest.raises(ValueError, match="Cannot write empty spool"): + dc.write(empty_spool, path, file_format="netcdf_cf") + + def test_multi_patch_write_error(self, multi_patch_spool, tmp_path): + """Test that multi-patch spool raises NotImplementedError.""" + path = tmp_path / "multi.nc" + + with pytest.raises( + NotImplementedError, match="Multi-patch spools not yet supported" + ): + dc.write(multi_patch_spool, path, file_format="netcdf_cf") + + def test_invalid_netcdf_file(self, invalid_hdf5_file): + """Test behavior with invalid NetCDF file.""" + with h5py.File(invalid_hdf5_file, "r") as h5file: + assert not is_netcdf4_file(h5file) + + def test_compression_options(self, compressed_netcdf_file): + """Test NetCDF file creation with compression options.""" + path, original_patch = compressed_netcdf_file + + spool = dc.read(path, file_format="netcdf_cf") + recovered_patch = spool[0] + + np.testing.assert_array_almost_equal( + original_patch.data, recovered_patch.data, decimal=6 + ) + +class TestNetCDFUtilsAdvanced: + """Additional tests for NetCDF utility functions.""" + + @pytest.fixture + def cf_compliant_file(self, tmp_path): + """Create a CF-compliant NetCDF file for testing.""" + _require_xarray_netcdf_engine() + path = tmp_path / "cf_compliant.nc" + patch = dc.get_example_patch("random_das") + patch = patch.update( + attrs={ + "station": "TEST_STATION", + "network": "TEST_NET", + "data_type": "strain_rate", + } + ) + dc.write(patch, path, file_format="netcdf_cf") + return path + + def test_extract_patch_attrs_from_netcdf(self, cf_compliant_file): + """Test extracting patch attributes from NetCDF file.""" + with h5py.File(cf_compliant_file, "r") as h5file: + attrs = extract_patch_attrs_from_netcdf(h5file) + assert isinstance(attrs, dict) + + def test_read_netcdf_coordinates(self, cf_compliant_file): + """Test reading coordinates from NetCDF file.""" + with h5py.File(cf_compliant_file, "r") as h5file: + coords = read_netcdf_coordinates(h5file) + + assert "time" in coords.coord_map + assert "distance" in coords.coord_map + + # Check that time coordinate is properly converted from CF format + time_coord = coords.coord_map["time"] + assert len(time_coord) > 0 + + def test_validate_cf_compliance(self, cf_compliant_file): + """Test CF compliance validation.""" + with h5py.File(cf_compliant_file, "r") as h5file: + issues = validate_cf_compliance(h5file) + + # Our implementation should produce CF-compliant files + assert len(issues) == 0, f"CF compliance issues found: {issues}" + + def test_validate_cf_compliance_with_issues(self, tmp_path, rng): + """Test CF compliance validation with non-compliant file.""" + path = tmp_path / "non_compliant.nc" + + # Create a file with CF issues + with h5py.File(path, "w") as h5file: + # Missing Conventions attribute + time_ds = h5file.create_dataset("time", data=np.arange(100)) + time_ds.make_scale("time") + # Missing units attribute + + h5file.create_dataset("data", data=rng.random((100, 50))) + # Missing long_name and units attributes + + with h5py.File(path, "r") as h5file: + issues = validate_cf_compliance(h5file) + + assert len(issues) > 0 + assert any("Conventions" in issue for issue in issues) + assert any("units" in issue for issue in issues) + + def test_coordinate_filtering_during_read(self, cf_compliant_file): + """Test coordinate filtering during NetCDF read.""" + # Read with time filtering + spool = dc.read(cf_compliant_file, file_format="netcdf_cf") + original_patch = spool[0] + + # Get time bounds for filtering + time_coord = original_patch.coords.get_array("time") + time_start = time_coord[10] + time_end = time_coord[50] + + # Read with filtering + filtered_spool = dc.read( + cf_compliant_file, file_format="netcdf_cf", time=(time_start, time_end) + ) + filtered_patch = filtered_spool[0] + + # Should have fewer time samples + assert filtered_patch.data.shape[1] < original_patch.data.shape[1] + + def test_different_data_types_cf_attrs(self): + """Test CF attributes for different data types.""" + # Test various data types + test_cases = [ + ("strain", "1", "Strain"), + ("velocity", "m/s", "Velocity"), + ("temperature", "K", "Temperature"), + ("pressure", "Pa", "Pressure"), + ("unknown_type", "1", "Distributed Acoustic Sensing data"), + ] + + for data_type, expected_units, expected_long_name in test_cases: + attrs = get_cf_data_attrs(data_type) + assert attrs["units"] == expected_units + assert attrs["long_name"] == expected_long_name + assert "_FillValue" in attrs + + def test_netcdf_format_detection_edge_cases(self, tmp_path, rng): + """Test NetCDF format detection edge cases.""" + _require_xarray_netcdf_engine() + # Test with newer CF version via file read + path1 = tmp_path / "newer_cf.nc" + with h5py.File(path1, "w") as h5file: + h5file.attrs["Conventions"] = b"CF-2.0" + h5file.attrs["_NCProperties"] = b"version=2,netcdf=test" + # Add required data for valid NetCDF + time_ds = h5file.create_dataset("time", data=np.arange(10)) + time_ds.make_scale("time") + dist_ds = h5file.create_dataset("distance", data=np.arange(50)) + dist_ds.make_scale("distance") + h5file.create_dataset("data", data=rng.random((50, 10))) + + # Test that it can be detected and read + spool = dc.read(path1, file_format="netcdf_cf") + assert len(spool) == 1 + + # Test with older supported CF version + path2 = tmp_path / "older_cf.nc" + with h5py.File(path2, "w") as h5file: + h5file.attrs["Conventions"] = b"CF-1.6" + h5file.attrs["_NCProperties"] = b"version=2,netcdf=test" + time_ds = h5file.create_dataset("time", data=np.arange(10)) + time_ds.make_scale("time") + dist_ds = h5file.create_dataset("distance", data=np.arange(50)) + dist_ds.make_scale("distance") + h5file.create_dataset("data", data=rng.random((50, 10))) + + # Should be able to read older CF versions + spool = dc.read(path2, file_format="netcdf_cf") + assert len(spool) == 1 + + def test_read_external_xdas_netcdf_file(self): + """External xdas NetCDF should keep its current readable structure.""" + _require_xarray_netcdf_engine() + path = fetch("xdas_netcdf.nc") + + patch = dc.read(path, file_format="netcdf_cf")[0] + + assert patch.dims == ("time", "distance") + assert patch.shape == (300, 401) + assert set(patch.coords.coord_map) == {"time", "distance"} + assert patch.attrs.tag == "" + assert patch.attrs["_source_patch_id"] == netcdf_utils.XDAS_PAYLOAD_VARIABLE + + def test_error_conditions(self, tmp_path): + """Test various error conditions.""" + _require_xarray_netcdf_engine() + # Test reading file with no data variables + path = tmp_path / "no_data.nc" + with h5py.File(path, "w") as h5file: + h5file.attrs["Conventions"] = b"CF-1.8" + time_ds = h5file.create_dataset("time", data=np.arange(10)) + time_ds.make_scale("time") + + # Use DASCore's standard reading interface to test error conditions + with pytest.raises(ValueError, match="No suitable data variable found"): + dc.read(path, file_format="netcdf_cf") + + def test_cf_time_edge_cases(self): + """Test CF time conversion edge cases.""" + # Test unsupported time unit + cf_times = np.array([0, 1, 2]) + invalid_units = "fortnights since 2023-01-01" + + with pytest.raises(ValueError, match="Unsupported time unit"): + cf_time_to_datetime64(cf_times, invalid_units) + + # Test various supported units + supported_units = [ + ("hours since 2023-01-01", "h"), + ("minutes since 2023-01-01", "min"), + ("milliseconds since 2023-01-01", "ms"), + ("microseconds since 2023-01-01", "us"), + ] + + for units, _ in supported_units: + cf_times = np.array([0, 1, 2]) + result = cf_time_to_datetime64(cf_times, units) + assert result.dtype == "datetime64[ns]" + assert len(result) == 3 From 866cfa8c4267c28861df0df8b127e7885ae05523 Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Thu, 9 Apr 2026 13:35:26 +0200 Subject: [PATCH 2/5] start refactor2 --- dascore/io/netcdf/cf.py | 134 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 134 insertions(+) create mode 100644 dascore/io/netcdf/cf.py diff --git a/dascore/io/netcdf/cf.py b/dascore/io/netcdf/cf.py new file mode 100644 index 000000000..6a114833c --- /dev/null +++ b/dascore/io/netcdf/cf.py @@ -0,0 +1,134 @@ +"""CF metadata helpers for NetCDF IO. + +See https://cfconventions.org/ for the Climate and Forecast metadata standard. +""" + +from __future__ import annotations + +import datetime +from typing import TYPE_CHECKING + +import numpy as np + +import dascore as dc +from dascore.utils.time import to_datetime64, to_float + +if TYPE_CHECKING: + from dascore.core.attrs import PatchAttrs + +# CF-compliant time reference (Unix epoch is standard) +CF_TIME_REFERENCE = "seconds since 1970-01-01 00:00:00" +CF_CALENDAR = "proleptic_gregorian" + +# CF standard names for DAS data types +CF_STANDARD_NAMES = { + "strain": "strain", + "strain_rate": "strain_rate", + "velocity": "velocity", + "acceleration": "acceleration", + "temperature": "air_temperature", + "pressure": "air_pressure", + "acoustic": "acoustic_signal", +} + + +def datetime64_to_cf_time(dt_array: np.ndarray) -> np.ndarray: + """Convert datetime64 array to CF time (seconds since 1970-01-01).""" + return to_float(dt_array) + + +def cf_time_to_datetime64(time_array: np.ndarray, units: str) -> np.ndarray: + """Convert CF time array with units string to datetime64[ns].""" + if " since " not in units: + msg = f"Invalid CF time units format: {units}" + raise ValueError(msg) + time_unit, ref_str = units.split(" since ", 1) + time_unit = time_unit.strip().lower() + unit_to_seconds = { + "days": 86400.0, + "hours": 3600.0, + "minutes": 60.0, + "seconds": 1.0, + "milliseconds": 1e-3, + "microseconds": 1e-6, + } + if time_unit not in unit_to_seconds: + msg = f"Unsupported time unit: {time_unit}" + raise ValueError(msg) + ref_epoch = to_float(to_datetime64(ref_str.strip())) + seconds_from_epoch = ( + np.asarray(time_array, dtype=np.float64) * unit_to_seconds[time_unit] + + ref_epoch + ) + return to_datetime64(seconds_from_epoch) + + +def get_cf_data_attrs(data_type: str = "acoustic_signal") -> dict[str, str | float]: + """Get CF-compliant attributes for a DAS data variable.""" + attrs = { + "long_name": "Distributed Acoustic Sensing data", + "_FillValue": np.nan, + } + data_type_lower = data_type.lower() if data_type else "acoustic" + if data_type_lower in CF_STANDARD_NAMES: + attrs["standard_name"] = CF_STANDARD_NAMES[data_type_lower] + else: + for key, std_name in CF_STANDARD_NAMES.items(): + if key in data_type_lower: + attrs["standard_name"] = std_name + break + else: + attrs["standard_name"] = "acoustic_signal" + if "strain_rate" in data_type_lower: + attrs["units"] = "1/s" + attrs["long_name"] = "Strain rate" + elif "strain" in data_type_lower: + attrs["units"] = "1" + attrs["long_name"] = "Strain" + elif "velocity" in data_type_lower: + attrs["units"] = "m/s" + attrs["long_name"] = "Velocity" + elif "temperature" in data_type_lower: + attrs["units"] = "K" + attrs["long_name"] = "Temperature" + elif "pressure" in data_type_lower: + attrs["units"] = "Pa" + attrs["long_name"] = "Pressure" + else: + attrs["units"] = "1" + return attrs + + +def get_cf_global_attrs( + patch_attrs: PatchAttrs, cf_version: str = "1.8" +) -> dict[str, str]: + """Get CF-compliant global attributes from PatchAttrs.""" + now = datetime.datetime.now(datetime.timezone.utc) + attrs = { + "Conventions": f"CF-{cf_version}", + "title": "DAS data from DASCore", + "source": f"DASCore v{dc.__version__}", + "history": f"{now.isoformat()}: Created by DASCore", + "references": "https://dascore.org", + "comment": "Distributed Acoustic Sensing data", + "date_created": now.isoformat(), + } + if patch_attrs.station: + attrs["station"] = patch_attrs.station + if patch_attrs.network: + attrs["network"] = patch_attrs.network + if patch_attrs.instrument_id: + attrs["instrument"] = patch_attrs.instrument_id + if patch_attrs.acquisition_id: + attrs["acquisition"] = patch_attrs.acquisition_id + if patch_attrs.tag: + attrs["tag"] = patch_attrs.tag + if patch_attrs.data_category: + attrs["data_category"] = patch_attrs.data_category + if hasattr(patch_attrs, "category") and patch_attrs.category: + attrs["category"] = patch_attrs.category + if patch_attrs.data_type: + attrs["data_type"] = patch_attrs.data_type + if hasattr(patch_attrs, "history") and patch_attrs.history: + attrs["processing_history"] = " | ".join(str(h) for h in patch_attrs.history) + return attrs From 71e23c9f81e312978f2e9d75391f3a93362737d7 Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Thu, 9 Apr 2026 13:44:44 +0200 Subject: [PATCH 3/5] pre-refactor3 --- dascore/io/netcdf/__init__.py | 6 + dascore/io/netcdf/cf.py | 134 ----------- dascore/io/netcdf/core.py | 8 +- dascore/io/netcdf/utils.py | 423 +++++++++------------------------- 4 files changed, 120 insertions(+), 451 deletions(-) delete mode 100644 dascore/io/netcdf/cf.py diff --git a/dascore/io/netcdf/__init__.py b/dascore/io/netcdf/__init__.py index f5c3a7a48..17e96d117 100644 --- a/dascore/io/netcdf/__init__.py +++ b/dascore/io/netcdf/__init__.py @@ -3,3 +3,9 @@ from __future__ import annotations from dascore.io.netcdf.core import NetCDFCFV18 +from dascore.io.netcdf.utils import ( + cf_time_to_datetime64, + datetime64_to_cf_time, + get_cf_data_attrs, + get_cf_global_attrs, +) diff --git a/dascore/io/netcdf/cf.py b/dascore/io/netcdf/cf.py deleted file mode 100644 index 6a114833c..000000000 --- a/dascore/io/netcdf/cf.py +++ /dev/null @@ -1,134 +0,0 @@ -"""CF metadata helpers for NetCDF IO. - -See https://cfconventions.org/ for the Climate and Forecast metadata standard. -""" - -from __future__ import annotations - -import datetime -from typing import TYPE_CHECKING - -import numpy as np - -import dascore as dc -from dascore.utils.time import to_datetime64, to_float - -if TYPE_CHECKING: - from dascore.core.attrs import PatchAttrs - -# CF-compliant time reference (Unix epoch is standard) -CF_TIME_REFERENCE = "seconds since 1970-01-01 00:00:00" -CF_CALENDAR = "proleptic_gregorian" - -# CF standard names for DAS data types -CF_STANDARD_NAMES = { - "strain": "strain", - "strain_rate": "strain_rate", - "velocity": "velocity", - "acceleration": "acceleration", - "temperature": "air_temperature", - "pressure": "air_pressure", - "acoustic": "acoustic_signal", -} - - -def datetime64_to_cf_time(dt_array: np.ndarray) -> np.ndarray: - """Convert datetime64 array to CF time (seconds since 1970-01-01).""" - return to_float(dt_array) - - -def cf_time_to_datetime64(time_array: np.ndarray, units: str) -> np.ndarray: - """Convert CF time array with units string to datetime64[ns].""" - if " since " not in units: - msg = f"Invalid CF time units format: {units}" - raise ValueError(msg) - time_unit, ref_str = units.split(" since ", 1) - time_unit = time_unit.strip().lower() - unit_to_seconds = { - "days": 86400.0, - "hours": 3600.0, - "minutes": 60.0, - "seconds": 1.0, - "milliseconds": 1e-3, - "microseconds": 1e-6, - } - if time_unit not in unit_to_seconds: - msg = f"Unsupported time unit: {time_unit}" - raise ValueError(msg) - ref_epoch = to_float(to_datetime64(ref_str.strip())) - seconds_from_epoch = ( - np.asarray(time_array, dtype=np.float64) * unit_to_seconds[time_unit] - + ref_epoch - ) - return to_datetime64(seconds_from_epoch) - - -def get_cf_data_attrs(data_type: str = "acoustic_signal") -> dict[str, str | float]: - """Get CF-compliant attributes for a DAS data variable.""" - attrs = { - "long_name": "Distributed Acoustic Sensing data", - "_FillValue": np.nan, - } - data_type_lower = data_type.lower() if data_type else "acoustic" - if data_type_lower in CF_STANDARD_NAMES: - attrs["standard_name"] = CF_STANDARD_NAMES[data_type_lower] - else: - for key, std_name in CF_STANDARD_NAMES.items(): - if key in data_type_lower: - attrs["standard_name"] = std_name - break - else: - attrs["standard_name"] = "acoustic_signal" - if "strain_rate" in data_type_lower: - attrs["units"] = "1/s" - attrs["long_name"] = "Strain rate" - elif "strain" in data_type_lower: - attrs["units"] = "1" - attrs["long_name"] = "Strain" - elif "velocity" in data_type_lower: - attrs["units"] = "m/s" - attrs["long_name"] = "Velocity" - elif "temperature" in data_type_lower: - attrs["units"] = "K" - attrs["long_name"] = "Temperature" - elif "pressure" in data_type_lower: - attrs["units"] = "Pa" - attrs["long_name"] = "Pressure" - else: - attrs["units"] = "1" - return attrs - - -def get_cf_global_attrs( - patch_attrs: PatchAttrs, cf_version: str = "1.8" -) -> dict[str, str]: - """Get CF-compliant global attributes from PatchAttrs.""" - now = datetime.datetime.now(datetime.timezone.utc) - attrs = { - "Conventions": f"CF-{cf_version}", - "title": "DAS data from DASCore", - "source": f"DASCore v{dc.__version__}", - "history": f"{now.isoformat()}: Created by DASCore", - "references": "https://dascore.org", - "comment": "Distributed Acoustic Sensing data", - "date_created": now.isoformat(), - } - if patch_attrs.station: - attrs["station"] = patch_attrs.station - if patch_attrs.network: - attrs["network"] = patch_attrs.network - if patch_attrs.instrument_id: - attrs["instrument"] = patch_attrs.instrument_id - if patch_attrs.acquisition_id: - attrs["acquisition"] = patch_attrs.acquisition_id - if patch_attrs.tag: - attrs["tag"] = patch_attrs.tag - if patch_attrs.data_category: - attrs["data_category"] = patch_attrs.data_category - if hasattr(patch_attrs, "category") and patch_attrs.category: - attrs["category"] = patch_attrs.category - if patch_attrs.data_type: - attrs["data_type"] = patch_attrs.data_type - if hasattr(patch_attrs, "history") and patch_attrs.history: - attrs["processing_history"] = " | ".join(str(h) for h in patch_attrs.history) - return attrs diff --git a/dascore/io/netcdf/core.py b/dascore/io/netcdf/core.py index 7da3f4e42..c68ccd405 100644 --- a/dascore/io/netcdf/core.py +++ b/dascore/io/netcdf/core.py @@ -15,19 +15,19 @@ from dascore.utils.misc import optional_import from .utils import ( + XDAS_PAYLOAD_VARIABLE, + coord_attrs, extract_patch_attrs_from_netcdf, find_main_data_variable, - get_xarray_data_var_name, - get_xarray_engine, get_cf_data_attrs, get_cf_global_attrs, get_cf_version, + get_xarray_data_var_name, + get_xarray_engine, is_netcdf4_file, iter_written_aux_coords, parse_cf_version, read_netcdf_coordinates, - XDAS_PAYLOAD_VARIABLE, - coord_attrs, ) diff --git a/dascore/io/netcdf/utils.py b/dascore/io/netcdf/utils.py index 6e1d5680d..6148d915e 100644 --- a/dascore/io/netcdf/utils.py +++ b/dascore/io/netcdf/utils.py @@ -1,4 +1,4 @@ -"""Utilities for NetCDF IO with CF conventions support. +"""NetCDF helper functions for DASCore IO. See https://cfconventions.org/ for the Climate and Forecast metadata standard. """ @@ -24,7 +24,6 @@ # CF-compliant time reference (Unix epoch is standard) CF_TIME_REFERENCE = "seconds since 1970-01-01 00:00:00" CF_CALENDAR = "proleptic_gregorian" -XDAS_PAYLOAD_VARIABLE = "__values__" # CF standard names for DAS data types CF_STANDARD_NAMES = { @@ -37,142 +36,7 @@ "acoustic": "acoustic_signal", } - -def coord_attrs(name: str, coord) -> dict[str, str]: - """Return CF-ish coordinate attrs for xarray-backed NetCDF output.""" - # Keep time explicit because CF consumers treat it as a special semantic - # axis, not just another coordinate with datetime-like values. - if name == "time": - return { - "standard_name": "time", - "long_name": "Time", - "axis": "T", - } - return { - "long_name": name.replace("_", " ").title(), - "standard_name": name.lower(), - "units": str(coord.units.units) if coord.units else ("m" if "depth" in name else "1"), - } - - -def get_xarray_data_var_name(dataset) -> str: - """Return the main xarray data variable name.""" - if "data" in dataset.data_vars: - return "data" - # XDAS-style files can surface the primary payload under a None key while - # exposing coord helper arrays (for example *_indices/*_values) as data vars. - if None in dataset.data_vars: - return None - if len(dataset.data_vars) == 1: - return next(iter(dataset.data_vars)) - msg = "No suitable data variable found in NetCDF file" - raise ValueError(msg) - - -def parse_cf_version(cf_version: str) -> tuple[int, int]: - """Parse a CF version string into comparable major/minor integers.""" - parts = cf_version.split(".") - major = int(parts[0]) - minor = int(parts[1]) if len(parts) > 1 else 0 - return major, minor - - -@cache -def get_xarray_engine(on_missing: str = "raise") -> str | None: - """Return the preferred xarray engine for NetCDF-4 files.""" - # Map importable module names onto the engine strings xarray expects. - module_to_engine = { - "netCDF4": "netcdf4", - "h5netcdf": "h5netcdf", - } - for module_name, engine_name in module_to_engine.items(): - mod = optional_import(module_name, on_missing="ignore") - if mod is not None: - return engine_name - # If no backend available ignore or raise. - if on_missing == "ignore": - return None - msg = ( - "Either netCDF4 or h5netcdf is required for NetCDF-4 read/write " - "functionality." - ) - raise MissingOptionalDependencyError(msg) - - -def iter_written_aux_coords(patch: dc.Patch): - """Yield names of auxiliary coordinates serialized to NetCDF.""" - for name, coord in patch.coords.coord_map.items(): - if coord._partial or name in patch.dims: - continue - yield name - - -def is_netcdf4_file(h5file: h5py.File) -> bool: - """ - Check if an HDF5 file follows NetCDF-4 conventions. - - Parameters - ---------- - h5file - Open h5py.File object - - Returns - ------- - bool - True if file appears to be NetCDF-4 format - - Notes - ----- - NetCDF-4 files are identified by: - - _NCProperties attribute (NetCDF-4 specific) - - Conventions attribute starting with "CF-" - """ - try: - # Check for NetCDF-specific markers - if "_NCProperties" in h5file.attrs: - return True - - # Check for Conventions attribute indicating CF compliance - conventions = h5file.attrs.get("Conventions", "") - if isinstance(conventions, bytes): - conventions = conventions.decode("utf-8", errors="ignore") - if conventions and "CF" in conventions: - return True - - return False - - except (AttributeError, KeyError): - return False - - -def get_cf_version(h5file: h5py.File) -> str | None: - """ - Extract CF convention version from NetCDF file. - - Parameters - ---------- - h5file - Open h5py.File object - - Returns - ------- - str | None - CF version string (e.g., "1.8") or None if not found - """ - conventions = h5file.attrs.get("Conventions", "") - if isinstance(conventions, bytes): - conventions = conventions.decode("utf-8", errors="ignore") - - # Handle various CF convention formats - if "CF-" in conventions: - # Format: "CF-1.8" or "CF-1.8, ACDD-1.3" - parts = conventions.split("CF-", 1)[1].split()[0].rstrip(",;") - return parts - elif conventions.startswith("CF "): - # Format: "CF 1.8" - return conventions.split()[1].rstrip(",;") - - return None +XDAS_PAYLOAD_VARIABLE = "__values__" def datetime64_to_cf_time(dt_array: np.ndarray) -> np.ndarray: @@ -181,15 +45,7 @@ def datetime64_to_cf_time(dt_array: np.ndarray) -> np.ndarray: def cf_time_to_datetime64(time_array: np.ndarray, units: str) -> np.ndarray: - """Convert CF time array with units string to datetime64[ns]. - - Parameters - ---------- - time_array - Numeric time values. - units - CF units string, e.g. "seconds since 1970-01-01 00:00:00". - """ + """Convert CF time array with units string to datetime64[ns].""" if " since " not in units: msg = f"Invalid CF time units format: {units}" raise ValueError(msg) @@ -206,7 +62,6 @@ def cf_time_to_datetime64(time_array: np.ndarray, units: str) -> np.ndarray: if time_unit not in unit_to_seconds: msg = f"Unsupported time unit: {time_unit}" raise ValueError(msg) - # Convert to seconds since 1970 then to datetime64 ref_epoch = to_float(to_datetime64(ref_str.strip())) seconds_from_epoch = ( np.asarray(time_array, dtype=np.float64) * unit_to_seconds[time_unit] @@ -216,41 +71,21 @@ def cf_time_to_datetime64(time_array: np.ndarray, units: str) -> np.ndarray: def get_cf_data_attrs(data_type: str = "acoustic_signal") -> dict[str, str | float]: - """ - Get CF-compliant attributes for DAS data variable. - - Parameters - ---------- - data_type - Type of DAS data - - Returns - ------- - dict - CF-compliant data attributes - """ + """Get CF-compliant attributes for a DAS data variable.""" attrs = { "long_name": "Distributed Acoustic Sensing data", "_FillValue": np.nan, } - - # Map data type to CF standard name data_type_lower = data_type.lower() if data_type else "acoustic" - - # Check for exact match first if data_type_lower in CF_STANDARD_NAMES: attrs["standard_name"] = CF_STANDARD_NAMES[data_type_lower] else: - # Check for partial matches for key, std_name in CF_STANDARD_NAMES.items(): if key in data_type_lower: attrs["standard_name"] = std_name break else: - # Default to generic acoustic signal attrs["standard_name"] = "acoustic_signal" - - # Add units based on data type if "strain_rate" in data_type_lower: attrs["units"] = "1/s" attrs["long_name"] = "Strain rate" @@ -267,30 +102,14 @@ def get_cf_data_attrs(data_type: str = "acoustic_signal") -> dict[str, str | flo attrs["units"] = "Pa" attrs["long_name"] = "Pressure" else: - # Generic units for acoustic data attrs["units"] = "1" - return attrs def get_cf_global_attrs( patch_attrs: PatchAttrs, cf_version: str = "1.8" ) -> dict[str, str]: - """ - Get CF-compliant global attributes from PatchAttrs. - - Parameters - ---------- - patch_attrs - DASCore patch attributes - cf_version - CF convention version - - Returns - ------- - dict - CF-compliant global attributes - """ + """Get CF-compliant global attributes from PatchAttrs.""" now = datetime.datetime.now(datetime.timezone.utc) attrs = { "Conventions": f"CF-{cf_version}", @@ -301,8 +120,6 @@ def get_cf_global_attrs( "comment": "Distributed Acoustic Sensing data", "date_created": now.isoformat(), } - - # Add optional attributes from patch if patch_attrs.station: attrs["station"] = patch_attrs.station if patch_attrs.network: @@ -313,40 +130,113 @@ def get_cf_global_attrs( attrs["acquisition"] = patch_attrs.acquisition_id if patch_attrs.tag: attrs["tag"] = patch_attrs.tag - - # Add data category/type info if patch_attrs.data_category: attrs["data_category"] = patch_attrs.data_category if hasattr(patch_attrs, "category") and patch_attrs.category: attrs["category"] = patch_attrs.category if patch_attrs.data_type: attrs["data_type"] = patch_attrs.data_type - - # Add processing history if available if hasattr(patch_attrs, "history") and patch_attrs.history: - history_str = " | ".join(str(h) for h in patch_attrs.history) - attrs["processing_history"] = history_str - + attrs["processing_history"] = " | ".join(str(h) for h in patch_attrs.history) return attrs +def coord_attrs(name: str, coord) -> dict[str, str]: + """Return CF-ish coordinate attrs for xarray-backed NetCDF output.""" + # Keep time explicit because CF consumers treat it as a special semantic + # axis, not just another coordinate with datetime-like values. + if name == "time": + return { + "standard_name": "time", + "long_name": "Time", + "axis": "T", + } + return { + "long_name": name.replace("_", " ").title(), + "standard_name": name.lower(), + "units": str(coord.units.units) if coord.units else ("m" if "depth" in name else "1"), + } + + +def get_xarray_data_var_name(dataset) -> str: + """Return the main xarray data variable name.""" + if "data" in dataset.data_vars: + return "data" + # XDAS-style files can surface the primary payload under a None key while + # exposing coord helper arrays (for example *_indices/*_values) as data vars. + if None in dataset.data_vars: + return None + if len(dataset.data_vars) == 1: + return next(iter(dataset.data_vars)) + msg = "No suitable data variable found in NetCDF file" + raise ValueError(msg) + + +@cache +def get_xarray_engine(on_missing: str = "raise") -> str | None: + """Return the preferred xarray engine for NetCDF-4 files.""" + # Map importable module names onto the engine strings xarray expects. + module_to_engine = { + "netCDF4": "netcdf4", + "h5netcdf": "h5netcdf", + } + for module_name, engine_name in module_to_engine.items(): + mod = optional_import(module_name, on_missing="ignore") + if mod is not None: + return engine_name + if on_missing == "ignore": + return None + msg = ( + "Either netCDF4 or h5netcdf is required for NetCDF-4 read/write " + "functionality." + ) + raise MissingOptionalDependencyError(msg) + + +def iter_written_aux_coords(patch: dc.Patch): + """Yield names of auxiliary coordinates serialized to NetCDF.""" + for name, coord in patch.coords.coord_map.items(): + if coord._partial or name in patch.dims: + continue + yield name + + +def parse_cf_version(cf_version: str) -> tuple[int, int]: + """Parse a CF version string into comparable major/minor integers.""" + parts = cf_version.split(".") + major = int(parts[0]) + minor = int(parts[1]) if len(parts) > 1 else 0 + return major, minor + + +def is_netcdf4_file(h5file: h5py.File) -> bool: + """Return True when an HDF5 file exposes strong NetCDF/CF markers.""" + try: + if "_NCProperties" in h5file.attrs: + return True + conventions = h5file.attrs.get("Conventions", "") + if isinstance(conventions, bytes): + conventions = conventions.decode("utf-8", errors="ignore") + return bool(conventions and "CF" in conventions) + except (AttributeError, KeyError): + return False + + +def get_cf_version(h5file: h5py.File) -> str | None: + """Extract the CF convention version string from a NetCDF file.""" + conventions = h5file.attrs.get("Conventions", "") + if isinstance(conventions, bytes): + conventions = conventions.decode("utf-8", errors="ignore") + if "CF-" in conventions: + return conventions.split("CF-", 1)[1].split()[0].rstrip(",;") + if conventions.startswith("CF "): + return conventions.split()[1].rstrip(",;") + return None + + def extract_patch_attrs_from_netcdf(h5file: h5py.File) -> dict: - """ - Extract patch attributes from NetCDF global attributes. - - Parameters - ---------- - h5file - Open h5py.File object - - Returns - ------- - dict - Dictionary of patch attributes - """ + """Extract DASCore patch attrs from NetCDF global attributes.""" attrs = {} - - # Map CF global attributes to patch attributes attr_mapping = { "station": "station", "network": "network", @@ -356,13 +246,12 @@ def extract_patch_attrs_from_netcdf(h5file: h5py.File) -> dict: "data_category": "data_category", "category": "category", } - for cf_name, patch_name in attr_mapping.items(): if cf_name in h5file.attrs: value = h5file.attrs[cf_name] if isinstance(value, bytes): value = value.decode("utf-8", errors="ignore") - if value: # Only add non-empty values + if value: attrs[patch_name] = value for cf_name in ("data_type", "source_data_type"): if cf_name not in h5file.attrs: @@ -372,64 +261,32 @@ def extract_patch_attrs_from_netcdf(h5file: h5py.File) -> dict: value = value.decode("utf-8", errors="ignore") if value: attrs.setdefault("data_type", value) - return attrs def _handle_time_interpolation( h5file: h5py.File, coord_name: str, coord_data: np.ndarray ) -> np.ndarray | None: - """ - Handle coordinate interpolation for time coordinates. - - Parameters - ---------- - h5file - Open h5py.File object - coord_name - Name of the coordinate (e.g., "time") - coord_data - Raw coordinate data array - - Returns - ------- - np.ndarray | None - Converted datetime64 array if interpolation data found, None otherwise - """ - # Compatibility path for xdas-produced NetCDF files that store time as CF - # tie points in `_values` plus optional `_indices`. This is - # not generic CF behavior and does not require an xdas dependency. - # Look for time_values with proper CF units + """Decode XDAS-style time tie points into a full datetime coordinate.""" time_values_name = f"{coord_name}_values" if time_values_name not in h5file: return None - time_values_var = h5file[time_values_name] units = time_values_var.attrs.get("units", "") if isinstance(units, bytes): units = units.decode("utf-8", errors="ignore") - if "since" not in units: return None - - # Get the time values and interpolate to full coordinate array time_values = time_values_var[:] time_indices_name = f"{coord_name}_indices" - if time_indices_name in h5file: - # Use indices for interpolation time_indices = h5file[time_indices_name][:] if len(time_values) >= 2 and len(time_indices) >= 2: - # Linear interpolation from tie points to full coordinate interpolated_values = np.interp(coord_data, time_indices, time_values) else: - # Fall back to using time_values directly interpolated_values = time_values else: - # Use time_values directly interpolated_values = time_values - - # Convert to datetime64 try: return cf_time_to_datetime64(interpolated_values, units) except (ValueError, KeyError): @@ -440,43 +297,31 @@ def read_netcdf_coordinates( h5file: h5py.File, data_var_name: str | None = None ) -> CoordManager: """Read coordinate information from a NetCDF file into a CoordManager.""" - coords = {} coord_order = [] - - # Use provided data variable name, or discover it main_data_var = data_var_name or find_main_data_variable(h5file) expected_shape = None if main_data_var and main_data_var in h5file: data_var = h5file[main_data_var] expected_shape = data_var.shape - - # Try to get dimension order from DIMENSION_LIST if "DIMENSION_LIST" in data_var.attrs: dim_list = data_var.attrs["DIMENSION_LIST"] for ref_array in dim_list: try: ref = ref_array[0] dim_scale = h5file[ref] - dim_name = dim_scale.name.strip("/") - coord_order.append(dim_name) + coord_order.append(dim_scale.name.strip("/")) except (IndexError, KeyError, TypeError): - # If we can't resolve reference, we'll fall back to discovery pass - # Collect dimension-scale coordinates dim_coords = {} for name, dataset in h5file.items(): if not isinstance(dataset, h5py.Dataset): continue if not (dataset.is_scale or "NAME" in dataset.attrs): continue - # Skip auxiliaries that don't match any data dimension size if expected_shape and dataset.shape[0] not in expected_shape: continue - if any( - skip in name.lower() - for skip in ("_points", "_indices", "_values", "_interpolation") - ): + if any(skip in name.lower() for skip in ("_points", "_indices", "_values", "_interpolation")): continue data = dataset[:] if name == "time" or dataset.attrs.get("axis") == "T": @@ -485,8 +330,7 @@ def read_netcdf_coordinates( units = units.decode("utf-8", errors="ignore") if "since" in units: try: - data = cf_time_to_datetime64(data, units) - dim_coords["time"] = data + dim_coords["time"] = cf_time_to_datetime64(data, units) continue except (ValueError, KeyError): pass @@ -496,7 +340,6 @@ def read_netcdf_coordinates( continue dim_coords[name] = data - # Build coords dict preserving dimension order from DIMENSION_LIST coords = {} if coord_order: for dim_name in coord_order: @@ -508,7 +351,6 @@ def read_netcdf_coordinates( else: coords = dim_coords - # Collect non-dimension coordinates stored with _DASCORE_DIMS attribute if main_data_var and main_data_var in h5file: coord_attr = h5file[main_data_var].attrs.get("coordinates", "") if isinstance(coord_attr, bytes): @@ -523,74 +365,35 @@ def read_netcdf_coordinates( if isinstance(dims_attr, bytes): dims_attr = dims_attr.decode("utf-8", errors="ignore") if dims_attr: - associated_dims = tuple(dims_attr.split(",")) - coords[name] = (associated_dims, dataset[:]) + coords[name] = (tuple(dims_attr.split(",")), dataset[:]) return dc.core.coordmanager.get_coord_manager(coords) def validate_cf_compliance(h5file: h5py.File) -> list[str]: - """ - Validate CF compliance and return list of issues. - - Parameters - ---------- - h5file - Open h5py.File object - - Returns - ------- - list[str] - List of CF compliance issues found - """ + """Validate CF compliance and return a list of issues.""" issues = [] - - # Check for required global attributes if "Conventions" not in h5file.attrs: issues.append("Missing required 'Conventions' global attribute") - - # Check coordinate variables for name, dataset in h5file.items(): if isinstance(dataset, h5py.Dataset) and dataset.is_scale: - # Check for required coordinate attributes if "units" not in dataset.attrs: issues.append(f"Coordinate '{name}' missing 'units' attribute") - - # Check time coordinate if name == "time" or dataset.attrs.get("standard_name") == "time": units = dataset.attrs.get("units", "") if not units or "since" not in str(units): issues.append(f"Time coordinate '{name}' has invalid units") - - # Check data variables for name, dataset in h5file.items(): if isinstance(dataset, h5py.Dataset) and not dataset.is_scale: - # Check for required data variable attributes if "units" not in dataset.attrs: issues.append(f"Data variable '{name}' missing 'units' attribute") if "long_name" not in dataset.attrs: issues.append(f"Data variable '{name}' missing 'long_name' attribute") - return issues def find_main_data_variable(h5file: h5py.File) -> str | None: - """ - Find the main data variable in the NetCDF file. - - Looks for 2D+ datasets that are not dimension scales. - Prioritizes variables with standard names suggesting DAS data. - - Parameters - ---------- - h5file - Open h5py.File object - - Returns - ------- - str | None - Name of the main data variable, or None if not found - """ + """Find the main data variable in a NetCDF file.""" priority_names = [ "data", "acoustic_data", @@ -600,28 +403,22 @@ def find_main_data_variable(h5file: h5py.File) -> str | None: "velocity", "amplitude", ] - candidates = [] - for name, item in h5file.items(): if not _is_data_variable_candidate(item): continue - - # Check for priority matches first if name in priority_names or _has_priority_standard_name(item, priority_names): return name - candidates.append(name) - return candidates[0] if candidates else None def _is_data_variable_candidate(item) -> bool: - """Check if item is a candidate for main data variable.""" + """Check if an item is a candidate for the main data variable.""" return isinstance(item, h5py.Dataset) and not item.is_scale and item.ndim >= 2 def _has_priority_standard_name(item, priority_names: list[str]) -> bool: - """Check if dataset has a standard_name matching priority names.""" + """Check if a dataset standard_name matches a known priority name.""" std_name = str(item.attrs.get("standard_name", "")).lower() return any(name in std_name for name in priority_names) From 6f53dee422c8cd039e466cc3061035b447c5b3e2 Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Thu, 9 Apr 2026 17:49:29 +0200 Subject: [PATCH 4/5] review --- dascore/io/netcdf/__init__.py | 8 +- dascore/io/netcdf/core.py | 207 +++---- dascore/io/netcdf/utils.py | 407 ++----------- dascore/utils/io.py | 3 +- tests/test_io/test_common_io.py | 8 +- tests/test_io/test_netcdf/test_netcdf.py | 739 ++--------------------- 6 files changed, 177 insertions(+), 1195 deletions(-) diff --git a/dascore/io/netcdf/__init__.py b/dascore/io/netcdf/__init__.py index 17e96d117..d726899cb 100644 --- a/dascore/io/netcdf/__init__.py +++ b/dascore/io/netcdf/__init__.py @@ -1,11 +1,5 @@ -"""NetCDF IO support for DASCore using CF (Climate and Forecast) conventions.""" +"""NetCDF IO support for DASCore.""" from __future__ import annotations from dascore.io.netcdf.core import NetCDFCFV18 -from dascore.io.netcdf.utils import ( - cf_time_to_datetime64, - datetime64_to_cf_time, - get_cf_data_attrs, - get_cf_global_attrs, -) diff --git a/dascore/io/netcdf/core.py b/dascore/io/netcdf/core.py index c68ccd405..82dc3425b 100644 --- a/dascore/io/netcdf/core.py +++ b/dascore/io/netcdf/core.py @@ -1,38 +1,29 @@ -"""Core NetCDF IO implementation with CF conventions.""" +"""Core NetCDF IO implementation built on xarray.""" from __future__ import annotations from pathlib import Path -import h5py - import dascore as dc from dascore.constants import SpoolType from dascore.io import FiberIO from dascore.io.core import ScanPayload, _make_scan_payload from dascore.utils.hdf5 import H5Reader -from dascore.utils.io import patch_to_xarray +from dascore.utils.io import patch_to_xarray, xarray_to_patch from dascore.utils.misc import optional_import from .utils import ( XDAS_PAYLOAD_VARIABLE, - coord_attrs, - extract_patch_attrs_from_netcdf, - find_main_data_variable, - get_cf_data_attrs, - get_cf_global_attrs, get_cf_version, + get_coord_manager_for_coordless_data_var, get_xarray_data_var_name, - get_xarray_engine, is_netcdf4_file, - iter_written_aux_coords, parse_cf_version, - read_netcdf_coordinates, ) class NetCDFCFV18(FiberIO): - """NetCDF-4 IO with CF-1.8 conventions, using xarray/netcdf4 for IO.""" + """NetCDF-4 IO using xarray for read/write and CF markers for detection.""" name = "NETCDF_CF" version = "1.8" @@ -55,102 +46,15 @@ def get_format(self, resource: H5Reader, **kwargs) -> tuple[str, str] | bool: def read(self, resource: Path, **kwargs) -> SpoolType: """Read a NetCDF-4 file into a Spool.""" xr = optional_import("xarray") - engine = get_xarray_engine() - # Read structural metadata first so xarray only has to resolve the final - # payload variable and dense coordinate values. - data_var_name, attrs_dict, coords = self._read_metadata(resource) - data_var_name, data_array, data = self._read_data_array( - resource, xr, engine, data_var_name - ) - patch = self._build_patch( - data=data, - data_array=data_array, - coords=coords, - attrs_dict=attrs_dict, - data_var_name=data_var_name, - ) - patch = self._apply_coordinate_filtering(patch, kwargs) + with xr.open_dataset(resource) as dataset: + data_var_name = get_xarray_data_var_name(dataset) + data_array = dataset[data_var_name].load() + patch = self._patch_from_dataset(dataset, data_var_name, data_array) + patch = self._select_from_kwargs(patch, kwargs) if not patch.data.size: return dc.spool([]) return dc.spool([patch]) - def _read_metadata(self, resource: Path): - """Read NetCDF metadata needed before opening through xarray.""" - with h5py.File(resource, "r") as h5file: - data_var_name = self._get_data_variable_name(h5file) - attrs_dict = self._get_patch_attrs(h5file) - coords = read_netcdf_coordinates(h5file, data_var_name) - return data_var_name, attrs_dict, coords - - def _read_data_array(self, resource: Path, xr, engine: str, data_var_name: str): - """Load the selected xarray data variable from disk.""" - # TODO: consider a lazy read path for large NetCDF arrays. - with xr.open_dataset(resource, engine=engine) as dataset: - data_array = dataset.get(data_var_name) - if data_array is None: - data_var_name = get_xarray_data_var_name(dataset) - data_array = dataset[data_var_name] - data = data_array.load().data - return data_var_name, data_array, data - - def _build_patch(self, *, data, data_array, coords, attrs_dict, data_var_name: str): - """Merge coordinate metadata and construct the output patch.""" - source_patch_id = ( - XDAS_PAYLOAD_VARIABLE if data_var_name is None else data_var_name - ) - # Start from the HDF-derived coordinates, then add any extra xarray-only - # coordinates that were materialized during decode. - coords_dict = { - name: ( - coord.values - if name in coords.dims - else (coords.dim_map[name], coord.values) - ) - for name, coord in coords.coord_map.items() - } - for name, coord in data_array.coords.items(): - if name not in coords_dict: - coords_dict[name] = (coord.dims, coord.values) - coords = dc.get_coord_manager(coords=coords_dict, dims=coords.dims) - return dc.Patch( - data=data, - coords=coords, - dims=coords.dims, - attrs=attrs_dict | {"_source_patch_id": source_patch_id}, - ) - - def _build_data_array(self, patch: dc.Patch): - """Convert a patch to an xarray DataArray with coordinate attrs.""" - data_array = patch_to_xarray(patch).rename("data") - for name, coord in patch.coords.coord_map.items(): - if coord._partial: - continue - data_array.coords[name].attrs.update(coord_attrs(name, coord)) - return data_array - - def _build_dataset(self, patch: dc.Patch, data_array): - """Create the xarray Dataset and attach CF metadata.""" - global_attrs = get_cf_global_attrs(patch.attrs, self.version) - if patch.attrs.data_type: - global_attrs["source_data_type"] = patch.attrs.data_type - dataset = data_array.to_dataset() - # NetCDF attrs cannot safely preserve DASCore's null-ish metadata, so - # filter those out before handing the dataset to xarray. - dataset.attrs.update( - { - key: value - for key, value in global_attrs.items() - if value not in (None, "") - } - ) - dataset["data"].attrs.update( - get_cf_data_attrs(patch.attrs.data_type or "acoustic_signal") - ) - aux_coord_names = tuple(iter_written_aux_coords(patch)) - if aux_coord_names: - dataset["data"].attrs["coordinates"] = " ".join(aux_coord_names) - return dataset - def _get_write_encoding(self, **kwargs): """Translate explicit write options into xarray encoding hints.""" compression = kwargs.get("compression") @@ -169,7 +73,7 @@ def _get_write_encoding(self, **kwargs): def write(self, spool: SpoolType, resource: Path, **kwargs) -> None: """ - Write a Spool to NetCDF-4 with CF-1.8 conventions. + Write a Spool to NetCDF-4 through xarray. Parameters ---------- @@ -180,50 +84,79 @@ def write(self, spool: SpoolType, resource: Path, **kwargs) -> None: explicit tuple of chunk sizes """ patch = self._validate_and_extract_patch(spool) - optional_import("xarray") - engine = get_xarray_engine() - # Build the xarray object first, then attach CF metadata and optional - # storage hints in separate steps to keep write policy isolated. - data_array = self._build_data_array(patch) - dataset = self._build_dataset(patch, data_array) + optional_import("xarray") # raises a helpful error if xarray is absent + dataset = patch_to_xarray(patch).rename("data").to_dataset() + dataset.attrs["Conventions"] = f"CF-{self.version}" encoding = self._get_write_encoding(**kwargs) dataset.to_netcdf( resource, - engine=engine, encoding={"data": encoding} if encoding else None, ) def scan(self, resource: H5Reader, **kwargs) -> list[ScanPayload]: - """Scan NetCDF file to extract metadata without loading data.""" - data_var_name = self._get_data_variable_name(resource) - data_var = resource[data_var_name] - coords = read_netcdf_coordinates(resource, data_var_name) - attrs_dict = self._get_patch_attrs(resource) + """Scan NetCDF file metadata without loading the full payload array.""" + xr = optional_import("xarray") + dataset_path = resource.filename + with xr.open_dataset(dataset_path) as dataset: + data_var_name = get_xarray_data_var_name(dataset) + # None is a valid xarray key for XDAS-style files whose primary + # payload is stored under a None variable name. + data_array = dataset[data_var_name] + coords = { + name: (coord.dims, coord.values) + for name, coord in data_array.coords.items() + } + attrs = dict(data_array.attrs) + dims = data_array.dims + shape = data_array.shape + dtype = str(data_array.dtype) + source_patch_id = self._get_source_patch_id(data_var_name) + coord_manager = self._coord_manager_from_data_array( + dataset, data_array, coords, dims, shape + ) return [ _make_scan_payload( - attrs=attrs_dict, - coords=coords, - dims=coords.dims, - shape=data_var.shape, - dtype=str(data_var.dtype), - source_patch_id=data_var_name, + attrs=attrs | {"_source_patch_id": source_patch_id}, + coords=coord_manager, + dims=dims, + shape=shape, + dtype=dtype, + source_patch_id=source_patch_id, ) ] - def _get_data_variable_name(self, resource: H5Reader) -> str: - """Return the main NetCDF data variable name or raise.""" - data_var_name = find_main_data_variable(resource) - if data_var_name is None: - msg = "No suitable data variable found in NetCDF file" - raise ValueError(msg) - return data_var_name - - def _get_patch_attrs(self, resource: H5Reader) -> dict: - """Extract patch attrs from a NetCDF resource.""" - return extract_patch_attrs_from_netcdf(resource) + def _get_source_patch_id(self, data_var_name): + """Normalize the selected xarray payload name to a patch id.""" + return XDAS_PAYLOAD_VARIABLE if data_var_name is None else data_var_name + + def _coord_manager_from_data_array(self, dataset, data_array, coords, dims, shape): + """Return coords from xarray when present or reconstruct dim coords.""" + if coords: + return dc.get_coord_manager(coords=coords, dims=dims) + return get_coord_manager_for_coordless_data_var(dataset, dims=dims, shape=shape) + + def _patch_from_dataset(self, dataset, data_var_name, data_array): + """Build one patch from an xarray dataset and selected data variable.""" + source_patch_id = self._get_source_patch_id(data_var_name) + attrs = dict(data_array.attrs) | {"_source_patch_id": source_patch_id} + if data_array.coords: + return xarray_to_patch(data_array).update(attrs=attrs) + coords = self._coord_manager_from_data_array( + dataset, + data_array, + coords={}, + dims=data_array.dims, + shape=data_array.shape, + ) + return dc.Patch( + data=data_array.data, + coords=coords, + dims=data_array.dims, + attrs=attrs, + ) - def _apply_coordinate_filtering(self, patch: dc.Patch, kwargs: dict) -> dc.Patch: - """Apply coordinate selection kwargs to a loaded patch.""" + def _select_from_kwargs(self, patch: dc.Patch, kwargs: dict) -> dc.Patch: + """Apply coordinate selection kwargs to one loaded patch.""" coord_kwargs = {k: v for k, v in kwargs.items() if k in patch.coords.coord_map} return patch.select(**coord_kwargs) if coord_kwargs else patch diff --git a/dascore/io/netcdf/utils.py b/dascore/io/netcdf/utils.py index 6148d915e..733086aa2 100644 --- a/dascore/io/netcdf/utils.py +++ b/dascore/io/netcdf/utils.py @@ -1,169 +1,21 @@ -"""NetCDF helper functions for DASCore IO. - -See https://cfconventions.org/ for the Climate and Forecast metadata standard. -""" +"""NetCDF helper functions for DASCore IO.""" from __future__ import annotations -import datetime -from functools import cache -from typing import TYPE_CHECKING - import h5py import numpy as np import dascore as dc -from dascore.exceptions import MissingOptionalDependencyError -from dascore.utils.misc import optional_import -from dascore.utils.time import to_datetime64, to_float - -if TYPE_CHECKING: - from dascore.core.attrs import PatchAttrs - from dascore.core.coordmanager import CoordManager - -# CF-compliant time reference (Unix epoch is standard) -CF_TIME_REFERENCE = "seconds since 1970-01-01 00:00:00" -CF_CALENDAR = "proleptic_gregorian" - -# CF standard names for DAS data types -CF_STANDARD_NAMES = { - "strain": "strain", - "strain_rate": "strain_rate", - "velocity": "velocity", - "acceleration": "acceleration", - "temperature": "air_temperature", - "pressure": "air_pressure", - "acoustic": "acoustic_signal", -} XDAS_PAYLOAD_VARIABLE = "__values__" -def datetime64_to_cf_time(dt_array: np.ndarray) -> np.ndarray: - """Convert datetime64 array to CF time (seconds since 1970-01-01).""" - return to_float(dt_array) - - -def cf_time_to_datetime64(time_array: np.ndarray, units: str) -> np.ndarray: - """Convert CF time array with units string to datetime64[ns].""" - if " since " not in units: - msg = f"Invalid CF time units format: {units}" - raise ValueError(msg) - time_unit, ref_str = units.split(" since ", 1) - time_unit = time_unit.strip().lower() - unit_to_seconds = { - "days": 86400.0, - "hours": 3600.0, - "minutes": 60.0, - "seconds": 1.0, - "milliseconds": 1e-3, - "microseconds": 1e-6, - } - if time_unit not in unit_to_seconds: - msg = f"Unsupported time unit: {time_unit}" - raise ValueError(msg) - ref_epoch = to_float(to_datetime64(ref_str.strip())) - seconds_from_epoch = ( - np.asarray(time_array, dtype=np.float64) * unit_to_seconds[time_unit] - + ref_epoch - ) - return to_datetime64(seconds_from_epoch) - - -def get_cf_data_attrs(data_type: str = "acoustic_signal") -> dict[str, str | float]: - """Get CF-compliant attributes for a DAS data variable.""" - attrs = { - "long_name": "Distributed Acoustic Sensing data", - "_FillValue": np.nan, - } - data_type_lower = data_type.lower() if data_type else "acoustic" - if data_type_lower in CF_STANDARD_NAMES: - attrs["standard_name"] = CF_STANDARD_NAMES[data_type_lower] - else: - for key, std_name in CF_STANDARD_NAMES.items(): - if key in data_type_lower: - attrs["standard_name"] = std_name - break - else: - attrs["standard_name"] = "acoustic_signal" - if "strain_rate" in data_type_lower: - attrs["units"] = "1/s" - attrs["long_name"] = "Strain rate" - elif "strain" in data_type_lower: - attrs["units"] = "1" - attrs["long_name"] = "Strain" - elif "velocity" in data_type_lower: - attrs["units"] = "m/s" - attrs["long_name"] = "Velocity" - elif "temperature" in data_type_lower: - attrs["units"] = "K" - attrs["long_name"] = "Temperature" - elif "pressure" in data_type_lower: - attrs["units"] = "Pa" - attrs["long_name"] = "Pressure" - else: - attrs["units"] = "1" - return attrs - - -def get_cf_global_attrs( - patch_attrs: PatchAttrs, cf_version: str = "1.8" -) -> dict[str, str]: - """Get CF-compliant global attributes from PatchAttrs.""" - now = datetime.datetime.now(datetime.timezone.utc) - attrs = { - "Conventions": f"CF-{cf_version}", - "title": "DAS data from DASCore", - "source": f"DASCore v{dc.__version__}", - "history": f"{now.isoformat()}: Created by DASCore", - "references": "https://dascore.org", - "comment": "Distributed Acoustic Sensing data", - "date_created": now.isoformat(), - } - if patch_attrs.station: - attrs["station"] = patch_attrs.station - if patch_attrs.network: - attrs["network"] = patch_attrs.network - if patch_attrs.instrument_id: - attrs["instrument"] = patch_attrs.instrument_id - if patch_attrs.acquisition_id: - attrs["acquisition"] = patch_attrs.acquisition_id - if patch_attrs.tag: - attrs["tag"] = patch_attrs.tag - if patch_attrs.data_category: - attrs["data_category"] = patch_attrs.data_category - if hasattr(patch_attrs, "category") and patch_attrs.category: - attrs["category"] = patch_attrs.category - if patch_attrs.data_type: - attrs["data_type"] = patch_attrs.data_type - if hasattr(patch_attrs, "history") and patch_attrs.history: - attrs["processing_history"] = " | ".join(str(h) for h in patch_attrs.history) - return attrs - - -def coord_attrs(name: str, coord) -> dict[str, str]: - """Return CF-ish coordinate attrs for xarray-backed NetCDF output.""" - # Keep time explicit because CF consumers treat it as a special semantic - # axis, not just another coordinate with datetime-like values. - if name == "time": - return { - "standard_name": "time", - "long_name": "Time", - "axis": "T", - } - return { - "long_name": name.replace("_", " ").title(), - "standard_name": name.lower(), - "units": str(coord.units.units) if coord.units else ("m" if "depth" in name else "1"), - } - - def get_xarray_data_var_name(dataset) -> str: """Return the main xarray data variable name.""" if "data" in dataset.data_vars: return "data" # XDAS-style files can surface the primary payload under a None key while - # exposing coord helper arrays (for example *_indices/*_values) as data vars. + # exposing coordinate helper arrays as additional data variables. if None in dataset.data_vars: return None if len(dataset.data_vars) == 1: @@ -172,35 +24,6 @@ def get_xarray_data_var_name(dataset) -> str: raise ValueError(msg) -@cache -def get_xarray_engine(on_missing: str = "raise") -> str | None: - """Return the preferred xarray engine for NetCDF-4 files.""" - # Map importable module names onto the engine strings xarray expects. - module_to_engine = { - "netCDF4": "netcdf4", - "h5netcdf": "h5netcdf", - } - for module_name, engine_name in module_to_engine.items(): - mod = optional_import(module_name, on_missing="ignore") - if mod is not None: - return engine_name - if on_missing == "ignore": - return None - msg = ( - "Either netCDF4 or h5netcdf is required for NetCDF-4 read/write " - "functionality." - ) - raise MissingOptionalDependencyError(msg) - - -def iter_written_aux_coords(patch: dc.Patch): - """Yield names of auxiliary coordinates serialized to NetCDF.""" - for name, coord in patch.coords.coord_map.items(): - if coord._partial or name in patch.dims: - continue - yield name - - def parse_cf_version(cf_version: str) -> tuple[int, int]: """Parse a CF version string into comparable major/minor integers.""" parts = cf_version.split(".") @@ -234,191 +57,43 @@ def get_cf_version(h5file: h5py.File) -> str | None: return None -def extract_patch_attrs_from_netcdf(h5file: h5py.File) -> dict: - """Extract DASCore patch attrs from NetCDF global attributes.""" - attrs = {} - attr_mapping = { - "station": "station", - "network": "network", - "instrument": "instrument_id", - "acquisition": "acquisition_id", - "tag": "tag", - "data_category": "data_category", - "category": "category", - } - for cf_name, patch_name in attr_mapping.items(): - if cf_name in h5file.attrs: - value = h5file.attrs[cf_name] - if isinstance(value, bytes): - value = value.decode("utf-8", errors="ignore") - if value: - attrs[patch_name] = value - for cf_name in ("data_type", "source_data_type"): - if cf_name not in h5file.attrs: - continue - value = h5file.attrs[cf_name] - if isinstance(value, bytes): - value = value.decode("utf-8", errors="ignore") - if value: - attrs.setdefault("data_type", value) - return attrs - - -def _handle_time_interpolation( - h5file: h5py.File, coord_name: str, coord_data: np.ndarray -) -> np.ndarray | None: - """Decode XDAS-style time tie points into a full datetime coordinate.""" - time_values_name = f"{coord_name}_values" - if time_values_name not in h5file: - return None - time_values_var = h5file[time_values_name] - units = time_values_var.attrs.get("units", "") - if isinstance(units, bytes): - units = units.decode("utf-8", errors="ignore") - if "since" not in units: +def _get_tie_point_coord(h5file, coord_name: str, coord_len: int) -> np.ndarray | None: + """Decode one XDAS-style tie-point coordinate array.""" + values_name = f"{coord_name}_values" + indices_name = f"{coord_name}_indices" + if values_name not in h5file: return None - time_values = time_values_var[:] - time_indices_name = f"{coord_name}_indices" - if time_indices_name in h5file: - time_indices = h5file[time_indices_name][:] - if len(time_values) >= 2 and len(time_indices) >= 2: - interpolated_values = np.interp(coord_data, time_indices, time_values) - else: - interpolated_values = time_values - else: - interpolated_values = time_values - try: - return cf_time_to_datetime64(interpolated_values, units) - except (ValueError, KeyError): - return None - - -def read_netcdf_coordinates( - h5file: h5py.File, data_var_name: str | None = None -) -> CoordManager: - """Read coordinate information from a NetCDF file into a CoordManager.""" - coord_order = [] - main_data_var = data_var_name or find_main_data_variable(h5file) - expected_shape = None - if main_data_var and main_data_var in h5file: - data_var = h5file[main_data_var] - expected_shape = data_var.shape - if "DIMENSION_LIST" in data_var.attrs: - dim_list = data_var.attrs["DIMENSION_LIST"] - for ref_array in dim_list: - try: - ref = ref_array[0] - dim_scale = h5file[ref] - coord_order.append(dim_scale.name.strip("/")) - except (IndexError, KeyError, TypeError): - pass - - dim_coords = {} - for name, dataset in h5file.items(): - if not isinstance(dataset, h5py.Dataset): - continue - if not (dataset.is_scale or "NAME" in dataset.attrs): - continue - if expected_shape and dataset.shape[0] not in expected_shape: - continue - if any(skip in name.lower() for skip in ("_points", "_indices", "_values", "_interpolation")): - continue - data = dataset[:] - if name == "time" or dataset.attrs.get("axis") == "T": - units = dataset.attrs.get("units", "") - if isinstance(units, bytes): - units = units.decode("utf-8", errors="ignore") - if "since" in units: - try: - dim_coords["time"] = cf_time_to_datetime64(data, units) - continue - except (ValueError, KeyError): - pass - time_coord = _handle_time_interpolation(h5file, name, data) - if time_coord is not None: - dim_coords["time"] = time_coord - continue - dim_coords[name] = data - - coords = {} - if coord_order: - for dim_name in coord_order: - if dim_name in dim_coords: - coords[dim_name] = dim_coords[dim_name] - for name, data in dim_coords.items(): - if name not in coords: - coords[name] = data - else: - coords = dim_coords - - if main_data_var and main_data_var in h5file: - coord_attr = h5file[main_data_var].attrs.get("coordinates", "") - if isinstance(coord_attr, bytes): - coord_attr = coord_attr.decode("utf-8", errors="ignore") - for name in coord_attr.split(): - if name in coords or name not in h5file: - continue - dataset = h5file[name] - if not isinstance(dataset, h5py.Dataset): - continue - dims_attr = dataset.attrs.get("_DASCORE_DIMS", "") - if isinstance(dims_attr, bytes): - dims_attr = dims_attr.decode("utf-8", errors="ignore") - if dims_attr: - coords[name] = (tuple(dims_attr.split(",")), dataset[:]) - - return dc.core.coordmanager.get_coord_manager(coords) - - -def validate_cf_compliance(h5file: h5py.File) -> list[str]: - """Validate CF compliance and return a list of issues.""" - issues = [] - if "Conventions" not in h5file.attrs: - issues.append("Missing required 'Conventions' global attribute") - for name, dataset in h5file.items(): - if isinstance(dataset, h5py.Dataset) and dataset.is_scale: - if "units" not in dataset.attrs: - issues.append(f"Coordinate '{name}' missing 'units' attribute") - if name == "time" or dataset.attrs.get("standard_name") == "time": - units = dataset.attrs.get("units", "") - if not units or "since" not in str(units): - issues.append(f"Time coordinate '{name}' has invalid units") - for name, dataset in h5file.items(): - if isinstance(dataset, h5py.Dataset) and not dataset.is_scale: - if "units" not in dataset.attrs: - issues.append(f"Data variable '{name}' missing 'units' attribute") - if "long_name" not in dataset.attrs: - issues.append(f"Data variable '{name}' missing 'long_name' attribute") - return issues - - -def find_main_data_variable(h5file: h5py.File) -> str | None: - """Find the main data variable in a NetCDF file.""" - priority_names = [ - "data", - "acoustic_data", - "das_data", - "strain", - "strain_rate", - "velocity", - "amplitude", - ] - candidates = [] - for name, item in h5file.items(): - if not _is_data_variable_candidate(item): - continue - if name in priority_names or _has_priority_standard_name(item, priority_names): - return name - candidates.append(name) - return candidates[0] if candidates else None - - -def _is_data_variable_candidate(item) -> bool: - """Check if an item is a candidate for the main data variable.""" - return isinstance(item, h5py.Dataset) and not item.is_scale and item.ndim >= 2 - - -def _has_priority_standard_name(item, priority_names: list[str]) -> bool: - """Check if a dataset standard_name matches a known priority name.""" - std_name = str(item.attrs.get("standard_name", "")).lower() - return any(name in std_name for name in priority_names) + values_var = h5file[values_name] + values = values_var[:] + if indices_name in h5file: + indices = h5file[indices_name][:] + if len(values) >= 2 and len(indices) >= 2: + sample_index = np.arange(coord_len, dtype=np.float64) + if np.issubdtype(np.asarray(values).dtype, np.datetime64): + value_ns = values.astype("datetime64[ns]").astype(np.int64) + values = np.interp(sample_index, indices, value_ns).astype(np.int64) + values = values.astype("datetime64[ns]") + else: + values = np.interp(sample_index, indices, values) + return values + + +def _get_dim_coord(h5file, coord_name: str, coord_len: int) -> np.ndarray: + """Return one dimension coordinate for a coord-less payload variable.""" + tied_values = _get_tie_point_coord(h5file, coord_name, coord_len) + if tied_values is not None: + return tied_values + if coord_name in h5file: + return h5file[coord_name][:] + return np.arange(coord_len) + + +def get_coord_manager_for_coordless_data_var( + h5file, dims: tuple[str, ...], shape: tuple[int, ...] +): + """Build dimension coordinates for payloads xarray exposes without coords.""" + coords = { + dim: _get_dim_coord(h5file, dim, size) + for dim, size in zip(dims, shape, strict=True) + } + return dc.get_coord_manager(coords=coords, dims=dims) diff --git a/dascore/utils/io.py b/dascore/utils/io.py index 8eb1abe31..53ff30a56 100644 --- a/dascore/utils/io.py +++ b/dascore/utils/io.py @@ -258,8 +258,7 @@ def patch_to_xarray(patch: PatchType): # Omit None-valued attrs because xarray backends may reject them during # NetCDF serialization, while a missing attr round-trips cleanly. attrs = { - key: value for key, value in dict(patch.attrs).items() - if value is not None + key: value for key, value in dict(patch.attrs).items() if value is not None } patch_dims = patch.dims coords = {} diff --git a/tests/test_io/test_common_io.py b/tests/test_io/test_common_io.py index 31ef65f92..d937c621a 100644 --- a/tests/test_io/test_common_io.py +++ b/tests/test_io/test_common_io.py @@ -9,6 +9,7 @@ """ from __future__ import annotations + from contextlib import suppress from functools import cache from io import BytesIO, UnsupportedOperation @@ -45,8 +46,11 @@ ) from dascore.utils.downloader import fetch, get_registry_df from dascore.utils.misc import all_close, iterate -from tests.test_io._common_io_test_utils import get_flat_io_test, skip_missing, skip_timeout - +from tests.test_io._common_io_test_utils import ( + get_flat_io_test, + skip_missing, + skip_timeout, +) # --- Fixtures diff --git a/tests/test_io/test_netcdf/test_netcdf.py b/tests/test_io/test_netcdf/test_netcdf.py index b13fb369d..92fcca20a 100644 --- a/tests/test_io/test_netcdf/test_netcdf.py +++ b/tests/test_io/test_netcdf/test_netcdf.py @@ -3,28 +3,19 @@ from __future__ import annotations import importlib.util +from typing import ClassVar import h5py import numpy as np import pytest import dascore as dc -from dascore.exceptions import MissingOptionalDependencyError from dascore.io.netcdf import core as netcdf_core from dascore.io.netcdf import utils as netcdf_utils from dascore.io.netcdf.utils import ( - cf_time_to_datetime64, - datetime64_to_cf_time, - extract_patch_attrs_from_netcdf, - find_main_data_variable, - get_cf_data_attrs, - get_cf_global_attrs, get_cf_version, is_netcdf4_file, - read_netcdf_coordinates, - validate_cf_compliance, ) -from dascore.utils.downloader import fetch def _get_xarray_netcdf_engine() -> str | None: @@ -131,138 +122,6 @@ def rng(): class TestNetCDFUtils: """Tests for NetCDF utility functions.""" - @pytest.fixture - def test_datetime_array(self): - """Create test datetime array for CF time conversion tests.""" - return np.array( - [ - "2023-01-01T00:00:00", - "2023-01-01T01:00:00", - "2023-01-01T02:00:00", - ], - dtype="datetime64[ns]", - ) - - @pytest.fixture - def expected_cf_times(self): - """Expected CF time values for test datetime array.""" - expected_base = 1672531200.0 # 2023-01-01T00:00:00 in epoch seconds - return np.array([expected_base, expected_base + 3600, expected_base + 7200]) - - @pytest.fixture - def test_patch_attrs(self): - """Create test patch attributes for global attrs tests.""" - return dc.PatchAttrs( - station="TEST_STATION", - network="TEST_NET", - instrument_id="TEST_INST", - ) - - def test_datetime64_to_cf_time(self, test_datetime_array, expected_cf_times): - """Test conversion of datetime64 to CF time format.""" - cf_times = datetime64_to_cf_time(test_datetime_array) - np.testing.assert_array_almost_equal(cf_times, expected_cf_times) - - def test_get_cf_data_attrs(self): - """Test CF data attributes generation.""" - attrs = get_cf_data_attrs("strain_rate") - assert attrs["standard_name"] == "strain_rate" - assert attrs["units"] == "1/s" - assert "_FillValue" in attrs - - def test_get_cf_global_attrs(self, test_patch_attrs): - """Test CF global attributes generation.""" - attrs = get_cf_global_attrs(test_patch_attrs, "1.8") - assert attrs["Conventions"] == "CF-1.8" - assert attrs["station"] == "TEST_STATION" - assert attrs["network"] == "TEST_NET" - assert attrs["instrument"] == "TEST_INST" - assert "institution" not in attrs - - def test_cf_time_to_datetime64(self): - """Test CF time to datetime64 conversion.""" - # Test with seconds since epoch - cf_times = np.array([0, 3600, 7200]) # 0, 1, 2 hours since epoch - units = "seconds since 1970-01-01 00:00:00" - - result = cf_time_to_datetime64(cf_times, units) - - expected = np.array( - ["1970-01-01T00:00:00", "1970-01-01T01:00:00", "1970-01-01T02:00:00"], - dtype="datetime64[ns]", - ) - - np.testing.assert_array_equal(result, expected) - - def test_cf_time_to_datetime64_different_units(self): - """Test CF time conversion with different time units.""" - # Test with days since epoch - cf_times = np.array([0, 1, 2]) - units = "days since 2023-01-01 00:00:00" - - result = cf_time_to_datetime64(cf_times, units) - - expected = np.array( - ["2023-01-01T00:00:00", "2023-01-02T00:00:00", "2023-01-03T00:00:00"], - dtype="datetime64[ns]", - ) - - np.testing.assert_array_equal(result, expected) - - def test_cf_time_to_datetime64_invalid_units(self): - """Test CF time conversion with invalid units.""" - cf_times = np.array([0, 1, 2]) - invalid_units = "invalid format" - - with pytest.raises(ValueError, match="Invalid CF time units format"): - cf_time_to_datetime64(cf_times, invalid_units) - - def test_find_main_data_variable(self, tmp_path, rng): - """Test finding main data variable in NetCDF file.""" - path = tmp_path / "test_data_var.nc" - - with h5py.File(path, "w") as h5file: - # Create dimension scales - time_ds = h5file.create_dataset("time", data=np.arange(100)) - time_ds.make_scale("time") - - # Create various datasets. - h5file.create_dataset("metadata", data=np.array([1, 2, 3])) - h5file.create_dataset("other_data", data=rng.random((50, 50))) - - # Create priority data variable - strain_data = h5file.create_dataset( - "strain_data", data=rng.random((100, 50)) - ) - strain_data.attrs["standard_name"] = b"strain" - - with h5py.File(path, "r") as h5file: - main_var = find_main_data_variable(h5file) - assert main_var == "strain_data" - - def test_find_main_data_variable_no_priority(self, tmp_path, rng): - """Test finding main data variable when no priority match.""" - path = tmp_path / "test_no_priority.nc" - - with h5py.File(path, "w") as h5file: - h5file.create_dataset("first_candidate", data=rng.random((50, 50))) - h5file.create_dataset("second_candidate", data=rng.random((60, 60))) - - with h5py.File(path, "r") as h5file: - main_var = find_main_data_variable(h5file) - assert main_var == "first_candidate" # Returns first candidate - - def test_find_main_data_variable_none_found(self, tmp_path): - """Test finding main data variable when none found.""" - path = tmp_path / "test_none_found.nc" - - with h5py.File(path, "w") as h5file: - h5file.create_dataset("only_1d", data=np.array([1, 2, 3])) # Only 1D data - - with h5py.File(path, "r") as h5file: - main_var = find_main_data_variable(h5file) - assert main_var is None - def test_get_cf_version(self, tmp_path): """Test extracting CF version from NetCDF file.""" path = tmp_path / "test_cf_version.nc" @@ -327,113 +186,6 @@ def test_is_netcdf4_file_from_dimension_scales(self, tmp_path): with h5py.File(path, "r") as h5file: assert not is_netcdf4_file(h5file) - def test_get_cf_data_attrs_partial_match(self): - """Partial data-type matches should map to known CF names.""" - attrs = get_cf_data_attrs("my_acceleration_trace") - assert attrs["standard_name"] == "acceleration" - - def test_get_cf_global_attrs_optional_fields(self): - """Optional patch attrs should be propagated to global metadata.""" - attrs = dc.PatchAttrs( - station="TEST_STATION", - network="TEST_NET", - instrument_id="TEST_INST", - acquisition_id="ACQ_01", - tag="taggy", - data_category="DAS", - data_type="strain_rate", - history=("step-1", "step-2"), - category="processed", - ) - - out = get_cf_global_attrs(attrs, "1.8") - - assert out["acquisition"] == "ACQ_01" - assert out["tag"] == "taggy" - assert out["data_category"] == "DAS" - assert out["data_type"] == "strain_rate" - assert out["category"] == "processed" - assert out["processing_history"] == "step-1 | step-2" - - def test_extract_patch_attrs_source_data_type_fallback(self, tmp_path): - """Source data type should populate data_type when primary attr is absent.""" - path = tmp_path / "source_data_type.nc" - with h5py.File(path, "w") as h5file: - h5file.attrs["network"] = np.bytes_("TEST_NET") - h5file.attrs["source_data_type"] = np.bytes_("strain_rate") - - with h5py.File(path, "r") as h5file: - attrs = extract_patch_attrs_from_netcdf(h5file) - - assert attrs["network"] == "TEST_NET" - assert attrs["data_type"] == "strain_rate" - - def test_handle_time_interpolation_without_indices(self, tmp_path): - """Time interpolation should fall back to raw time_values when needed.""" - path = tmp_path / "time_interp.nc" - with h5py.File(path, "w") as h5file: - h5file.create_dataset("time", data=np.arange(3)) - time_values = h5file.create_dataset("time_values", data=np.arange(3)) - time_values.attrs["units"] = "seconds since 1970-01-01 00:00:00" - - with h5py.File(path, "r") as h5file: - out = netcdf_utils._handle_time_interpolation(h5file, "time", np.arange(3)) - - expected = np.array( - ["1970-01-01T00:00:00", "1970-01-01T00:00:01", "1970-01-01T00:00:02"], - dtype="datetime64[ns]", - ) - np.testing.assert_array_equal(out, expected) - - def test_handle_time_interpolation_invalid_units_returns_none(self, tmp_path): - """Interpolation should fail quietly when auxiliary units are invalid.""" - path = tmp_path / "time_interp_invalid.nc" - with h5py.File(path, "w") as h5file: - h5file.create_dataset("time", data=np.arange(3)) - time_values = h5file.create_dataset("time_values", data=np.arange(3)) - time_values.attrs["units"] = "invalid" - - with h5py.File(path, "r") as h5file: - out = netcdf_utils._handle_time_interpolation(h5file, "time", np.arange(3)) - - assert out is None - - def test_read_netcdf_coordinates_with_non_dim_coord_attrs(self, tmp_path): - """Non-dimensional coordinates should honor _DASCORE_DIMS metadata.""" - path = tmp_path / "non_dim_coords.nc" - with h5py.File(path, "w") as h5file: - time_ds = h5file.create_dataset("time", data=np.arange(4)) - time_ds.make_scale("time") - dist_ds = h5file.create_dataset("distance", data=np.arange(3)) - dist_ds.make_scale("distance") - data_ds = h5file.create_dataset("data", data=np.ones((3, 4))) - data_ds.dims[0].attach_scale(dist_ds) - data_ds.dims[1].attach_scale(time_ds) - data_ds.attrs["coordinates"] = np.bytes_("latitude") - lat_ds = h5file.create_dataset("latitude", data=np.linspace(1.0, 2.0, 3)) - lat_ds.attrs["_DASCORE_DIMS"] = np.bytes_("distance") - - with h5py.File(path, "r") as h5file: - coords = read_netcdf_coordinates(h5file) - - assert coords.dims == ("distance", "time") - assert "latitude" in coords.coord_map - assert coords.dim_map["latitude"] == ("distance",) - - def test_validate_cf_compliance_invalid_time_units(self, tmp_path): - """Invalid time units should be reported as a CF issue.""" - path = tmp_path / "invalid_time_units.nc" - with h5py.File(path, "w") as h5file: - h5file.attrs["Conventions"] = "CF-1.8" - time_ds = h5file.create_dataset("time", data=np.arange(4)) - time_ds.make_scale("time") - time_ds.attrs["units"] = "seconds" - - with h5py.File(path, "r") as h5file: - issues = validate_cf_compliance(h5file) - - assert any("invalid units" in issue for issue in issues) - def test_is_netcdf4_file_handles_attribute_error(self): """Attribute errors during detection should return False.""" @@ -475,128 +227,10 @@ def test_get_cf_version_decodes_bytes(self, tmp_path): with h5py.File(path, "r") as h5file: assert get_cf_version(h5file) == "1.9" - def test_handle_time_interpolation_decodes_units_and_handles_bad_units( - self, tmp_path - ): - """Interpolation should decode bytes and handle unsupported units.""" - decoded_path = tmp_path / "time_interp_bytes.nc" - with h5py.File(decoded_path, "w") as h5file: - h5file.create_dataset("time", data=np.arange(3)) - time_values = h5file.create_dataset("time_values", data=np.arange(3)) - time_values.attrs["units"] = np.bytes_("seconds since 1970-01-01 00:00:00") - time_indices = h5file.create_dataset("time_indices", data=np.array([0])) - time_indices[...] = np.array([0]) - - with h5py.File(decoded_path, "r") as h5file: - out = netcdf_utils._handle_time_interpolation(h5file, "time", np.arange(3)) - - assert out.dtype == "datetime64[ns]" - - invalid_path = tmp_path / "time_interp_unsupported.nc" - with h5py.File(invalid_path, "w") as h5file: - h5file.create_dataset("time", data=np.arange(3)) - time_values = h5file.create_dataset("time_values", data=np.arange(3)) - time_values.attrs["units"] = "fortnights since 2023-01-01" - - with h5py.File(invalid_path, "r") as h5file: - out = netcdf_utils._handle_time_interpolation(h5file, "time", np.arange(3)) - - assert out is None - - def test_read_netcdf_coordinates_handles_bad_refs_and_auxiliaries(self, tmp_path): - """Coordinate reader should tolerate bad refs and skip auxiliaries.""" - path = tmp_path / "coords_bad_refs.nc" - with h5py.File(path, "w") as h5file: - h5file.create_group("aux_group") - h5file.create_dataset("time_values", data=np.arange(3)) - data_ds = h5file.create_dataset("data", data=np.ones((2, 3))) - data_ds.attrs["DIMENSION_LIST"] = np.array([[b"/missing"]], dtype="S8") - data_ds.attrs["coordinates"] = np.bytes_("missing aux_group latitude") - lat_ds = h5file.create_dataset("latitude", data=np.arange(2)) - lat_ds.attrs["_DASCORE_DIMS"] = np.bytes_("distance") - dist_ds = h5file.create_dataset("distance", data=np.arange(2)) - dist_ds.make_scale("distance") - time_ds = h5file.create_dataset("time", data=np.arange(3)) - time_ds.make_scale("time") - time_ds.attrs["units"] = np.bytes_("fortnights since 2023-01-01") - - with h5py.File(path, "r") as h5file: - coords = read_netcdf_coordinates(h5file) - - assert coords.dims == ("distance", "time") - assert "latitude" in coords.coord_map - - def test_read_netcdf_coordinates_without_dimension_list(self, tmp_path): - """Coordinate reader should fall back to discovered order when needed.""" - path = tmp_path / "coords_no_dimension_list.nc" - with h5py.File(path, "w") as h5file: - dist_ds = h5file.create_dataset("distance", data=np.arange(2)) - dist_ds.make_scale("distance") - time_ds = h5file.create_dataset("time", data=np.arange(3)) - time_ds.make_scale("time") - time_ds.attrs["units"] = "seconds since 1970-01-01 00:00:00" - h5file.create_dataset("data", data=np.ones((2, 3))) - - with h5py.File(path, "r") as h5file: - coords = read_netcdf_coordinates(h5file) - - assert coords.dims == ("distance", "time") - np.testing.assert_array_equal(coords.coord_map["distance"].values, np.arange(2)) - - def test_read_netcdf_coordinates_adds_unordered_dim_coords(self, tmp_path): - """Coords discovered outside DIMENSION_LIST should still be retained.""" - path = tmp_path / "coords_extra_scale.nc" - with h5py.File(path, "w") as h5file: - dist_ds = h5file.create_dataset("distance", data=np.arange(2)) - dist_ds.make_scale("distance") - time_ds = h5file.create_dataset("time", data=np.arange(3)) - time_ds.make_scale("time") - time_ds.attrs["units"] = np.bytes_("seconds since 1970-01-01 00:00:00") - extra_ds = h5file.create_dataset("channel", data=np.arange(2)) - extra_ds.make_scale("channel") - data_ds = h5file.create_dataset("data", data=np.ones((2, 3))) - data_ds.dims[0].attach_scale(dist_ds) - data_ds.dims[1].attach_scale(time_ds) - - with h5py.File(path, "r") as h5file: - coords = read_netcdf_coordinates(h5file) - - assert "channel" in coords.coord_map - - def test_read_netcdf_coordinates_skips_auxiliary_named_scales(self, tmp_path): - """Auxiliary scale names should be skipped before coord collection.""" - path = tmp_path / "coords_aux_skip.nc" - with h5py.File(path, "w") as h5file: - dist_ds = h5file.create_dataset("distance", data=np.arange(2)) - dist_ds.make_scale("distance") - time_ds = h5file.create_dataset("time", data=np.arange(3)) - time_ds.make_scale("time") - time_ds.attrs["units"] = "seconds since 1970-01-01 00:00:00" - aux_ds = h5file.create_dataset("time_values", data=np.arange(3)) - aux_ds.make_scale("time_values") - h5file.create_dataset("data", data=np.ones((2, 3))) - - with h5py.File(path, "r") as h5file: - coords = read_netcdf_coordinates(h5file) - - assert "time_values" not in coords.coord_map - class TestNetCDFCoreHelpers: """Direct tests for lightweight NetCDF core helpers.""" - def test_coord_attrs_cover_distance_and_depth(self, example_patch): - """Coordinate helper should populate distance and depth CF attrs.""" - distance_attrs = netcdf_utils.coord_attrs( - "distance", example_patch.coords.coord_map["distance"] - ) - depth_coord = dc.get_coord(data=np.arange(3), units="ft") - depth_attrs = netcdf_utils.coord_attrs("sensor_depth", depth_coord) - - assert distance_attrs["standard_name"] == "distance" - assert distance_attrs["units"] == "m" - assert depth_attrs["units"] == "ft" - def test_get_xarray_data_var_name(self): """Dataset helper should find expected data variables or raise.""" xr = pytest.importorskip("xarray") @@ -605,7 +239,7 @@ def test_get_xarray_data_var_name(self): ds_multi = xr.Dataset({"signal": (("x",), [1, 2]), "other": (("x",), [3, 4])}) class _DatasetWithNone: - data_vars = {None: object(), "distance_indices": object()} + data_vars: ClassVar = {None: object(), "distance_indices": object()} assert netcdf_utils.get_xarray_data_var_name(ds_with_data) == "data" assert netcdf_utils.get_xarray_data_var_name(_DatasetWithNone()) is None @@ -650,173 +284,20 @@ def test_get_format_accepts_comma_separated_conventions(self, tmp_path): with h5py.File(path, "r") as h5file: assert formatter.get_format(h5file) == ("NETCDF_CF", "1.8") - def test_get_data_variable_name_raises_for_missing_data(self, tmp_path): - """Formatter should raise when no main data variable exists.""" - formatter = netcdf_core.NetCDFCFV18() - path = tmp_path / "missing_data.nc" - with h5py.File(path, "w") as h5file: - h5file.create_dataset("time", data=np.arange(3)) - - with h5py.File(path, "r") as h5file: - with pytest.raises(ValueError, match="No suitable data variable found"): - formatter._get_data_variable_name(h5file) - - def test_apply_coordinate_filtering_handles_selected_and_unrelated_kwargs( - self, example_patch - ): - """Filtering should only use known coordinate kwargs.""" - formatter = netcdf_core.NetCDFCFV18() - - unfiltered = formatter._apply_coordinate_filtering(example_patch, {"tag": "x"}) - filtered = formatter._apply_coordinate_filtering( - example_patch, {"distance": (10, 20)} - ) - - assert unfiltered.equals(example_patch) - assert filtered.data.shape[0] < example_patch.data.shape[0] - - def test_read_uses_xarray_dataset_and_merges_missing_coords( - self, minimal_cf_netcdf_path, monkeypatch - ): - """Read should accept xarray-backed data and merge extra coords.""" - - class FakeCoord: - def __init__(self, dims, values): - self.dims = dims - self.values = values - - class FakeDataArray: - def __init__(self, data): - self.data = data - self.coords = { - "distance": FakeCoord(("distance",), np.arange(data.shape[0])), - "time": FakeCoord(("time",), np.arange(data.shape[1])), - "latitude": FakeCoord( - ("distance",), np.linspace(1.0, 2.0, data.shape[0]) - ), - } - - def load(self): - return self - - class FakeDataset: - def __init__(self, data_array): - self.data_vars = {None: data_array} - self._data_array = data_array - - def __enter__(self): - return self - - def __exit__(self, *args): - return False - - def get(self, name): - return None - - def __getitem__(self, item): - assert item is None - return self._data_array - - fake_coords = dc.get_coord_manager( - coords={"distance": np.arange(2), "time": np.arange(3)}, - dims=("distance", "time"), - ) - fake_data_array = FakeDataArray(np.arange(6).reshape(2, 3)) - fake_dataset = FakeDataset(fake_data_array) - fake_xarray = type( - "FakeXarray", - (), - {"open_dataset": staticmethod(lambda *args, **kwargs: fake_dataset)}, - ) + def test_get_write_encoding_invalid_compression_raises(self): + """Write encoding should reject unsupported compression values.""" formatter = netcdf_core.NetCDFCFV18() - def _optional_import(name, on_missing="raise"): - if name == "xarray": - return fake_xarray - if name == "netCDF4": - return object() - return None + with pytest.raises(ValueError, match="only gzip compression"): + formatter._get_write_encoding(compression="szip") - monkeypatch.setattr(netcdf_core, "optional_import", _optional_import) - monkeypatch.setattr( - netcdf_core, - "read_netcdf_coordinates", - lambda *args: fake_coords, - ) - monkeypatch.setattr( - formatter, "_get_data_variable_name", lambda resource: "__values__" - ) - monkeypatch.setattr( - formatter, "_get_patch_attrs", lambda resource: {"tag": "fake"} - ) - - spool = formatter.read(minimal_cf_netcdf_path) - patch = spool[0] - - np.testing.assert_array_equal(patch.data, fake_data_array.data) - assert patch.attrs.tag == "fake" - assert "latitude" in patch.coords.coord_map - - def test_read_falls_back_to_single_xarray_data_var( - self, minimal_cf_netcdf_path, monkeypatch - ): - """Read should use the only xarray data var when HDF-derived name misses.""" - - class FakeDataArray: - def __init__(self): - self.data = np.ones((2, 2)) - self.coords = {} - - def load(self): - return self - - class FakeDataset: - def __init__(self): - self.data_vars = {"signal": FakeDataArray()} - - def __enter__(self): - return self - - def __exit__(self, *args): - return False - - def get(self, name): - return None - - def __getitem__(self, item): - return self.data_vars[item] - - fake_coords = dc.get_coord_manager( - coords={"distance": np.arange(2), "time": np.arange(2)}, - dims=("distance", "time"), - ) - fake_xarray = type( - "FakeXarray", - (), - {"open_dataset": staticmethod(lambda *args, **kwargs: FakeDataset())}, - ) + def test_get_write_encoding_explicit_chunks(self): + """Write encoding should pass explicit chunk sizes through.""" formatter = netcdf_core.NetCDFCFV18() - def _optional_import(name, on_missing="raise"): - if name == "xarray": - return fake_xarray - if name == "netCDF4": - return object() - return None + out = formatter._get_write_encoding(chunks=(10, 20)) - monkeypatch.setattr(netcdf_core, "optional_import", _optional_import) - monkeypatch.setattr( - netcdf_core, "read_netcdf_coordinates", lambda *args: fake_coords - ) - monkeypatch.setattr( - formatter, "_get_data_variable_name", lambda resource: "missing" - ) - monkeypatch.setattr( - formatter, "_get_patch_attrs", lambda resource: {"tag": "fallback"} - ) - - spool = formatter.read(minimal_cf_netcdf_path) - assert spool[0].attrs.tag == "fallback" + assert out["chunksizes"] == (10, 20) def test_read_returns_empty_spool_for_empty_filtered_patch( self, minimal_cf_netcdf_path, monkeypatch @@ -826,7 +307,17 @@ def test_read_returns_empty_spool_for_empty_filtered_patch( class FakeDataArray: def __init__(self): self.data = np.ones((1, 1)) - self.coords = {} + + def _make_coord(dims, vals): + return type("Coord", (), {"dims": dims, "values": vals})() + + self.coords = { + "distance": _make_coord(("distance",), np.array([0])), + "time": _make_coord(("time",), np.array([0])), + } + self.attrs = {} + self.dims = ("distance", "time") + self.shape = self.data.shape def load(self): return self @@ -834,6 +325,7 @@ def load(self): class FakeDataset: def __init__(self): self.data_vars = {"data": FakeDataArray()} + self.attrs = {} def __enter__(self): return self @@ -844,16 +336,15 @@ def __exit__(self, *args): def get(self, name): return self.data_vars["data"] + def __getitem__(self, item): + return self.data_vars[item] + formatter = netcdf_core.NetCDFCFV18() fake_xarray = type( "FakeXarray", (), {"open_dataset": staticmethod(lambda *args, **kwargs: FakeDataset())}, ) - fake_coords = dc.get_coord_manager( - coords={"distance": np.arange(1), "time": np.arange(1)}, - dims=("distance", "time"), - ) empty_patch = dc.Patch( data=np.empty((0, 0)), coords=dc.get_coord_manager( @@ -872,24 +363,13 @@ def _optional_import(name, on_missing="raise"): return None monkeypatch.setattr(netcdf_core, "optional_import", _optional_import) - monkeypatch.setattr( - netcdf_core, "read_netcdf_coordinates", lambda *args: fake_coords - ) - monkeypatch.setattr( - formatter, "_get_data_variable_name", lambda resource: "data" - ) - monkeypatch.setattr( - formatter, "_get_patch_attrs", lambda resource: {"tag": "empty"} - ) - monkeypatch.setattr( - formatter, "_apply_coordinate_filtering", lambda patch, kwargs: empty_patch - ) + monkeypatch.setattr(dc.Patch, "select", lambda self, **kwargs: empty_patch) - spool = formatter.read(minimal_cf_netcdf_path) + spool = formatter.read(minimal_cf_netcdf_path, time=(0, 1)) assert len(spool) == 0 def test_write_uses_xarray_dataset_path(self, example_patch, tmp_path, monkeypatch): - """Write should pass CF attrs and encoding through the xarray path.""" + """Write should pass the xarray dataset and encoding through unchanged.""" class FakeCoord: def __init__(self): @@ -949,22 +429,12 @@ def _optional_import(name, on_missing="raise"): lambda patch: fake_data_array, ) - patch_with_partial = example_patch.new( - coords=example_patch.coords.update( - latitude=("distance", np.linspace(0.0, 1.0, example_patch.shape[0])), - quality=("distance", dc.get_coord(shape=(example_patch.shape[0],))), - ) - ) out_path = tmp_path / "write_stub.nc" - formatter.write( - patch_with_partial.update(attrs={"data_type": "strain_rate"}), out_path - ) + formatter.write(example_patch, out_path) dataset = fake_data_array.dataset - assert dataset.attrs["Conventions"] == "CF-1.8" - assert dataset.attrs["source_data_type"] == "strain_rate" - assert dataset.data.attrs["standard_name"] == "strain_rate" - assert dataset.data.attrs["coordinates"] == "latitude" + assert dataset.attrs == {"Conventions": "CF-1.8"} + assert dataset.data.attrs == {} assert dataset.to_netcdf_calls _args, kwargs = dataset.to_netcdf_calls[0] assert kwargs["encoding"] is None @@ -1003,13 +473,6 @@ def test_write_netcdf(self, example_patch, tmp_path): # Check it's detected as NetCDF assert is_netcdf4_file(h5file) - # Check global attributes - assert "Conventions" in h5file.attrs - conventions = h5file.attrs["Conventions"] - if isinstance(conventions, bytes): - conventions = conventions.decode() - assert "CF-" in conventions - # Check coordinates exist assert "time" in h5file assert "distance" in h5file @@ -1138,7 +601,7 @@ def test_current_netcdf_output_is_readable_by_xarray(self, patch_variant, tmp_pa def test_written_file_exposes_expected_metadata_to_xarray( self, patch_with_attrs, tmp_path ): - """Xarray should see the expected CF metadata on DASCore output.""" + """Xarray should see the basic dataset structure on DASCore output.""" xr = pytest.importorskip("xarray") engine = _require_xarray_netcdf_engine() @@ -1146,21 +609,12 @@ def test_written_file_exposes_expected_metadata_to_xarray( dc.write(patch_with_attrs, path, file_format="netcdf_cf") with xr.open_dataset(path, engine=engine) as dataset: - assert dataset.attrs["Conventions"] == "CF-1.8" - assert dataset.attrs["source_data_type"] == "strain_rate" - assert dataset.attrs["station"] == "TEST_STATION" - assert dataset.attrs["network"] == "TEST_NET" - assert list(dataset.data_vars) == ["data"] - assert list(dataset.coords) == ["distance", "time"] - assert dataset["data"].attrs["standard_name"] == "strain_rate" - assert dataset["data"].attrs["units"] == "1/s" - assert dataset["time"].attrs["axis"] == "T" - assert dataset["time"].attrs["standard_name"] == "time" - assert dataset["distance"].attrs["standard_name"] == "distance" - assert dataset["distance"].attrs["units"] == "m" - assert "featureType" not in dataset.attrs - assert "coordinates" not in dataset["data"].attrs + assert dataset.attrs["Conventions"] == "CF-1.8" + assert set(dataset.coords) == {"distance", "time"} + assert dataset["data"].attrs["station"] == "TEST_STATION" + assert dataset["data"].attrs["network"] == "TEST_NET" + assert dataset["data"].attrs["data_type"] == "strain_rate" def test_dascore_and_xarray_roundtrip_agree_for_non_dim_coords( self, patch_with_non_dim_coords, tmp_path @@ -1252,9 +706,28 @@ def test_compression_options(self, compressed_netcdf_file): original_patch.data, recovered_patch.data, decimal=6 ) + class TestNetCDFUtilsAdvanced: """Additional tests for NetCDF utility functions.""" + def test_coordless_coord_manager_falls_back_without_tie_points(self): + """Coord fallback should use direct vars or arange without tie points.""" + xr = pytest.importorskip("xarray") + dataset = xr.Dataset( + data_vars={"distance": (("distance",), np.array([0.0, 2.0, 4.0]))} + ) + + coords = netcdf_utils.get_coord_manager_for_coordless_data_var( + dataset, + dims=("distance", "time"), + shape=(3, 4), + ) + + np.testing.assert_array_equal( + coords.coord_map["distance"].values, [0.0, 2.0, 4.0] + ) + np.testing.assert_array_equal(coords.coord_map["time"].values, np.arange(4)) + @pytest.fixture def cf_compliant_file(self, tmp_path): """Create a CF-compliant NetCDF file for testing.""" @@ -1271,53 +744,6 @@ def cf_compliant_file(self, tmp_path): dc.write(patch, path, file_format="netcdf_cf") return path - def test_extract_patch_attrs_from_netcdf(self, cf_compliant_file): - """Test extracting patch attributes from NetCDF file.""" - with h5py.File(cf_compliant_file, "r") as h5file: - attrs = extract_patch_attrs_from_netcdf(h5file) - assert isinstance(attrs, dict) - - def test_read_netcdf_coordinates(self, cf_compliant_file): - """Test reading coordinates from NetCDF file.""" - with h5py.File(cf_compliant_file, "r") as h5file: - coords = read_netcdf_coordinates(h5file) - - assert "time" in coords.coord_map - assert "distance" in coords.coord_map - - # Check that time coordinate is properly converted from CF format - time_coord = coords.coord_map["time"] - assert len(time_coord) > 0 - - def test_validate_cf_compliance(self, cf_compliant_file): - """Test CF compliance validation.""" - with h5py.File(cf_compliant_file, "r") as h5file: - issues = validate_cf_compliance(h5file) - - # Our implementation should produce CF-compliant files - assert len(issues) == 0, f"CF compliance issues found: {issues}" - - def test_validate_cf_compliance_with_issues(self, tmp_path, rng): - """Test CF compliance validation with non-compliant file.""" - path = tmp_path / "non_compliant.nc" - - # Create a file with CF issues - with h5py.File(path, "w") as h5file: - # Missing Conventions attribute - time_ds = h5file.create_dataset("time", data=np.arange(100)) - time_ds.make_scale("time") - # Missing units attribute - - h5file.create_dataset("data", data=rng.random((100, 50))) - # Missing long_name and units attributes - - with h5py.File(path, "r") as h5file: - issues = validate_cf_compliance(h5file) - - assert len(issues) > 0 - assert any("Conventions" in issue for issue in issues) - assert any("units" in issue for issue in issues) - def test_coordinate_filtering_during_read(self, cf_compliant_file): """Test coordinate filtering during NetCDF read.""" # Read with time filtering @@ -1338,23 +764,6 @@ def test_coordinate_filtering_during_read(self, cf_compliant_file): # Should have fewer time samples assert filtered_patch.data.shape[1] < original_patch.data.shape[1] - def test_different_data_types_cf_attrs(self): - """Test CF attributes for different data types.""" - # Test various data types - test_cases = [ - ("strain", "1", "Strain"), - ("velocity", "m/s", "Velocity"), - ("temperature", "K", "Temperature"), - ("pressure", "Pa", "Pressure"), - ("unknown_type", "1", "Distributed Acoustic Sensing data"), - ] - - for data_type, expected_units, expected_long_name in test_cases: - attrs = get_cf_data_attrs(data_type) - assert attrs["units"] == expected_units - assert attrs["long_name"] == expected_long_name - assert "_FillValue" in attrs - def test_netcdf_format_detection_edge_cases(self, tmp_path, rng): """Test NetCDF format detection edge cases.""" _require_xarray_netcdf_engine() @@ -1388,19 +797,10 @@ def test_netcdf_format_detection_edge_cases(self, tmp_path, rng): # Should be able to read older CF versions spool = dc.read(path2, file_format="netcdf_cf") assert len(spool) == 1 - - def test_read_external_xdas_netcdf_file(self): - """External xdas NetCDF should keep its current readable structure.""" - _require_xarray_netcdf_engine() - path = fetch("xdas_netcdf.nc") - - patch = dc.read(path, file_format="netcdf_cf")[0] - - assert patch.dims == ("time", "distance") - assert patch.shape == (300, 401) + patch = spool[0] assert set(patch.coords.coord_map) == {"time", "distance"} assert patch.attrs.tag == "" - assert patch.attrs["_source_patch_id"] == netcdf_utils.XDAS_PAYLOAD_VARIABLE + assert patch.attrs["_source_patch_id"] == "data" def test_error_conditions(self, tmp_path): """Test various error conditions.""" @@ -1415,26 +815,3 @@ def test_error_conditions(self, tmp_path): # Use DASCore's standard reading interface to test error conditions with pytest.raises(ValueError, match="No suitable data variable found"): dc.read(path, file_format="netcdf_cf") - - def test_cf_time_edge_cases(self): - """Test CF time conversion edge cases.""" - # Test unsupported time unit - cf_times = np.array([0, 1, 2]) - invalid_units = "fortnights since 2023-01-01" - - with pytest.raises(ValueError, match="Unsupported time unit"): - cf_time_to_datetime64(cf_times, invalid_units) - - # Test various supported units - supported_units = [ - ("hours since 2023-01-01", "h"), - ("minutes since 2023-01-01", "min"), - ("milliseconds since 2023-01-01", "ms"), - ("microseconds since 2023-01-01", "us"), - ] - - for units, _ in supported_units: - cf_times = np.array([0, 1, 2]) - result = cf_time_to_datetime64(cf_times, units) - assert result.dtype == "datetime64[ns]" - assert len(result) == 3 From f4475a3843f4aa664e6c49f4e8dcb2b8eb00827d Mon Sep 17 00:00:00 2001 From: Derrick Chambers Date: Thu, 9 Apr 2026 17:55:59 +0200 Subject: [PATCH 5/5] skip netcdf when needed --- tests/test_io/test_netcdf/test_netcdf.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/tests/test_io/test_netcdf/test_netcdf.py b/tests/test_io/test_netcdf/test_netcdf.py index 92fcca20a..ef3f419e5 100644 --- a/tests/test_io/test_netcdf/test_netcdf.py +++ b/tests/test_io/test_netcdf/test_netcdf.py @@ -17,6 +17,8 @@ is_netcdf4_file, ) +pytest.importorskip("xarray") + def _get_xarray_netcdf_engine() -> str | None: """Return an xarray engine that can open NetCDF-4 files, if available.""" @@ -363,6 +365,22 @@ def _optional_import(name, on_missing="raise"): return None monkeypatch.setattr(netcdf_core, "optional_import", _optional_import) + monkeypatch.setattr( + netcdf_core, + "xarray_to_patch", + lambda data_array: dc.Patch( + data=data_array.data, + coords=dc.get_coord_manager( + coords={ + name: (coord.dims, coord.values) + for name, coord in data_array.coords.items() + }, + dims=data_array.dims, + ), + dims=data_array.dims, + attrs=dict(data_array.attrs), + ), + ) monkeypatch.setattr(dc.Patch, "select", lambda self, **kwargs: empty_patch) spool = formatter.read(minimal_cf_netcdf_path, time=(0, 1))