gri-terrain 0.2.6

Terrain elevation data access for geolocation applications

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") 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.

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 (-500m to 9000m)
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)

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.