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

* Add design spec for geotiff performance and memory controls

Covers dtype on read, compression_level on write, and VRT tiled
output from to_geotiff for streaming dask writes.

* Add implementation plan for geotiff performance controls (#1083)

* Add compression_level parameter to to_geotiff and write

Thread a compression_level: int | None = None parameter from the
public to_geotiff() API down to the compress() call in the writer.
Validation rejects out-of-range levels for deflate (1-9), zstd (1-22),
and lz4 (0-16); codecs without a level concept (lzw, packbits, jpeg)
accept any value and ignore it.

* Add compression_level to write_geotiff_gpu and forward from to_geotiff

write_geotiff_gpu now accepts compression_level: int | None = None so
callers that pass the parameter through to_geotiff on the GPU path no
longer get an unexpected keyword argument error. The value is accepted
silently since nvCOMP does not currently expose level control.

* Clean up compression_level plumbing: keyword args, docstring, dead import

- Use keyword form compression_level=compression_level at the two
  _write_tiled call sites in write() for clarity and refactor safety.
- Update compress() docstring to cover all codecs with level support
  (deflate, zstd, lz4) instead of only mentioning deflate.
- Remove unused tempfile import from test_compression_level.py.

* Add dtype parameter to open_geotiff for on-read dtype casting

Adds `dtype=None` to `open_geotiff`, `read_geotiff_dask`, `read_geotiff_gpu`,
and `read_vrt`. When specified, each path casts the array to the target dtype
after nodata masking. Float-to-int casts raise ValueError to prevent
accidental fractional data loss. Resolves part of #1083.

* Fix dtype cast skipped on GPU non-tiled fallback early return

The stripped-file fallback in read_geotiff_gpu returned before
reaching the dtype cast block. Add the cast before the early return
so the dtype parameter is honored for non-tiled GeoTIFFs on the
GPU path.

* Fix dask dtype validation against post-nodata-masking dtype

The dask path validated the dtype cast against the raw file dtype
(e.g. uint16), but nodata masking promotes integer arrays to float64
inside the delayed worker. Requesting dtype='int32' on a uint16 file
with nodata would pass validation then silently cast float64+NaN to
int32, producing garbage. Now compute the effective post-masking dtype
and validate against that.

Also remove the dead `dtype` positional parameter from
_delayed_read_window and add tests for the nodata+dtype interaction
on both eager and dask paths.

* Add VRT tiled output to to_geotiff (#1083)

When to_geotiff receives a path ending in .vrt, write a directory of
tiled GeoTIFFs with a VRT index file instead of a monolithic TIFF.
This lets dask arrays stream to disk without materializing the full
array in RAM.

- Dask inputs get one tile per chunk, written in parallel via
  dask.delayed
- Numpy inputs are sliced by tile_size
- cog=True and overview_levels raise ValueError with VRT output
- Non-empty tiles directory raises FileExistsError
- 11 new tests cover numpy/dask round-trips, tile naming, relative
  paths, compression forwarding, and edge cases

* Fix OOM risk, 3D guard, and unused param in _write_vrt_tiled

- Use scheduler='synchronous' in dask.compute so tiles write one at a
  time, keeping peak memory to a single chunk instead of N_cores chunks
- Reject 3D arrays in both the dask and numpy VRT paths with a clear
  error message, since raw.chunks[0] would read the band dimension
  instead of rows for (bands, y, x) shaped arrays
- Remove unused gpu parameter from _write_vrt_tiled signature and
  call site; CuPy detection already happens inside _write_single_tile

* Add dtype, compression_level, and VRT output examples to README GeoTIFF I/O section

* Add user guide notebook for GeoTIFF performance controls

Demonstrates dtype on read, compression_level on write, and VRT tiled
output using self-contained tempdir examples.

* Fix VRT relative paths on Windows: use forward slashes in XML (#1083)

os.path.relpath produces backslashes on Windows, but VRT XML
expects forward slashes for cross-platform portability.

Author: brendancol

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/plans/2026-03-30-geotiff-perf-controls.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.

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

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