Add morphological raster operators (#949)

README.md

@@ -171,6 +171,17 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e
 
 -------
 
+### **Morphological**
+
+| Name | Description | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |
+|:----------:|:------------|:----------------------:|:--------------------:|:-------------------:|:------:|
+| [Erode](xrspatial/morphology.py) | Morphological erosion (local minimum over structuring element) | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Dilate](xrspatial/morphology.py) | Morphological dilation (local maximum over structuring element) | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Opening](xrspatial/morphology.py) | Erosion then dilation (removes small bright features) | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Closing](xrspatial/morphology.py) | Dilation then erosion (fills small dark gaps) | ✅️ | ✅️ | ✅️ | ✅️ |
+
+-------
+
 ### **Fire**
 
 | Name | Description | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |

docs/source/reference/index.rst

@@ -15,6 +15,7 @@ Reference
    focal
    hydrology
    interpolation
+   morphology
    multispectral
    pathfinding
    proximity

docs/source/reference/morphology.rst

@@ -0,0 +1,40 @@
+..  _reference.morphology:
+
+**********
+Morphology
+**********
+
+Erode
+=====
+.. autosummary::
+   :toctree: _autosummary
+
+   xrspatial.morphology.morph_erode
+
+Dilate
+======
+.. autosummary::
+   :toctree: _autosummary
+
+   xrspatial.morphology.morph_dilate
+
+Opening
+=======
+.. autosummary::
+   :toctree: _autosummary
+
+   xrspatial.morphology.morph_opening
+
+Closing
+=======
+.. autosummary::
+   :toctree: _autosummary
+
+   xrspatial.morphology.morph_closing
+
+Kernel Construction
+===================
+.. autosummary::
+   :toctree: _autosummary
+
+   xrspatial.morphology._circle_kernel

examples/user_guide/17_Morphological_Operators.ipynb

@@ -0,0 +1,217 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Morphological Operators\n",
+    "\n",
+    "This notebook demonstrates the four morphological raster operators\n",
+    "in `xrspatial.morphology`:\n",
+    "\n",
+    "- **`morph_erode`** -- local minimum (shrinks bright regions)\n",
+    "- **`morph_dilate`** -- local maximum (expands bright regions)\n",
+    "- **`morph_opening`** -- erosion then dilation (removes small bright features)\n",
+    "- **`morph_closing`** -- dilation then erosion (fills small dark gaps)\n",
+    "\n",
+    "These operations work on any numeric raster, including grayscale\n",
+    "elevation data and binary classification masks.  All four backends\n",
+    "(numpy, cupy, dask+numpy, dask+cupy) are supported."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import xarray as xr\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "from xrspatial.morphology import (\n",
+    "    _circle_kernel,\n",
+    "    morph_closing,\n",
+    "    morph_dilate,\n",
+    "    morph_erode,\n",
+    "    morph_opening,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 1. Synthetic elevation data\n",
+    "\n",
+    "Create a small raster with a peak and a pit to show how erosion and\n",
+    "dilation behave."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rng = np.random.default_rng(42)\n",
+    "y, x = np.mgrid[-3:3:0.05, -3:3:0.05]\n",
+    "z = np.exp(-(x**2 + y**2)) - 0.5 * np.exp(-((x - 1.5)**2 + (y - 1)**2) / 0.3)\n",
+    "z += 0.05 * rng.standard_normal(z.shape)  # add noise\n",
+    "\n",
+    "raster = xr.DataArray(z, dims=['y', 'x'], name='elevation')\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(6, 5))\n",
+    "raster.plot(ax=ax, cmap='terrain')\n",
+    "ax.set_title('Original raster')\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 2. Erosion and dilation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "kernel = np.ones((5, 5), dtype=np.uint8)\n",
+    "\n",
+    "eroded = morph_erode(raster, kernel=kernel, boundary='nearest')\n",
+    "dilated = morph_dilate(raster, kernel=kernel, boundary='nearest')\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
+    "for ax, data, title in zip(\n",
+    "    axes,\n",
+    "    [raster, eroded, dilated],\n",
+    "    ['Original', 'Eroded (local min)', 'Dilated (local max)'],\n",
+    "):\n",
+    "    data.plot(ax=ax, cmap='terrain', vmin=z.min(), vmax=z.max())\n",
+    "    ax.set_title(title)\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 3. Opening and closing\n",
+    "\n",
+    "Opening removes small bright spikes; closing fills small dark pits."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "opened = morph_opening(raster, kernel=kernel, boundary='nearest')\n",
+    "closed = morph_closing(raster, kernel=kernel, boundary='nearest')\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
+    "for ax, data, title in zip(\n",
+    "    axes,\n",
+    "    [raster, opened, closed],\n",
+    "    ['Original', 'Opening', 'Closing'],\n",
+    "):\n",
+    "    data.plot(ax=ax, cmap='terrain', vmin=z.min(), vmax=z.max())\n",
+    "    ax.set_title(title)\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 4. Binary mask cleanup\n",
+    "\n",
+    "A common use case: clean up noisy classification results.  Opening\n",
+    "removes salt noise; closing removes pepper noise."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Generate a binary mask with noise\n",
+    "mask = np.zeros((100, 100), dtype=np.float64)\n",
+    "mask[20:80, 20:80] = 1.0  # large square\n",
+    "\n",
+    "# Sprinkle salt-and-pepper noise\n",
+    "noise = rng.random(mask.shape)\n",
+    "mask[noise < 0.02] = 1.0   # salt\n",
+    "mask[noise > 0.98] = 0.0   # pepper\n",
+    "\n",
+    "noisy = xr.DataArray(mask, dims=['y', 'x'], name='mask')\n",
+    "cleaned = morph_opening(\n",
+    "    morph_closing(noisy, kernel=kernel, boundary='nearest'),\n",
+    "    kernel=kernel,\n",
+    "    boundary='nearest',\n",
+    ")\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n",
+    "noisy.plot(ax=axes[0], cmap='gray')\n",
+    "axes[0].set_title('Noisy mask')\n",
+    "cleaned.plot(ax=axes[1], cmap='gray')\n",
+    "axes[1].set_title('After closing + opening')\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 5. Circular structuring element\n",
+    "\n",
+    "Use `_circle_kernel(radius)` to create a disk-shaped kernel\n",
+    "instead of a square."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "disk = _circle_kernel(3)\n",
+    "print('Circular kernel (radius=3):')\n",
+    "print(disk)\n",
+    "\n",
+    "eroded_disk = morph_erode(raster, kernel=disk, boundary='nearest')\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n",
+    "raster.plot(ax=axes[0], cmap='terrain')\n",
+    "axes[0].set_title('Original')\n",
+    "eroded_disk.plot(ax=axes[1], cmap='terrain')\n",
+    "axes[1].set_title('Eroded with circular kernel (r=3)')\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

xrspatial/__init__.py

@@ -41,6 +41,10 @@
 from xrspatial.flow_length import flow_length  # noqa
 from xrspatial.flow_path import flow_path  # noqa
 from xrspatial.focal import mean  # noqa
+from xrspatial.morphology import morph_closing  # noqa
+from xrspatial.morphology import morph_dilate  # noqa
+from xrspatial.morphology import morph_erode  # noqa
+from xrspatial.morphology import morph_opening  # noqa
 from xrspatial.hand import hand  # noqa
 from xrspatial.hillshade import hillshade  # noqa
 from xrspatial.mahalanobis import mahalanobis  # noqa

xrspatial/accessor.py

@@ -173,6 +173,24 @@ def focal_mean(self, **kwargs):
         from .focal import mean
         return mean(self._obj, **kwargs)
 
+    # ---- Morphological ----
+
+    def morph_erode(self, **kwargs):
+        from .morphology import morph_erode
+        return morph_erode(self._obj, **kwargs)
+
+    def morph_dilate(self, **kwargs):
+        from .morphology import morph_dilate
+        return morph_dilate(self._obj, **kwargs)
+
+    def morph_opening(self, **kwargs):
+        from .morphology import morph_opening
+        return morph_opening(self._obj, **kwargs)
+
+    def morph_closing(self, **kwargs):
+        from .morphology import morph_closing
+        return morph_closing(self._obj, **kwargs)
+
     # ---- Proximity / Distance ----
 
     def proximity(self, **kwargs):
@@ -501,6 +519,24 @@ def focal_mean(self, **kwargs):
         from .focal import mean
         return mean(self._obj, **kwargs)
 
+    # ---- Morphological ----
+
+    def morph_erode(self, **kwargs):
+        from .morphology import morph_erode
+        return morph_erode(self._obj, **kwargs)
+
+    def morph_dilate(self, **kwargs):
+        from .morphology import morph_dilate
+        return morph_dilate(self._obj, **kwargs)
+
+    def morph_opening(self, **kwargs):
+        from .morphology import morph_opening
+        return morph_opening(self._obj, **kwargs)
+
+    def morph_closing(self, **kwargs):
+        from .morphology import morph_closing
+        return morph_closing(self._obj, **kwargs)
+
     # ---- Diffusion ----
 
     def diffuse(self, **kwargs):