Is your feature request related to a problem?

Complex numbers are not natively handled by netCDF, and currently only the h5netcdf engine in xarray can handle them. However, it produces files unreadable by current versions of netCDF. There are also many different conventions, and h5netcdf can only read a few of them.

Describe the solution you'd like

I've written a drop-in extension to netCDF that seamlessly handles complex numbers, and can read many different conventions.

The Python API is built on top of netCDF4-python, so only a very simple change to xarray is needed to use it:

```diff modified xarray/backends/api.py @@ -115,7 +115,10 @@ def _get_default_engine_gz() -> Literal["scipy"]: def _get_default_engine_netcdf() -> Literal["netcdf4", "scipy"]: engine: Literal["netcdf4", "scipy"] try: - import netCDF4 # noqa: F401 + try: + import nc_complex # noqa: F401 + except ImportError: + import netCDF4 # noqa: F401

     engine = "netcdf4"
 except ImportError:  # pragma: no cover

modified xarray/backends/netCDF4_.py @@ -327,7 +327,10 @@ class NetCDF4DataStore(WritableCFDataStore): def init( self, manager, group=None, mode=None, lock=NETCDF4_PYTHON_LOCK, autoclose=False ): - import netCDF4 + try: + import nc_complex as netCDF4 + except ImportError: + import netCDF4

     if isinstance(manager, netCDF4.Dataset):
         if group is None:

@@ -364,7 +367,10 @@ class NetCDF4DataStore(WritableCFDataStore): lock_maker=None, autoclose=False, ): - import netCDF4 + try: + import nc_complex as netCDF4 + except ImportError: + import netCDF4

     if isinstance(filename, os.PathLike):
         filename = os.fspath(filename)

```

and then it works out of the box:

```python import numpy as np import xarray as xr

complex_array = np.array([0 + 0j, 1 + 0j, 0 + 1j, 1 + 1j, 0.25 + 0.75j], dtype="c16")

ds = xr.Dataset({"complex_data": ("x", complex_array)}) ds.to_netcdf("test.nc")

with xr.open_dataset("test.nc") as df: print(df.complex_data) assert ds == df

<xarray.DataArray 'complex_data' (x: 5)>

[5 values with dtype=complex128]

Dimensions without coordinates: x

```

If the library is built against the latest main of netCDF-C, then it also supports files written by h5netcdf. And while the package isn't on PyPI yet, built wheels with this feature are available from CI jobs.

nc_complex is not quite ready yet, but I wanted to gauge interest for adding it as an optional replacement for netCDF4 in xarray.

Describe alternatives you've considered

No response

Additional context

No response

See #8288

• [ ] Closes #xxxx
• [ ] Tests added
• [ ] User visible changes (including notable bug fixes) are documented in whats-new.rst
• [ ] New functions/methods are listed in api.rst

These changes allow me to make a pretty minimal subclass that just swaps out the netCDF4.Dataset class for my version:

```python from xarray.backends.netCDF4_ import NetCDF4BackendEntrypoint, NetCDF4DataStore from ._wrapper import Dataset # This is my netCDF extension

class NcComplexDataStore(NetCDF4DataStore): @classmethod def DatasetClass(self) -> type: return Dataset

class NcComplexBackendEntrypoint(NetCDF4BackendEntrypoint): DataStore = NcComplexDataStore ```

Uses sphinx-inline-tabs for pretty and compact options.

Also reword surrounding text a little for clarity.

This is more consistent with the other exceptions raised in this function, and helps to more quickly diagnose issues when e.g. creating datasets.

• [x] User visible changes (including notable bug fixes) are documented in whats-new.rst

Is this a large enough change to warrant an entry in whats-new?

What happened?

Calling np.divide() passing a bool expression involving a DataArray results in a RecursionError starting at apply_ufunc. Passing DataArray.data works fine.

What did you expect to happen?

No response

Minimal Complete Verifiable Example

```Python import numpy as np import xarray as xr

N = 10 array = np.random.random((N,)) da = xr.DataArray(np.random.random((N,)), dims=("x",), coords={"x": np.arange(N)})

This works

print(np.divide(array, da, where=(da > 0.2).data))

This gives recursion error

print(np.divide(array, da, where=da > 0.2)) ```

MVCE confirmation

• [X] Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• [X] Complete example — the example is self-contained, including all data and the text of any traceback.
• [X] Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• [X] New issue — a search of GitHub Issues suggests this is not a duplicate.

Relevant log output

```Python Traceback (most recent call last): File "divide_recursion_mvce.py", line 7, in <module> print(np.divide(array, da, where=da > 0.2)) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "venv/lib64/python3.11/site-packages/xarray/core/arithmetic.py", line 86, in array_ufunc return apply_ufunc( ^^^^^^^^^^^^ File "venv/lib64/python3.11/site-packages/xarray/core/computation.py", line 1197, in apply_ufunc return apply_dataarray_vfunc( ^^^^^^^^^^^^^^^^^^^^^^ File "venv/lib64/python3.11/site-packages/xarray/core/computation.py", line 304, in apply_dataarray_vfunc result_var = func(data_vars) ^^^^^^^^^^^^^^^^ File "venv/lib64/python3.11/site-packages/xarray/core/computation.py", line 761, in apply_variable_ufunc result_data = func(input_data) ^^^^^^^^^^^^^^^^^ File "venv/lib64/python3.11/site-packages/xarray/core/arithmetic.py", line 86, in array_ufunc return apply_ufunc( ^^^^^^^^^^^^ File "venv/lib64/python3.11/site-packages/xarray/core/computation.py", line 1210, in apply_ufunc return apply_array_ufunc(func, *args, dask=dask) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

...

File "venv/lib64/python3.11/site-packages/xarray/core/computation.py", line 1210, in apply_ufunc return apply_array_ufunc(func, *args, dask=dask) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "venv/lib64/python3.11/site-packages/xarray/core/computation.py", line 815, in apply_array_ufunc if any(is_chunked_array(arg) for arg in args): ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "venv/lib64/python3.11/site-packages/xarray/core/computation.py", line 815, in <genexpr> if any(is_chunked_array(arg) for arg in args): ^^^^^^^^^^^^^^^^^^^^^ File "venv/lib64/python3.11/site-packages/xarray/core/pycompat.py", line 99, in is_chunked_array return is_duck_dask_array(x) or (is_duck_array(x) and hasattr(x, "chunks")) ^^^^^^^^^^^^^^^^^^^^^ File "venv/lib64/python3.11/site-packages/xarray/core/pycompat.py", line 95, in is_duck_dask_array return is_duck_array(x) and is_dask_collection(x) ^^^^^^^^^^^^^^^^^^^^^ File "venv/lib64/python3.11/site-packages/xarray/core/pycompat.py", line 87, in is_dask_collection if module_available("dask"): ^^^^^^^^^^^^^^^^^^^^^^^^ File "venv/lib64/python3.11/site-packages/xarray/core/utils.py", line 1160, in module_available return importlib.util.find_spec(module) is not None ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "<frozen importlib.util>", line 103, in find_spec File "<frozen importlib._bootstrap>", line 1080, in _find_spec File "<frozen importlib._bootstrap>", line 924, in find_spec File "<frozen importlib._bootstrap>", line 241, in _call_with_frames_removed RecursionError: maximum recursion depth exceeded while calling a Python object ```

Anything else we need to know?

No response

Environment