Do not crash on arrays of strings

[bouweandela]

I'm trying to use ncdata to load Pangeo data into iris, but encountering issues with a dimension that contains string variables. Complete example code:

import intake
import ncdata.iris_xarray

cat = intake.open_esm_datastore(
    "https://storage.googleapis.com/cmip6/pangeo-cmip6.json",
).search(
    table_id="Amon",
    source_id="GFDL-CM4",
    experiment_id="historical",
    variable_id="tas",
    member_id="r1i1p1f1",
)

ds = cat.to_dask()
# ds = ds.drop_vars(["member_id"])
print(ds)
cubes = ncdata.iris_xarray.cubes_from_xarray(ds)
print(cubes)

A similar issue was reported in #111. Here is a minimal reproducer:

In [1]: import xarray as xr

In [2]: import ncdata.iris_xarray

In [3]: ds = xr.DataArray.from_dict({
   ...:     "coords": {
   ...:         "member_id": {"dims": "member_id", "data": ["r1i1p1f1"]}
   ...:     },
   ...:     "dims": "member_id",
   ...:     "data": [1],
   ...:     "name": "a",
   ...: }).to_dataset()
   ...: 

In [4]: ncdata.iris_xarray.cubes_from_xarray(ds)
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[4], line 1
----> 1 ncdata.iris_xarray.cubes_from_xarray(ds)

File ~/mambaforge/envs/esmvaltool/lib/python3.13/site-packages/ncdata/iris_xarray.py:54, in cubes_from_xarray(xrds, xr_save_kwargs, iris_load_kwargs)
     52 xr_save_kwargs = xr_save_kwargs or {}
     53 ncdata = from_xarray(xrds, **xr_save_kwargs)
---> 54 cubes = to_iris(ncdata, **iris_load_kwargs)
     55 return cubes

File ~/mambaforge/envs/esmvaltool/lib/python3.13/site-packages/ncdata/iris.py:52, in to_iris(ncdata, **iris_load_kwargs)
     50     dslikes = [Nc4DatasetLike(data) for data in ncdata]
     51 else:
---> 52     dslikes = Nc4DatasetLike(ncdata)
     53 cubes = CubeList(iris.load(dslikes, **iris_load_kwargs))
     54 return cubes

File ~/mambaforge/envs/esmvaltool/lib/python3.13/site-packages/ncdata/dataset_like.py:112, in Nc4DatasetLike.__init__(self, ncdata)
    108 self._ncdata = ncdata
    109 # N.B. we need to create + store our OWN variables, as they are wrappers for
    110 #  the underlying NcVariable objects, with different properties.
    111 self.variables = {
--> 112     name: Nc4VariableLike._from_ncvariable(ncvar)
    113     for name, ncvar in self._ncdata.variables.items()
    114 }
    115 self.dimensions = {
    116     name: Nc4DimensionLike(dim) for name, dim in self._ncdata.dimensions.items()
    117 }

File ~/mambaforge/envs/esmvaltool/lib/python3.13/site-packages/ncdata/dataset_like.py:253, in Nc4VariableLike._from_ncvariable(cls, ncvar, dtype)
    251 if dtype is None:
    252     dtype = ncvar.dtype
--> 253 self = cls(
    254     ncvar=ncvar,
    255     datatype=dtype,
    256 )
    257 return self

File ~/mambaforge/envs/esmvaltool/lib/python3.13/site-packages/ncdata/dataset_like.py:246, in Nc4VariableLike.__init__(self, ncvar, datatype)
    240     array = da.ma.masked_array(da.zeros(self.shape, self.datatype), mask=True)
    242 # Convert from the "inner" raw form to masked-and-scaled, and then back again
    243 # (by assigning self._data_array).
    244 # This ensures that our inner data is stored with the correct fill-value, as
    245 # determined from our attributes and the netCDF4 default fill-values.
--> 246 array = self._maskandscale_inner_to_apparent(array)
    247 self._data_array = array

File ~/mambaforge/envs/esmvaltool/lib/python3.13/site-packages/ncdata/dataset_like.py:303, in Nc4VariableLike._maskandscale_inner_to_apparent(self, array)
    297 """
    298 Convert raw data values to masked+scaled.
    299 
    300 This replicates the netCDF4 default auto-scaling procedure.
    301 """
    302 # N.B. fill-value matches the internal raw (unscaled) values and dtype
--> 303 fv = self._get_fillvalue()
    304 if fv is not None:
    305     array = da.ma.masked_equal(array, fv)

File ~/mambaforge/envs/esmvaltool/lib/python3.13/site-packages/ncdata/dataset_like.py:293, in Nc4VariableLike._get_fillvalue(self)
    290     if self.dtype.itemsize != 1:
    291         # NOTE: single-byte types have NO default fill-value
    292         dtype_code = self.dtype.str[1:]
--> 293         fv = netCDF4.default_fillvals[dtype_code]
    294 return fv

KeyError: 'U8'

The code changes here result at least being able to load the cube, even if iris still not accepts string type DimCoords.

By the way, are you aware that xarray does not support using the default fill-value at all? pydata/xarray#2742

@pp-mo Do you think this pull request could be a step in the right direction to making ncdata work with string arrays? I would be happy to continue working on this if you think it's worth it or close if you think a different approach is needed.

[pp-mo]

Do you think this pull request could be a step in the right direction to making ncdata work with string arrays?

Sorry @bouweandela I've been tied up elsewhere, and forgot to respond here.
I will get round to considering this shortly!

[pp-mo]

Thanks for your MRE. but I think it shows that there are further things going wrong here, that more effort is needed to fix...

So, there's already a "way around" the problem that datasets like this won't load into ncdata.
As in this "how-to" example

But when we apply that, we see that there are other problems here associated with the Dask representation of the data :

>>> ds = xr.DataArray.from_dict({
...      "coords": {
...          "member_id": {"dims": "member_id", "data": ["r1i1p1f1"]}
...      },
...      "dims": "member_id",
...      "data": [1],
...      "name": "a",
... }).to_dataset()
>>> ds.to_netcdf("tmp.nc")
>>> ncds = ncdata.netcdf4.from_nc4("tmp.nc", dim_chunks={"member_id": 1})
>>> print(ncds)
<NcData: /
    dimensions:
        member_id = 1

    variables:
        <NcVariable(<U0): member_id(member_id)>
        <NcVariable(int64): a(member_id)>
>
>>> v = ncds.variables['member_id']
>>> print(v)
<NcVariable(<U0): member_id(member_id)>
>>> print(v.data)
dask.array<array, shape=(1,), dtype=<U1, chunksize=(1,), chunktype=numpy.ndarray>
>>> print(repr(v.data.compute()))
array(['r1i1p1f1'], dtype=object)
>>> 

So, the result when loaded has the wrong datatype, and the realisation of the dask array doesn't deliver the expected dtype either!
I think the "<U0" variable dtype possibly implies variable-length strings, but I'm not immediately sure how it got there.
The fact that the dask array dtype doesn't match what you get when you compute it is, unfortunately, just a known problem with dask -- it means that the proxy code is not delivering the expected result data.

It is nevertheless possible this way to load the file, get the variable content and re-present it in a form acceptable to ncdata -- and Iris. For instance, as in the example.

To be clear, what parts of this are important to you ?
- e.g. the behaviour in ncdata, loading from xarray, or just loading the data into Iris ?

[pp-mo]

BTW it is still work-in-progress to properly handle the ordinary character arrays in iris data,
so, unfortunately, a more thorough handling of variable-length strings is still somewhat further behind that.

[pp-mo]

Ping @bouweandela : I've done nothing more on this yet, but thought about it.
Can you explain what you most immediately need to achieve (e.g. for a workaround) ?

[bouweandela]

What we're trying to achieve is to add a way to load data from any intake-esm catalogue into ESMValTool in ESMValGroup/ESMValCore#2690. So the data comes in as an Xarray dataset and I need to convert it to Iris cubes. It would be rather convenient if that just 'worked' (provided the datasets follow the CF conventions) and there would be no need to manipulate the Xarray datasets before conversion because then we need to write special case code, for example ESMValGroup/ESMValCore#2690 (comment), for every new catalogue we add.

Our analyses do not even use the "member_id" coordinate, so even if it gets loaded as an AuxCoord instead of a DimCoord, they would run through just fine with the current behaviour of Iris if the code didn't stumble on the conversion step. Even models with many ensemble members in CMIP6 had less 100 of these, so if I would concatenate this data along the member_id coordinate I would end up with an array with fewer than 100 elements, so computational efficiency does not seem important at all.