Repository: pydata/xarray

Issue:
Saving and loading an array of strings changes datatype to object
### What is your issue?

See the code below
```python
import xarray as xr
# Make an xarray with an array of strings
da = xr.DataArray(
    data=[["a", "b"], ["c", "d"]],
    dims=["x", "y"],
    coords={"x": [0, 1], "y": [0, 1]},
)
da.to_netcdf("test.nc", mode='w')
# Load the xarray back in
da_loaded = xr.load_dataarray("test.nc")
assert da.dtype == da_loaded.dtype, "Dtypes don't match"
```

Now `da_loaded.dtype` is `dtype('O')`.

Same happens with `engine="h5netcdf"`.

```
import xarray
xarray.show_versions()
```
```
INSTALLED VERSIONS
------------------
commit: None
python: 3.11.0 | packaged by conda-forge | (main, Jan 14 2023, 12:26:40) [Clang 14.0.6 ]
python-bits: 64
OS: Darwin
OS-release: 22.3.0
machine: arm64
processor: arm
byteorder: little
LC_ALL: en_US.UTF-8
LANG: en_US.UTF-8
LOCALE: ('en_US', 'UTF-8')
libhdf5: 1.12.2
libnetcdf: None

xarray: 2023.2.0
pandas: 1.5.3
numpy: 1.24.2
scipy: 1.10.0
netCDF4: None
pydap: None
h5netcdf: 1.1.0
h5py: 3.8.0
Nio: None
zarr: None
cftime: None
nc_time_axis: None
PseudoNetCDF: None
rasterio: None
cfgrib: None
iris: None
bottleneck: None
dask: None
distributed: None
matplotlib: 3.7.0
cartopy: None
seaborn: 0.12.2
numbagg: None
fsspec: 2023.1.0
cupy: None
pint: None
sparse: None
flox: None
numpy_groupies: None
setuptools: 67.3.2
pip: 23.0
conda: None
pytest: 7.2.1
mypy: 1.0.1
IPython: 8.10.0
sphinx: 5.3.0
```

Discussion:
kmuehlbauer: Resurrected this PR, rebased on latest main and force-pushed. 
kmuehlbauer: So, the h5netcdf-part does only work for versions h5py >= 3.7 and h5netcdf >= 1.1.0 (this fails for the 3.9 min deps). How can I account for that?  Is there some best-practice?
kmuehlbauer: > So, the h5netcdf-part does only work for versions h5py >= 3.7 and h5netcdf >= 1.1.0 (this fails for the 3.9 min deps). How can I account for that? Is there some best-practice?

We are already unblocked for the h5py-version as the minimum dependency version is already at 3.7. h5netcdf get's unblocked on November, 23rd (in two month's time).

Should we just wait it out and make vlen strings available in one block?
dcherian: > h5netcdf get's unblocked on November, 23rd

I'm OK with bumping it early for this. If you want we could wait till after the next release to bring it down to one month :)
kmuehlbauer: > I'm OK with bumping it early for this. If you want we could wait till after the next release to bring it down to one month :)

Yes, wait until after next release and bump early.
kmuehlbauer: @dcherian Finally, I think this is ready.
dcherian: Thanks @kmuehlbauer 

Patch:
diff --git a/ci/requirements/min-all-deps.yml b/ci/requirements/min-all-deps.yml
index 7d0f29c0..22076f88 100644
--- a/ci/requirements/min-all-deps.yml
+++ b/ci/requirements/min-all-deps.yml
@@ -16,11 +16,11 @@ dependencies:
   - dask-core=2022.7
   - distributed=2022.7
   - flox=0.5
-  - h5netcdf=1.0
+  - h5netcdf=1.1
   # h5py and hdf5 tend to cause conflicts
   # for e.g. hdf5 1.12 conflicts with h5py=3.1
   # prioritize bumping other packages instead
-  - h5py=3.6
+  - h5py=3.7
   - hdf5=1.12
   - hypothesis
   - iris=3.2
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 2fb76cfe..928a2c95 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -32,13 +32,14 @@ New Features
   By `Sam Levang <https://github.com/slevang>`_.
 - Allow the usage of h5py drivers (eg: ros3) via h5netcdf (:pull:`8360`).
   By `Ezequiel Cimadevilla <https://github.com/zequihg50>`_.
+- Enable VLEN string fill_values, preserve VLEN string dtypes (:issue:`1647`, :issue:`7652`, :issue:`7868`, :pull:`7869`).
+  By `Kai Mühlbauer <https://github.com/kmuehlbauer>`_.
 
 Breaking changes
 ~~~~~~~~~~~~~~~~
 - drop support for `cdms2 <https://github.com/CDAT/cdms>`_. Please use
   `xcdat <https://github.com/xCDAT/xcdat>`_ instead (:pull:`8441`).
   By `Justus Magin <https://github.com/keewis>`_.
-
 - Following pandas, :py:meth:`infer_freq` will return ``"Y"``, ``"YS"``,
   ``"QE"``, ``"ME"``, ``"h"``, ``"min"``, ``"s"``, ``"ms"``, ``"us"``, or
   ``"ns"`` instead of ``"A"``, ``"AS"``, ``"Q"``, ``"M"``, ``"H"``, ``"T"``,
@@ -46,6 +47,8 @@ Breaking changes
   deprecation of the latter frequency strings (:issue:`8394`, :pull:`8415`). By
   `Spencer Clark <https://github.com/spencerkclark>`_.
 - Bump minimum tested pint version to ``>=0.22``. By `Deepak Cherian <https://github.com/dcherian>`_.
+- Minimum supported versions for the following packages have changed: ``h5py >=3.7``, ``h5netcdf>=1.1``.
+  By `Kai Mühlbauer <https://github.com/kmuehlbauer>`_.
 
 Deprecations
 ~~~~~~~~~~~~
diff --git a/xarray/backends/h5netcdf_.py b/xarray/backends/h5netcdf_.py
index a68a44b5..d9385fc6 100644
--- a/xarray/backends/h5netcdf_.py
+++ b/xarray/backends/h5netcdf_.py
@@ -271,15 +271,6 @@ class H5NetCDFStore(WritableCFDataStore):
         dtype = _get_datatype(variable, raise_on_invalid_encoding=check_encoding)
 
         fillvalue = attrs.pop("_FillValue", None)
-        if dtype is str and fillvalue is not None:
-            raise NotImplementedError(
-                "h5netcdf does not yet support setting a fill value for "
-                "variable-length strings "
-                "(https://github.com/h5netcdf/h5netcdf/issues/37). "
-                f"Either remove '_FillValue' from encoding on variable {name!r} "
-                "or set {'dtype': 'S1'} in encoding to use the fixed width "
-                "NC_CHAR type."
-            )
 
         if dtype is str:
             dtype = h5py.special_dtype(vlen=str)
diff --git a/xarray/backends/netCDF4_.py b/xarray/backends/netCDF4_.py
index f21f15bf..1aee4c1c 100644
--- a/xarray/backends/netCDF4_.py
+++ b/xarray/backends/netCDF4_.py
@@ -490,16 +490,6 @@ class NetCDF4DataStore(WritableCFDataStore):
 
         fill_value = attrs.pop("_FillValue", None)
 
-        if datatype is str and fill_value is not None:
-            raise NotImplementedError(
-                "netCDF4 does not yet support setting a fill value for "
-                "variable-length strings "
-                "(https://github.com/Unidata/netcdf4-python/issues/730). "
-                f"Either remove '_FillValue' from encoding on variable {name!r} "
-                "or set {'dtype': 'S1'} in encoding to use the fixed width "
-                "NC_CHAR type."
-            )
-
         encoding = _extract_nc4_variable_encoding(
             variable, raise_on_invalid=check_encoding, unlimited_dims=unlimited_dims
         )
diff --git a/xarray/coding/variables.py b/xarray/coding/variables.py
index df660f90..48719760 100644
--- a/xarray/coding/variables.py
+++ b/xarray/coding/variables.py
@@ -562,3 +562,15 @@ class NonStringCoder(VariableCoder):
 
     def decode(self):
         raise NotImplementedError()
+
+
+class ObjectVLenStringCoder(VariableCoder):
+    def encode(self):
+        return NotImplementedError
+
+    def decode(self, variable: Variable, name: T_Name = None) -> Variable:
+        if variable.dtype == object and variable.encoding.get("dtype", False) == str:
+            variable = variable.astype(variable.encoding["dtype"])
+            return variable
+        else:
+            return variable
diff --git a/xarray/conventions.py b/xarray/conventions.py
index cf207f0c..75f816e6 100644
--- a/xarray/conventions.py
+++ b/xarray/conventions.py
@@ -265,6 +265,10 @@ def decode_cf_variable(
             var = strings.CharacterArrayCoder().decode(var, name=name)
         var = strings.EncodedStringCoder().decode(var)
 
+    if original_dtype == object:
+        var = variables.ObjectVLenStringCoder().decode(var)
+        original_dtype = var.dtype
+
     if mask_and_scale:
         for coder in [
             variables.UnsignedIntegerCoder(),


Test patch:
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index 80b6951d..85248b5c 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -864,12 +864,13 @@ class CFEncodedBase(DatasetIOBase):
         assert check_vlen_dtype(original["a"].dtype) == str
         with self.roundtrip(original) as actual:
             assert_identical(original, actual)
-            assert object == actual["a"].dtype
-            assert actual["a"].dtype == original["a"].dtype
-            # only check metadata for capable backends
-            # eg. NETCDF3 based backends do not roundtrip metadata
-            if actual["a"].dtype.metadata is not None:
-                assert check_vlen_dtype(actual["a"].dtype) == str
+            if np.issubdtype(actual["a"].dtype, object):
+                # only check metadata for capable backends
+                # eg. NETCDF3 based backends do not roundtrip metadata
+                if actual["a"].dtype.metadata is not None:
+                    assert check_vlen_dtype(actual["a"].dtype) == str
+            else:
+                assert actual["a"].dtype == np.dtype("<U1")
 
     @pytest.mark.parametrize(
         "decoded_fn, encoded_fn",
@@ -1374,32 +1375,39 @@ class NetCDF4Base(NetCDFBase):
             with self.open(tmp_file, group="data/2") as actual2:
                 assert_identical(data2, actual2)
 
-    def test_encoding_kwarg_vlen_string(self) -> None:
-        for input_strings in [[b"foo", b"bar", b"baz"], ["foo", "bar", "baz"]]:
-            original = Dataset({"x": input_strings})
-            expected = Dataset({"x": ["foo", "bar", "baz"]})
-            kwargs = dict(encoding={"x": {"dtype": str}})
-            with self.roundtrip(original, save_kwargs=kwargs) as actual:
-                assert actual["x"].encoding["dtype"] is str
-                assert_identical(actual, expected)
-
-    def test_roundtrip_string_with_fill_value_vlen(self) -> None:
+    @pytest.mark.parametrize(
+        "input_strings, is_bytes",
+        [
+            ([b"foo", b"bar", b"baz"], True),
+            (["foo", "bar", "baz"], False),
+            (["foó", "bár", "baź"], False),
+        ],
+    )
+    def test_encoding_kwarg_vlen_string(
+        self, input_strings: list[str], is_bytes: bool
+    ) -> None:
+        original = Dataset({"x": input_strings})
+
+        expected_string = ["foo", "bar", "baz"] if is_bytes else input_strings
+        expected = Dataset({"x": expected_string})
+        kwargs = dict(encoding={"x": {"dtype": str}})
+        with self.roundtrip(original, save_kwargs=kwargs) as actual:
+            assert actual["x"].encoding["dtype"] == "<U3"
+            assert actual["x"].dtype == "<U3"
+            assert_identical(actual, expected)
+
+    @pytest.mark.parametrize("fill_value", ["XXX", "", "bár"])
+    def test_roundtrip_string_with_fill_value_vlen(self, fill_value: str) -> None:
         values = np.array(["ab", "cdef", np.nan], dtype=object)
         expected = Dataset({"x": ("t", values)})
 
-        # netCDF4-based backends don't support an explicit fillvalue
-        # for variable length strings yet.
-        # https://github.com/Unidata/netcdf4-python/issues/730
-        # https://github.com/h5netcdf/h5netcdf/issues/37
-        original = Dataset({"x": ("t", values, {}, {"_FillValue": "XXX"})})
-        with pytest.raises(NotImplementedError):
-            with self.roundtrip(original) as actual:
-                assert_identical(expected, actual)
+        original = Dataset({"x": ("t", values, {}, {"_FillValue": fill_value})})
+        with self.roundtrip(original) as actual:
+            assert_identical(expected, actual)
 
         original = Dataset({"x": ("t", values, {}, {"_FillValue": ""})})
-        with pytest.raises(NotImplementedError):
-            with self.roundtrip(original) as actual:
-                assert_identical(expected, actual)
+        with self.roundtrip(original) as actual:
+            assert_identical(expected, actual)
 
     def test_roundtrip_character_array(self) -> None:
         with create_tmp_file() as tmp_file: