NetCDF-4 support

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

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."""

dascore/io/netcdf/__init__.py

@@ -0,0 +1,5 @@
+"""NetCDF IO support for DASCore."""
+
+from __future__ import annotations
+
+from dascore.io.netcdf.core import NetCDFCFV18

dascore/io/netcdf/core.py

@@ -0,0 +1,172 @@
+"""Core NetCDF IO implementation built on xarray."""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+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, xarray_to_patch
+from dascore.utils.misc import optional_import
+
+from .utils import (
+    XDAS_PAYLOAD_VARIABLE,
+    get_cf_version,
+    get_coord_manager_for_coordless_data_var,
+    get_xarray_data_var_name,
+    is_netcdf4_file,
+    parse_cf_version,
+)
+
+
+class NetCDFCFV18(FiberIO):
+    """NetCDF-4 IO using xarray for read/write and CF markers for detection."""
+
+    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")
+        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 _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 through xarray.
+
+        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")  # 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,
+            encoding={"data": encoding} if encoding else None,
+        )
+
+    def scan(self, resource: H5Reader, **kwargs) -> list[ScanPayload]:
+        """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 | {"_source_patch_id": source_patch_id},
+                coords=coord_manager,
+                dims=dims,
+                shape=shape,
+                dtype=dtype,
+                source_patch_id=source_patch_id,
+            )
+        ]
+
+    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 _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
+
+    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]

dascore/io/netcdf/utils.py

@@ -0,0 +1,99 @@
+"""NetCDF helper functions for DASCore IO."""
+
+from __future__ import annotations
+
+import h5py
+import numpy as np
+
+import dascore as dc
+
+XDAS_PAYLOAD_VARIABLE = "__values__"
+
+
+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 coordinate helper arrays as additional data variables.
+    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
+
+
+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 _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
+    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)

dascore/utils/io.py

@@ -255,7 +255,11 @@ 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():

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"

tests/test_io/test_common_io.py

@@ -29,6 +29,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
@@ -87,6 +88,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