Add rescale and standardize normalization utilities

README.md

@@ -352,6 +352,8 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e
 | Name | Description | Source | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |
 |:----------:|:------------|:------:|:----------------------:|:--------------------:|:-------------------:|:------:|
 | [Preview](xrspatial/preview.py) | Downsamples a raster to target pixel dimensions for visualization | Custom | ✅️ | ✅️ | ✅️ | 🔄 |
+| [Rescale](xrspatial/normalize.py) | Min-max normalization to a target range (default [0, 1]) | Standard | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Standardize](xrspatial/normalize.py) | Z-score normalization (subtract mean, divide by std) | Standard | ✅️ | ✅️ | ✅️ | ✅️ |
 
 #### Usage
 

docs/source/reference/utilities.rst

@@ -46,6 +46,14 @@ Preview
 
     xrspatial.preview.preview
 
+Normalization
+=============
+.. autosummary::
+    :toctree: _autosummary
+
+    xrspatial.normalize.rescale
+    xrspatial.normalize.standardize
+
 Diagnostics
 ===========
 .. autosummary::

examples/user_guide/34_Normalization.ipynb

@@ -0,0 +1,101 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "x2762y1yuan",
+   "source": "# Normalization: rescale and standardize\n\nTwo common preprocessing steps before combining rasters or feeding them into models:\n\n- **rescale** maps values to a target range (default [0, 1]) using min-max normalization.\n- **standardize** centers values at zero with unit variance (z-score normalization).\n\nBoth functions handle NaN and infinite values (they pass through unchanged) and work on all four xarray-spatial backends: NumPy, CuPy, Dask+NumPy, and Dask+CuPy.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "wzy5yzkbde",
+   "source": "%matplotlib inline\nimport numpy as np\nimport xarray as xr\nimport matplotlib.pyplot as plt\n\nfrom xrspatial.normalize import rescale, standardize\nfrom xrspatial import generate_terrain",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5zey0oy2cji",
+   "source": "## Synthetic terrain\n\nGenerate a 500x500 elevation raster with values roughly in the 0-1200 range. We'll sprinkle in a few NaN cells to show how they're preserved.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "1n68h492hdu",
+   "source": "terrain = generate_terrain(canvas=xr.DataArray(np.zeros((500, 500)), dims=['y', 'x']))\n\n# Add some NaN holes\nrng = np.random.default_rng(42)\nmask = rng.random(terrain.shape) < 0.005\nterrain.values[mask] = np.nan\n\nprint(f\"Shape: {terrain.shape}\")\nprint(f\"Range: {float(np.nanmin(terrain)):.1f} to {float(np.nanmax(terrain)):.1f}\")\nprint(f\"NaN cells: {int(np.isnan(terrain.values).sum())}\")\n\nfig, ax = plt.subplots(figsize=(7, 6))\nterrain.plot.imshow(ax=ax, cmap='terrain', add_colorbar=True,\n                    cbar_kwargs={'label': 'Elevation'})\nax.set_title('Raw elevation')\nax.set_axis_off()\nplt.tight_layout()",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "h7fv68zd0dq",
+   "source": "## rescale: min-max normalization\n\nBy default, `rescale()` maps finite values to [0, 1]. You can supply a custom range with `new_min` and `new_max`.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "z9j191ufbs",
+   "source": "scaled_01 = rescale(terrain)\nscaled_byte = rescale(terrain, new_min=0, new_max=255)\n\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\n\nterrain.plot.imshow(ax=axes[0], cmap='terrain', add_colorbar=True)\naxes[0].set_title('Original')\naxes[0].set_axis_off()\n\nscaled_01.plot.imshow(ax=axes[1], cmap='terrain', add_colorbar=True)\naxes[1].set_title('rescale() -> [0, 1]')\naxes[1].set_axis_off()\n\nscaled_byte.plot.imshow(ax=axes[2], cmap='terrain', add_colorbar=True)\naxes[2].set_title('rescale(0, 255)')\naxes[2].set_axis_off()\n\nplt.tight_layout()\n\nprint(f\"[0,1] range:   {float(np.nanmin(scaled_01)):.4f} to {float(np.nanmax(scaled_01)):.4f}\")\nprint(f\"[0,255] range: {float(np.nanmin(scaled_byte)):.1f} to {float(np.nanmax(scaled_byte)):.1f}\")\nprint(f\"NaN preserved: {int(np.isnan(scaled_01.values).sum())} cells\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "u8clrxl4r2i",
+   "source": "## standardize: z-score normalization\n\n`standardize()` subtracts the mean and divides by the standard deviation of finite values. The result has mean ~0 and std ~1. Use `ddof=1` for sample standard deviation.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "e64naydwlya",
+   "source": "zscored = standardize(terrain)\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\n\nterrain.plot.imshow(ax=axes[0], cmap='terrain', add_colorbar=True)\naxes[0].set_title('Original')\naxes[0].set_axis_off()\n\nzscored.plot.imshow(ax=axes[1], cmap='RdBu_r', add_colorbar=True,\n                    cbar_kwargs={'label': 'Z-score'})\naxes[1].set_title('standardize()')\naxes[1].set_axis_off()\n\nplt.tight_layout()\n\nfinite = zscored.values[np.isfinite(zscored.values)]\nprint(f\"Mean:  {finite.mean():.2e}\")\nprint(f\"Std:   {finite.std():.6f}\")\nprint(f\"Range: {finite.min():.3f} to {finite.max():.3f}\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0ql3g8hvphj",
+   "source": "## Practical use case: combining layers with different scales\n\nWhen combining elevation and slope into a composite index, the raw values live on different scales. Rescaling both to [0, 1] puts them on equal footing.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "9hhldlalk7f",
+   "source": "from xrspatial import slope\n\nslp = slope(terrain)\n\n# Raw values are on very different scales\nprint(f\"Elevation range: {float(np.nanmin(terrain)):.0f} to {float(np.nanmax(terrain)):.0f}\")\nprint(f\"Slope range:     {float(np.nanmin(slp)):.1f} to {float(np.nanmax(slp)):.1f}\")\n\n# Rescale both to [0, 1] and combine\nelev_norm = rescale(terrain)\nslope_norm = rescale(slp)\n\n# Simple composite: high elevation + steep slope = high risk\ncomposite = 0.6 * elev_norm + 0.4 * slope_norm\n\nfig, axes = plt.subplots(1, 3, figsize=(18, 5))\n\nelev_norm.plot.imshow(ax=axes[0], cmap='terrain', add_colorbar=True)\naxes[0].set_title('Elevation [0, 1]')\naxes[0].set_axis_off()\n\nslope_norm.plot.imshow(ax=axes[1], cmap='YlOrRd', add_colorbar=True)\naxes[1].set_title('Slope [0, 1]')\naxes[1].set_axis_off()\n\ncomposite.plot.imshow(ax=axes[2], cmap='inferno', add_colorbar=True)\naxes[2].set_title('Weighted composite (0.6 elev + 0.4 slope)')\naxes[2].set_axis_off()\n\nplt.tight_layout()",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c6w0w124c4c",
+   "source": "## Accessor syntax\n\nBoth functions are available through the `.xrs` accessor on DataArrays.\n\n```python\nimport xrspatial\n\nterrain.xrs.rescale()\nterrain.xrs.standardize(ddof=1)\n```",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "gcc8pi1a7u",
+   "source": "import xrspatial  # registers .xrs accessor\n\naccessor_result = terrain.xrs.rescale()\nnp.testing.assert_array_equal(accessor_result.values, scaled_01.values)\nprint(\"Accessor output matches function output.\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

xrspatial/__init__.py

@@ -67,6 +67,8 @@
 from xrspatial.multispectral import sipi  # noqa
 from xrspatial.pathfinding import a_star_search  # noqa
 from xrspatial.pathfinding import multi_stop_search  # noqa
+from xrspatial.normalize import rescale  # noqa
+from xrspatial.normalize import standardize  # noqa
 from xrspatial.perlin import perlin  # noqa
 from xrspatial.preview import preview  # noqa
 from xrspatial.proximity import allocation  # noqa

xrspatial/accessor.py

@@ -321,6 +321,16 @@ def preview(self, **kwargs):
         from .preview import preview
         return preview(self._obj, **kwargs)
 
+    # ---- Normalization ----
+
+    def rescale(self, **kwargs):
+        from .normalize import rescale
+        return rescale(self._obj, **kwargs)
+
+    def standardize(self, **kwargs):
+        from .normalize import standardize
+        return standardize(self._obj, **kwargs)
+
     # ---- Raster to vector ----
 
     def polygonize(self, **kwargs):
@@ -637,6 +647,16 @@ def preview(self, **kwargs):
         from .preview import preview
         return preview(self._obj, **kwargs)
 
+    # ---- Normalization ----
+
+    def rescale(self, **kwargs):
+        from .normalize import rescale
+        return rescale(self._obj, **kwargs)
+
+    def standardize(self, **kwargs):
+        from .normalize import standardize
+        return standardize(self._obj, **kwargs)
+
     # ---- Fire ----
 
     def burn_severity_class(self, **kwargs):