Add kriging user guide notebook (#951)

Covers basic interpolation, variance output, variogram model
comparison, and a soil pH mapping example.

Author: brendancol

examples/user_guide/17_Kriging.ipynb

@@ -0,0 +1,230 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "cell-0",
+   "metadata": {},
+   "source": [
+    "# Kriging Interpolation\n",
+    "\n",
+    "The `kriging()` function performs Ordinary Kriging, a geostatistical interpolation method that produces optimal, unbiased predictions from scattered point observations. Unlike IDW, kriging accounts for the spatial correlation structure of the data through a variogram model.\n",
+    "\n",
+    "Key features:\n",
+    "- Automatic experimental variogram computation and model fitting\n",
+    "- Three variogram models: spherical, exponential, gaussian\n",
+    "- Optional kriging variance (prediction uncertainty) output\n",
+    "- All four backends: NumPy, Dask, CuPy, Dask+CuPy"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cell-1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import xarray as xr\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "from xrspatial.interpolate import kriging"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cell-2",
+   "metadata": {},
+   "source": [
+    "## 1. Basic interpolation from point observations\n",
+    "\n",
+    "Generate scattered sample points from a known surface and use kriging to reconstruct the full field."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cell-3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# True surface: z = sin(x) * cos(y)\n",
+    "rng = np.random.RandomState(42)\n",
+    "n_pts = 40\n",
+    "x_pts = rng.uniform(0, 6, n_pts)\n",
+    "y_pts = rng.uniform(0, 6, n_pts)\n",
+    "z_pts = np.sin(x_pts) * np.cos(y_pts) + rng.normal(0, 0.05, n_pts)\n",
+    "\n",
+    "# Output grid\n",
+    "x_grid = np.linspace(0, 6, 60)\n",
+    "y_grid = np.linspace(0, 6, 60)\n",
+    "template = xr.DataArray(\n",
+    "    np.zeros((len(y_grid), len(x_grid))),\n",
+    "    dims=['y', 'x'],\n",
+    "    coords={'y': y_grid, 'x': x_grid},\n",
+    ")\n",
+    "\n",
+    "result = kriging(x_pts, y_pts, z_pts, template)\n",
+    "\n",
+    "# True surface for comparison\n",
+    "gx, gy = np.meshgrid(x_grid, y_grid)\n",
+    "true_surface = np.sin(gx) * np.cos(gy)\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(12, 5))\n",
+    "im0 = axes[0].imshow(true_surface, extent=[0, 6, 6, 0], cmap='viridis')\n",
+    "axes[0].scatter(x_pts, y_pts, c='red', s=15, edgecolors='k', linewidth=0.5)\n",
+    "axes[0].set_title('True surface + sample points')\n",
+    "fig.colorbar(im0, ax=axes[0], shrink=0.7)\n",
+    "\n",
+    "im1 = axes[1].imshow(result.values, extent=[0, 6, 6, 0], cmap='viridis')\n",
+    "axes[1].set_title('Kriging prediction')\n",
+    "fig.colorbar(im1, ax=axes[1], shrink=0.7)\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cell-4",
+   "metadata": {},
+   "source": [
+    "## 2. Kriging variance\n",
+    "\n",
+    "Set `return_variance=True` to get prediction uncertainty. Variance is low near observed points and higher in data-sparse regions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cell-5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pred, var = kriging(x_pts, y_pts, z_pts, template, return_variance=True)\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(12, 5))\n",
+    "im0 = axes[0].imshow(pred.values, extent=[0, 6, 6, 0], cmap='viridis')\n",
+    "axes[0].scatter(x_pts, y_pts, c='red', s=15, edgecolors='k', linewidth=0.5)\n",
+    "axes[0].set_title('Prediction')\n",
+    "fig.colorbar(im0, ax=axes[0], shrink=0.7)\n",
+    "\n",
+    "im1 = axes[1].imshow(var.values, extent=[0, 6, 6, 0], cmap='magma')\n",
+    "axes[1].scatter(x_pts, y_pts, c='cyan', s=15, edgecolors='k', linewidth=0.5)\n",
+    "axes[1].set_title('Kriging variance')\n",
+    "fig.colorbar(im1, ax=axes[1], shrink=0.7)\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cell-6",
+   "metadata": {},
+   "source": [
+    "## 3. Variogram model comparison\n",
+    "\n",
+    "The `variogram_model` parameter controls the spatial correlation model. Different models produce subtly different predictions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cell-7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "models = ['spherical', 'exponential', 'gaussian']\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
+    "\n",
+    "for ax, model in zip(axes, models):\n",
+    "    r = kriging(x_pts, y_pts, z_pts, template, variogram_model=model)\n",
+    "    im = ax.imshow(r.values, extent=[0, 6, 6, 0], cmap='viridis')\n",
+    "    ax.set_title(f'{model}')\n",
+    "    fig.colorbar(im, ax=ax, shrink=0.7)\n",
+    "\n",
+    "plt.suptitle('Variogram model comparison', y=1.02)\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cell-8",
+   "metadata": {},
+   "source": [
+    "## 4. Practical example: soil property mapping\n",
+    "\n",
+    "Simulate soil pH measurements at random field locations and produce a continuous map with uncertainty."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cell-9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Simulated soil pH: smooth trend + spatially correlated noise\n",
+    "rng = np.random.RandomState(7)\n",
+    "n_samples = 50\n",
+    "x_soil = rng.uniform(0, 100, n_samples)  # meters\n",
+    "y_soil = rng.uniform(0, 100, n_samples)\n",
+    "\n",
+    "# Trend: pH increases toward the northeast\n",
+    "ph = 5.5 + 0.015 * x_soil + 0.010 * y_soil + rng.normal(0, 0.3, n_samples)\n",
+    "\n",
+    "# Dense prediction grid\n",
+    "xg = np.linspace(0, 100, 80)\n",
+    "yg = np.linspace(0, 100, 80)\n",
+    "template_soil = xr.DataArray(\n",
+    "    np.zeros((len(yg), len(xg))),\n",
+    "    dims=['y', 'x'],\n",
+    "    coords={'y': yg, 'x': xg},\n",
+    ")\n",
+    "\n",
+    "ph_pred, ph_var = kriging(\n",
+    "    x_soil, y_soil, ph, template_soil,\n",
+    "    variogram_model='spherical', return_variance=True,\n",
+    ")\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(13, 5))\n",
+    "\n",
+    "im0 = axes[0].imshow(\n",
+    "    ph_pred.values, extent=[0, 100, 100, 0],\n",
+    "    cmap='RdYlGn', vmin=5, vmax=8,\n",
+    ")\n",
+    "axes[0].scatter(x_soil, y_soil, c=ph, cmap='RdYlGn', vmin=5, vmax=8,\n",
+    "                s=30, edgecolors='k', linewidth=0.5)\n",
+    "axes[0].set_title('Predicted soil pH')\n",
+    "axes[0].set_xlabel('East (m)')\n",
+    "axes[0].set_ylabel('North (m)')\n",
+    "fig.colorbar(im0, ax=axes[0], shrink=0.7, label='pH')\n",
+    "\n",
+    "im1 = axes[1].imshow(\n",
+    "    ph_var.values, extent=[0, 100, 100, 0], cmap='magma',\n",
+    ")\n",
+    "axes[1].scatter(x_soil, y_soil, c='cyan', s=15, edgecolors='k', linewidth=0.5)\n",
+    "axes[1].set_title('Prediction variance')\n",
+    "axes[1].set_xlabel('East (m)')\n",
+    "fig.colorbar(im1, ax=axes[1], shrink=0.7, label='Variance')\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.10.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}