# VRTs


# What a VRT is

A [GDAL virtual raster](https://gdal.org/en/stable/drivers/raster/vrt.html) (VRT) is a small XML file that describes a raster built from other rasters. It holds no pixels of its own. Instead, it lists the source files and says how to combine them: how to mosaic them, which bands to use and how to reproject them. GDAL and rasterio open a VRT like any other raster and read from the sources only when you ask for pixels.

This VRT mosaics two tiles side by side (shortened):

``` xml
<VRTDataset rasterXSize="2000" rasterYSize="1000">
  <SRS>PROJCS["WGS 84 / UTM zone 32N", ...]</SRS>
  <GeoTransform>600000.0, 10.0, 0.0, 5200000.0, 0.0, -10.0</GeoTransform>
  <VRTRasterBand dataType="Float32" band="1">
    <ComplexSource>
      <SourceFilename relativeToVRT="1">a.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <DstRect xOff="0" yOff="0" xSize="1000" ySize="1000" />
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">b.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <DstRect xOff="1000" yOff="0" xSize="1000" ySize="1000" />
    </ComplexSource>
  </VRTRasterBand>
</VRTDataset>
```

`DstRect` places each tile in the output, in pixels: `a.tif` fills the left half and `b.tif` the right half.


# When a VRT helps

A VRT is useful when you need a combined or reprojected raster but don't want to write a copy of it:

- Mosaicking tiles into one raster.
- Reprojecting a raster onto a target grid, such as an [MGRS square](mgrs-grid.md).
- Stacking bands from separate files, e.g. an RGB raster from three single-band files.
- Converting values, such as backscatter from linear power to decibels.

Every read goes back to the sources and recalculates the output. If you read the same pixels often, or chain many steps, write the result to a COG once and read that.


# What freezebase adds

rasterio opens VRTs, but the only one it can build is a `WarpedVRT`, which exists only in memory. The others take GDAL's command-line tools, such as `gdalbuildvrt` and `gdalwarp -of VRT`, or the GDAL Python bindings. freezebase builds them in Python:

- It needs only rasterio, not the GDAL Python bindings.
- A VRT is a Python object, so you can inspect and edit its bands, nodata values and sources before writing it.
- Steps chain in memory, so you can combine the workflows above, for example mosaicking and then reprojecting, without writing files in between.
- Sources can be on S3, with the same credentials as the rest of freezebase.
- [write_cog](../reference/raster.write_cog.md#freezebase.raster.write_cog) turns any VRT into a COG in one call.


# A first example

Mosaic two tiles and save the result as a COG:

``` python
from freezebase import vrt
from freezebase.raster import write_cog

mosaic = vrt.mosaic(["a.tif", "b.tif"])
write_cog(mosaic, "mosaic.tif")
```

[vrt.mosaic](../reference/vrt.mosaic.md#freezebase.vrt.mosaic) reads only the headers of the two files. [write_cog](../reference/raster.write_cog.md#freezebase.raster.write_cog) then reads the pixels through GDAL and writes the COG without loading the whole raster into memory.


# Common tasks

Each function below takes file paths or VRTs and returns a new VRT. Nothing is written until you call [write](../reference/vrt.VRT.md#freezebase.vrt.VRT.write) or [write_cog](../reference/raster.write_cog.md#freezebase.raster.write_cog).


## Mosaic tiles

[vrt.mosaic](../reference/vrt.mosaic.md#freezebase.vrt.mosaic) places tiles on one grid. Where tiles overlap, the later tile in the list wins.

``` python
mosaic = vrt.mosaic(["a.tif", "b.tif"], band_names=["VV"])
```

All tiles need the same CRS, band count, data types and nodata values. Tiles with different pixel sizes need `resolution`: `"finest"` takes the smallest pixel size, and a tuple `(x, y)` sets it. `bounds` fixes the extent, which otherwise covers every tile.


## Reproject onto a grid

[vrt.warp](../reference/vrt.warp.md#freezebase.vrt.warp) reprojects a raster. Give it a target CRS and pixel size, and optionally an extent. This puts a scene onto an MGRS square:

``` python
from freezebase.mgrs import MGRSGrid

square = MGRSGrid(aoi)[0]
on_square = vrt.warp("scene.tif", crs=square.crs, resolution=10.0, bounds=square.boundingbox)
```

Without `resolution`, [warp](../reference/vrt.warp.md#freezebase.vrt.warp) picks a pixel size close to the source's. `like="other.tif"` copies the CRS, extent and pixel size of another raster instead.

[warp](../reference/vrt.warp.md#freezebase.vrt.warp) keeps the source's data type, so averaging `uint16` values rounds the result. Pass `dtype="float32"` to keep the fractions.


## Stack bands from several files

[vrt.stack](../reference/vrt.stack.md#freezebase.vrt.stack) puts the bands of several rasters on the same grid into one raster, in order:

``` python
rgb = vrt.stack(["red.tif", "green.tif", "blue.tif"], band_names=["red", "green", "blue"])
```

`bands` picks which bands to take from each file, for example `bands=[[1], [2, 3]]`. Files on different grids need [warp](../reference/vrt.warp.md#freezebase.vrt.warp) first.


## Convert values

[vrt.derive](../reference/vrt.derive.md#freezebase.vrt.derive) converts values with one of GDAL's built-in pixel functions. It can also combine several inputs, for example into a ratio:

``` python
db = vrt.derive("sigma0.tif", vrt.PixelFunction.db())
ratio = vrt.derive(["vv.tif", "vh.tif"], vrt.PixelFunction.div())
```

[PixelFunction](../reference/vrt.PixelFunction.md#freezebase.vrt.PixelFunction) covers decibels, ratios, differences, the mean, the minimum and the maximum, among others. With one input file, [derive](../reference/vrt.derive.md#freezebase.vrt.derive) converts each of its bands. With several, it computes one band from band 1 of each. The result is `float32` with NaN as nodata unless you pass `dtype` and `nodata`, and a nodata pixel in any input gives a nodata pixel in the output.


## Chain steps

Each function accepts a VRT wherever it accepts a file, so steps chain without writing files:

``` python
web = vrt.warp(vrt.mosaic(["a.tif", "b.tif"]), crs="EPSG:3857", resolution=20.0)
db = vrt.derive(web, vrt.PixelFunction.db())

with db.open() as ds:  # nothing is written to disk
    data = ds.read(1)
```

The final VRT contains the earlier steps, and reading it needs GDAL 3.11 or newer. To keep a step as its own file, write it and pass the returned path on:

``` python
zone = vrt.mosaic(["a.tif", "b.tif"]).write("zone.vrt")
web = vrt.warp(zone, crs="EPSG:3857", resolution=20.0).write("web.vrt")
```

Write the steps to files when other programs need to open them, or when the chain is long.


## Inspect and edit a VRT

[vrt.read](../reference/vrt.read.md#freezebase.vrt.read) loads a VRT file, including one written by `gdalbuildvrt`, and [vrt.from_raster](../reference/vrt.from_raster.md#freezebase.vrt.from_raster) wraps a single raster file. Change what you need, then write it:

``` python
v = vrt.from_raster("scene.tif")
v.bands[0].description = "VH"
v.bands[0].nodata = 0.0
v.write("scene.vrt")
```

A [VRT](../reference/vrt.VRT.md#freezebase.vrt.VRT) has a CRS, a transform, a size and a list of bands. Each [Band](../reference/vrt.Band.md#freezebase.vrt.Band) lists the sources it reads from, and [VRT.sources](../reference/vrt.VRT.md#freezebase.vrt.VRT.sources) lists them for the whole VRT.


# Saving the result

You can save a VRT as a `.vrt` file, which stays small and keeps reading from the sources, or as a raster file that holds the pixels. [VRT.write](../reference/vrt.VRT.md#freezebase.vrt.VRT.write) writes the `.vrt` file and returns its path. [write_cog](../reference/raster.write_cog.md#freezebase.raster.write_cog) writes a Cloud Optimized GeoTIFF, with an optional profile:

``` python
write_cog(mosaic, "mosaic.tif", {"compress": "zstd"})
```

[rewrite_tiff](../reference/raster.rewrite_tiff.md#freezebase.raster.rewrite_tiff) also takes a VRT, for a format other than COG. It ignores `dtype` and `nodata` in the profile; use [write_cog](../reference/raster.write_cog.md#freezebase.raster.write_cog) to change them.

A raster file has one data type and one nodata value for all its bands. If the VRT's bands differ, pass `dtype` or `nodata` in the profile, or [write_cog](../reference/raster.write_cog.md#freezebase.raster.write_cog) raises `ValueError`.

- `nodata` replaces each band's nodata value, and pixels that held the old value get the new one. `nodata: None` keeps the VRT's values.
- `dtype` converts every band. Fractions are rounded and values outside the new type's range are clipped, so a valid value can end up equal to the nodata value: 300 written as `uint8` with nodata 255 becomes 255 and reads as nodata. [write_cog](../reference/raster.write_cog.md#freezebase.raster.write_cog) raises `ValueError` if a band's nodata value does not fit the new type.

``` python
dem = vrt.mosaic(["a.tif", "b.tif"])  # float32, nodata NaN
write_cog(dem, "dem.tif", {"dtype": "int16", "nodata": -9999})
```


# VRTs on S3

An S3 source is a `UPath` built with [make_s3_upath()](../reference/s3.make_s3_upath.md#freezebase.s3.make_s3_upath), and keeps that path's credentials and endpoint. [Accessing S3](accessing-s3.md) describes the options.


## How source paths are saved

In Python, every source path is absolute. When [VRT.write](../reference/vrt.VRT.md#freezebase.vrt.VRT.write) saves a VRT, it stores a source path relative to the VRT if both are local or both are in the same S3 bucket. Other sources are stored as absolute GDAL paths such as `/vsis3/bucket/tile.tif` or `/vsicurl/https://...`. `relative=False` stores every path as absolute.

[vrt.read](../reference/vrt.read.md#freezebase.vrt.read) reverses this: relative paths are resolved against the VRT's location, and `/vsis3/` paths become `s3://` paths. These get the storage options (credentials and endpoint) of the VRT if it is on S3, or the ones you pass:

``` python
from freezebase import vrt
from freezebase.s3 import make_s3_upath

root = make_s3_upath("s3://bucket/run", profile="research")
local = vrt.read("warped.vrt", storage_options=root.storage_options)
tiles = local.source_files(storage_options=root.storage_options)
```

[source_files()](../reference/vrt.VRT.md#freezebase.vrt.VRT.source_files) lists every file the VRT reads in the end, including those behind nested VRTs and behind `.vrt` files it points to. [VRT.sources](../reference/vrt.VRT.md#freezebase.vrt.VRT.sources) lists only the direct sources and reads nothing. Following a `.vrt` file over HTTPS needs aiohttp, which the `s3` extra installs.


## Credentials

[VRT.open()](../reference/vrt.VRT.md#freezebase.vrt.VRT.open) opens a VRT with the credentials its sources need, whether they are local, on S3 or public over HTTP(S):

``` python
with mosaic.open() as ds:
    data = ds.read(1)
```

To open a written `.vrt` file with rasterio yourself, wrap it in [VRT.env()](../reference/vrt.VRT.md#freezebase.vrt.VRT.env):

``` python
import rasterio

path = mosaic.write("mosaic.vrt")
with mosaic.env(), rasterio.open(path) as ds:
    data = ds.read(1)
```

Both raise `ValueError` if the sources need different credentials. If the VRT points to other `.vrt` files with sources on S3, pass those sources' `storage_options` to [open()](../reference/vrt.VRT.md#freezebase.vrt.VRT.open) or [env()](../reference/vrt.VRT.md#freezebase.vrt.VRT.env).


# Converting heights between vertical datums

A height is measured from a reference surface, its vertical datum. Many DEMs give heights above a geoid; the Copernicus DEM, for example, uses EGM2008. GNSS receivers give heights above the ellipsoid. The two differ by tens of metres.

[warp](../reference/vrt.warp.md#freezebase.vrt.warp) converts heights when both the source CRS and the target `crs` have a vertical datum (a compound CRS such as `EPSG:4326+3855`) or are 3D (such as `EPSG:4979`). There is no extra option:

``` python
ellipsoidal = vrt.warp("copdem.tif", crs="EPSG:4979", src_crs="EPSG:4326+3855")
```

`src_crs` adds the vertical datum when the file does not record it.


## Getting the geoid grid

PROJ needs each geoid model as a grid file, and neither rasterio nor pyproj ships any. Either set `PROJ_NETWORK=ON` before the process starts, so PROJ downloads the parts it needs, or download the grid once into a folder that also holds rasterio's PROJ data:

``` python
import shutil

from rasterio.env import PROJDataFinder

shutil.copytree(PROJDataFinder().search(), "proj_data", dirs_exist_ok=True)
```

``` bash
python -m pyproj sync --file us_nga_egm08_25 --target-directory proj_data  # 80 MB
```

Then point rasterio at the folder:

``` python
import rasterio

with rasterio.Env(PROJ_DATA="proj_data"):
    ellipsoidal = vrt.warp("copdem.tif", crs="EPSG:4979", src_crs="EPSG:4326+3855")
```

PROJ remembers for the rest of the process whether it found a grid. Set `PROJ_DATA` or `PROJ_NETWORK` before the first reprojection; setting it after a failed lookup has no effect.


## When the conversion fails

- If PROJ cannot find the grid, [warp](../reference/vrt.warp.md#freezebase.vrt.warp) raises `ValueError` instead of ignoring the geoid. `allow_ballpark=True` accepts PROJ's fallback, which leaves the heights unchanged.
- Either both CRSs need a vertical datum or a third axis, or neither does; a mix raises `ValueError`. To keep the heights as they are, pass a 2D `src_crs` or a 2D `crs`.
- Heights in other units, such as feet, are converted as `gdalwarp` converts them.
- Only single-band rasters without a scale or offset are converted. Warp the height band on its own, convert the stored values to heights first, or pass a 2D `src_crs` to skip the conversion.