Repository: corteva/rioxarray

Issue:
Incorrect missing value handling when _FillValue is float
I am working with [this raster file](https://github.com/corteva/rioxarray/files/5226354/no2_2019_01.tif.zip) which was generated (I believe) with R. 

```python
>>> import rioxarray
>>> raster = rioxarray.open_rasterio('no2_2019_01.tif')
>>> raster.attrs
{'STATISTICS_MAXIMUM': 0.00022005502250977,
 'STATISTICS_MEAN': nan,
 'STATISTICS_MINIMUM': 7.5976657853971e-07,
 'STATISTICS_STDDEV': nan,
 '_FillValue': -3.4e+38,
 'scale_factor': 1.0,
 'add_offset': 0.0,
 'grid_mapping': 'spatial_ref'}
```
The missing values are thus specified with the value `-3.4e+38`. When invoking `.plot()`, the following figure is generated:

```python
>>> raster.plot()
```
![Capture d’écran 2020-09-15 à 18 24 04](https://user-images.githubusercontent.com/16692099/93237635-ad160f00-f780-11ea-8399-57767b262e20.png)

The large negative value is not being ignored when plotting. I ended up manually fixing this issue by masking the array:

```python
>>> import numpy as np
>>> raster = raster.rio.write_nodata(np.float32(raster.rio.nodata))
>>> raster = raster.where(raster != raster.rio.nodata)
>>> raster = raster.rio.write_nodata(np.nan)
>>> raster.plot()
```
![Capture d’écran 2020-09-15 à 18 32 44](https://user-images.githubusercontent.com/16692099/93238539-ddaa7880-f781-11ea-8423-4b07534d548b.png)

I expected the latter figure to be produced without needing to manually change the array nor the attributes of the array. 


#### Environment Information

```python
>>> rioxarray.show_versions()
rioxarray (0.0.31) deps:
  rasterio: 1.1.5
    xarray: 0.16.0
      GDAL: 3.0.4

Other python deps:
     scipy: 1.5.2
    pyproj: 2.6.1.post1

System:
    python: 3.8.5 | packaged by conda-forge | (default, Aug 29 2020, 01:18:42)  [Clang 10.0.1 ]
executable: /Users/bilelomrani/opt/miniconda3/envs/raster/bin/python
   machine: macOS-10.15.6-x86_64-i386-64bit
```

#### Installation method
 - conda

#### Conda environment information (if you installed with conda):

<br/>
Environment (<code>conda list</code>):
<details>

```
$ conda list | grep -E "rasterio|xarray|gdal"
libgdal                   3.0.4               h242383b_10    conda-forge
rasterio                  1.1.5            py38hac06332_1    conda-forge
rioxarray                 0.0.31                     py_0    conda-forge
xarray                    0.16.0                     py_0    conda-forge
```
</details>

<br/>
Details about  <code>conda</code> and system ( <code>conda info</code> ):
<details>

```
$ conda info
active environment : raster
    active env location : /Users/bilelomrani/opt/miniconda3/envs/raster
            shell level : 2
       user config file : /Users/bilelomrani/.condarc
 populated config files : /Users/bilelomrani/.condarc
          conda version : 4.8.5
    conda-build version : not installed
         python version : 3.8.5.final.0
       virtual packages : __osx=10.15.6
       base environment : /Users/bilelomrani/opt/miniconda3  (writable)
           channel URLs : https://conda.anaconda.org/conda-forge/osx-64
                          https://conda.anaconda.org/conda-forge/noarch
                          https://repo.anaconda.com/pkgs/main/osx-64
                          https://repo.anaconda.com/pkgs/main/noarch
                          https://repo.anaconda.com/pkgs/r/osx-64
                          https://repo.anaconda.com/pkgs/r/noarch
          package cache : /Users/bilelomrani/opt/miniconda3/pkgs
                          /Users/bilelomrani/.conda/pkgs
       envs directories : /Users/bilelomrani/opt/miniconda3/envs
                          /Users/bilelomrani/.conda/envs
               platform : osx-64
             user-agent : conda/4.8.5 requests/2.24.0 CPython/3.8.5 Darwin/19.6.0 OSX/10.15.6
                UID:GID : 501:20
             netrc file : /Users/bilelomrani/.netrc
           offline mode : False
```
</details>


Discussion:

Patch:
diff --git a/docs/history.rst b/docs/history.rst
index 2a483959..fc6bac04 100644
--- a/docs/history.rst
+++ b/docs/history.rst
@@ -4,6 +4,7 @@ History
 Latest
 ------
 - BUG: Ensure transform correct in rio.clip without coords (pull #165)
+- BUG: Ensure the nodata value matches the dtype (pull #166)
 
 0.0.31
 ------
diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py
index 02f7cb8c..837b92ea 100644
--- a/rioxarray/rioxarray.py
+++ b/rioxarray/rioxarray.py
@@ -288,6 +288,21 @@ def _get_data_var_message(obj):
         return ""
 
 
+def _ensure_nodata_dtype(original_nodata, new_dtype):
+    """
+    Convert the nodata to the new datatype and raise warning
+    if the value of the nodata value changed.
+    """
+    original_nodata = float(original_nodata)
+    nodata = np.dtype(new_dtype).type(original_nodata)
+    if not np.isnan(nodata) and original_nodata != nodata:
+        warnings.warn(
+            f"The nodata value ({original_nodata}) has been automatically "
+            f"changed to ({nodata}) to match the dtype of the data."
+        )
+    return nodata
+
+
 class XRasterBase(object):
     """This is the base class for the GIS extensions for xarray"""
 
@@ -1073,6 +1088,7 @@ class RasterArray(XRasterBase):
         data_obj = self._get_obj(inplace=inplace)
         input_nodata = False if input_nodata is None else input_nodata
         if input_nodata is not False:
+            input_nodata = _ensure_nodata_dtype(input_nodata, self._obj.dtype)
             data_obj.rio.update_attrs(dict(_FillValue=input_nodata), inplace=True)
         else:
             new_vars = dict(data_obj.attrs)
@@ -1084,7 +1100,10 @@ class RasterArray(XRasterBase):
     @property
     def encoded_nodata(self):
         """Return the encoded nodata value for the dataset if encoded."""
-        return self._obj.encoding.get("_FillValue")
+        encoded_nodata = self._obj.encoding.get("_FillValue")
+        if encoded_nodata is None:
+            return None
+        return _ensure_nodata_dtype(encoded_nodata, self._obj.dtype)
 
     @property
     def nodata(self):
@@ -1117,6 +1136,7 @@ class RasterArray(XRasterBase):
             self._nodata = False
             return None
 
+        self._nodata = _ensure_nodata_dtype(self._nodata, self._obj.dtype)
         return self._nodata
 
     def reproject(
@@ -1696,6 +1716,15 @@ class RasterArray(XRasterBase):
                 "dtype",
             )
         }
+        rio_nodata = (
+            self.encoded_nodata if self.encoded_nodata is not None else self.nodata
+        )
+        if rio_nodata is not None:
+            # Ensure dtype of output data matches the expected dtype.
+            # This check is added here as the dtype of the data is
+            # converted right before writing.
+            rio_nodata = _ensure_nodata_dtype(rio_nodata, dtype)
+
         with rasterio.open(
             raster_path,
             "w",
@@ -1706,9 +1735,7 @@ class RasterArray(XRasterBase):
             dtype=dtype,
             crs=self.crs,
             transform=self.transform(recalc=recalc_transform),
-            nodata=(
-                self.encoded_nodata if self.encoded_nodata is not None else self.nodata
-            ),
+            nodata=rio_nodata,
             **out_profile,
         ) as dst:
 
@@ -1726,7 +1753,7 @@ class RasterArray(XRasterBase):
                     out_data = self._obj
                 if self.encoded_nodata is not None:
                     out_data = out_data.fillna(self.encoded_nodata)
-                data = out_data.astype(dtype).load().data
+                data = out_data.values.astype(dtype)
                 if data.ndim == 2:
                     dst.write(data, 1, window=window)
                 else:


Test patch:
diff --git a/test/integration/test_integration_merge.py b/test/integration/test_integration_merge.py
index 9889ff3e..8efcf805 100644
--- a/test/integration/test_integration_merge.py
+++ b/test/integration/test_integration_merge.py
@@ -81,6 +81,7 @@ def test_merge_arrays__res():
     assert merged.coords["band"].values == [1]
     assert sorted(merged.coords) == ["band", "spatial_ref", "x", "y"]
     assert merged.rio.crs == rds.rio.crs
+    assert_almost_equal(merged.attrs.pop("_FillValue"), rds.attrs.pop("_FillValue"))
     assert merged.attrs == rds.attrs
     assert_almost_equal(nansum(merged), 13556564)
 
diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py
index 717bef75..fce56e66 100644
--- a/test/integration/test_integration_rioxarray.py
+++ b/test/integration/test_integration_rioxarray.py
@@ -1271,6 +1271,35 @@ def test_to_raster__dataset__different_nodata(tmpdir):
             rds.rio.to_raster(str(tmp_raster))
 
 
+@pytest.mark.parametrize("windowed", [True, False])
+def test_to_raster__different_dtype(tmp_path, windowed):
+    test_da = xarray.DataArray(
+        numpy.zeros((5, 5)),
+        dims=("y", "x"),
+        coords={"y": numpy.arange(1, 6), "x": numpy.arange(2, 7)},
+    )
+    test_da.values[1, 1] = -1.1
+    test_nd = test_da.rio.write_nodata(-1.1)
+    test_nd.rio.write_transform(
+        Affine.from_gdal(425047, 3.0, 0.0, 4615780, 0.0, -3.0), inplace=True
+    )
+    test_nd.rio.write_crs("EPSG:4326", inplace=True)
+    tmpfile = tmp_path / "dtype.tif"
+    with pytest.warns(
+        UserWarning,
+        match=(
+            r"The nodata value \(-1.1\) has been automatically changed to "
+            r"\(255\) to match the dtype of the data."
+        ),
+    ):
+        test_nd.rio.to_raster(tmpfile, dtype=numpy.uint8, windowed=windowed)
+    xds = rioxarray.open_rasterio(tmpfile)
+    assert str(xds.dtype) == "uint8"
+    assert xds.attrs["_FillValue"] == 255
+    assert xds.rio.nodata == 255
+    assert xds.squeeze().values[1, 1] == 255
+
+
 def test_missing_spatial_dimensions():
     test_ds = xarray.Dataset()
     with pytest.raises(DimensionError):
@@ -1647,6 +1676,45 @@ def test_nodata_writer__remove():
     assert not test_nd.attrs
 
 
+@pytest.mark.parametrize("nodata", [-1.1, "-1.1"])
+def test_nodata_writer__different_dtype(nodata):
+    test_da = xarray.DataArray(
+        numpy.zeros((5, 5), dtype=int),
+        dims=("y", "x"),
+        coords={"y": numpy.arange(1, 6), "x": numpy.arange(2, 7)},
+    )
+    with pytest.warns(
+        UserWarning,
+        match=(
+            r"The nodata value \(-1.1\) has been automatically changed to "
+            r"\(-1\) to match the dtype of the data."
+        ),
+    ):
+        test_nd = test_da.rio.write_nodata(nodata)
+    assert not test_da.attrs
+    assert test_nd.attrs["_FillValue"] == -1
+    assert test_nd.rio.nodata == -1
+
+
+@pytest.mark.parametrize("nodata", [-1.1, "-1.1"])
+def test_nodata_reader__different_dtype(nodata):
+    test_da = xarray.DataArray(
+        numpy.zeros((5, 5), dtype=numpy.uint8),
+        dims=("y", "x"),
+        coords={"y": numpy.arange(1, 6), "x": numpy.arange(2, 7)},
+        attrs={"_FillValue": nodata},
+    )
+    assert test_da.attrs["_FillValue"] == nodata
+    with pytest.warns(
+        UserWarning,
+        match=(
+            r"The nodata value \(-1.1\) has been automatically changed to "
+            r"\(255\) to match the dtype of the data."
+        ),
+    ):
+        assert test_da.rio.nodata == 255
+
+
 def test_isel_window():
     with rioxarray.open_rasterio(
         os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc")