Add COG overview generation user guide notebook (#1150)

Covers automatic and explicit overview levels, resampling methods,
and round-trip verification. Self-contained with synthetic data.

Author: brendancol

examples/user_guide/50_COG_Overview_Generation.ipynb

@@ -0,0 +1,230 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Cloud Optimized GeoTIFF: Overview (Pyramid) Generation\n",
+    "\n",
+    "Cloud Optimized GeoTIFFs (COGs) store internal overview images at reduced\n",
+    "resolution, so map viewers and tiling servers can fetch the right zoom\n",
+    "level without downloading the full file.\n",
+    "\n",
+    "xarray-spatial generates these overviews natively during `to_geotiff()`\n",
+    "with `cog=True`. No GDAL or `gdaladdo` post-processing required.\n",
+    "\n",
+    "This notebook covers:\n",
+    "- Writing a COG with automatic overviews\n",
+    "- Choosing resampling methods\n",
+    "- Specifying explicit overview levels\n",
+    "- Verifying the result"
+   ]
+  },
+  {
+   "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.geotiff import open_geotiff, to_geotiff\n",
+    "from xrspatial.geotiff._header import parse_header, parse_all_ifds"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Create synthetic terrain data\n",
+    "\n",
+    "A 256x256 elevation surface with some peaks, to give the overviews\n",
+    "something visually meaningful to downsample."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rng = np.random.RandomState(42)\n",
+    "y = np.linspace(45.0, 44.0, 256)\n",
+    "x = np.linspace(-120.0, -119.0, 256)\n",
+    "\n",
+    "# Smooth terrain from Gaussian peaks\n",
+    "xx, yy = np.meshgrid(x, y)\n",
+    "terrain = (\n",
+    "    800 * np.exp(-((xx + 119.5)**2 + (yy - 44.5)**2) / 0.02)\n",
+    "    + 600 * np.exp(-((xx + 119.3)**2 + (yy - 44.7)**2) / 0.05)\n",
+    "    + 100 * rng.rand(256, 256)\n",
+    ").astype(np.float32)\n",
+    "\n",
+    "da = xr.DataArray(\n",
+    "    terrain, dims=['y', 'x'],\n",
+    "    coords={'y': y, 'x': x},\n",
+    "    attrs={'crs': 4326},\n",
+    "    name='elevation',\n",
+    ")\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(6, 5))\n",
+    "da.plot(ax=ax, cmap='terrain')\n",
+    "ax.set_title('Synthetic elevation')\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Write a COG with automatic overviews\n",
+    "\n",
+    "Pass `cog=True` and xarray-spatial handles the rest: it halves the\n",
+    "dimensions until the smallest overview fits within a single tile,\n",
+    "then writes all levels into the file with IFDs at the start (the COG\n",
+    "layout requirement for efficient HTTP range requests)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import tempfile, os\n",
+    "\n",
+    "tmpdir = tempfile.mkdtemp(prefix='cog_demo_')\n",
+    "cog_path = os.path.join(tmpdir, 'auto_overviews.tif')\n",
+    "\n",
+    "to_geotiff(da, cog_path, cog=True, compression='deflate')\n",
+    "\n",
+    "# Check how many IFDs (resolution levels) were written\n",
+    "with open(cog_path, 'rb') as f:\n",
+    "    raw = f.read()\n",
+    "\n",
+    "header = parse_header(raw)\n",
+    "ifds = parse_all_ifds(raw, header)\n",
+    "for i, ifd in enumerate(ifds):\n",
+    "    label = 'Full resolution' if i == 0 else f'Overview {i}'\n",
+    "    print(f'{label}: {ifd.width} x {ifd.height}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Resampling methods\n",
+    "\n",
+    "Different data types need different resampling:\n",
+    "- **mean** (default): continuous data like elevation or temperature\n",
+    "- **nearest**: categorical or index data (land cover classes)\n",
+    "- **mode**: majority-class downsampling for classified rasters\n",
+    "- **min** / **max** / **median**: when you need conservative bounds"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "methods = ['mean', 'nearest', 'min', 'max']\n",
+    "fig, axes = plt.subplots(1, len(methods), figsize=(14, 3))\n",
+    "\n",
+    "for ax, method in zip(axes, methods):\n",
+    "    path = os.path.join(tmpdir, f'cog_{method}.tif')\n",
+    "    to_geotiff(da, path, cog=True, compression='deflate',\n",
+    "               overview_resampling=method)\n",
+    "\n",
+    "    # Read back the first overview level\n",
+    "    result = open_geotiff(path, overview_level=1)\n",
+    "    result.plot(ax=ax, cmap='terrain', add_colorbar=False)\n",
+    "    ax.set_title(f'{method} ({result.shape[0]}x{result.shape[1]})')\n",
+    "    ax.set_xlabel('')\n",
+    "    ax.set_ylabel('')\n",
+    "\n",
+    "plt.suptitle('Overview level 1 by resampling method', y=1.02)\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Explicit overview levels\n",
+    "\n",
+    "For fine-grained control, pass `overview_levels` as a list of\n",
+    "decimation factors. Each factor halves the previous level once."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "explicit_path = os.path.join(tmpdir, 'explicit_levels.tif')\n",
+    "to_geotiff(da, explicit_path, cog=True, compression='deflate',\n",
+    "           overview_levels=[2, 4, 8])\n",
+    "\n",
+    "with open(explicit_path, 'rb') as f:\n",
+    "    raw = f.read()\n",
+    "\n",
+    "header = parse_header(raw)\n",
+    "ifds = parse_all_ifds(raw, header)\n",
+    "for i, ifd in enumerate(ifds):\n",
+    "    label = 'Full resolution' if i == 0 else f'Overview {i}'\n",
+    "    print(f'{label}: {ifd.width} x {ifd.height}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Verify round-trip\n",
+    "\n",
+    "Full-resolution pixel values are preserved through the COG write."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "result = open_geotiff(cog_path)\n",
+    "max_diff = float(np.max(np.abs(result.values - terrain)))\n",
+    "print(f'Max pixel difference: {max_diff}')\n",
+    "assert max_diff < 1e-5, 'Values should round-trip exactly'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Clean up\n",
+    "import shutil\n",
+    "shutil.rmtree(tmpdir)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}