Repository: dcs4cop/xcube

Issue:
`fill_value` and `dtype` not used correctly in `rasterize_features()` 
**Describe the bug**

Function `xcube.core.geom.rasterize_features()` allows specifying how a given feature shall be rasterized as an xarray variable by  passing target `fill_value` and `dtype` among other things for each feature. 

Currently `fill_value` and `dtype` are not used correctly to produce the desired result: If `dtype` is given and integer, the variable will have that dtype, although a given `fill _value` should actually be encoded as NaN in the resulting data. However the latter is impossible for integer data.

**Expected behavior**

If `fill_value` is given, the resulting variable must be a floating point type with NaN values where the feature value was `fill_value`. In addition, the _encoding_ of the variable must specify `_FillValue` and the given `dtype`. 


Discussion:

Patch:
diff --git a/CHANGES.md b/CHANGES.md
index a240011c..23cddf6d 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -4,6 +4,13 @@
 
 ### Fixes
 
+* Implementation of function `xcube.core.geom.rasterize_features()` 
+  has been changed to account for consistent use of a target variable's
+  `fill_value` and `dtype` for a given feature.
+  In-memory (decoded) variables now always use dtype `float64` and use 
+  `np.nan` to represent missing values. Persisted (encoded) variable data
+  will make use of the target `fill_value` and `dtype`. (#778)
+
 ## Changes in 0.13.0.dev7
 
 ### Other
diff --git a/xcube/core/geom.py b/xcube/core/geom.py
index 079c1619..1130c4a8 100644
--- a/xcube/core/geom.py
+++ b/xcube/core/geom.py
@@ -111,6 +111,10 @@ def rasterize_features(
                                # Deprecated, no longer used.
     }
 
+    Note that newly created variables will have data type `np.float64`
+    because `np.nan` is used to encode missing values. `fill_value` and
+    `dtype` are used to encode the variables when persisting the data.
+
     Currently, the coordinates of the geometries in the given
     *features* must use the same CRS as the given *dataset*.
 
@@ -235,20 +239,18 @@ def rasterize_features(
             'name',
             feature_prop_name.replace(' ', '_')
         )
-        var_dtype = var_prop_mapping.get('dtype', np.float64)
+        var_dtype = np.dtype(var_prop_mapping.get('dtype', np.float64))
         var_fill_value = var_prop_mapping.get('fill_value', np.nan)
         var_attrs = var_prop_mapping.get('attrs', {})
 
         feature_image = rasterized_features[feature_index]
-        if feature_image.dtype != var_dtype:
-            feature_image = da.Array.astype(feature_image,
-                                            dtype=var_dtype)
 
         feature_var = xr.DataArray(feature_image,
                                    coords=yx_coords,
                                    dims=yx_dims,
                                    attrs=var_attrs)
-        feature_var.encoding.update(fill_value=var_fill_value)
+        feature_var.encoding.update(_FillValue=var_fill_value,
+                                    dtype=var_dtype)
         dataset[var_name] = feature_var
 
     return dataset


Test patch:
diff --git a/test/core/test_geom.py b/test/core/test_geom.py
index c5fb1615..4b58bcf8 100644
--- a/test/core/test_geom.py
+++ b/test/core/test_geom.py
@@ -110,11 +110,14 @@ class RasterizeFeaturesIntoDataset(unittest.TestCase):
                                      features,
                                      ['a', 'b', 'c'],
                                      var_props=dict(
-                                         b=dict(name='b', dtype=np.float32,
+                                         b=dict(name='b',
+                                                dtype=np.float32,
                                                 fill_value=np.nan,
                                                 attrs=dict(units='meters')),
-                                         c=dict(name='c2', dtype=np.uint8,
-                                                fill_value=0, converter=int)
+                                         c=dict(name='c2',
+                                                dtype=np.uint8,
+                                                fill_value=0,
+                                                converter=int)
                                      ),
                                      tile_size=tile_size,
                                      in_place=False)
@@ -140,15 +143,32 @@ class RasterizeFeaturesIntoDataset(unittest.TestCase):
         self.assertEquals((10, 10), dataset.a.shape)
         self.assertEquals((10, 10), dataset.b.shape)
         self.assertEquals((10, 10), dataset.c2.shape)
+
+        # Assert in-memory (decoding) data types are correct.
         self.assertEquals(np.float64, dataset.a.dtype)
-        self.assertEquals(np.float32, dataset.b.dtype)
-        self.assertEquals(np.uint8, dataset.c2.dtype)
+        self.assertEquals(np.float64, dataset.b.dtype)
+        self.assertEquals(np.float64, dataset.c2.dtype)
+
+        # Assert external representation (encoding) information
+        # is correctly set up.
+        # See also test.core.test_xarray.XarrayEncodingTest
+        self.assertIs(nan, dataset.a.encoding.get('_FillValue'))
+        self.assertIs(nan, dataset.b.encoding.get('_FillValue'))
+        self.assertEqual(0, dataset.c2.encoding.get('_FillValue'))
+        self.assertEqual(np.dtype('float64'), dataset.a.encoding.get('dtype'))
+        self.assertEqual(np.dtype('float32'), dataset.b.encoding.get('dtype'))
+        self.assertEqual(np.dtype('uint8'), dataset.c2.encoding.get('dtype'))
+
+        # Other metadata
         self.assertEquals({}, dataset.a.attrs)
         self.assertEquals({'units': 'meters'}, dataset.b.attrs)
         self.assertEquals({}, dataset.c2.attrs)
         self.assertEquals((y_name, x_name), dataset.a.dims)
         self.assertEquals((y_name, x_name), dataset.b.dims)
         self.assertEquals((y_name, x_name), dataset.c2.dims)
+
+        # Assert actual data is correct
+
         actual_a_values = dataset.a.values
         expected_a_values = np.array(
             [[0.6, 0.6, 0.6, 0.6, 0.6, nan, nan, 0.8, 0.8, 0.8],
@@ -185,17 +205,16 @@ class RasterizeFeaturesIntoDataset(unittest.TestCase):
                                        actual_b_values)
         actual_c_values = dataset.c2.values
         expected_c_values = np.array(
-            [[8, 8, 8, 8, 8, 0, 0, 6, 6, 6],
-             [8, 8, 8, 8, 8, 0, 0, 6, 6, 6],
-             [8, 8, 8, 8, 8, 0, 0, 6, 6, 6],
-             [8, 8, 8, 8, 8, 0, 0, 6, 6, 6],
-             [8, 8, 8, 8, 8, 0, 0, 6, 6, 6],
-             [9, 9, 9, 9, 9, 0, 0, 7, 7, 7],
-             [9, 9, 9, 9, 9, 0, 0, 7, 7, 7],
-             [9, 9, 9, 9, 9, 0, 0, 7, 7, 7],
-             [9, 9, 9, 9, 9, 0, 0, 7, 7, 7],
-             [9, 9, 9, 9, 9, 0, 0, 7, 7, 7]],
-            dtype=np.uint8
+            [[8, 8, 8, 8, 8, nan, nan, 6, 6, 6],
+             [8, 8, 8, 8, 8, nan, nan, 6, 6, 6],
+             [8, 8, 8, 8, 8, nan, nan, 6, 6, 6],
+             [8, 8, 8, 8, 8, nan, nan, 6, 6, 6],
+             [8, 8, 8, 8, 8, nan, nan, 6, 6, 6],
+             [9, 9, 9, 9, 9, nan, nan, 7, 7, 7],
+             [9, 9, 9, 9, 9, nan, nan, 7, 7, 7],
+             [9, 9, 9, 9, 9, nan, nan, 7, 7, 7],
+             [9, 9, 9, 9, 9, nan, nan, 7, 7, 7],
+             [9, 9, 9, 9, 9, nan, nan, 7, 7, 7]]
         )
         if inverse_y:
             expected_c_values = expected_c_values[::-1, :]
diff --git a/test/core/test_xarray.py b/test/core/test_xarray.py
index 839e0c49..e44a0fd0 100644
--- a/test/core/test_xarray.py
+++ b/test/core/test_xarray.py
@@ -1,8 +1,10 @@
 import unittest
 
+import numpy as np
 import xarray as xr
 
 from test.sampledata import new_test_dataset
+from xcube.core.dsio import rimraf
 from xcube.core.gridmapping import GridMapping
 from xcube.core.new import new_cube
 from xcube.core.xarray import DatasetAccessor
@@ -65,17 +67,112 @@ class XCubeDatasetAccessorTest(unittest.TestCase):
         self.assertIsInstance(dataset.xcube.vars_to_dim(), xr.Dataset)
 
     def test_levels(self):
-        dataset = new_test_dataset(["2010-01-01", "2010-01-02", "2010-01-03", "2010-01-04", "2010-01-05"],
-                                   precipitation=0.4, temperature=275.2)
+        dataset = new_test_dataset(
+            ["2010-01-01", "2010-01-02", "2010-01-03", "2010-01-04",
+             "2010-01-05"],
+            precipitation=0.4, temperature=275.2)
         levels = dataset.xcube.levels(spatial_tile_shape=(45, 45))
         self.assertIsInstance(levels, list)
         self.assertEqual(3, len(levels))
-        self.assertTrue(all(isinstance(level, xr.Dataset) for level in levels))
+        self.assertTrue(
+            all(isinstance(level, xr.Dataset) for level in levels))
         self.assertTrue(all("precipitation" in level for level in levels))
         self.assertTrue(all("temperature" in level for level in levels))
         self.assertEqual([(5, 180, 360), (5, 90, 180), (5, 45, 90)],
                          [level.precipitation.shape for level in levels])
-        self.assertEqual([((1, 1, 1, 1, 1), (45, 45, 45, 45), (45, 45, 45, 45, 45, 45, 45, 45)),
+        self.assertEqual([((1, 1, 1, 1, 1), (45, 45, 45, 45),
+                           (45, 45, 45, 45, 45, 45, 45, 45)),
                           ((1, 1, 1, 1, 1), (45, 45), (45, 45, 45, 45)),
                           ((1, 1, 1, 1, 1), (45,), (45, 45))],
                          [level.precipitation.chunks for level in levels])
+
+
+class XarrayEncodingTest(unittest.TestCase):
+    """This test demonstrates how xarray uses a variable's encoding dict
+    when persisting data.
+    """
+
+    zarr_path = "test.zarr"
+    nc_path = "test.nc"
+
+    def tearDown(self) -> None:
+        rimraf(self.zarr_path)
+        rimraf(self.nc_path)
+
+    def test_zarr_uint8_encoding(self):
+        self.assert_persistence_ok(engine="zarr",
+                                   encoded_dtype=np.dtype("uint8"),
+                                   decoded_dtype=np.dtype("float32"))
+
+    def test_zarr_int16_encoding(self):
+        self.assert_persistence_ok(engine="zarr",
+                                   encoded_dtype=np.dtype("int16"),
+                                   decoded_dtype=np.dtype("float32"))
+
+    def test_zarr_uint32_encoding(self):
+        self.assert_persistence_ok(engine="zarr",
+                                   encoded_dtype=np.dtype("uint32"),
+                                   decoded_dtype=np.dtype("float64"))
+
+    def test_zarr_int64_encoding(self):
+        self.assert_persistence_ok(engine="zarr",
+                                   encoded_dtype=np.dtype("int64"),
+                                   decoded_dtype=np.dtype("float64"))
+
+    def test_nc_uint8_encoding(self):
+        self.assert_persistence_ok(engine="netcdf4",
+                                   encoded_dtype=np.dtype("uint8"),
+                                   decoded_dtype=np.dtype("float32"))
+
+    def test_nc_int16_encoding(self):
+        self.assert_persistence_ok(engine="netcdf4",
+                                   encoded_dtype=np.dtype("int16"),
+                                   decoded_dtype=np.dtype("float32"))
+
+    def test_nc_uint32_encoding(self):
+        self.assert_persistence_ok(engine="netcdf4",
+                                   encoded_dtype=np.dtype("uint32"),
+                                   decoded_dtype=np.dtype("float64"))
+
+    def test_nc_int64_encoding(self):
+        self.assert_persistence_ok(engine="netcdf4",
+                                   encoded_dtype=np.dtype("int64"),
+                                   decoded_dtype=np.dtype("float64"))
+
+    def assert_persistence_ok(self,
+                              engine: str,
+                              encoded_dtype: np.dtype,
+                              decoded_dtype: np.dtype,
+                              encoded_fill_value: int = 2,
+                              decoded_fill_value: float = np.nan):
+        a = xr.DataArray(np.array([1, np.nan, 3],
+                                  dtype=np.float64),
+                         dims="x")
+        a.encoding.update(_FillValue=encoded_fill_value,
+                          dtype=encoded_dtype)
+        ds = xr.Dataset(dict(A=a))
+
+        if engine == "zarr":
+            path = self.zarr_path
+            ds.to_zarr(path, mode="w")
+        else:
+            path = self.nc_path
+            ds.to_netcdf(path, mode="w")
+
+        # No kwargs
+        ds = xr.open_dataset(path, engine=engine)
+        np.testing.assert_equal(ds.A.values,
+                                np.array([1, decoded_fill_value, 3]))
+        self.assertEqual(decoded_dtype, ds.A.dtype)
+
+        # decode_cf=True
+        ds = xr.open_dataset(path, engine=engine, decode_cf=True)
+        np.testing.assert_equal(ds.A.values,
+                                np.array([1, decoded_fill_value, 3]))
+        self.assertEqual(decoded_dtype, ds.A.dtype)
+
+        # decode_cf=False
+        ds = xr.open_dataset(path, engine=engine, decode_cf=False)
+        np.testing.assert_equal(ds.A.values,
+                                np.array([1, encoded_fill_value, 3]))
+        self.assertEqual(encoded_dtype, ds.A.dtype)