gri-terrain

Terrain elevation data access for geolocation applications.

Status: Alpha -- The three data sources (DTED, GeoTIFF/COG, Copernicus DEM), tile caching, interpolation, ray-terrain intersection, and spline surface normals are implemented and tested. A higher-level multi-source Terrain facade with automatic NaN-triggered fallback is planned but not yet available; query each source directly for now.

Overview

gri-terrain provides access to terrain elevation data from multiple sources with caching and interpolation. It is part of the GRI FOSS (GeoSol Research Free and Open Source Software) ecosystem. Requires Python 3.12+.

Features

• Multiple data sources: DTED, GeoTIFF/COG, Copernicus DEM
• Vectorized lookups: Efficient batch elevation queries
• Interpolation options: Nearest neighbor, bilinear, bicubic
• Tile caching: In-memory LRU caching, plus an on-disk cache for downloaded Copernicus tiles
• Pre-caching: Load a region's tiles up front for fast repeat lookups
• Ray-terrain intersection: Adaptive-step ray tracing against a source's terrain
• gri-utils interop: Sources yield ECEF sheets for gri-utils terrain math (interpolation, stitching, visibility/shadow, observable intersection, spline surface normals)

Installation

pip install gri-terrain

Or for development:

git clone https://gitlab.com/geosol-foss/python/gri-terrain.git
cd gri-terrain
uv sync

Quick Start

Query a source directly. Every source exposes the same vectorized get_altitude(lat, lon, *, interpolation="bicubic", datum="hae") method.

import numpy as np
from gri_terrain import DTEDSource

# Point a source at a directory of DTED tiles (level auto-detected)
source = DTEDSource("/path/to/dted")

# Vectorized lookup (lat/lon in WGS84 degrees, elevation in meters)
lats = np.array([40.0, 41.0, 42.0])
lons = np.array([-105.0, -106.0, -107.0])
altitudes = source.get_altitude(lats, lons)

# A single point works too (returns a length-1 array)
altitude = source.get_altitude(np.array([40.0]), np.array([-105.0]))[0]

Locations outside a source's coverage return NaN rather than raising.

Vertical Datum

Elevations are returned as height above the WGS84 ellipsoid by default, not height above sea level. The two differ by the geoid undulation, up to about -107 m to +85 m depending on location, so the distinction is never cosmetic.

# Height above the WGS84 ellipsoid: what lla_to_xyz and every other
# coordinate transform expects. This is the default.
alt = source.get_altitude(lats, lons)

# Height above mean sea level: what the source file holds, and what a map
# or a GPS receiver's "elevation" shows.
alt = source.get_altitude(lats, lons, datum="source")

Every common terrain product publishes orthometric height, measured from a geoid rather than from the ellipsoid: DTED and SRTM are EGM96, Copernicus DEM is EGM2008. Handing those values straight to a coordinate transform displaces every point radially by the undulation, and for a shallow ray that vertical error becomes a horizontal one an order of magnitude larger.

Sources determine their own datum: DTEDSource reads the vertical datum from each file's DSI record, GeoTiffSource reads a compound CRS when the file has one and otherwise assumes EGM96, and CopernicusSource declares EGM2008. Pass vertical_datum= to any of them to override a file whose metadata is wrong:

# DTED resampled from Copernicus still reports "MSL" in its header
source = DTEDSource("/path/to/dted", level=2, vertical_datum="egm2008")

The conversion uses the EGM96 grid shipped with gri-utils. See the gri-utils geoid README for the model and its accuracy.

Data Sources

Copernicus DEM (Default)

High-quality global elevation data from ESA, freely available on AWS:

• GLO-30: 30m resolution (default)
• GLO-90: 90m resolution

import numpy as np
from gri_terrain import CopernicusSource

source = CopernicusSource(resolution="30m")  # or "90m"
alt = source.get_altitude(np.array([40.0]), np.array([-105.0]))

DTED

Digital Terrain Elevation Data (military standard):

• Level 0: ~1km resolution
• Level 1: ~100m resolution
• Level 2: ~30m resolution

from gri_terrain.sources import DTEDSource

# Load a specific DTED level
source = DTEDSource("/path/to/dted", level=1)

# Load best available resolution per tile (requires best/ symlinks)
source = DTEDSource("/path/to/dted", level="best")

terrain = Terrain(sources=[source])

GeoTIFF

Local elevation data in GeoTIFF format:

import numpy as np
from gri_terrain import GeoTiffSource

# Accepts a single file or a directory of tiles
source = GeoTiffSource("/path/to/elevation.tif")
alt = source.get_altitude(np.array([40.0]), np.array([-105.0]))

Interpolation

Three interpolation methods are supported:

# Nearest neighbor (fastest)
alt = source.get_altitude(lats, lons, interpolation="nearest")

# Bilinear (good balance)
alt = source.get_altitude(lats, lons, interpolation="bilinear")

# Bicubic (smoothest, default)
alt = source.get_altitude(lats, lons, interpolation="bicubic")

Pre-caching

Load all tiles covering a bounding box up front, for fast repeat lookups in a known operating area. Returns the number of tiles loaded and raises MemoryError if the region exceeds the in-memory cache budget.

source.precache_region(
    lat_min=39.0, lat_max=41.0,
    lon_min=-106.0, lon_max=-104.0,
)

For CopernicusSource this also downloads the covering tiles to the on-disk cache, making the region available offline afterward.

Dependencies

• numpy
• scipy
• rasterio (GeoTIFF support)
• dted (DTED file parsing)
• gri-utils (coordinate conversions)

Ray-Terrain Intersection

Find where a ray intersects the terrain surface:

import numpy as np
from gri_terrain import DTEDSource, ray_terrain
from gri_utils.conversion import lla_to_xyz, xyz_to_lla

source = DTEDSource("/path/to/dted")

# Observer at 45N, 0E, 10km altitude looking down
origin_lla = np.array([45.0, 0.0, 10000.0])
origin_xyz = lla_to_xyz(origin_lla)

# Direction toward Earth center (descending)
direction_xyz = -origin_xyz / np.linalg.norm(origin_xyz)

# Find intersection
hit_xyz = ray_terrain(source, origin_xyz, direction_xyz)

if hit_xyz is not None:
    hit_lla = xyz_to_lla(hit_xyz)
    print(f"Hit at {hit_lla[0]:.4f}N, {hit_lla[1]:.4f}E, {hit_lla[2]:.1f}m")

The algorithm uses adaptive step sizes based on:

1. Altitude band skip: Fast-forward to the terrain altitude band (-800m to 9200m above the ellipsoid, wide enough to cover the geoid undulation on top of real terrain)
2. Tile skip: When above a tile's maximum elevation, skip to tile boundary
3. Slope-based skip: Use 45-degree max terrain slope assumption for safe step sizes

The altitude_epsilon_m parameter offsets the terrain surface (useful for vegetation canopy or safety margins):

# Find where ray passes within 10m of terrain
hit = ray_terrain(terrain, origin, direction, altitude_epsilon_m=10.0)

Specular Reflection

Where a signal from tx_xyz bounces off terrain to reach rx_xyz:

from gri_terrain import DTEDSource, specular_point

source = DTEDSource("/path/to/dted")
hit = specular_point(source, tx_xyz, rx_xyz, wavelength_m=0.1903)

hit.point_xyz        # ECEF bounce point
hit.delay_s          # excess delay over the direct path
hit.delay_var_s2     # from the tile's own declared vertical accuracy
hit.d_delay_d_tx     # exact sensitivity to transmitter motion
hit.point_cov        # first Fresnel zone, as a covariance

The delay sensitivities are exact and cost nothing extra: because the path length is stationary in the bounce point, its derivative with respect to an endpoint has no term in how that point moves. delay_var_s2 follows the same way from the surface height uncertainty, which DTEDSource reads from each tile's accuracy record rather than assuming.

Ocean and gaps in coverage fall back to the geoid, so a bounce off water lands on the sea surface rather than on the ellipsoid tens of metres below it.

For a moving emitter or collector, pass the previous epoch's point as guess_lla. The solve then starts beside its own answer, and branch_changed tells you when the bounce has jumped to a different facet -- across which the delay stays continuous but its gradient does not.

Surface Normals from Gridded Elevation

Spline-smoothed unit surface normals (ECEF) from a regular (lat, lon) elevation grid are pure array math and live in gri-utils, not here:

import numpy as np
from gri_utils.terrain import grid_normals_spline

lats = np.linspace(39.5, 40.5, 121)   # deg
lons = np.linspace(-105.5, -104.5, 121)
alt_grid = load_elevation(lats, lons)  # shape (121, 121), meters

normals_ecef = grid_normals_spline(lats, lons, alt_grid)
# shape (121, 121, 3), unit vectors in ECEF

To compute normals from a gri-terrain source, get a sheet first (source.get_covering_sheets(...)), recover its altitude grid, or feed the sheet to the other gri_utils.terrain routines (visibility, shadow, intersection).

Development Status

• Source abstraction (TerrainSource ABC)
• DTED source elevation lookup
• GeoTIFF/COG source elevation lookup
• Copernicus DEM source (on-demand AWS download + disk cache)
• Vectorized lookups with nearest/bilinear/bicubic interpolation
• Tile caching (in-memory LRU; Copernicus adds an on-disk cache)
• Pre-caching a region
• Ray-terrain intersection (over a source; algorithm in gri-utils)
• Sheets for gri-utils terrain math (get_covering_sheets / covers)
• Higher-level multi-source Terrain facade with automatic NaN-triggered fallback

Attribution

When using Copernicus DEM data, include:

(c) DLR e.V. 2010-2014 and (c) Airbus Defence and Space GmbH 2014-2018 provided under COPERNICUS by the European Union and ESA; all rights reserved

Other Projects

Current list of other GRI FOSS Projects we are building and maintaining.

License

MIT License. See LICENSE for details.