Add GPU COG overview support (#1150)

README.md

@@ -145,6 +145,7 @@ Native GeoTIFF and Cloud Optimized GeoTIFF reader/writer. No GDAL required.
 |:-----|:------------|:-----:|:----:|:--------:|:-------------:|:-----:|
 | [open_geotiff](xrspatial/geotiff/__init__.py) | Read GeoTIFF / COG / VRT | ✅️ | ✅️ | ✅️ | ✅️ | ✅️ |
 | [to_geotiff](xrspatial/geotiff/__init__.py) | Write DataArray as GeoTIFF / COG | ✅️ | ✅️ | ✅️ | ✅️ | ✅️ |
+| [write_geotiff_gpu](xrspatial/geotiff/__init__.py) | GPU-accelerated GeoTIFF / COG write | | | ✅️ | | |
 | [write_vrt](xrspatial/geotiff/__init__.py) | Generate VRT mosaic from GeoTIFFs | ✅️ | | | | |
 
 `open_geotiff` and `to_geotiff` auto-dispatch to the correct backend:
@@ -163,6 +164,10 @@ open_geotiff('mosaic.vrt')                           # VRT mosaic (auto-detected
 to_geotiff(cupy_array, 'out.tif')                    # auto-detects GPU
 to_geotiff(data, 'out.tif', gpu=True)                # force GPU compress
 to_geotiff(data, 'ortho.tif', compression='jpeg')    # JPEG for orthophotos
+to_geotiff(data, 'cog.tif', cog=True)                # COG with auto overviews
+to_geotiff(data, 'cog.tif', cog=True,                # COG with explicit levels
+           overview_levels=[2, 4, 8],
+           overview_resampling='nearest')
 write_vrt('mosaic.vrt', ['tile1.tif', 'tile2.tif'])  # generate VRT
 
 open_geotiff('dem.tif', dtype='float32')             # half memory
@@ -189,7 +194,7 @@ ds.xrs.open_geotiff('large_dem.tif')                 # read windowed to Dataset
 - Cloud storage: S3 (`s3://`), GCS (`gs://`), Azure (`az://`) via fsspec
 - GPUDirect Storage: SSD→GPU direct DMA via KvikIO (optional)
 - Thread-safe mmap reads, atomic writes, HTTP connection reuse (urllib3)
-- Overview generation: mean, nearest, min, max, median, mode, cubic
+- Overview generation (CPU and GPU): mean, nearest, min, max, median, mode, cubic
 - Planar config, big-endian byte swap, PixelIsArea/PixelIsPoint
 
 **Read performance** (real-world files, A6000 GPU):

docs/source/reference/geotiff.rst

@@ -0,0 +1,22 @@
+..  _reference.geotiff:
+
+***************
+GeoTIFF / COG
+***************
+
+Reading
+=======
+.. autosummary::
+    :toctree: _autosummary
+
+    xrspatial.geotiff.open_geotiff
+    xrspatial.geotiff.read_vrt
+
+Writing
+=======
+.. autosummary::
+    :toctree: _autosummary
+
+    xrspatial.geotiff.to_geotiff
+    xrspatial.geotiff.write_geotiff_gpu
+    xrspatial.geotiff.write_vrt

docs/source/reference/index.rst

@@ -14,6 +14,7 @@ Reference
    fire
    flood
    focal
+   geotiff
    hydrology
    interpolation
    kde

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
+}