Repository: corteva/rioxarray

Issue:
Add support for `to_raster()` with 2D datasets.
It is often times useful to write a dataset to a raster. This is only possible if the data is 2D as the band dimension is taken by the different variables.

It would be nice to support metadata and writing band descriptions.

Discussion:

Patch:
diff --git a/docs/examples/convert_to_raster.ipynb b/docs/examples/convert_to_raster.ipynb
index ba053fc7..d67fabbf 100644
--- a/docs/examples/convert_to_raster.ipynb
+++ b/docs/examples/convert_to_raster.ipynb
@@ -7,10 +7,15 @@
     "# Example - Convert dataset to raster (GeoTiff)\n",
     "\n",
     "Often, it is desirable to take a variable (band) out of your dataset and export it to a raster.\n",
-    "This is possible with the [rio.to_raster()](../rioxarray.rst#rioxarray.rioxarray.RasterArray.to_raster) method. It does most of the work for you so you don't\n",
+    "This is possible with the `rio.to_raster()`method. It does most of the work for you so you don't\n",
     "have to.\n",
     "\n",
-    "Note: The [rio.to_raster()](../rioxarray.rst#rioxarray.rioxarray.RasterArray.to_raster) method only works on a 2-dimensional or 3-dimensional `xarray.DataArray`"
+    "Note: The `rio.to_raster()` method only works on a 2-dimensional or 3-dimensional `xarray.DataArray` or a 2-dimensional `xarray.Dataset`.\n",
+    "\n",
+    "API Reference:\n",
+    "\n",
+    "- DataArray: [rio.to_raster()](../rioxarray.rst#rioxarray.rioxarray.RasterArray.to_raster)\n",
+    "- Dataset: [rio.to_raster()](../rioxarray.rst#rioxarray.rioxarray.RasterDataset.to_raster)\n"
    ]
   },
   {
@@ -19,10 +24,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import rioxarray\n",
-    "import rasterio\n",
-    "import xarray\n",
-    "import re"
+    "import rioxarray"
    ]
   },
   {
@@ -43,11 +45,12 @@
        "<pre>&lt;xarray.Dataset&gt;\n",
        "Dimensions:      (time: 2, x: 10, y: 10)\n",
        "Coordinates:\n",
-       "  * time         (time) object 2016-12-19 10:27:29 2016-12-29 12:52:41.659696\n",
        "  * y            (y) float64 8.085e+06 8.085e+06 ... 8.085e+06 8.085e+06\n",
        "  * x            (x) float64 4.663e+05 4.663e+05 ... 4.663e+05 4.663e+05\n",
+       "  * time         (time) object 2016-12-19 10:27:29 2016-12-29 12:52:41.659696\n",
        "    spatial_ref  int64 0\n",
        "Data variables:\n",
+       "    blue         (time, y, x) float64 ...\n",
        "    green        (time, y, x) float64 ...\n",
        "Attributes:\n",
        "    coordinates:  spatial_ref</pre>"
@@ -56,11 +59,12 @@
        "<xarray.Dataset>\n",
        "Dimensions:      (time: 2, x: 10, y: 10)\n",
        "Coordinates:\n",
-       "  * time         (time) object 2016-12-19 10:27:29 2016-12-29 12:52:41.659696\n",
        "  * y            (y) float64 8.085e+06 8.085e+06 ... 8.085e+06 8.085e+06\n",
        "  * x            (x) float64 4.663e+05 4.663e+05 ... 4.663e+05 4.663e+05\n",
+       "  * time         (time) object 2016-12-19 10:27:29 2016-12-29 12:52:41.659696\n",
        "    spatial_ref  int64 0\n",
        "Data variables:\n",
+       "    blue         (time, y, x) float64 ...\n",
        "    green        (time, y, x) float64 ...\n",
        "Attributes:\n",
        "    coordinates:  spatial_ref"
@@ -74,7 +78,6 @@
    "source": [
     "rds = rioxarray.open_rasterio(\n",
     "    \"../../test/test_data/input/PLANET_SCOPE_3D.nc\",\n",
-    "    variable=[\"green\"],\n",
     ")\n",
     "rds"
    ]
@@ -83,110 +86,60 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Managing Information Loss with xarray operations\n",
-    "\n",
-    "Sometimes, you can lose important information from your dataset when performing operations.\n",
-    "You will likely want to keep track of your `nodata` and your `CRS`.\n",
-    "\n",
-    "API Reference:\n",
+    "## Converting Dataset to raster\n",
     "\n",
-    "- [rio.to_raster()](../rioxarray.rst#rioxarray.rioxarray.RasterArray.to_raster)\n",
-    "- [rio.write_crs()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.write_crs)\n",
-    "- [rio.update_attrs()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.update_attrs)\n",
-    "- [rio.crs](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.crs)\n",
-    "- [rio.nodata](../rioxarray.rst#rioxarray.rioxarray.RasterArray.nodata)"
+    "Dataset: [rio.to_raster()](../rioxarray.rst#rioxarray.rioxarray.RasterDataset.to_raster)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 3,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "({}, CRS.from_epsg(32722), None)"
-      ]
-     },
-     "execution_count": 3,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "new_ds = rds.green + rds.green\n",
-    "new_ds.attrs, new_ds.rio.crs, new_ds.rio.nodata"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 4,
-   "metadata": {},
    "outputs": [],
    "source": [
-    "new_ds.rio.to_raster(\"combination.tif\")"
+    "# note how one time slice was selected on export to make the dataset 2D\n",
+    "rds.isel(time=0).rio.to_raster(\"planet_scope.tif\")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "{\"bounds\": [466266.0, 8084670.0, 466296.0, 8084700.0], \"colorinterp\": [\"gray\", \"undefined\"], \"count\": 2, \"crs\": \"EPSG:32722\", \"descriptions\": [\"green\", \"green\"], \"driver\": \"GTiff\", \"dtype\": \"float64\", \"height\": 10, \"indexes\": [1, 2], \"interleave\": \"pixel\", \"lnglat\": [-51.31732641226951, -17.322997474192466], \"mask_flags\": [[\"all_valid\"], [\"all_valid\"]], \"nodata\": null, \"res\": [3.0, 3.0], \"shape\": [10, 10], \"tiled\": false, \"transform\": [3.0, 0.0, 466266.0, 0.0, -3.0, 8084700.0, 0.0, 0.0, 1.0], \"units\": [null, null], \"width\": 10}\n"
+      "{\"bounds\": [466266.0, 8084670.0, 466296.0, 8084700.0], \"colorinterp\": [\"gray\", \"undefined\"], \"count\": 2, \"crs\": \"EPSG:32722\", \"descriptions\": [\"blue\", \"green\"], \"driver\": \"GTiff\", \"dtype\": \"float64\", \"height\": 10, \"indexes\": [1, 2], \"interleave\": \"pixel\", \"lnglat\": [-51.31732641226951, -17.322997474192466], \"mask_flags\": [[\"nodata\"], [\"nodata\"]], \"nodata\": NaN, \"res\": [3.0, 3.0], \"shape\": [10, 10], \"tiled\": false, \"transform\": [3.0, 0.0, 466266.0, 0.0, -3.0, 8084700.0, 0.0, 0.0, 1.0], \"units\": [null, null], \"width\": 10}\n"
      ]
     }
    ],
    "source": [
-    "!rio info combination.tif"
+    "!rio info planet_scope.tif"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": 6,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "({'grid_mapping': 'spatial_ref',\n",
-       "  'nodata': 0,\n",
-       "  'units': ('DN', 'DN'),\n",
-       "  '_FillValue': nan,\n",
-       "  'transform': (3.0, 0.0, 466266.0, 0.0, -3.0, 8084700.0),\n",
-       "  'scale_factor': 1.0,\n",
-       "  'add_offset': 0.0},\n",
-       " CRS.from_epsg(32722),\n",
-       " nan)"
-      ]
-     },
-     "execution_count": 6,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
    "source": [
-    "new_ds = rds.green + rds.green\n",
-    "new_ds.rio.write_crs(rds.green.rio.crs, inplace=True)\n",
-    "new_ds.rio.update_attrs(rds.green.attrs, inplace=True)\n",
-    "new_ds.attrs, new_ds.rio.crs, new_ds.rio.nodata"
+    "## Converting DataArray to raster\n",
+    "\n",
+    "DataArray: [rio.to_raster()](../rioxarray.rst#rioxarray.rioxarray.RasterArray.to_raster)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [],
    "source": [
-    "new_ds.rio.to_raster(\"combination_keep_attrs.tif\")"
+    "# note how selecting one variable allowed for multiple time steps in a single raster\n",
+    "rds.green.rio.to_raster(\"planet_scope_green.tif\")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [
     {
@@ -198,7 +151,7 @@
     }
    ],
    "source": [
-    "!rio info combination_keep_attrs.tif"
+    "!rio info planet_scope_green.tif"
    ]
   }
  ],
@@ -218,7 +171,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.7"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,
diff --git a/docs/examples/examples.rst b/docs/examples/examples.rst
index 9f4858ed..f349c5cc 100644
--- a/docs/examples/examples.rst
+++ b/docs/examples/examples.rst
@@ -10,6 +10,7 @@ This page contains links to a collection of examples of how to use rioxarray.
    :caption: Notebooks:
 
    crs_management.ipynb
+   manage_information_loss.ipynb
    clip_geom.ipynb
    clip_box.ipynb
    reproject.ipynb
diff --git a/docs/examples/manage_information_loss.ipynb b/docs/examples/manage_information_loss.ipynb
new file mode 100644
index 00000000..cae1fb55
--- /dev/null
+++ b/docs/examples/manage_information_loss.ipynb
@@ -0,0 +1,199 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Managing Information Loss with xarray operations\n",
+    "\n",
+    "Sometimes, you can lose important information from your dataset when performing operations.\n",
+    "You will likely want to keep track of the attributes, `nodata`, and `CRS`.\n",
+    "\n",
+    "API Reference:\n",
+    "\n",
+    "- [rio.to_raster()](../rioxarray.rst#rioxarray.rioxarray.RasterDataset.to_raster)\n",
+    "- [rio.write_crs()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.write_crs)\n",
+    "- [rio.update_attrs()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.update_attrs)\n",
+    "- [rio.crs](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.crs)\n",
+    "- [rio.nodata](../rioxarray.rst#rioxarray.rioxarray.RasterArray.nodata)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import rioxarray"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "See docs for [rioxarray.open_rasterio](../rioxarray.rst#rioxarray-open-rasterio)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rds = rioxarray.open_rasterio(\n",
+    "    \"../../test/test_data/input/PLANET_SCOPE_3D.nc\",\n",
+    "    variable=[\"green\"],\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Notice the original data:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "({'grid_mapping': 'spatial_ref',\n",
+       "  'nodata': 0,\n",
+       "  'units': ('DN', 'DN'),\n",
+       "  '_FillValue': nan,\n",
+       "  'transform': (3.0, 0.0, 466266.0, 0.0, -3.0, 8084700.0),\n",
+       "  'scale_factor': 1.0,\n",
+       "  'add_offset': 0.0},\n",
+       " CRS.from_epsg(32722),\n",
+       " nan)"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "rds.green.attrs, rds.green.rio.crs, rds.green.rio.nodata"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Notice how information is lost in the operation:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "({}, CRS.from_epsg(32722), None)"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "new_ds = rds.green + rds.green\n",
+    "new_ds.attrs, new_ds.rio.crs, new_ds.rio.nodata"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The solution is to save the original attributes and then copy them over\n",
+    "once the operation is complete:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "({'grid_mapping': 'spatial_ref',\n",
+       "  'nodata': 0,\n",
+       "  'units': ('DN', 'DN'),\n",
+       "  '_FillValue': nan,\n",
+       "  'transform': (3.0, 0.0, 466266.0, 0.0, -3.0, 8084700.0),\n",
+       "  'scale_factor': 1.0,\n",
+       "  'add_offset': 0.0},\n",
+       " CRS.from_epsg(32722),\n",
+       " nan)"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "new_ds = rds.green + rds.green\n",
+    "new_ds.rio.write_crs(rds.green.rio.crs, inplace=True)\n",
+    "new_ds.rio.update_attrs(rds.green.attrs, inplace=True)\n",
+    "new_ds.attrs, new_ds.rio.crs, new_ds.rio.nodata"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "new_ds.rio.to_raster(\"combination_keep_attrs.tif\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "{\"bounds\": [466266.0, 8084670.0, 466296.0, 8084700.0], \"colorinterp\": [\"gray\", \"undefined\"], \"count\": 2, \"crs\": \"EPSG:32722\", \"descriptions\": [\"green\", \"green\"], \"driver\": \"GTiff\", \"dtype\": \"float64\", \"height\": 10, \"indexes\": [1, 2], \"interleave\": \"pixel\", \"lnglat\": [-51.31732641226951, -17.322997474192466], \"mask_flags\": [[\"nodata\"], [\"nodata\"]], \"nodata\": NaN, \"res\": [3.0, 3.0], \"shape\": [10, 10], \"tiled\": false, \"transform\": [3.0, 0.0, 466266.0, 0.0, -3.0, 8084700.0, 0.0, 0.0, 1.0], \"units\": [null, null], \"width\": 10}\n"
+     ]
+    }
+   ],
+   "source": [
+    "!rio info combination_keep_attrs.tif"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/docs/history.rst b/docs/history.rst
index 9b2c6889..96e0dfc5 100644
--- a/docs/history.rst
+++ b/docs/history.rst
@@ -1,6 +1,10 @@
 History
 =======
 
+0.0.21
+-------
+- ENH: Added to_raster method for Datasets (issue #76)
+
 0.0.20
 ------
 - BUG: ensure band_key is list when iterating over bands for mask and scale (pull #87)
diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py
index 68d84244..c1ab22e3 100644
--- a/rioxarray/rioxarray.py
+++ b/rioxarray/rioxarray.py
@@ -13,6 +13,7 @@ datacube is licensed under the Apache License, Version 2.0:
 import copy
 import math
 import warnings
+from uuid import uuid4
 
 import numpy as np
 import rasterio.warp
@@ -233,6 +234,10 @@ def _write_metatata_to_raster(raster_handle, xarray_dataset, tags):
         + FILL_VALUE_NAMES
         + ("transform", "scales", "scale_factor", "add_offset", "offsets")
     )
+    # this is for when multiple values are used
+    # in this case, it will be stored in the raster description
+    if not isinstance(tags.get("long_name"), str):
+        skip_tags += ("long_name",)
     tags = {key: value for key, value in tags.items() if key not in skip_tags}
     raster_handle.update_tags(**tags)
 
@@ -1475,3 +1480,84 @@ class RasterDataset(XRasterBase):
         for var in self.vars:
             interpolated_dataset[var] = self._obj[var].rio.interpolate_na(method=method)
         return interpolated_dataset
+
+    def to_raster(
+        self,
+        raster_path,
+        driver="GTiff",
+        dtype=None,
+        tags=None,
+        windowed=False,
+        recalc_transform=True,
+        **profile_kwargs,
+    ):
+        """
+        Export the Dataset to a raster file. Only works with 2D data.
+
+        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.
+        windowed: bool, optional
+            If True, it will write using the windows of the output raster.
+            This only works if the output raster is tiled. As such, if you
+            set this to True, the output raster will be tiled.
+            Default is False.
+        **profile_kwargs
+            Additional keyword arguments to pass into writing the raster. The
+            nodata, transform, crs, count, width, and height attributes
+            are ignored.
+
+        """
+        variable_dim = "band_{}".format(uuid4())
+        data_array = self._obj.to_array(dim=variable_dim)
+        # write data array names to raster
+        data_array.attrs["long_name"] = data_array[variable_dim].values.tolist()
+        # ensure raster metadata preserved
+        scales = []
+        offsets = []
+        nodatavals = []
+        crs_list = []
+        for data_var in data_array[variable_dim].values:
+            scales.append(self._obj[data_var].attrs.get("scale_factor", 1.0))
+            offsets.append(self._obj[data_var].attrs.get("add_offset", 0.0))
+            nodatavals.append(self._obj[data_var].rio.nodata)
+            crs_list.append(self._obj[data_var].rio.crs)
+        data_array.attrs["scales"] = scales
+        data_array.attrs["offsets"] = offsets
+        nodata = nodatavals[0]
+        if (
+            all(nodataval == nodata for nodataval in nodatavals)
+            or np.isnan(nodatavals).all()
+        ):
+            data_array.rio.write_nodata(nodata, inplace=True)
+        else:
+            raise RioXarrayError(
+                "All nodata values must be the same when exporting to raster. "
+                "Current values: {}".format(nodatavals)
+            )
+        crs = crs_list[0]
+        if all(crs_i == crs for crs_i in crs_list):
+            data_array.rio.write_crs(crs, inplace=True)
+        else:
+            raise RioXarrayError(
+                "All CRS must be the same when exporting to raster. "
+                "Current values: {}".format(crs_list)
+            )
+        # write it to a raster
+        data_array.rio.to_raster(
+            raster_path=raster_path,
+            driver=driver,
+            dtype=dtype,
+            tags=tags,
+            windowed=windowed,
+            recalc_transform=recalc_transform,
+            **profile_kwargs,
+        )


Test patch:
diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py
index 53f0404b..4526589e 100644
--- a/test/integration/test_integration_rioxarray.py
+++ b/test/integration/test_integration_rioxarray.py
@@ -810,7 +810,9 @@ def test_to_raster__custom_description(tmpdir):
         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() if key != "nodata"
+            key: str(value)
+            for key, value in xds.attrs.items()
+            if key not in ("nodata", "long_name")
         }
 
     with rasterio.open(str(tmp_raster)) as rds:
@@ -918,6 +920,72 @@ def test_to_raster__preserve_profile__none_nodata(windowed, tmpdir):
         assert rds.nodata is None
 
 
+def test_to_raster__dataset(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:
+        mda.isel(time=0).rio.to_raster(str(tmp_raster))
+
+    with rioxarray.open_rasterio(str(tmp_raster)) as rdscompare:
+        assert rdscompare.scale_factor == 1.0
+        assert rdscompare.add_offset == 0.0
+        assert rdscompare.long_name == ("blue", "green")
+        assert rdscompare.rio.crs == mda.rio.crs
+        assert numpy.isnan(rdscompare.rio.nodata)
+
+
+def test_to_raster__dataset__mask_and_scale(tmpdir):
+    output_raster = tmpdir.join("tmmx_20190121.tif")
+    with rioxarray.open_rasterio(
+        os.path.join(TEST_INPUT_DATA_DIR, "tmmx_20190121.nc")
+    ) as rds:
+        rds.isel(band=0).rio.to_raster(str(output_raster))
+
+    with rioxarray.open_rasterio(str(output_raster)) as rdscompare:
+        assert rdscompare.scale_factor == 0.1
+        assert rdscompare.add_offset == 220.0
+        assert rdscompare.long_name == "air_temperature"
+        assert rdscompare.rio.crs == rds.rio.crs
+        assert rdscompare.rio.nodata == rds.air_temperature.rio.nodata
+
+
+def test_to_raster__dataset__different_crs(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:
+        rds = mda.isel(time=0)
+        attrs = rds.green.attrs
+        attrs["crs"] = "EPSG:4326"
+        attrs.pop("grid_mapping")
+        rds.green.attrs = attrs
+        attrs = rds.blue.attrs
+        attrs["crs"] = "EPSG:32722"
+        attrs.pop("grid_mapping")
+        rds.blue.attrs = attrs
+        rds = rds.drop_vars("spatial_ref")
+        with pytest.raises(
+            RioXarrayError, match="All CRS must be the same when exporting to raster."
+        ):
+            rds.rio.to_raster(str(tmp_raster))
+
+
+def test_to_raster__dataset__different_nodata(tmpdir):
+    tmp_raster = tmpdir.join("planet_3d_raster.tif")
+    with xarray.open_dataset(
+        os.path.join(TEST_INPUT_DATA_DIR, "PLANET_SCOPE_3D.nc"), mask_and_scale=False,
+    ) as mda:
+        rds = mda.isel(time=0)
+        rds.green.rio.write_nodata(1234, inplace=True)
+        rds.blue.rio.write_nodata(2345, inplace=True)
+        with pytest.raises(
+            RioXarrayError,
+            match="All nodata values must be the same when exporting to raster.",
+        ):
+            rds.rio.to_raster(str(tmp_raster))
+
+
 def test_missing_spatial_dimensions():
     test_ds = xarray.Dataset()
     with pytest.raises(DimensionError):