xarray-contrib · brendancol · Feb 20, 2026 · Feb 20, 2026 · Feb 20, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,17 @@
 -----------
 
 
+### Version 0.5.3 (unreleased)
+- Fixes #134: Add `xr.Dataset` as input type for appropriate modules.
+  Most public functions now transparently accept `xr.Dataset` in addition
+  to `xr.DataArray`. Single-input functions (slope, aspect, curvature,
+  hillshade, focal.mean, all classification functions, proximity/allocation/
+  direction) iterate over data variables and return a Dataset. Multi-input
+  functions (all multispectral indices) accept a Dataset with band-name
+  keyword arguments. `zonal.stats` computes per-variable statistics and
+  returns a merged DataFrame with prefixed columns.
+
+
 ### Version 0.5.2 - 2025-12-18
 - Make dask optional (#835)
 - Fixes 832 update citation info in readme (#834)

diff --git a/README.md b/README.md
@@ -258,6 +258,21 @@ my_dataarray = xr.DataArray(...)
 hillshaded_dataarray = hillshade(my_dataarray)
 ```
 
+##### Dataset Support
+
+Most functions also accept an `xr.Dataset`. Single-input functions (surface, classification, focal, proximity) apply the operation to each data variable and return a Dataset. Multi-input functions (multispectral indices) accept a Dataset with band-name keyword arguments.
+
+```python
+# Single-input: returns a Dataset with slope computed for each variable
+slope_ds = slope(my_dataset)
+
+# Multi-input: map Dataset variables to band parameters
+ndvi_result = ndvi(my_dataset, nir='band_5', red='band_4')
+
+# Zonal stats: columns prefixed by variable name
+stats_df = zonal.stats(zones, my_dataset)  # → elevation_mean, temperature_mean, ...
+```
+
 Check out the user guide [here](/examples/user_guide/).
 
 ------

diff --git a/docs/source/getting_started/usage.rst b/docs/source/getting_started/usage.rst
@@ -15,6 +15,27 @@ Basic Pattern
    my_dataarray = xr.DataArray(...)
    hillshaded_dataarray = hillshade(my_dataarray)
 
+
+Dataset Support
+================
+
+Most functions also accept an ``xr.Dataset``. Single-input functions apply
+the operation to each data variable and return a Dataset. Multi-input
+functions (multispectral indices) accept a Dataset with band-name keyword
+arguments.
+
+.. code-block:: python
+
+   from xrspatial import slope
+   from xrspatial.multispectral import ndvi
+
+   # Single-input: returns a Dataset with slope for each variable
+   slope_ds = slope(my_dataset)
+
+   # Multi-input: map Dataset variables to band parameters
+   ndvi_result = ndvi(my_dataset, nir='band_5', red='band_4')
+
+
 Check out the user guide `here <https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide>`_.
 
 

diff --git a/docs/source/user_guide/data_types.rst b/docs/source/user_guide/data_types.rst
@@ -227,11 +227,49 @@ Best Practices
    combined = (ndvi_result + savi_result) / 2
 
 
+Dataset Input Support
+=====================
+
+Most functions accept an ``xr.Dataset`` in addition to ``xr.DataArray``.
+When a Dataset is passed, the operation is applied to each data variable
+independently and the result is returned as a new Dataset.
+
+Single-input functions (surface, classification, focal, proximity):
+
+.. code-block:: python
+
+   from xrspatial import slope
+
+   # Apply slope to every variable in the Dataset
+   slope_ds = slope(my_dataset)
+   # Returns an xr.Dataset with the same variable names
+
+Multi-input functions (multispectral indices) accept a Dataset with keyword
+arguments that map band aliases to variable names:
+
+.. code-block:: python
+
+   from xrspatial.multispectral import ndvi
+
+   # Map Dataset variables to band parameters
+   ndvi_result = ndvi(my_dataset, nir='band_5', red='band_4')
+
+``zonal.stats`` also accepts a Dataset for the ``values`` parameter, returning
+a merged DataFrame with columns prefixed by variable name:
+
+.. code-block:: python
+
+   from xrspatial.zonal import stats
+
+   df = stats(zones, my_dataset)
+   # Columns: zone, elevation_mean, elevation_max, ..., temperature_mean, ...
+
+
 Summary
 =======
 
-- **Input**: xarray-spatial accepts any numeric data type (int or float)
+- **Input**: xarray-spatial accepts any numeric data type (int or float), as either ``xr.DataArray`` or ``xr.Dataset``
 - **Processing**: All calculations are performed in float32 precision
-- **Output**: Results are returned as float32 DataArrays
+- **Output**: Results are returned as float32 DataArrays (or a Dataset of float32 DataArrays when a Dataset is passed)
 - **Consistency**: This behavior is consistent across NumPy, Dask, and CuPy backends
 - **Rationale**: Float32 provides adequate precision for geospatial analysis while using half the memory of float64
diff --git a/docs/source/user_guide/multispectral.ipynb b/docs/source/user_guide/multispectral.ipynb
@@ -10,13 +10,7 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": [
-    "Xarray-spatial's Multispectral tools provide a range of functions pertaining to remote sensing data such as satellite imagery. A range of functions are available to calculate various vegetation and environmental parameters from the range of band data available for an area. These functions accept and output data in the form of xarray.DataArray rasters.\n",
-    "\n",
-    "- [Generate terrain](#Generate-Terrain-Data)  \n",
-    "- [Bump](#Bump)  \n",
-    "- [NDVI](#NDVI)  "
-   ]
+   "source": "Xarray-spatial's Multispectral tools provide a range of functions pertaining to remote sensing data such as satellite imagery. A range of functions are available to calculate various vegetation and environmental parameters from the range of band data available for an area. These functions accept and output data in the form of xarray.DataArray rasters. They also accept an xr.Dataset as the first argument with band-name keyword arguments to map variables to bands (e.g. `ndvi(ds, nir='B5', red='B4')`).\n\n- [Generate terrain](#Generate-Terrain-Data)  \n- [Bump](#Bump)  \n- [NDVI](#NDVI)  "
   },
   {
    "cell_type": "markdown",
@@ -64,11 +58,7 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": [
-    "The following functions apply to image data with bands in different parts of the UV/Visible/IR spectrum (multispectral), so we'll bring in some multispectral satellite image data to work with.\n",
-    "\n",
-    "Below, we loaded all of the images and transformed them into the form of an xarray DataArray to use in the Xarray-spatial functions."
-   ]
+   "source": "The following functions apply to image data with bands in different parts of the UV/Visible/IR spectrum (multispectral), so we'll bring in some multispectral satellite image data to work with.\n\nBelow, we loaded all of the images and transformed them into the form of an xarray DataArray to use in the Xarray-spatial functions. Note: you can also load bands into an `xr.Dataset` and pass it directly to multispectral functions with band-name keyword arguments (e.g. `ndvi(ds, nir='nir', red='red')`)."
   },
   {
    "cell_type": "code",
@@ -1408,4 +1398,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 2
-}
+}
diff --git a/docs/source/user_guide/surface.ipynb b/docs/source/user_guide/surface.ipynb
@@ -10,17 +10,7 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": [
-    "With the Surface tools, you can quantify and visualize a terrain landform represented by a digital elevation model.\n",
-    "\n",
-    "Starting with a raster elevation surface that represented as an Xarray DataArray, these tools help you in identifying some specific patterns that were not readily apparent in the original surface. Return of each function is also an Xarray DataArray.\n",
-    "\n",
-    "- [Hillshade](#Hillshade): Creates a shaded relief from a surface raster by considering the illumination source angle and shadows.\n",
-    "- [Slope](#Slope): Identifies the slope from each cell of a raster.\n",
-    "- [Curvature](#Curvature): Calculates the curvature of a raster surface.\n",
-    "- [Aspect](#Aspect): Derives the aspect from each cell of a raster surface.\n",
-    "- [Viewshed](#Viewshed): Determines visible locations in the input raster surface from a viewpoint with some optional observer features."
-   ]
+   "source": "With the Surface tools, you can quantify and visualize a terrain landform represented by a digital elevation model.\n\nStarting with a raster elevation surface that represented as an Xarray DataArray (or an Xarray Dataset containing multiple elevation variables), these tools help you in identifying some specific patterns that were not readily apparent in the original surface. When a DataArray is passed, the return is a DataArray. When a Dataset is passed, the function is applied to each variable independently and the return is a Dataset.\n\n- [Hillshade](#Hillshade): Creates a shaded relief from a surface raster by considering the illumination source angle and shadows.\n- [Slope](#Slope): Identifies the slope from each cell of a raster.\n- [Curvature](#Curvature): Calculates the curvature of a raster surface.\n- [Aspect](#Aspect): Derives the aspect from each cell of a raster surface.\n- [Viewshed](#Viewshed): Determines visible locations in the input raster surface from a viewpoint with some optional observer features."
   },
   {
    "cell_type": "markdown",
@@ -932,4 +922,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 2
-}
+}
diff --git a/examples/user_guide/1_Surface.ipynb b/examples/user_guide/1_Surface.ipynb
@@ -3,26 +3,7 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": [
-    "# Xarray-spatial\n",
-    "### User Guide: Surface tools\n",
-    "-----\n",
-    "With the Surface tools, you can quantify and visualize a terrain landform represented by a digital elevation model.\n",
-    "\n",
-    "Starting with a raster elevation surface, represented as an Xarray DataArray, these tools can help you identify some specific patterns that may not be readily apparent in the original surface. The return of each function is also an Xarray DataArray.\n",
-    "\n",
-    "[Hillshade](#Hillshade): Creates a shaded relief from a surface raster by considering the illumination source angle and shadows.\n",
-    "\n",
-    "[Slope](#Slope): Identifies the slope for each cell of a raster.\n",
-    "\n",
-    "[Curvature](#Curvature): Calculates the curvature of a raster surface.\n",
-    "\n",
-    "[Aspect](#Aspect): Derives the aspect for each cell of a raster surface.\n",
-    "\n",
-    "[Viewshed](#Viewshed): Determines visible locations in the input raster surface from a viewpoint with an optional observer height.\n",
-    "\n",
-    "-----------\n"
-   ]
+   "source": "# Xarray-spatial\n### User Guide: Surface tools\n-----\nWith the Surface tools, you can quantify and visualize a terrain landform represented by a digital elevation model.\n\nStarting with a raster elevation surface, represented as an Xarray DataArray (or an Xarray Dataset containing multiple elevation variables), these tools can help you identify some specific patterns that may not be readily apparent in the original surface. When a DataArray is passed, the return is a DataArray. When a Dataset is passed, the function is applied to each variable independently and the return is a Dataset.\n\n[Hillshade](#Hillshade): Creates a shaded relief from a surface raster by considering the illumination source angle and shadows.\n\n[Slope](#Slope): Identifies the slope for each cell of a raster.\n\n[Curvature](#Curvature): Calculates the curvature of a raster surface.\n\n[Aspect](#Aspect): Derives the aspect for each cell of a raster surface.\n\n[Viewshed](#Viewshed): Determines visible locations in the input raster surface from a viewpoint with an optional observer height.\n\n-----------\n"
   },
   {
    "cell_type": "markdown",
@@ -524,4 +505,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 4
-}
+}
diff --git a/examples/user_guide/6_Remote_Sensing.ipynb b/examples/user_guide/6_Remote_Sensing.ipynb
@@ -3,24 +3,7 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": [
-    "# Xarray-spatial\n",
-    "### User Guide: Remote Sensing tools\n",
-    "-----\n",
-    "\n",
-    "Xarray-spatial's Remote Sensing tools provide a range of functions pertaining to remote sensing data such as satellite imagery. A range of functions are available to calculate various vegetation and environmental parameters from the range of band data available for an area. These functions accept and output data in the form of xarray.DataArray rasters.\n",
-    "\n",
-    "[True Color](#True-Color)   \n",
-    "[Vegetation Index](#Vegetation-Index): [NDVI](#NDVI), [SAVI](#SAVI), [ARVI](#ARVI), [EVI](#EVI)   \n",
-    "[Green Chlorophyll Index - GCI](#Green-Chlorophyll-Index-(GCI))  \n",
-    "[Normalized Burn Ratio](#Normalized-Burn-Ratio): [NBR](#NBR), [NBR2](#NBR2)    \n",
-    "[Normalized Difference Moisture Index - NDMI](#Normalized-Difference-Moisture-Index-(NDMI))  \n",
-    "[Structure Insensitive Pigment Index - SIPI](#Structure-Insensitive-Pigment-Index-(SIPI))  \n",
-    "[Enhanced Built-Up and Bareness Index - EBBI](#Enhanced-Built-Up-and-Bareness-Index-(EBBI))  \n",
-    "[Bump Mapping](#Bump-Mapping)  \n",
-    "\n",
-    "-----------\n"
-   ]
+   "source": "# Xarray-spatial\n### User Guide: Remote Sensing tools\n-----\n\nXarray-spatial's Remote Sensing tools provide a range of functions pertaining to remote sensing data such as satellite imagery. A range of functions are available to calculate various vegetation and environmental parameters from the range of band data available for an area. These functions accept and output data in the form of xarray.DataArray rasters. They also accept an xr.Dataset as the first argument with band-name keyword arguments to map variables to bands (e.g. `ndvi(ds, nir='B5', red='B4')`).\n\n[True Color](#True-Color)   \n[Vegetation Index](#Vegetation-Index): [NDVI](#NDVI), [SAVI](#SAVI), [ARVI](#ARVI), [EVI](#EVI)   \n[Green Chlorophyll Index - GCI](#Green-Chlorophyll-Index-(GCI))  \n[Normalized Burn Ratio](#Normalized-Burn-Ratio): [NBR](#NBR), [NBR2](#NBR2)    \n[Normalized Difference Moisture Index - NDMI](#Normalized-Difference-Moisture-Index-(NDMI))  \n[Structure Insensitive Pigment Index - SIPI](#Structure-Insensitive-Pigment-Index-(SIPI))  \n[Enhanced Built-Up and Bareness Index - EBBI](#Enhanced-Built-Up-and-Bareness-Index-(EBBI))  \n[Bump Mapping](#Bump-Mapping)  \n\n-----------\n"
   },
   {
    "cell_type": "markdown",
@@ -68,11 +51,7 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": [
-    "The following functions apply to image data with bands in different parts of the UV/Visible/IR spectrum (multispectral), so we'll bring in some multispectral satellite image data to work with.\n",
-    "\n",
-    "Below, we loaded all of the images and transformed them into the form of an xarray DataArray to use in the Xarray-spatial functions."
-   ]
+   "source": "The following functions apply to image data with bands in different parts of the UV/Visible/IR spectrum (multispectral), so we'll bring in some multispectral satellite image data to work with.\n\nBelow, we loaded all of the images and transformed them into the form of an xarray DataArray to use in the Xarray-spatial functions. Note: you can also load bands into an `xr.Dataset` and pass it directly to multispectral functions with band-name keyword arguments (e.g. `ndvi(ds, nir='nir', red='red')`)."
   },
   {
    "cell_type": "code",
@@ -695,4 +674,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 4
-}
+}
diff --git a/xrspatial/aspect.py b/xrspatial/aspect.py
@@ -19,6 +19,7 @@
 from xrspatial.utils import _extract_latlon_coords
 from xrspatial.utils import cuda_args
 from xrspatial.utils import ngjit
+from xrspatial.dataset_support import supports_dataset
 
 
 def _geodesic_cuda_dims(shape):
@@ -270,6 +271,7 @@ def _run_dask_cupy_geodesic(data, lat_2d, lon_2d, a2, b2, z_factor):
 # Public API
 # =====================================================================
 
+@supports_dataset
 def aspect(agg: xr.DataArray,
            name: Optional[str] = 'aspect',
            method: str = 'planar',
@@ -296,9 +298,11 @@ def aspect(agg: xr.DataArray,
 
     Parameters
     ----------
-    agg : xarray.DataArray
+    agg : xarray.DataArray or xr.Dataset
         2D NumPy, CuPy, or Dask with NumPy-backed xarray DataArray
         of elevation values.
+        If a Dataset is passed, the operation is applied to each
+        data variable independently.
     name : str, default='aspect'
         Name of ouput DataArray.
     method : str, default='planar'
@@ -313,7 +317,10 @@ def aspect(agg: xr.DataArray,
 
     Returns
     -------
-    aspect_agg : xarray.DataArray of the same type as `agg`
+    aspect_agg : xarray.DataArray or xr.Dataset
+        If `agg` is a DataArray, returns a DataArray of the same type.
+        If `agg` is a Dataset, returns a Dataset with aspect computed
+        for each data variable.
         2D aggregate array of calculated aspect values.
         All other input attributes are preserved.