Repository: corteva/rioxarray

Issue:
Better nodata handling
Currently the nodata property searches for nodata in `_FillValue`. But it is not always the case (e.g. `open_rasterio` uses `nodatavals`). Also, the logic for getting/setting nodata in `reproject` is needs to be better generalized. Additionally, when writing `to_raster` is is probably to set the nodata value to the one in `encoding` if it exists and to use `fillna` to unset NaN values with the `encoding` value. 

Discussion:

Patch:
diff --git a/appveyor.yml b/appveyor.yml
index c19f8d52..ce67f833 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -12,7 +12,6 @@ install:
   - conda config --set always_yes yes
   - conda config --add channels conda-forge
   - conda config --add channels conda-forge
-  - conda update -q conda python
 
 #-----------------------------------------------------------------------------
 # Create conda environment
diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py
index 797b2010..acefb4b1 100644
--- a/rioxarray/rioxarray.py
+++ b/rioxarray/rioxarray.py
@@ -31,7 +31,7 @@ from rioxarray.exceptions import (
 )
 from rioxarray.crs import crs_to_wkt
 
-FILL_VALUE_NAMES = ("_FillValue", "missing_value", "fill_value")
+FILL_VALUE_NAMES = ("_FillValue", "missing_value", "fill_value", "nodata")
 UNWANTED_RIO_ATTRS = ("nodata", "nodatavals", "crs", "is_tiled", "res")
 DEFAULT_GRID_MAP = "spatial_ref"
 
@@ -73,18 +73,20 @@ def _get_grid_map_name(src_data_array):
 
 
 def _generate_attrs(src_data_array, dst_affine, dst_nodata):
-
     # add original attributes
     new_attrs = copy.deepcopy(src_data_array.attrs)
     # remove all nodata information
     for unwanted_attr in FILL_VALUE_NAMES + UNWANTED_RIO_ATTRS:
-        try:
-            del new_attrs[unwanted_attr]
-        except KeyError:
-            pass
+        new_attrs.pop(unwanted_attr, None)
 
     # add nodata information
-    new_attrs["_FillValue"] = src_data_array.encoding.get("_FillValue", dst_nodata)
+    fill_value = (
+        src_data_array.rio.nodata
+        if src_data_array.rio.nodata is not None
+        else dst_nodata
+    )
+    if src_data_array.rio.encoded_nodata is None and fill_value is not None:
+        new_attrs["_FillValue"] = fill_value
 
     # add raster spatial information
     new_attrs["transform"] = tuple(dst_affine)
@@ -149,20 +151,11 @@ def _add_attrs_proj(new_data_array, src_data_array):
         new_data_array.rio.transform(recalc=True),
         new_data_array.rio.nodata,
     )
-    if new_data_array.rio.nodata is None:
-        try:
-            new_attrs["_FillValue"] = src_data_array.attrs["nodata"]
-        except KeyError:
-            pass
-        try:
-            new_attrs["_FillValue"] = src_data_array.attrs["nodatavals"][0]
-        except KeyError:
-            pass
     # remove fill value if it already exists in the encoding
     # this is for data arrays pulling the encoding from a
     # source data array instead of being generated anew.
     if "_FillValue" in new_data_array.encoding:
-        del new_attrs["_FillValue"]
+        new_attrs.pop("_FillValue", None)
 
     new_data_array.attrs = new_attrs
 
@@ -275,19 +268,37 @@ class RasterArray(XRasterBase):
             self._crs = CRS.from_user_input(self._obj.attrs["crs"])
         return self._crs
 
+    @property
+    def encoded_nodata(self):
+        """Return the encoded nodata value for the dataset if encoded."""
+        return self._obj.encoding.get("_FillValue")
+
     @property
     def nodata(self):
         """Get the nodata value for the dataset."""
         if self._nodata is not None:
-            return self._nodata
+            return None if self._nodata is False else self._nodata
 
-        if self._obj.encoding.get("_FillValue") is not None:
+        if self.encoded_nodata is not None:
             self._nodata = np.nan
         else:
             self._nodata = self._obj.attrs.get(
                 "_FillValue",
-                self._obj.attrs.get("missing_value", self._obj.attrs.get("fill_value")),
+                self._obj.attrs.get(
+                    "missing_value",
+                    self._obj.attrs.get("fill_value", self._obj.attrs.get("nodata")),
+                ),
             )
+        if self._nodata is None:
+            try:
+                self._nodata = self._obj.attrs["nodatavals"][0]
+            except (KeyError, IndexError):
+                pass
+
+        if self._nodata is None:
+            self._nodata = False
+            return None
+
         return self._nodata
 
     def resolution(self, recalc=False):
@@ -481,12 +492,16 @@ class RasterArray(XRasterBase):
             dst_data = np.zeros((dst_height, dst_width), dtype=self._obj.dtype.type)
 
         try:
-            dst_nodata = self._obj.dtype.type(self.nodata or -9999)
+            dst_nodata = self._obj.dtype.type(
+                self.nodata if self.nodata is not None else -9999
+            )
         except ValueError:
             # if integer, set nodata to -9999
             dst_nodata = self._obj.dtype.type(-9999)
 
-        src_nodata = self._obj.dtype.type(self._obj.attrs.get("nodata", dst_nodata))
+        src_nodata = self._obj.dtype.type(
+            self.nodata if self.nodata is not None else dst_nodata
+        )
         rasterio.warp.reproject(
             source=src_data,
             destination=dst_data,
@@ -515,6 +530,7 @@ class RasterArray(XRasterBase):
             coords=_make_coords(self._obj, dst_affine, dst_width, dst_height, dst_crs),
             dims=tuple(dst_dims),
             attrs=new_attrs,
+            encoding=self._obj.encoding,
         )
         return add_spatial_ref(xda, dst_crs, DEFAULT_GRID_MAP)
 
@@ -788,6 +804,7 @@ class RasterArray(XRasterBase):
             coords=self._obj.coords,
             dims=self._obj.dims,
             attrs=self._obj.attrs,
+            encoding=self._obj.encoding,
         )
 
         # make sure correct attributes preserved & projection added
@@ -834,10 +851,17 @@ class RasterArray(XRasterBase):
             dtype=dtype,
             crs=self.crs,
             transform=self.transform(recalc=True),
-            nodata=self.nodata,
+            nodata=(
+                self.encoded_nodata if self.encoded_nodata is not None else self.nodata
+            ),
             **profile_kwargs,
         ) as dst:
-            data = self._obj.values.astype(dtype)
+
+            if self.encoded_nodata is not None:
+                out_data = self._obj.fillna(self.encoded_nodata)
+            else:
+                out_data = self._obj
+            data = out_data.astype(dtype).load().data
             if data.ndim == 2:
                 dst.write(data, 1)
             else:
diff --git a/sphinx/history.rst b/sphinx/history.rst
index 19b0abda..c2d61e7b 100644
--- a/sphinx/history.rst
+++ b/sphinx/history.rst
@@ -1,6 +1,12 @@
 History
 =======
 
+0.0.5
+-----
+- Find nodata and nodatavals in 'nodata' property (pull #12)
+- Added 'encoded_nodata' property to DataArray (pull #12)
+- Write the raster with encoded_nodata instead of NaN for nodata (pull #12)
+
 0.0.4
 ------
 - Added ability to export data array to raster (pull #8)


Test patch:
diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py
index 6db69b27..1708a019 100644
--- a/test/integration/test_integration_rioxarray.py
+++ b/test/integration/test_integration_rioxarray.py
@@ -659,6 +659,24 @@ def test_geographic_reproject():
         _assert_xarrays_equal(mds_repr, mdc)
 
 
+def test_geographic_reproject__missing_nodata():
+    sentinel_2_geographic = os.path.join(
+        TEST_INPUT_DATA_DIR, "sentinel_2_L1C_geographic.nc"
+    )
+    sentinel_2_utm = os.path.join(
+        TEST_COMPARE_DATA_DIR, "sentinel_2_L1C_utm__auto_nodata.nc"
+    )
+    with xarray.open_dataset(
+        sentinel_2_geographic, autoclose=True
+    ) as mda, xarray.open_dataset(sentinel_2_utm, autoclose=True) as mdc:
+        mda.red.attrs.pop("nodata")
+        mda.nir.attrs.pop("nodata")
+        mds_repr = mda.rio.reproject("+init=epsg:32721")
+        # mds_repr.to_netcdf(sentinel_2_utm)
+        # test
+        _assert_xarrays_equal(mds_repr, mdc)
+
+
 def test_geographic_resample_integer():
     sentinel_2_geographic = os.path.join(
         TEST_INPUT_DATA_DIR, "sentinel_2_L1C_geographic.nc"
@@ -676,7 +694,7 @@ def test_geographic_resample_integer():
 
 
 def test_to_raster(tmpdir):
-    tmp_raster = tmpdir.join("modis_raster.tiff")
+    tmp_raster = tmpdir.join("modis_raster.tif")
     test_tags = {"test": "1"}
     with xarray.open_dataarray(
         os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc"), autoclose=True
@@ -685,21 +703,23 @@ def test_to_raster(tmpdir):
         xds = mda.copy()
 
     with rasterio.open(str(tmp_raster)) as rds:
+        assert rds.count == 1
         assert rds.crs == xds.rio.crs
         assert_array_equal(rds.transform, xds.rio.transform())
-        assert_array_equal(rds.nodata, xds.rio.nodata)
-        assert_array_equal(rds.read(1), xds.values)
+        assert_array_equal(rds.nodata, xds.rio.encoded_nodata)
+        assert_array_equal(rds.read(1), xds.fillna(xds.rio.encoded_nodata).values)
         assert rds.count == 1
         assert rds.tags() == {"AREA_OR_POINT": "Area", **test_tags}
 
 
 def test_to_raster_3d(tmpdir):
-    tmp_raster = tmpdir.join("planet_3d_raster.tiff")
+    tmp_raster = tmpdir.join("planet_3d_raster.tif")
     with xarray.open_dataset(
         os.path.join(TEST_INPUT_DATA_DIR, "PLANET_SCOPE_3D.nc"), autoclose=True
     ) as mda:
-        mda.green.rio.to_raster(str(tmp_raster))
-        xds = mda.green.copy()
+        xds = mda.green.fillna(mda.green.rio.encoded_nodata)
+        xds.rio._nodata = mda.green.rio.encoded_nodata
+        xds.rio.to_raster(str(tmp_raster))
 
     with rasterio.open(str(tmp_raster)) as rds:
         assert rds.crs == xds.rio.crs