Add dtype, compression_level, and VRT output to geotiff I/O (#1083)

README.md

@@ -165,6 +165,12 @@ to_geotiff(data, 'out.tif', gpu=True)                # force GPU compress
 to_geotiff(data, 'ortho.tif', compression='jpeg')    # JPEG for orthophotos
 write_vrt('mosaic.vrt', ['tile1.tif', 'tile2.tif'])  # generate VRT
 
+open_geotiff('dem.tif', dtype='float32')             # half memory
+open_geotiff('dem.tif', dtype='float32', chunks=512) # Dask + half memory
+to_geotiff(data, 'out.tif', compression_level=1)     # fast scratch write
+to_geotiff(data, 'out.tif', compression_level=22)    # max compression
+to_geotiff(dask_da, 'mosaic.vrt')                    # stream Dask to VRT
+
 # Accessor methods
 da.xrs.to_geotiff('out.tif', compression='lzw')     # write from DataArray
 ds.xrs.open_geotiff('large_dem.tif')                 # read windowed to Dataset extent

docs/superpowers/specs/2026-03-30-geotiff-perf-controls-design.md

@@ -0,0 +1,147 @@
+# GeoTIFF Performance and Memory Controls
+
+Adds three parameters to `open_geotiff` and `to_geotiff` that let callers
+control memory usage, compression speed, and large-raster write strategy.
+All three are opt-in; default behaviour is unchanged.
+
+## 1. `dtype` parameter on `open_geotiff`
+
+### API
+
+```python
+open_geotiff(source, *, dtype=None, ...)
+```
+
+`dtype` accepts any numpy dtype string or object (`np.float32`, `'float32'`,
+etc.). `None` preserves the file's native dtype (current behaviour).
+
+### Read paths
+
+| Path | Behaviour |
+|------|-----------|
+| Eager (numpy) | Output array allocated at target dtype. Each decoded tile/strip cast before copy-in. Peak overhead: one tile at native dtype. |
+| Dask | Each delayed chunk function casts after decode. Output chunks are target dtype. Same per-tile overhead. |
+| GPU (CuPy) | Cast on device after decode. |
+| Dask + CuPy | Combination of dask and GPU paths. |
+
+### Numba LZW fast path
+
+The LZW decoder is a numba JIT function that emits values one at a time into a
+byte buffer. A variant will decode each value and cast inline to the target
+dtype so the per-tile buffer is never allocated at native dtype. Other codecs
+(deflate, zstd) return byte buffers from C libraries where per-value
+interception isn't possible, so those fall back to the tile-level cast.
+
+### Validation
+
+- Narrowing float casts (float64 to float32): allowed.
+- Narrowing int casts (int64 to int16): allowed (user asked for it explicitly).
+- Widening casts (float32 to float64, uint8 to int32): allowed.
+- Float to int: `ValueError` (lossy in a way users often don't intend).
+- Unsupported casts (e.g. complex128 to uint8): `ValueError`.
+
+## 2. `compression_level` parameter on `to_geotiff`
+
+### API
+
+```python
+to_geotiff(data, path, *, compression='zstd', compression_level=None, ...)
+```
+
+`compression_level` is `int | None`. `None` uses the codec's existing default.
+
+### Ranges
+
+| Codec | Range | Default | Direction |
+|-------|-------|---------|-----------|
+| deflate | 1 -- 9 | 6 | 1 = fastest, 9 = smallest |
+| zstd | 1 -- 22 | 3 | 1 = fastest, 22 = smallest |
+| lz4 | 0 -- 16 | 0 | 0 = fastest |
+| lzw | n/a | n/a | No level support; ignored silently |
+| jpeg | n/a | n/a | Quality is a separate axis; ignored |
+| packbits | n/a | n/a | Ignored |
+| none | n/a | n/a | Ignored |
+
+### Plumbing
+
+`to_geotiff` passes `compression_level` to `write()`, which passes it to
+`compress()`. The internal `compress()` already accepts a `level` argument; we
+just thread it through the two intermediate call sites that currently hardcode
+it.
+
+### Validation
+
+- Out-of-range level for a codec that supports levels: `ValueError`.
+- Level set for a codec without level support: silently ignored.
+
+### GPU path
+
+`write_geotiff_gpu` also accepts and forwards the level to nvCOMP batch
+compression, which supports levels for zstd and deflate.
+
+## 3. VRT output from `to_geotiff` via `.vrt` extension
+
+### Trigger
+
+When `path` ends in `.vrt`, `to_geotiff` writes a tiled VRT instead of a
+monolithic TIFF. No new parameter needed.
+
+### Output layout
+
+```
+output.vrt
+output_tiles/
+  tile_0000_0000.tif   # row_col, zero-padded
+  tile_0000_0001.tif
+  ...
+```
+
+Directory name derived from the VRT stem (`foo.vrt` -> `foo_tiles/`).
+Zero-padding width scales to the grid dimensions.
+
+### Behaviour per input type
+
+| Input | Tiling strategy | Memory profile |
+|-------|----------------|----------------|
+| Dask DataArray | One tile per dask chunk. Each task computes its chunk and writes one `.tif`. | One chunk in RAM at a time (scheduler controlled). |
+| Dask + CuPy | Same, GPU compress per tile. | One chunk in GPU memory at a time. |
+| Numpy / ndarray | Slice into `tile_size`-sized pieces, write each. | Source array already in RAM; tile slices are views (no duplication). |
+| CuPy | Same as numpy but GPU compress. | Source on GPU; tiles are views. |
+
+### Per-tile properties
+
+- Same `compression`, `compression_level`, `predictor`, `nodata`, `crs` as the
+  parent call.
+- `tiled=True` with the caller's `tile_size` (internal TIFF tiling within each
+  chunk-file).
+- GeoTransform adjusted to each tile's spatial position (row/col offset from
+  the full raster origin).
+- No COG overviews on individual tiles.
+
+### VRT generation
+
+After all tiles are written, call `write_vrt()` with relative paths. The VRT
+XML references each tile by its spatial extent and band mapping.
+
+### Edge cases and validation
+
+- `cog=True` with a `.vrt` path: `ValueError` (mutually exclusive).
+- Tiles directory exists and is non-empty: `FileExistsError` to prevent silent
+  overwrites.
+- Tiles directory doesn't exist: created automatically.
+- `overview_levels` with `.vrt` path: `ValueError` (overviews don't apply).
+
+### Dask scheduling
+
+For dask inputs, all delayed tile-write tasks are submitted to
+`dask.compute()` at once. The scheduler manages parallelism and memory. Each
+task is: compute chunk, compress, write tile file. No coordination between
+tasks.
+
+## Out of scope
+
+- Streaming write of a monolithic `.tif` from dask input (tracked as a separate
+  issue). Users who need a single file from a large dask array can write to VRT
+  and convert externally, or ensure sufficient RAM.
+- JPEG quality parameter (separate concern from compression level).
+- Automatic chunk-size recommendation based on available memory.

examples/user_guide/46_GeoTIFF_Performance.ipynb

@@ -0,0 +1,188 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "# GeoTIFF Performance Controls\n",
+        "\n",
+        "Three opt-in controls for memory, compression speed, and large-raster writes:\n",
+        "\n",
+        "- **`dtype` on read** — cast to a narrower type at load time to reduce memory use\n",
+        "- **`compression_level` on write** — trade write speed for file size (or vice versa)\n",
+        "- **VRT tiled output** — stream a chunked dask array to a directory of tiles without loading it all at once"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "import numpy as np\n",
+        "import xarray as xr\n",
+        "import os\n",
+        "import tempfile\n",
+        "from xrspatial.geotiff import open_geotiff, to_geotiff"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "# Create a 1000x1000 float64 DEM-like raster (~8MB)\n",
+        "rng = np.random.default_rng(42)\n",
+        "elevation = rng.normal(loc=500, scale=100, size=(1000, 1000)).astype(np.float64)\n",
+        "y = np.linspace(40.0, 41.0, 1000)\n",
+        "x = np.linspace(-106.0, -105.0, 1000)\n",
+        "dem = xr.DataArray(elevation, dims=['y', 'x'],\n",
+        "                   coords={'y': y, 'x': x},\n",
+        "                   attrs={'crs': 4326})\n",
+        "print(f\"DEM shape: {dem.shape}, dtype: {dem.dtype}, size: {dem.nbytes / 1e6:.1f} MB\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## 1. dtype on read\n",
+        "\n",
+        "Pass `dtype` to `open_geotiff` to cast the raster to a narrower type at load time.\n",
+        "Reading a float64 file as float32 halves memory use without any extra copy.\n",
+        "The cast happens inside rasterio before the array reaches Python, so it works on\n",
+        "all read paths: eager, dask, and GPU."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "with tempfile.TemporaryDirectory() as tmpdir:\n",
+        "    path = os.path.join(tmpdir, 'dem_f64.tif')\n",
+        "    to_geotiff(dem, path)\n",
+        "\n",
+        "    # Read at native dtype (float64)\n",
+        "    native = open_geotiff(path)\n",
+        "    print(f\"Native: dtype={native.dtype}, size={native.nbytes / 1e6:.1f} MB\")\n",
+        "\n",
+        "    # Read as float32 -- half the memory\n",
+        "    downcast = open_geotiff(path, dtype='float32')\n",
+        "    print(f\"Downcast: dtype={downcast.dtype}, size={downcast.nbytes / 1e6:.1f} MB\")\n",
+        "\n",
+        "    # Works with dask too\n",
+        "    dask_f32 = open_geotiff(path, dtype='float32', chunks=256)\n",
+        "    print(f\"Dask chunks dtype: {dask_f32.dtype}\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## 2. compression_level\n",
+        "\n",
+        "`to_geotiff` accepts a `compression_level` argument alongside `compression`.\n",
+        "For zstd the range is 1–22 (1 = fastest, 22 = smallest file).\n",
+        "For deflate the range is 1–9.\n",
+        "The default is the codec's own default when `compression_level` is omitted.\n",
+        "\n",
+        "Use a low level when write speed matters (streaming pipelines, scratch files).\n",
+        "Use a high level for archival or network transfer where file size dominates."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "import time\n",
+        "\n",
+        "with tempfile.TemporaryDirectory() as tmpdir:\n",
+        "    results = []\n",
+        "    for level in [1, 3, 10, 22]:\n",
+        "        path = os.path.join(tmpdir, f'dem_zstd_l{level}.tif')\n",
+        "        t0 = time.perf_counter()\n",
+        "        to_geotiff(dem, path, compression='zstd', compression_level=level)\n",
+        "        elapsed = time.perf_counter() - t0\n",
+        "        size_kb = os.path.getsize(path) / 1024\n",
+        "        results.append((level, elapsed, size_kb))\n",
+        "        print(f\"  level={level:2d}  time={elapsed:.3f}s  size={size_kb:.0f} KB\")\n",
+        "\n",
+        "    print(f\"\\nLevel 1 vs 22: {results[0][2]/results[-1][2]:.1f}x size difference\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## 3. VRT tiled output\n",
+        "\n",
+        "Pass a `.vrt` path to `to_geotiff` and it writes a directory of GeoTIFF tiles\n",
+        "plus a VRT index file that GDAL treats as a single dataset.\n",
+        "\n",
+        "Each tile corresponds to one dask chunk and is written independently, so only\n",
+        "one chunk is in memory at a time. This makes it practical to write arrays that\n",
+        "are larger than RAM.\n",
+        "\n",
+        "The VRT uses relative paths, so the whole output directory is portable."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "with tempfile.TemporaryDirectory() as tmpdir:\n",
+        "    # Chunk the DEM for dask processing\n",
+        "    dask_dem = dem.chunk({'y': 500, 'x': 500})\n",
+        "\n",
+        "    vrt_path = os.path.join(tmpdir, 'tiled_dem.vrt')\n",
+        "    to_geotiff(dask_dem, vrt_path, compression='zstd')\n",
+        "\n",
+        "    # Show what was created\n",
+        "    tiles_dir = os.path.join(tmpdir, 'tiled_dem_tiles')\n",
+        "    print(\"Files created:\")\n",
+        "    print(f\"  {os.path.basename(vrt_path)}\")\n",
+        "    for f in sorted(os.listdir(tiles_dir)):\n",
+        "        size = os.path.getsize(os.path.join(tiles_dir, f)) / 1024\n",
+        "        print(f\"  tiled_dem_tiles/{f}  ({size:.0f} KB)\")\n",
+        "\n",
+        "    # Read it back via VRT\n",
+        "    result = open_geotiff(vrt_path)\n",
+        "    print(f\"\\nRound-trip: shape={result.shape}, dtype={result.dtype}\")\n",
+        "    print(f\"Max difference: {float(np.abs(result.values - dem.values).max()):.2e}\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Summary\n",
+        "\n",
+        "| Feature | Parameter | Where to use |\n",
+        "|---|---|---|\n",
+        "| dtype cast | `open_geotiff(..., dtype='float32')` | Reduce read memory by half |\n",
+        "| compression level | `to_geotiff(..., compression_level=1)` | Fast scratch writes; set high for archival |\n",
+        "| VRT tiled output | `to_geotiff(..., 'out.vrt')` | Stream large dask arrays to disk without OOM |"
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "name": "python",
+      "version": "3.10.0"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 5
+}