Repository: pydata/xarray

Issue:
DataArray read from netcdf with unexpected type
#### Code Sample, a copy-pastable example if possible

```python
import numpy as np
import h5netcdf

filename = "mask_and_scale_float32.nc"

with h5netcdf.File(filename, 'w') as f:
    f.dimensions = {'x': 5}
    v = f.create_variable('hello', ('x',), dtype=np.uint16)
    v[:] = np.ones(5, dtype=np.uint16)
    v[0] = np.uint16(65535)
    v.attrs['_FillValue'] = np.uint16(65535)
    v.attrs['scale_factor'] = np.float32(2)
    v.attrs['add_offset'] = np.float32(0.5)


import xarray as xr

v = xr.open_dataset(filename, mask_and_scale=True)['hello']
print(v.dtype)
```

#### Problem description

The `scale_factor` and `add_offset` being float32, I would expect the result from loading to be a float32 array. However, we get a float64 array instead. A float32 array for a very large dataset is better for faster computations.

#### Expected Output

float32

#### Output of ``xr.show_versions()``

<details>
INSTALLED VERSIONS
------------------
commit: None
python: 2.7.5.final.0
python-bits: 64
OS: Linux
OS-release: 3.10.0-693.11.6.el7.x86_64
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_GB.UTF-8
LOCALE: None.None

xarray: 0.10.0
pandas: 0.21.0
numpy: 1.13.3
scipy: 0.18.1
netCDF4: 1.1.8
h5netcdf: 0.4.2
Nio: None
bottleneck: None
cyordereddict: None
dask: 0.16.0+37.g1fef002
matplotlib: 2.1.0
cartopy: None
seaborn: None
setuptools: 38.2.4
pip: 9.0.1
conda: None
pytest: 3.1.3
IPython: 5.5.0
sphinx: 1.3.6

</details>


Discussion:
Zac-HD: Thanks - I was actually writing up an issue and decided it would be easier to demonstrate the proposed fix in a PR, but I'll open an issue first next time.

The checkbox about flake8 could be removed from the issue template now - since #1824 we run flake8 on everything in CI so if tests pass flake8 is passing too.

Re: tests: what do you (and @shoyer) think about using [Hypothesis](https://hypothesis.readthedocs.io) for some property-based tests of variable coding?  "encoding then decoding is a no-op" is a classic property 😄  Upside, more powerful and better at finding edge cases; downside slower simply because it checks more cases (a configurable number).
shoyer: As mentioned in https://github.com/pydata/xarray/issues/1842, maybe we should also make a point not to upcast float32 input?

One possible concern with changing precision from float64 -> float32 is that some reduce operations like mean could become due to lower precision. So it's a good think @fujiisoup wrote https://github.com/pydata/xarray/pull/1841 so we can specify `dtype` in reductions :).

With regards to Hypothesis: I haven't used it, but it does seem very intriguing. I'm sure that it could turn up quite a few bugs in xarray. Would it make sense to add it as an optional dependency to the test suite? Even if we use Hypothesis to cover more edge cases, I think we will still want normal test coverage for most behavior.
Zac-HD: Added tests; float32 not upcast, float16 (["intended for storage of many floating-point values where higher precision is not needed, not for performing arithmetic"](https://en.wikipedia.org/wiki/Half-precision_floating-point_format)) *is* upcast but only to float32.

I'll open a new issue to add a basic suite of property-based tests to Xarray 😄
jhamman: This all looks good now. I merged in master and resolved a small conflict in the docs. I'll merge after the tests pass. 

Patch:
diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md
index 3a324810..5e9aa06f 100644
--- a/.github/PULL_REQUEST_TEMPLATE.md
+++ b/.github/PULL_REQUEST_TEMPLATE.md
@@ -1,5 +1,4 @@
  - [ ] Closes #xxxx (remove if there is no corresponding issue, which should only be the case for minor changes)
  - [ ] Tests added (for all bug fixes or enhancements)
  - [ ] Tests passed (for all non-documentation changes)
- - [ ] Passes ``git diff upstream/master **/*py | flake8 --diff`` (remove if you did not edit any Python files)
  - [ ] Fully documented, including `whats-new.rst` for all changes and `api.rst` for new API (remove if this change should not be visible to users, e.g., if it is an internal clean-up, or if this is part of a larger project that will be documented later)
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 495b088c..48de04ed 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -57,6 +57,11 @@ Enhancements
 - :py:func:`~plot.line()` learned to draw multiple lines if provided with a
   2D variable.
   By `Deepak Cherian <https://github.com/dcherian>`_.
+- Reduce memory usage when decoding a variable with a scale_factor, by
+  converting 8-bit and 16-bit integers to float32 instead of float64
+  (:pull:`1840`), and keeping float16 and float32 as float32 (:issue:`1842`).
+  Correspondingly, encoded variables may also be saved with a smaller dtype.
+  By `Zac Hatfield-Dodds <https://github.com/Zac-HD>`_.
 
 .. _Zarr: http://zarr.readthedocs.io/
 
@@ -76,14 +81,12 @@ Bug fixes
 - Fixed encoding of multi-dimensional coordinates in
   :py:meth:`~Dataset.to_netcdf` (:issue:`1763`).
   By `Mike Neish <https://github.com/neishm>`_.
-
 - Fixed chunking with non-file-based rasterio datasets (:issue:`1816`) and
   refactored rasterio test suite.
   By `Ryan Abernathey <https://github.com/rabernat>`_
 - Bug fix in open_dataset(engine='pydap') (:issue:`1775`)
   By `Keisuke Fujii <https://github.com/fujiisoup>`_.
-
-- Bug fix in vectorized assignment  (:issue:`1743`, `1744`).
+- Bug fix in vectorized assignment  (:issue:`1743`, :issue:`1744`).
   Now item assignment to :py:meth:`~DataArray.__setitem__` checks
 - Bug fix in vectorized assignment  (:issue:`1743`, :issue:`1744`).
   Now item assignment to :py:meth:`DataArray.__setitem__` checks
diff --git a/xarray/coding/variables.py b/xarray/coding/variables.py
index bf2ded8b..5d32970e 100644
--- a/xarray/coding/variables.py
+++ b/xarray/coding/variables.py
@@ -205,6 +205,25 @@ def _scale_offset_decoding(data, scale_factor, add_offset, dtype):
     return data
 
 
+def _choose_float_dtype(dtype, has_offset):
+    """Return a float dtype that can losslessly represent `dtype` values."""
+    # Keep float32 as-is.  Upcast half-precision to single-precision,
+    # because float16 is "intended for storage but not computation"
+    if dtype.itemsize <= 4 and np.issubdtype(dtype, np.floating):
+        return np.float32
+    # float32 can exactly represent all integers up to 24 bits
+    if dtype.itemsize <= 2 and np.issubdtype(dtype, np.integer):
+        # A scale factor is entirely safe (vanishing into the mantissa),
+        # but a large integer offset could lead to loss of precision.
+        # Sensitivity analysis can be tricky, so we just use a float64
+        # if there's any offset at all - better unoptimised than wrong!
+        if not has_offset:
+            return np.float32
+    # For all other types and circumstances, we just use float64.
+    # (safe because eg. complex numbers are not supported in NetCDF)
+    return np.float64
+
+
 class CFScaleOffsetCoder(VariableCoder):
     """Scale and offset variables according to CF conventions.
 
@@ -216,7 +235,8 @@ class CFScaleOffsetCoder(VariableCoder):
         dims, data, attrs, encoding = unpack_for_encoding(variable)
 
         if 'scale_factor' in encoding or 'add_offset' in encoding:
-            data = data.astype(dtype=np.float64, copy=True)
+            dtype = _choose_float_dtype(data.dtype, 'add_offset' in encoding)
+            data = data.astype(dtype=dtype, copy=True)
             if 'add_offset' in encoding:
                 data -= pop_to(encoding, attrs, 'add_offset', name=name)
             if 'scale_factor' in encoding:
@@ -230,7 +250,7 @@ class CFScaleOffsetCoder(VariableCoder):
         if 'scale_factor' in attrs or 'add_offset' in attrs:
             scale_factor = pop_to(attrs, encoding, 'scale_factor', name=name)
             add_offset = pop_to(attrs, encoding, 'add_offset', name=name)
-            dtype = np.float64
+            dtype = _choose_float_dtype(data.dtype, 'add_offset' in attrs)
             transform = partial(_scale_offset_decoding,
                                 scale_factor=scale_factor,
                                 add_offset=add_offset,


Test patch:
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index 013131a1..9c4db8b4 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -65,7 +65,7 @@ def open_example_dataset(name, *args, **kwargs):
 
 
 def create_masked_and_scaled_data():
-    x = np.array([np.nan, np.nan, 10, 10.1, 10.2])
+    x = np.array([np.nan, np.nan, 10, 10.1, 10.2], dtype=np.float32)
     encoding = {'_FillValue': -1, 'add_offset': 10,
                 'scale_factor': np.float32(0.1), 'dtype': 'i2'}
     return Dataset({'x': ('t', x, {}, encoding)})
@@ -80,7 +80,7 @@ def create_encoded_masked_and_scaled_data():
 def create_unsigned_masked_scaled_data():
     encoding = {'_FillValue': 255, '_Unsigned': 'true', 'dtype': 'i1',
                 'add_offset': 10, 'scale_factor': np.float32(0.1)}
-    x = np.array([10.0, 10.1, 22.7, 22.8, np.nan])
+    x = np.array([10.0, 10.1, 22.7, 22.8, np.nan], dtype=np.float32)
     return Dataset({'x': ('t', x, {}, encoding)})
 
 
diff --git a/xarray/tests/test_coding.py b/xarray/tests/test_coding.py
index d1b54fca..a6faea87 100644
--- a/xarray/tests/test_coding.py
+++ b/xarray/tests/test_coding.py
@@ -1,5 +1,7 @@
 import numpy as np
 
+import pytest
+
 import xarray as xr
 from xarray.core.pycompat import suppress
 from xarray.coding import variables
@@ -36,3 +38,15 @@ def test_coder_roundtrip():
     coder = variables.CFMaskCoder()
     roundtripped = coder.decode(coder.encode(original))
     assert_identical(original, roundtripped)
+
+
+@pytest.mark.parametrize('dtype', 'u1 u2 i1 i2 f2 f4'.split())
+def test_scaling_converts_to_float32(dtype):
+    original = xr.Variable(('x',), np.arange(10, dtype=dtype),
+                           encoding=dict(scale_factor=10))
+    coder = variables.CFScaleOffsetCoder()
+    encoded = coder.encode(original)
+    assert encoded.dtype == np.float32
+    roundtripped = coder.decode(encoded)
+    assert_identical(original, roundtripped)
+    assert roundtripped.dtype == np.float32