Ability to force float64 instead of float32 issue #2304 #2751

DevDaoud · 2019-02-07T15:10:58Z

Ability to promote to float64 instead of float32 when dealing with int variables with scale_factor.
added parameter force_promote_float64 (False by default) to open_dataset and open_dataarray that enables this behaviour when True.

Closes float32 instead of float64 when decoding int16 with scale_factor netcdf var using xarray #2304
Tests added
Fully documented, including whats-new.rst for all changes and api.rst for new API

When dealing with int variables with scale factor xarray transform them to float32 by default. This enhancement aims to give the possibity to force that to float64 when needed. The behaviour is unchanged when not specifying the parameter force_promote_float64.

DevDaoud · 2019-02-08T10:41:34Z

@shoyer did you have a look at this ?

shoyer · 2019-02-08T16:15:11Z

@daoudjahdou sorry I missed your comment over in #2304. Threading the new parameter down though other functions is definitely the preferred approach here -- mutating global variables makes it much harder to follow control flow. Though I would probably stop short of modifying maybe_promote and would put this logic inside CFMaskCoder and/or CFScaleOffsetCoder in xarray.coding.variables.

…oder

DevDaoud · 2019-02-12T09:21:59Z

@shoyer the logic is now propagated down to _choose_float_dtype inside CFScaleOffsetCoder, please let me know what you think.

shoyer · 2019-02-12T18:24:51Z

@daoudjahdou what do you think of @magau's proposed solution over in #2304 (comment)? Checking the dtypes of of add_offset and scale_factor might be a better alternative than requiring setting a custom option.

DevDaoud · 2019-02-13T08:01:39Z

@shoyer yes sure i'll update the pull request with the mentioned modifications.

http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/ch08.html

DevDaoud · 2019-02-18T14:25:00Z

@shoyer now scale_factor and add_offset are taken into account when encoding and decoding data.
if none of them is present or if they're not a subtype of np.generic the old behaviour takes place.

shoyer · 2019-02-19T05:03:29Z

xarray/coding/variables.py

-            dtype = _choose_float_dtype(data.dtype, 'add_offset' in encoding)
+            scale_factor = pop_to(attrs, encoding, 'scale_factor', name=name)
+            add_offset = pop_to(attrs, encoding, 'add_offset', name=name)
+            dtype = _choose_float_dtype(data.dtype, scale_factor, add_offset)
            data = data.astype(dtype=dtype, copy=True)
            if 'add_offset' in encoding:


This line (and the line below for scale_factor) needs to be updated to use the add_offset variable if it isn't None rather than checking for it in encoding again. (pop_to removes an item from the dict in which it's found.)

shoyer · 2019-02-19T05:09:10Z

xarray/coding/variables.py

+    # from the variable (containing the packed data)
+    # then the unpacked data should match
+    # the type of these attributes, which must both be of type float or both
+    # be of type double. An additional restriction in this case is that


What should we do if add_offset/scale_factor are not floats?

Right now in your implementation, if scale_factor is mistakenly np.int64(10) then the data would get decoded as int64, which seems like a bad idea (e.g., if the data is some float type). I think we would should probably ignore the type of scale_factor/add_offset in these cases instead.

Hi @shoyer! Sorry for my delayed intervenience... First of all I'd like to thank and congratulate you for your efficiency and say that I'm very glad to see that you've taken into account my considerations about this topic.

About the comment above I'd say that, since the CFScaleOffsetCoder will always multiply the raw data by the scale_factor, in case the scale_factor is an integer the decoded values will always result into integers too (unless the raw data dtype is float, as you've mentioned).
Even in these cases I'd return the decoded values with the same dtype as the scale_factor / add_offset, which is the expected behavior (data providers must be aware...). Users can always use open_dataset with decode_cf=False and call the decode_cf_variable function later, after fixing the variable.attrs['scale_factor / add_offset'] dtypes.
A possible solution is add the scale_factor and add_offset as arguments of the open_dataset function, to overwrite the original ones, which wold receive labeled values like: scale_factor={'var1': whatever, 'var2': whatever}, add_offset: {'var1': whatever, 'var2': whatever}. This wold also resolve the #2304 issue (with a better approach, in my opinion).

By the way, I've also noted that the CFScaleOffsetCoder.encode changes the values.dtype before applying the inverse operation (i.e the dividing by the scale_factor), this will result in truncated values. I've already faced this situation, with my own netcdf encoding processes, and would advice you to only change the dtype after applying the transformation and use the numpy.round function (for integers encoding) to avoid the truncation artifacts.

Once again, sorry for my extensive comments and not opening an issue, as I should, but my intention is that you evaluate my considerations and take your own decisions with your colleagues, and keep the nice work you all have being done. 👍

Cheers!

shoyer · 2019-02-19T05:11:17Z

xarray/coding/variables.py

+    # We first return scale_factor type
+    # as multiplying takes priority over
+    # adding values
+    if dtype is not type(scale_factor) and \


Do we actually need to check whether the type of scale_factor matches dtype? I think the logic works the same either way.

shoyer · 2019-02-19T05:12:28Z

xarray/coding/variables.py

+    # as multiplying takes priority over
+    # adding values
+    if dtype is not type(scale_factor) and \
+            isinstance(scale_factor, np.generic):


nit: per PEP8, please prefer using parentheses () to split across multiple lines instead of \

shoyer · 2019-02-19T05:15:29Z

xarray/tests/test_backends.py

@@ -1156,6 +1157,25 @@ def test_mask_and_scale(self):
                expected = create_masked_and_scaled_data()
                assert_identical(expected, ds)

+    def test_mask_and_scale_with_float64_scale_factor(self):


It would be nice to test all of the edge cases just using CFScaleOffsetCoder, e.g.,

what happens if add_offset and scale_factor have different types?

what if they're integers instead of floats?

See the tests in tests/test_coding.py for examples.

shoyer · 2019-02-19T05:22:57Z

xarray/coding/variables.py

+
+    # We first return scale_factor type
+    # as multiplying takes priority over
+    # adding values


I'm a little nervous about making this assumption. The NumPy casting rules for a * x + b would give the result the largest dtype of any of a, x and b. Maybe we should similarly pick whichever of scale_factor, data and add_offset has the largest float dtype?

…hdou/xarray into force_float64_issue_2304

DevDaoud · 2019-02-22T06:38:32Z

@shoyer i changed the implementation, and took into consideration your comments.
now returning largest type takes place only when decoding.
i added a test with all types of scale factor, add_offset and variable.

DevDaoud · 2019-02-28T08:09:58Z

@shoyer did you have a look at this ?

shoyer · 2019-03-04T05:48:18Z

xarray/coding/variables.py

+                 np.dtype(type(add_offset)),
+                 np.dtype(dtype)]
+        # scaled_type should be the largest type we find
+        scaled_dtype = max(types)


Use xarray.core.dtypes.result_type() here instead of max().

I don't know how NumPy defines comparison between dtypes and I can't find it documented anywhere, so I'm a little nervous about relying on it.

@shoyer i've tested it with all numpy dtypes pair combinations and it turns out that it always gives the largest, i can for sure modify it

shoyer · 2019-03-04T05:52:16Z

xarray/coding/variables.py

+    # It is not advised to unpack an int into a float as there
+    # is a potential precision loss.
+
+    if mode is 'decoding':


Rather than a mode argument, could you write two separate functions, _choose_encoding_float_dtype and _choose_decoding_float_dtype? It would be OK for one to call the other, but mode arguments are mild code smell.

shoyer · 2019-03-04T05:57:31Z

xarray/tests/test_backends.py

@@ -685,7 +685,9 @@ def test_roundtrip_mask_and_scale(self, decoded_fn, encoded_fn):
        with self.roundtrip(decoded) as actual:
            for k in decoded.variables:
                assert (decoded.variables[k].dtype
-                        == actual.variables[k].dtype)
+                        == actual.variables[k].dtype or
+                        (decoded.variables[k].dtype == np.float32 and


It's not clear to me why this test needs to change. Can you give an example of what behavior changed?

@shoyer
Here is my understanding of what xarray does :
In all cases encoding returns either np.float32 or np.float64
Encoding only cares about existence of add_offset when the variable is an integer
If the variable is a float then the encoded dtype is np.float32
If the variable is an integer with add_offset then the encoded dtype np.float64
If the variable is an integer with no add_affset the encoded dtype np.float32
In all other cases the encoded dtype is np.float64

Here is my understanding of what cf-conventions imply:
decoding is the equivalent of unpacking mentioned in the cf-convention
in all cases decoding takes the encoded dtype which is either np.float32 or np.float64
then decoded type is the largest between scale_factor, add_offset and encoded type and not original variable

That being said in that test there are cases where
original type, scale_factor, add_offset are on the following configuration:
np.float32, np.int64, np.float32
then the encoded type is np.float32 as the variable is np.float32
when decoding the largest (dtypes.result_type()) between the three (encoded, sf, ao) is np.float64

and that's why "decoded" (orignal) in the test is np.float32 and actual (roundtripped) is np.float64

In my next commit i took into consideration your remark and i wrote every possible case and what to expect as an encoded or decoded dtype.

shoyer · 2019-03-04T05:58:49Z

xarray/tests/test_coding.py

+    # We make sure that roundtripped is larger than
+    # the original
+    assert roundtripped.dtype.itemsize >= original.dtype.itemsize
+    assert (roundtripped.dtype is np.dtype(np.float64)


nit: it's safer to use == instead of is for comparing dtype objects

shoyer · 2019-03-04T06:07:32Z

xarray/tests/test_coding.py

+                                        np.int64(10), np.uint8(10),
+                                        np.uint16(10), np.uint32(10),
+                                        np.uint64(10), np.uint64(10)])
+def test_scaling_according_to_cf_convention(dtype, scale_factor, add_offset):


This is a great set of invariant checks! But it would still be nice to see a more specific list of unit tests verifying that the decoded dtype is exactly as expected for a set of known inputs, e.g.,

@pytest.mark.parametrize( ('input_dtype', 'scale_factor', 'add_offset', 'expected_dtype'), [ (np.float32, 1, 0, np.float32), (np.float32, 1, 10, np.float64), ... ] ) def test_cfscaleoffset_decoding_dtype(input_dtype, scale_factor, add_offset, expected_dtype): original = xr.Variable(...) decoded = variables.CFScaleOffsetCoder().decode(original) assert decoded.dtype == expected_dtype

This would us confidence that it's handling all the edge cases properly, e.g., only giving a bigger dtype when truly warranted.

DevDaoud · 2019-03-08T10:29:32Z

@shoyer tests are failing but it does'nt seem to be coming from this PR, i saw the same error on other PRs as well, my tests were working fine until i made a git pull.

shoyer · 2019-03-15T05:42:18Z

@daoudjahdou thanks for sticking with this.

I'm concerned about how lose the "roundtrip" invariant in some edge cases, where the user did not specifically indicate the desired dtype but rather used a Python type. Right now xarray seems to use inconsistent rules for casting attribute types, based on the netCDF backend:

In [8]: for engine in ['scipy', 'netcdf4', 'h5netcdf']:
   ...:     ds = xarray.Dataset({'x': np.float32(0)})
   ...:     encoding = {'x': {'scale_factor': 1, 'add_offset': 1}}
   ...:     ds.to_netcdf(f'test-{engine}.nc', engine=engine, encoding=encoding)
   ...:

In [9]: ! ncdump -h test-scipy.nc
netcdf test-scipy {
variables:
	float x ;
		x:add_offset = 1 ;
		x:scale_factor = 1 ;
		x:_FillValue = NaNf ;
}

In [10]: ! ncdump -h test-netcdf4.nc
netcdf test-netcdf4 {
variables:
	float x ;
		x:_FillValue = NaNf ;
		x:add_offset = 1LL ;
		x:scale_factor = 1LL ;
}

In [11]: ! ncdump -h test-h5netcdf.nc
netcdf test-h5netcdf {
variables:
	float x ;
		x:_FillValue = NaNf ;
		x:add_offset = 1LL ;
		x:scale_factor = 1LL ;
}

At least for add_offset and scale_factor, I think we should be more careful when determining the attribute dtype. Python's int and float don't have any intrinsic dtype. If the attributes are specified as native Python types (not NumPy types) and can be safely cast into the dtype of the data, we should probably use that instead.

DevDaoud · 2019-03-15T14:00:34Z

@shoyer do you mean that we consider that by default when we deal with Python's int and float we cast them to np.int64 and np.float64 ?

shoyer · 2019-03-16T03:50:00Z

@shoyer do you mean that we consider that by default when we deal with Python's int and float we cast them to np.int64 and np.float64 ?

Currently we cast to int32/float64 for netCDF3 files (which don't support int64) and int64/float64 for netCDF4 files.

That default behavior doesn't really make sense if the user provided builtin Python numbers for these attributes.

For example, if I write something like the following, the data for x should probably decode to float32:

ds = xarray.Dataset({'x': np.float32(0)})
encoding = {'x': {'scale_factor': 1, 'add_offset': 1, 'dtype': 'int8'}}
ds.to_netcdf('test.nc', encoding=encoding)

With the current implementation in this PR, I think the data will come back as float64.

To get something that will decode as float64 without the original data being float64, the user should need to intentionally choose dtypes for scale_factor and add_offsert, e.g.,

ds = xarray.Dataset({'x': np.float32(0)})
encoding = {'x': {'scale_factor': np.float64(1), 'add_offset': np.float64(1), 'dtype': 'int8'}}
ds.to_netcdf('test.nc', encoding=encoding)

Does this make sense to you? From my perspective, two goals for this change should be:

Adhere to CF conventions for deciding how to decode data from disk.
Provide a consistent interface for xarray users, so that by default data is restored in the same form that it was written.

DevDaoud · 2019-04-04T09:55:14Z

@shoyer sorry for the delayed response.
dtypes.result_type(1, np.float32(1)) returns dtype('float64')
That's what makes this behaviour for Python's int and float
Keeping the consistency would then require testing if scale_factor*var_dtype + add_offset fit in var_dtype in case of Python's int and float.
Correct me if i'm wrong but this is a bit hard to do without knowing max,min values to avoid overlapping
Evaluating those could break performance if only used for encoding and decoding.
Do have any idea how this could be acheived ? or is it simpler to keep it as it is ?

shoyer · 2019-04-08T05:19:47Z

Correct me if i'm wrong but this is a bit hard to do without knowing max,min values to avoid overlapping
Evaluating those could break performance if only used for encoding and decoding.

Right, we definitely don't want to inspect the data to verify if it fits.

I think we will will need explicitly checks for isinstance(x, (int, float)) (for scale_factor and add_offset) to identify these cases in the encode() method. Let me try to come up with a concrete suggestion... EDIT: See my suggestion below. I don't think we want these explicit checks, after all.

shoyer · 2019-04-08T05:41:35Z

xarray/coding/variables.py

@@ -217,9 +258,9 @@ class CFScaleOffsetCoder(VariableCoder):

    def encode(self, variable, name=None):
        dims, data, attrs, encoding = unpack_for_encoding(variable)
-
        if 'scale_factor' in encoding or 'add_offset' in encoding:


Alternative logic for this clause that I think might work:

if 'scale_factor' in encoding or 'add_offset' in encoding: scale_factor = encoding.pop('scale_factor', 1) add_offset = encoding.pop('add_offset, 0) if 'dtype' in encoding: dtype = pop_to(encoding, attrs, 'dtype', name=name) else: dtype = _choose_encoding_float_dtype(data.dtype, add_offset != 0) # save scale_factor and add_offset with dtype matching the decoded # data, per CF conventions safe_setitem(attrs, 'scale_factor', data.dtype.type(scale_factor), name) safe_setitem(attrs, 'add_offset', data.dtype.type(add_offset), name) # apply scale/offset encoding data = data.astype(dtype=dtype, copy=True) if add_offset: data -= add_offset if scale_factor != 1: # multiplication is faster than division data *= 1/scale_factor

The idea here is to always write the scale_factor and add_offset attributes according to CF conventions, which ensures that data always gets read out the right way. This does remove some user flexibility, but xarray is OK with being opinionated about how it writes metadata -- it is not a general purpose tool for writing completely flexible netCDF files.

shoyer · 2019-04-08T05:48:28Z

xarray/coding/variables.py

+
+        if add_offset is None:
+            types = (np.dtype(type(scale_factor)),
+                     np.dtype(dtype))


It would be nice if these five lines creating types handled add_offset and scale_factor in a symmetric way, e.g.,

types = [np.dtype(dtype)] if scale_factor is not None: types.append(np.dtype(type(scale_offset)) if add_offset is not None: types.append(np.dtype(type(add_offset))

shoyer · 2019-04-08T05:50:08Z

xarray/coding/variables.py

+    # It is not advised to unpack an int into a float as there
+    # is a potential precision loss.
+
+    if scale_factor or add_offset:


It's safer to use is not None for checking default values:

Suggested change

if scale_factor or add_offset:

if scale_factor is not None or add_offset is not None:

shoyer · 2019-04-08T05:58:01Z

xarray/coding/variables.py

+                and np.issubdtype(scaled_dtype, np.floating)):
+            return scaled_dtype
+
+    return _choose_encoding_float_dtype(dtype, add_offset is not None)


Potentially we could use the add_offset value here -- 0 and None are equivalent as far as picking an automatic float dtype is concerned, indicating "no offset":

Suggested change

return _choose_encoding_float_dtype(dtype, add_offset is not None)

return _choose_encoding_float_dtype(dtype, add_offset)

DevDaoud · 2019-05-23T07:48:39Z

@shoyer i've tested the solution provided, it works like a charm with my tests however many tests are broken on test_backends.py with cases where we lose precision i'll give you more detail.

fmaussion · 2019-11-26T12:23:44Z

I've noted that the very recent ERA5 data is also cast to float32 at decoding, although the scale and offset params are float64. I think it would be a very valuable addition to xarray to do this properly: @daoudjahdou are you still willing to work on this PR? Otherwise I'll try to pick this up when time permits.

Daoud Jahdou added 5 commits February 7, 2019 15:10

updated documentation pydata#2304

093a338

Fixed typos doc pydata#2304

47d4207

Minor modification to documentation

3804d56

Updated doc

0142546

Propagated the force_promote_float64 parameter down to CFScaleOffsetC…

5b421c8

…oder

removing dtypes constant mutation

c1a2a03

implementing cf convention regarding scale_factor and add_offset types

b05193d

http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/ch08.html

shoyer reviewed Feb 19, 2019

View reviewed changes

DevDaoud and others added 4 commits February 21, 2019 17:44

Merge branch 'master' into force_float64_issue_2304

b42c744

Taking into account comments and adding edge case tests

a94f2af

Merge branch 'force_float64_issue_2304' of https://github.com/daoudja…

c394509

…hdou/xarray into force_float64_issue_2304

removed unused code

99599ed

shoyer reviewed Mar 4, 2019

View reviewed changes

DevDaoud and others added 4 commits March 7, 2019 18:45

Merge branch 'master' into force_float64_issue_2304

1de0884

Adding more tests with expected outputs

638332c

Renaming test_coding methods

41da205

Replacing array by tuple

29d191a

Merge branch 'master' into force_float64_issue_2304

c535df2

shoyer reviewed Apr 8, 2019

View reviewed changes

keewis added topic-backends and removed topic-backends labels Nov 24, 2020

mankoff mentioned this pull request Oct 6, 2022

Improved CF decoding #6812

Open

2 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Ability to force float64 instead of float32 issue #2304 #2751

Ability to force float64 instead of float32 issue #2304 #2751

DevDaoud commented Feb 7, 2019 •

edited

DevDaoud commented Feb 8, 2019

shoyer commented Feb 8, 2019

DevDaoud commented Feb 12, 2019 •

edited

shoyer commented Feb 12, 2019

DevDaoud commented Feb 13, 2019

DevDaoud commented Feb 18, 2019 •

edited

shoyer Feb 19, 2019

shoyer Feb 19, 2019

magau Feb 21, 2019 •

edited

shoyer Feb 19, 2019

shoyer Feb 19, 2019

shoyer Feb 19, 2019

shoyer Feb 19, 2019

DevDaoud commented Feb 22, 2019 •

edited

DevDaoud commented Feb 28, 2019

shoyer Mar 4, 2019

DevDaoud Mar 4, 2019

shoyer Mar 4, 2019

shoyer Mar 4, 2019

DevDaoud Mar 8, 2019 •

edited

shoyer Mar 4, 2019

shoyer Mar 4, 2019

DevDaoud commented Mar 8, 2019

shoyer commented Mar 15, 2019

DevDaoud commented Mar 15, 2019

shoyer commented Mar 16, 2019

DevDaoud commented Apr 4, 2019

shoyer commented Apr 8, 2019 •

edited

shoyer Apr 8, 2019

shoyer Apr 8, 2019

shoyer Apr 8, 2019

shoyer Apr 8, 2019

DevDaoud commented May 23, 2019

fmaussion commented Nov 26, 2019

	if scale_factor or add_offset:
	if scale_factor is not None or add_offset is not None:

	return _choose_encoding_float_dtype(dtype, add_offset is not None)
	return _choose_encoding_float_dtype(dtype, add_offset)

Ability to force float64 instead of float32 issue #2304 #2751

Are you sure you want to change the base?

Ability to force float64 instead of float32 issue #2304 #2751

Conversation

DevDaoud commented Feb 7, 2019 • edited

DevDaoud commented Feb 8, 2019

shoyer commented Feb 8, 2019

DevDaoud commented Feb 12, 2019 • edited

shoyer commented Feb 12, 2019

DevDaoud commented Feb 13, 2019

DevDaoud commented Feb 18, 2019 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

magau Feb 21, 2019 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

DevDaoud commented Feb 22, 2019 • edited

DevDaoud commented Feb 28, 2019

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

DevDaoud Mar 8, 2019 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

DevDaoud commented Mar 8, 2019

shoyer commented Mar 15, 2019

DevDaoud commented Mar 15, 2019

shoyer commented Mar 16, 2019

DevDaoud commented Apr 4, 2019

shoyer commented Apr 8, 2019 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

DevDaoud commented May 23, 2019

fmaussion commented Nov 26, 2019

DevDaoud commented Feb 7, 2019 •

edited

DevDaoud commented Feb 12, 2019 •

edited

DevDaoud commented Feb 18, 2019 •

edited

magau Feb 21, 2019 •

edited

DevDaoud commented Feb 22, 2019 •

edited

DevDaoud Mar 8, 2019 •

edited

shoyer commented Apr 8, 2019 •

edited