Repository: corteva/rioxarray

Issue:
Write out metadata to tif files on export
1. Preserve attributes in the tags
2. Band names: https://gis.stackexchange.com/a/284242/144357

Discussion:

Patch:
diff --git a/rioxarray/_io.py b/rioxarray/_io.py
index 7931571e..aac5d51c 100644
--- a/rioxarray/_io.py
+++ b/rioxarray/_io.py
@@ -272,7 +272,7 @@ def _get_rasterio_attrs(riods, masked):
         # The offset values for the raster bands
         attrs["offsets"] = riods.offsets
     if hasattr(riods, "descriptions") and any(riods.descriptions):
-        if len(riods.descriptions) == 1:
+        if len(set(riods.descriptions)) == 1:
             attrs["long_name"] = riods.descriptions[0]
         else:
             # Descriptions for each dataset band
diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py
index 5f665ebf..21e3ed52 100644
--- a/rioxarray/rioxarray.py
+++ b/rioxarray/rioxarray.py
@@ -12,7 +12,6 @@ datacube is licensed under the Apache License, Version 2.0:
 """
 import copy
 import warnings
-from datetime import datetime
 
 import numpy as np
 import rasterio.warp
@@ -30,6 +29,7 @@ from rioxarray.exceptions import (
     MissingCRS,
     NoDataInBounds,
     OneDimensionalRaster,
+    RioXarrayError,
     TooManyDimensions,
 )
 
@@ -200,6 +200,31 @@ def _make_dst_affine(src_data_array, src_crs, dst_crs, dst_resolution=None):
     return dst_affine, dst_width, dst_height
 
 
+def _write_metatata_to_raster(raster_handle, xarray_dataset, tags):
+    """
+    Write the metadata stored in the xarray object to raster metadata
+    """
+    # write metadata to raster
+    tags = {} if tags is None else tags
+    raster_handle.update_tags(**xarray_dataset.attrs, **tags)
+
+    # write band name information
+    long_name = xarray_dataset.attrs.get("long_name")
+    if isinstance(long_name, (tuple, list)):
+        if len(long_name) != raster_handle.count:
+            raise RioXarrayError(
+                "Number of names in the 'long_name' attribute does not equal "
+                "the number of bands."
+            )
+        for iii, band_description in enumerate(long_name):
+            raster_handle.set_band_description(iii + 1, band_description)
+    else:
+        band_description = long_name or xarray_dataset.name
+        if band_description:
+            for iii in range(raster_handle.count):
+                raster_handle.set_band_description(iii + 1, band_description)
+
+
 class XRasterBase(object):
     """This is the base class for the GIS extensions for xarray"""
 
@@ -1177,6 +1202,10 @@ class RasterArray(XRasterBase):
             ),
             **out_profile,
         ) as dst:
+
+            _write_metatata_to_raster(dst, self._obj, tags)
+
+            # write data to raster
             if windowed:
                 window_iter = dst.block_windows(1)
             else:
@@ -1193,8 +1222,6 @@ class RasterArray(XRasterBase):
                     dst.write(data, 1, window=window)
                 else:
                     dst.write(data, window=window)
-            if tags is not None:
-                dst.update_tags(**tags)
 
 
 @xarray.register_dataset_accessor("rio")
@@ -1268,7 +1295,6 @@ class RasterDataset(XRasterBase):
                 dst_affine_width_height=dst_affine_width_height,
                 resampling=resampling,
             )
-        resampled_dataset.attrs["creation_date"] = str(datetime.utcnow())
         return resampled_dataset
 
     def reproject_match(self, match_data_array, resampling=Resampling.nearest):
@@ -1301,7 +1327,6 @@ class RasterDataset(XRasterBase):
             resampled_dataset[var] = self._obj[var].rio.reproject_match(
                 match_data_array, resampling=resampling
             )
-        resampled_dataset.attrs["creation_date"] = str(datetime.utcnow())
         return resampled_dataset
 
     def clip_box(self, minx, miny, maxx, maxy, auto_expand=False, auto_expand_limit=3):
@@ -1341,7 +1366,6 @@ class RasterDataset(XRasterBase):
                 auto_expand=auto_expand,
                 auto_expand_limit=auto_expand_limit,
             )
-        clipped_dataset.attrs["creation_date"] = str(datetime.utcnow())
         return clipped_dataset
 
     def clip(self, geometries, crs, all_touched=False, drop=True, invert=False):
@@ -1395,7 +1419,6 @@ class RasterDataset(XRasterBase):
             clipped_dataset[var] = self._obj[var].rio.clip(
                 geometries, crs=crs, all_touched=all_touched, drop=drop, invert=invert
             )
-        clipped_dataset.attrs["creation_date"] = str(datetime.utcnow())
         return clipped_dataset
 
     def interpolate_na(self, method="nearest"):
@@ -1415,5 +1438,4 @@ class RasterDataset(XRasterBase):
         interpolated_dataset = xarray.Dataset(attrs=self._obj.attrs)
         for var in self.vars:
             interpolated_dataset[var] = self._obj[var].rio.interpolate_na(method=method)
-        interpolated_dataset.attrs["creation_date"] = str(datetime.utcnow())
         return interpolated_dataset
diff --git a/sphinx/history.rst b/sphinx/history.rst
index f8334a1f..b6390983 100644
--- a/sphinx/history.rst
+++ b/sphinx/history.rst
@@ -4,6 +4,8 @@ History
 0.0.18
 ------
 - Fixed windowed writing to require tiled output raster (pull #66)
+- Write data array attributes using `rio.to_raster` (issue #64)
+- Write variable name to descriptions if possible in `rio.to_raster` (issue #64)
 
 0.0.17
 ------


Test patch:
diff --git a/test/conftest.py b/test/conftest.py
index 8f1ab378..419b748f 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -1,6 +1,5 @@
 import os
 
-import xarray
 from numpy.testing import assert_almost_equal, assert_array_equal
 
 from rioxarray.rioxarray import UNWANTED_RIO_ATTRS
@@ -14,9 +13,6 @@ def _assert_xarrays_equal(input_xarray, compare_xarray, precision=7):
     # xarray.testing.assert_equal(input_xarray, compare_xarray)
     def assert_attrs_equal(input_xr, compare_xr):
         """check attrubutes that matter"""
-        if isinstance(input_xr, xarray.Dataset):
-            assert "creation_date" in input_xr.attrs
-
         for attr in compare_xr.attrs:
             if attr == "transform":
                 assert_almost_equal(input_xr.attrs[attr], compare_xr.attrs[attr][:6])
diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py
index 4477afd0..2a48acdb 100644
--- a/test/integration/test_integration_rioxarray.py
+++ b/test/integration/test_integration_rioxarray.py
@@ -1,5 +1,4 @@
 import os
-import sys
 from functools import partial
 
 import numpy
@@ -17,6 +16,7 @@ from rioxarray.exceptions import (
     MissingCRS,
     NoDataInBounds,
     OneDimensionalRaster,
+    RioXarrayError,
 )
 from rioxarray.rioxarray import _make_coords
 from test.conftest import (
@@ -147,10 +147,7 @@ def test_clip_box(modis_clip):
             maxx=xdi.x[6].values,
             maxy=xdi.y[4].values,
         )
-        if isinstance(xdc, xarray.Dataset):
-            xdc.attrs["creation_date"] = clipped_ds.attrs["creation_date"]
         _assert_xarrays_equal(clipped_ds, xdc)
-
         # make sure it safely writes to netcdf
         clipped_ds.to_netcdf(modis_clip["output"])
 
@@ -167,9 +164,6 @@ def test_clip_box__auto_expand(modis_clip):
             auto_expand=True,
         )
 
-        if isinstance(xdc, xarray.Dataset):
-            xdc.attrs["creation_date"] = clipped_ds.attrs["creation_date"]
-
         _assert_xarrays_equal(clipped_ds, xdc)
         # make sure it safely writes to netcdf
         clipped_ds.to_netcdf(modis_clip["output"])
@@ -747,6 +741,7 @@ def test_to_raster(open_method, windowed, recalc_transform, tmpdir):
             tags=test_tags,
         )
         xds = mda.copy().squeeze()
+        xds_attrs = {key: str(value) for key, value in mda.attrs.items()}
 
     with rasterio.open(str(tmp_raster)) as rds:
         assert rds.count == 1
@@ -755,7 +750,7 @@ def test_to_raster(open_method, windowed, recalc_transform, tmpdir):
         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}
+        assert rds.tags() == {"AREA_OR_POINT": "Area", **test_tags, **xds_attrs}
 
 
 @pytest.mark.parametrize(
@@ -773,12 +768,49 @@ def test_to_raster_3d(open_method, windowed, tmpdir):
         xds = mda.green.fillna(mda.green.rio.encoded_nodata)
         xds.rio._nodata = mda.green.rio.encoded_nodata
         xds.rio.to_raster(str(tmp_raster), windowed=windowed)
+        xds_attrs = {key: str(value) for key, value in xds.attrs.items()}
 
     with rasterio.open(str(tmp_raster)) as rds:
         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(), xds.values)
+        assert rds.tags() == {"AREA_OR_POINT": "Area", **xds_attrs}
+        assert rds.descriptions == ("green", "green")
+
+    # test roundtrip
+    with rioxarray.open_rasterio(str(tmp_raster)) as rds:
+        assert rds.attrs["long_name"] == "green"
+
+
+def test_to_raster__custom_description(tmpdir):
+    tmp_raster = tmpdir.join("planet_3d_raster.tif")
+    with xarray.open_dataset(
+        os.path.join(TEST_INPUT_DATA_DIR, "PLANET_SCOPE_3D.nc")
+    ) as mda:
+        xds = mda.green.fillna(mda.green.rio.encoded_nodata)
+        xds.attrs["long_name"] = ("one", "two")
+        xds.rio.to_raster(str(tmp_raster))
+        xds_attrs = {key: str(value) for key, value in xds.attrs.items()}
+
+    with rasterio.open(str(tmp_raster)) as rds:
+        assert rds.tags() == {"AREA_OR_POINT": "Area", **xds_attrs}
+        assert rds.descriptions == ("one", "two")
+
+    # test roundtrip
+    with rioxarray.open_rasterio(str(tmp_raster)) as rds:
+        assert rds.attrs["long_name"] == ("one", "two")
+
+
+def test_to_raster__custom_description__wrong(tmpdir):
+    tmp_raster = tmpdir.join("planet_3d_raster.tif")
+    with xarray.open_dataset(
+        os.path.join(TEST_INPUT_DATA_DIR, "PLANET_SCOPE_3D.nc")
+    ) as mda:
+        xds = mda.green.fillna(mda.green.rio.encoded_nodata)
+        xds.attrs["long_name"] = ("one", "two", "three")
+        with pytest.raises(RioXarrayError):
+            xds.rio.to_raster(str(tmp_raster))
 
 
 @pytest.mark.xfail(reason="Precision issues with windowed writing on python 3.6")