Module geotiff 

A launch site’s height from a user’s elevation file: a GeoTIFF digital elevation model (DEM) on a geographic (latitude and longitude) grid.

A GeoTIFF is a TIFF image whose pixels are heights and whose tags say where on Earth they lie. The caller reads the whole file into memory and passes its bytes; ElevationRaster::parse reads the tags, and ElevationRaster::height_at finds the pixel a latitude and longitude fall in and decodes only the tile or strip that holds it, so a lookup adds about one tile’s memory to the file’s. The TIFF itself (codecs, predictors, tiles, strips, byte order, BigTIFF) is decoded by image-rs’s tiff crate (docs.rs); this module reads the geographic tags as the OGC GeoTIFF Standard 1.1 (OGC 19-008r4, 2019) defines them, and where GDAL reads a file differently from the standard it either follows GDAL or refuses the file, so a height read here is the one GDAL’s readers give for the same point.

Where a pixel lies. The file ties raster coordinates to longitude and latitude with a tiepoint (I, J) ↦ (X, Y) and a pixel size (S_x, S_y) (ModelTiepointTag, ModelPixelScaleTag; §7.3), or with an affine matrix (ModelTransformationTag; its terms a, b, d and e, f, h give X = a·I + b·J + d, Y = e·I + f·J + h). A positive S_y means latitude falls as rows go down (Requirement 10.4). As GDAL does, both become the longitude and latitude of the outer corner of pixel (0, 0) and a signed pixel size:

λ₀ = X − I·S_x, φ₀ = Y − J·(−S_y), Δλ = S_x, Δφ = −S_y, or from the matrix λ₀ = d, Δλ = a, φ₀ = h, Δφ = f.

The raster type (GTRasterTypeGeoKey, §7.2.1) says what a raster coordinate names. In the usual pixel is area, (0, 0) is the outer corner of the first pixel; in pixel is point it is that pixel’s center, so the corner is half a pixel back: λ₀ −= Δλ/2, φ₀ −= Δφ/2. A point (φ, λ) then lies in column ⌊(λ − λ₀)/Δλ⌋ and row ⌊(φ − φ₀)/Δφ⌋: the pixel whose area holds it, with no interpolation between pixels. On a 1-arc-second grid (about 31 m by 26 m at 33° N) that is the height of the ground within about 20 m of the site. A point within about 10⁻¹³ of a pixel’s width of an edge can fall on either side, here and in GDAL, by rounding.

Refused where the standard and GDAL disagree, or where GDAL needs more than this reads: a negative S_y (GDAL reads it as north-up, the standard as south-up), both a pixel scale and a matrix (the standard forbids it; GDAL takes the scale), a rotated matrix, several tiepoints (ground control points) and an internal nodata mask.

What is read. One band of unsigned or signed 8-, 16- or 32-bit integers, or 32- or 64-bit floats; uncompressed, LZW, Deflate or PackBits, with or without the horizontal or floating-point predictor; tiles or strips; little- or big-endian; classic TIFF or BigTIFF. The first image in the file is the one read (later ones are a cloud-optimized GeoTIFF’s overviews). The CRS must be geographic (GTModelTypeGeoKey 2, or a geodetic CRS key with no model type, as GeoTIFF 1.0 writers leave it; GDAL reads the latter as a local CRS, at the same pixel positions) in degrees from Greenwich, and one of NEAR_WGS84: datums within a few meters of WGS 84, where a point’s WGS 84 latitude and longitude read the right pixel to within a few meters (more near the rupture of a large earthquake since the datum was fixed; the list gives examples). Its EPSG code is reported. Any other, and a projected file (UTM, say), is refused with GeoTiffError::Unsupported naming its code; gdalwarp -t_srs EPSG:4326 in.tif out.tif turns it into one this reads.

Heights. A pixel’s raw value v becomes a height in meters as (v·scale + offset)·unit, with GDAL’s scale and offset:

• S_z from ModelPixelScaleTag and Z₀ − z₀·S_z from the tiepoint’s heights, where GDAL certainly applies them: a GeoTIFF 1.1 directory of model type 2 naming a vertical CRS from VERTICAL_CRS_UNITS, with no VerticalDatumGeoKey, beside a geographic CRS other than WGS 84 3D. With no vertical key, or in a GeoTIFF 1.0 directory (where GDAL drops the vertical CRS), they are ignored, as GDAL ignores them. Between those, whether GDAL applies them turns on how it resolves the keys, and the file is refused unless they give GDAL’s own scale 1 and offset 0 and GDAL_METADATA gives no other scale.
• Otherwise the scale and offset items of GDAL’s GDAL_METADATA tag; otherwise 1 and 0. GDAL matches that tag’s items with quirks, so an item with one of the roles read here is refused if it has a namespace, a capital in an attribute’s name, a sample that isn’t plain digits, a value that isn’t a single text node (CDATA included), or the IMAGE_STRUCTURE domain. GDAL’s own skips (no name, no sample, another band) are skipped.

The unit is the one the file states by VerticalUnitsGeoKey (meters, international feet or US survey feet), by a vertical CRS from VERTICAL_CRS_UNITS, or by the unittype item of GDAL_METADATA; two that disagree are refused, as is a vertical CRS off the list (GDAL takes its unit from EPSG’s registry, whatever the key says). Vertical keys GDAL drops with their unit, or reads by rules of its own, are refused: a private value (above 32767) in any of them (dropped with a model type, read without one), any beside WGS 84 3D, VerticalDatumGeoKey 6030 beside WGS 84 with model type 2 (GDAL makes it WGS 84 3D), and any with no model type and no unit key. In a GeoTIFF 1.0 directory GDAL drops the vertical CRS but keeps its unit; this module reports both; with no model type GDAL’s local CRS ignores VerticalGeoKey, which this module reports too. A unit name is read after trimming ASCII blanks; GDAL drops leading blanks typed as they are and keeps the rest ("ft " is feet here). A file that states no unit is read as meters, flagged by RasterInfo::vertical_unit_stated; GDAL reports no unit there, except beside a vertical datum key alone, where it assumes meters too; a file in feet that states none reads 3.28 times too high. The vertical datum (VerticalGeoKey, NAVD88 or EGM2008, say) is reported, not applied. A value equal to the file’s nodata value (GDAL’s GDAL_NODATA tag) or a NaN reads as no height; a nodata value the sample type can’t hold exactly, such as 12.5 on integers, matches nothing.

Bounds on a hostile file. A tile or strip larger than MAX_CHUNK_BYTES decoded is refused at ElevationRaster::parse, and ElevationRaster::values grows the raster fallibly as rows of tiles decode. GDAL_METADATA nested more than 16 deep is refused before it is parsed. The tiff crate prints one debug line to standard error when a tag’s value passes its 1 MiB limit (its own dbg!).

Guide: A launch site’s elevation walks through an example and says how the reader is checked: against GDAL’s reading, through rasterio, of seven files and a whole USGS tile. The choices are in ADR-128, a site’s height from a user’s GeoTIFF.

use hpr_io::geotiff::ElevationRaster;

// A real program reads its file: `let bytes = std::fs::read(path)?;`.
let bytes = include_bytes!("../../tests/fixtures/geotiff/usgs-f32-lzw-fp-tiles.tif");
let raster = ElevationRaster::parse(bytes)?;
// Spaceport America's runway: 1,400.691 m in the USGS's terrain model.
let height_m = raster.height_at(32.99, -106.97)?;
assert_eq!(height_m.map(|h| (h * 1000.0).round() / 1000.0), Some(1400.691));

Structs

Bounds

A raster’s edges, degrees.

ElevationRaster

A GeoTIFF elevation file, its tags read and its pixels left packed in the borrowed bytes.

Pixel

A pixel’s place in the raster, counted from 0 at the first row and column.

RasterInfo

What an elevation file says about itself.

Enums

GeoTiffError

Why a GeoTIFF elevation file could not be read, or a height not taken from it.

RasterType

What a raster coordinate names (GTRasterTypeGeoKey, OGC 19-008r4 §7.2.1).

SampleType

The sample type of the raster’s pixels.

VerticalUnit

The unit a pixel’s value is in, converted to meters by VerticalUnit::meters.

Constants

MAX_CHUNK_BYTES

The largest tile or strip, decoded, a file may have: 256 MiB, the tiff crate’s own limit on a decoded chunk, which its padding of a floating-point tile would otherwise bypass.

MAX_VALUES_PIXELS

The most pixels ElevationRaster::values decodes into one vector: 2²⁸, 2 GiB of f64.

NEAR_WGS84

The geographic CRSs read, by EPSG code: datums whose latitude and longitude lie within a few meters of WGS 84’s, so a point given in WGS 84 reads the right pixel, or its neighbour on a grid finer than a few meters. They part by plate motion since each was fixed, and by earthquakes: near the rupture of a large one since a datum was fixed, such as Chile’s in 2010 for SIRGAS 2000 or Wenchuan’s in 2008 for CGCS2000, the ground moved several meters. JGD2000 is left out: Japan’s 2011 earthquake moved its north-east by more than 5 m, and JGD2011 replaced it.

VERTICAL_CRS_UNITS

The vertical CRSs read, by EPSG code, with their unit from EPSG’s registry: (code, unit). All are gravity-related heights, positive up. A file naming another is refused, as GDAL would take its unit from the registry, which this reader doesn’t hold.