Raster transformation
=====================

Raster inputs are handled by
:meth:`~vyperdatum.transformer.Transformer.transform_raster`. Horizontal and
vertical datum transformations are applied, and the output is written to the
local disk. The behavior described here applies to GDAL-supported raster
inputs reached through ``transform`` or by calling ``transform_raster``
directly.

Same horizontal CRS requirement
-------------------------------

As of the current release, raster transformation requires the input and
output horizontal CRSs to match. The pipeline anchors the output to the
input's pixel grid (origin, resolution, and dimensions), which is
well-defined only when both ends share a horizontal CRS. A strict
authority-code comparison is performed at the top of ``transform_raster``;
when the horizontal CRSs differ, ``NotImplementedError`` is raised and the
point-cloud drivers are indicated as the supported alternative for cross-CRS
work. The check is positioned before the internal error-handling block so
that the exception propagates rather than being captured. A development-only
bypass parameter exists but is not forwarded by ``transform`` and is not
intended for production use.

The cutline two-pass pipeline
-----------------------------

For vertical shifts that rely on NOAA partial-coverage grids, a single
combined warp was found to be unreliable, and a two-pass approach is used
instead:

1. The input is clipped to the cutline polygon's bounding envelope,
   producing a cut intermediate that is pixel-aligned to the input.
2. Pass 1 warps the clipped intermediate through the PROJ pipeline and
   materializes the result eagerly to a GeoTIFF, so that no coordinate
   operation is carried lazily into the next stage.
3. Pass 2 expands the result back to the input's extent and applies the
   cutline polygon as a mask, with the coordinate operation suppressed so
   that the step affects extent only.
4. Non-elevation bands are restored from the input.

The cut intermediate, the Pass 1 output, the cutline file, and any tiled
intermediates are removed after Pass 2 completes.

Large inputs: tiled Pass 1
--------------------------

Very large rasters were observed to fail during the setup phase of a single
Pass 1 warp. For inputs whose post-clip pixel count exceeds 250 million,
Pass 1 no longer runs as a single warp. Instead, the cut intermediate is
processed in tiles, each tile is warped independently to a GeoTIFF in a
subdirectory beside the output file, and a VRT mosaic over the tile files is
built and used as the input to Pass 2.

The tile dimensions default to 4096 by 4096 and are configurable through the
``VYPER_PASS1_TILE_SIZE`` environment variable. The 250-million-pixel
routing threshold is a module-level constant in ``transformer.py``. Inputs
below the threshold continue to use the single-warp path unchanged. The tiled
path is reachable only on the same-horizontal-CRS path; the horizontal-CRS
check is applied before the tiling decision is reached.

The tile mosaic is safe with respect to the lazy-evaluation concern noted
above. Each tile GeoTIFF is materialized eagerly during its warp, so the VRT
built over the tiles is a pure routing index and Pass 2 reads computed pixel
values directly, with no transformation re-evaluated.

If any tile's warp fails, the tiles written so far and the subdirectory are
removed and the exception is re-raised, so that no partial output is left in
place. Tiles are processed sequentially in the current release.

Output renderability
--------------------

Two output characteristics are aligned with the input so that downstream
tools render the output efficiently:

- The GeoTIFF block size of the input is propagated to the output, mirroring
  the existing BAG block-size handling.
- Overview pyramids are generated by default. Pyramid levels are derived
  automatically by doubling until the smallest level fits within
  approximately 256 pixels on the longer axis. Average resampling is used by
  default for continuous data, with an override for categorical bands, and a
  horizontal predictor is set automatically for float rasters under DEFLATE
  or LZW compression. Overview generation is controlled by the ``overview``
  parameter of ``transform_raster``, which defaults to ``True``.