gri-terrain 0.2.4

Terrain elevation data access for geolocation applications

gri-terrain

Terrain elevation data access for geolocation applications.

Status: Early Development -- Class interfaces are defined but source-level elevation lookups are not yet implemented. Ray-terrain intersection is functional.

Overview

gri-terrain provides unified access to terrain elevation data from multiple sources with automatic fallback, 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
• Automatic fallback: Query sources in priority order
• Vectorized lookups: Efficient batch elevation queries
• Interpolation options: Nearest neighbor, bilinear, bicubic
• Tile caching: Memory and disk caching for performance
• Pre-caching: Download regions for offline use
• Ray-terrain intersection: Adaptive-step ray tracing against gridded elevation
• Surface normals: Bicubic-spline normals in ECEF from gridded elevation, geodetically correct across wide latitude ranges

Installation

pip install gri-terrain

Or for development:

git clone https://gitlab.com/geosol-foss/python/gri-terrain.git
cd gri-terrain
. .init_venv.sh

Quick Start

from gri_terrain import Terrain

# Use default sources (Copernicus DEM with DTED-0 fallback)
terrain = Terrain()

# Get elevation at a single point (not yet implemented)
altitude = terrain.get_altitude(lat=40.0, lon=-105.0)

# Vectorized lookup (not yet implemented)
import numpy as np
lats = np.array([40.0, 41.0, 42.0])
lons = np.array([-105.0, -106.0, -107.0])
altitudes = terrain.get_altitude(lats, lons)

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

from gri_terrain.sources import CopernicusSource

source = CopernicusSource(resolution="30m")  # not yet implemented
terrain = Terrain(sources=[source])

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:

from gri_terrain.sources import GeoTiffSource

source = GeoTiffSource("/path/to/elevation.tif")  # not yet implemented
terrain = Terrain(sources=[source])

Interpolation

Three interpolation methods are supported:

# Nearest neighbor (fastest)
alt = terrain.get_altitude(lat, lon, interpolation="nearest")

# Bilinear (default, good balance)
alt = terrain.get_altitude(lat, lon, interpolation="bilinear")

# Bicubic (smoothest)
alt = terrain.get_altitude(lat, lon, interpolation="bicubic")

Pre-caching

Download tiles for offline use:

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

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:

from gri_terrain.intersect import ray_terrain
from gri_utils.conversion import lla_to_xyz
import numpy as np

# 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(terrain, 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

Compute unit surface normals in ECEF from a regular (lat, lon) elevation grid. Fits a bicubic spline to the elevation array and evaluates analytic partial derivatives, which removes the step-noise that finite-difference gradients produce on integer-meter DEM quantization. Slopes are converted from degrees to ENU meters using the ellipsoid's meridional and prime-vertical radii of curvature at each grid point, so results are correct across wide latitude ranges (no flat-Earth approximation).

import numpy as np
from gri_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

Pass a non-WGS84 ellipsoid via the ellipsoid= keyword (accepts any gri_utils.constants.GeoEllipsoid).

Development Status

• Source abstraction (TerrainSource ABC)
• Terrain class with fallback chain
• Ray-terrain intersection
• Surface normals from gridded elevation (spline-based)
• DTED source elevation lookup
• GeoTIFF/COG source elevation lookup
• Copernicus DEM source (AWS access)
• Tile caching (memory and disk)
• Pre-caching for offline use

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.