Repository: corteva/rioxarray

Issue:
Add ability to export to geotiff
It would be useful to write the dataframe out to a geotiff file.

Discussion:
snowman2: Appveyor builds not working:
```
Traceback (most recent call last):
  File "C:\Miniconda3-x64\lib\site-packages\conda\exceptions.py", line 1043, in __call__
    return func(*args, **kwargs)
  File "C:\Miniconda3-x64\lib\site-packages\conda\cli\main.py", line 84, in _main
    exit_code = do_call(args, p)
  File "C:\Miniconda3-x64\lib\site-packages\conda\cli\conda_argparse.py", line 80, in do_call
    module = import_module(relative_mod, __name__.rsplit('.', 1)[0])
  File "C:\Miniconda3-x64\lib\importlib\__init__.py", line 127, in import_module
    return _bootstrap._gcd_import(name[level:], package, level)
  File "<frozen importlib._bootstrap>", line 1006, in _gcd_import
  File "<frozen importlib._bootstrap>", line 983, in _find_and_load
  File "<frozen importlib._bootstrap>", line 967, in _find_and_load_unlocked
  File "<frozen importlib._bootstrap>", line 677, in _load_unlocked
  File "<frozen importlib._bootstrap_external>", line 728, in exec_module
  File "<frozen importlib._bootstrap>", line 219, in _call_with_frames_removed
  File "C:\Miniconda3-x64\lib\site-packages\conda\cli\main_create.py", line 10, in <module>
    from .install import install
  File "C:\Miniconda3-x64\lib\site-packages\conda\cli\install.py", line 19, in <module>
    from ..core.index import calculate_channel_urls, get_index
  File "C:\Miniconda3-x64\lib\site-packages\conda\core\index.py", line 9, in <module>
    from .package_cache_data import PackageCacheData
  File "C:\Miniconda3-x64\lib\site-packages\conda\core\package_cache_data.py", line 15, in <module>
    from conda_package_handling.api import InvalidArchiveError
  File "C:\Miniconda3-x64\lib\site-packages\conda_package_handling\api.py", line 3, in <module>
    from libarchive.exception import ArchiveError as _LibarchiveArchiveError
  File "C:\Miniconda3-x64\lib\site-packages\libarchive\__init__.py", line 1, in <module>
    from .entry import ArchiveEntry
  File "C:\Miniconda3-x64\lib\site-packages\libarchive\entry.py", line 6, in <module>
    from . import ffi
  File "C:\Miniconda3-x64\lib\site-packages\libarchive\ffi.py", line 27, in <module>
    libarchive = ctypes.cdll.LoadLibrary(libarchive_path)
  File "C:\Miniconda3-x64\lib\ctypes\__init__.py", line 434, in LoadLibrary
    return self._dlltype(name)
  File "C:\Miniconda3-x64\lib\ctypes\__init__.py", line 356, in __init__
    self._handle = _dlopen(self._name, mode)
TypeError: LoadLibrary() argument 1 must be str, not None
During handling of the above exception, another exception occurred:
```
djhoese: Would it be possible to also add a "tags" keyword argument?
snowman2: > Would it be possible to also add a "tags" keyword argument?

That could be useful. Do you have an example of how you would like this to work?
djhoese: In `trollimage` we use it for the image DATETIME by default: https://github.com/pytroll/trollimage/blob/master/trollimage/xrimage.py#L340

Otherwise, we have a lot of users who have custom visualization clients or other tools where they've historically had to throw extra metadata in to tiff tags/metadata (@mraspaud may have examples). But a pretty common one is people adding a `scale_factor` and `add_offset` to their metadata so the visualization client can use them to get "real" data values back.
snowman2: > Oh and what about a colorinterp?

Maybe a default `**kwargs` you can pass into writing the raster would be the solution here. I have been thinking about this and there are many options that would be nice to add.
djhoese: > Maybe a default **kwargs you can pass into writing the raster would be the solution here. 

Yeah I didn't want to suggest it since most people get annoyed of "catch all" kwargs if they can be avoided.
snowman2: > Yeah I didn't want to suggest it since most people get annoyed of "catch all" kwargs if they can be avoided.

I definitely understand that. I have been debating it in my head a bit. The catch all is annoying, but if something is added in another version of rasterio or a driver specific option is added and it is not explicitly set here, the function won't be able to handle it. Another option is to add an argument for a rasterio profile dictionary.
snowman2: Looks like rasterio.open uses kwargs, so I don't think it would be terrible to use kwargs here as well:
https://github.com/mapbox/rasterio/blob/fc9861cef7a70f5136dd4b47920d05ea88312137/rasterio/__init__.py#L58
snowman2: I think adding a `tags` dict option that would be called in `update_tags` would solve your problems, correct?
https://github.com/pytroll/trollimage/blob/63fa32f2d40bb65ebc39c4be1fb1baf8f163db98/trollimage/xrimage.py#L357
djhoese: It would make it more flexible for users, yes. I forgot about the `update_tags` call being needed.
snowman2: @djhoese I am planning to merge these changes in soon as the tests pass (for non-windows builds at least). See any problems with the code in its current state?
snowman2: > The only other helpful long term solution would be to use windowed writing like we do in trollimage to achieve better performance with dask (saving chunk by chunk), but this has its own issues we've discovered with GDAL's in-memory caching and making large output files if you don't have a lot of memory available.

That is a good idea. With the issues you are mentioning, I think I will create a separate issue/PR for that to better understand the issues.
snowman2: Thanks for the review @djhoese!

Patch:
diff --git a/rioxarray/exceptions.py b/rioxarray/exceptions.py
index acf629e6..723a3915 100644
--- a/rioxarray/exceptions.py
+++ b/rioxarray/exceptions.py
@@ -18,3 +18,11 @@ class OneDimensionalRaster(RioXarrayError):
 
 class SingleVariableDataset(RioXarrayError):
     """This is for when you have a dataset with a single variable."""
+
+
+class TooManyDimensions(RioXarrayError):
+    """This is raised when there are more dimensions than is supported by the method"""
+
+
+class InvalidDimensionOrder(RioXarrayError):
+    """This is raised when there the dimensions are not ordered correctly."""
diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py
index 68c10e76..797b2010 100644
--- a/rioxarray/rioxarray.py
+++ b/rioxarray/rioxarray.py
@@ -23,7 +23,12 @@ from rasterio.enums import Resampling
 from rasterio.features import geometry_mask
 from scipy.interpolate import griddata
 
-from rioxarray.exceptions import NoDataInBounds, OneDimensionalRaster
+from rioxarray.exceptions import (
+    InvalidDimensionOrder,
+    NoDataInBounds,
+    OneDimensionalRaster,
+    TooManyDimensions,
+)
 from rioxarray.crs import crs_to_wkt
 
 FILL_VALUE_NAMES = ("_FillValue", "missing_value", "fill_value")
@@ -327,6 +332,34 @@ class RasterArray(XRasterBase):
         bottom = float(self._obj[self.y_dim][-1])
         return left, bottom, right, top
 
+    def _check_dimensions(self):
+        """
+        This function validates that the dimensions 2D/3D and
+        they are are in the proper order.
+
+        Returns:
+        --------
+        str or None: Name extra dimension.
+        """
+        extra_dims = list(set(list(self._obj.dims)) - set([self.x_dim, self.y_dim]))
+        if len(extra_dims) > 1:
+            raise TooManyDimensions("Only 2D and 3D data arrays supported.")
+        elif extra_dims and self._obj.dims != (extra_dims[0], self.y_dim, self.x_dim):
+            raise InvalidDimensionOrder(
+                "Invalid dimension order. Expected order: {0}. "
+                "You can use `DataArray.transpose{0}`"
+                " to reorder your dimensions.".format(
+                    (extra_dims[0], self.y_dim, self.x_dim)
+                )
+            )
+        elif not extra_dims and self._obj.dims != (self.y_dim, self.x_dim):
+            raise InvalidDimensionOrder(
+                "Invalid dimension order. Expected order: {0}"
+                "You can use `DataArray.transpose{0}` "
+                "to reorder your dimensions.".format((self.y_dim, self.x_dim))
+            )
+        return extra_dims[0] if extra_dims else None
+
     def bounds(self, recalc=False):
         """Determine the bounds of the `xarray.DataArray`
 
@@ -438,17 +471,13 @@ class RasterArray(XRasterBase):
             dst_affine, dst_width, dst_height = _make_dst_affine(
                 self._obj, self.crs, dst_crs, resolution
             )
-        extra_dims = list(set(list(self._obj.dims)) - set([self.x_dim, self.y_dim]))
-        if len(extra_dims) > 1:
-            raise RuntimeError("Reproject only supports 2D and 3D datasets.")
-        if extra_dims:
-            assert self._obj.dims == (extra_dims[0], self.y_dim, self.x_dim)
+        extra_dim = self._check_dimensions()
+        if extra_dim:
             dst_data = np.zeros(
-                (self._obj[extra_dims[0]].size, dst_height, dst_width),
+                (self._obj[extra_dim].size, dst_height, dst_width),
                 dtype=self._obj.dtype.type,
             )
         else:
-            assert self._obj.dims == (self.y_dim, self.x_dim)
             dst_data = np.zeros((dst_height, dst_width), dtype=self._obj.dtype.type)
 
         try:
@@ -742,13 +771,10 @@ class RasterArray(XRasterBase):
         :class:`xarray.DataArray`: An interpolated :class:`xarray.DataArray` object.
 
         """
-        extra_dims = list(set(list(self._obj.dims)) - set([self.x_dim, self.y_dim]))
-        if len(extra_dims) > 1:
-            raise RuntimeError("Interpolate only supports 2D and 3D datasets.")
-        if extra_dims:
-            assert self._obj.dims == (extra_dims[0], self.y_dim, self.x_dim)
+        extra_dim = self._check_dimensions()
+        if extra_dim:
             interp_data = []
-            for _, sub_xds in self._obj.groupby(extra_dims[0]):
+            for _, sub_xds in self._obj.groupby(extra_dim):
                 interp_data.append(
                     self._interpolate_na(sub_xds.load().data, method=method)
                 )
@@ -769,6 +795,56 @@ class RasterArray(XRasterBase):
 
         return interp_array
 
+    def to_raster(
+        self, raster_path, driver="GTiff", dtype=None, tags=None, **profile_kwargs
+    ):
+        """
+        Export the DataArray to a raster file.
+
+        Parameters
+        ----------
+        raster_path: str
+            The path to output the raster to.
+        driver: str, optional
+            The name of the GDAL/rasterio driver to use to export the raster.
+            Default is "GTiff".
+        dtype: str, optional
+            The data type to write the raster to. Default is the datasets dtype.
+        tags: dict, optional
+            A dictionary of tags to write to the raster.
+        **profile_kwargs
+            Additional keyword arguments to pass into writing the raster. The
+            nodata, transform, crs, count, width, and height attributes
+            are automatically added.
+
+        """
+        width, height = self.shape
+        dtype = str(self._obj.dtype) if dtype is None else dtype
+        extra_dim = self._check_dimensions()
+        count = 1
+        if extra_dim is not None:
+            count = self._obj[extra_dim].size
+        with rasterio.open(
+            raster_path,
+            "w",
+            driver=driver,
+            height=int(height),
+            width=int(width),
+            count=count,
+            dtype=dtype,
+            crs=self.crs,
+            transform=self.transform(recalc=True),
+            nodata=self.nodata,
+            **profile_kwargs,
+        ) as dst:
+            data = self._obj.values.astype(dtype)
+            if data.ndim == 2:
+                dst.write(data, 1)
+            else:
+                dst.write(data)
+            if tags is not None:
+                dst.update_tags(**tags)
+
 
 @xarray.register_dataset_accessor("rio")
 class RasterDataset(XRasterBase):
diff --git a/sphinx/history.rst b/sphinx/history.rst
index ff7fe9eb..19b0abda 100644
--- a/sphinx/history.rst
+++ b/sphinx/history.rst
@@ -1,4 +1,6 @@
 History
 =======
 
-In pre-release stages now.
+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 0d0cf491..6db69b27 100644
--- a/test/integration/test_integration_rioxarray.py
+++ b/test/integration/test_integration_rioxarray.py
@@ -2,6 +2,7 @@ import os
 
 import numpy
 import pytest
+import rasterio
 import xarray
 from numpy.testing import assert_almost_equal, assert_array_equal
 from rasterio.crs import CRS
@@ -672,3 +673,36 @@ def test_geographic_resample_integer():
         # mds_interp.to_netcdf(sentinel_2_interp)
         # test
         _assert_xarrays_equal(mds_interp, mdc)
+
+
+def test_to_raster(tmpdir):
+    tmp_raster = tmpdir.join("modis_raster.tiff")
+    test_tags = {"test": "1"}
+    with xarray.open_dataarray(
+        os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc"), autoclose=True
+    ) as mda:
+        mda.rio.to_raster(str(tmp_raster), tags=test_tags)
+        xds = mda.copy()
+
+    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(1), xds.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")
+    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()
+
+    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)