Raster transformation

Raster inputs are handled by 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.