Metadata-Version: 2.4
Name: gri-terrain
Version: 0.2.4
Summary: Terrain elevation data access for geolocation applications
Project-URL: Homepage, https://geosolresearch.com
Project-URL: Repository, https://gitlab.com/geosol-foss/python/gri-terrain
Project-URL: Issues, https://gitlab.com/geosol-foss/python/gri-terrain/-/issues
Project-URL: Changelog, https://gitlab.com/geosol-foss/python/gri-terrain/-/releases
Author-email: GeoSol Research Inc <contact@geosolresearch.com>
License-Expression: MIT
License-File: LICENSE
Keywords: DEM,DTED,GeoTIFF,elevation,terrain
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering
Classifier: Topic :: Scientific/Engineering :: GIS
Requires-Python: >=3.12
Requires-Dist: dted>=1.0.3
Requires-Dist: gri-utils>=0.3.4
Requires-Dist: numpy>=2.3.3
Requires-Dist: rasterio>=1.4.3
Requires-Dist: scipy>=1.16.3
Requires-Dist: tqdm>=4.67.0
Description-Content-Type: text/markdown

[![GeoSol Research Logo](https://geosolresearch.com/logos/foss_logo.png "GeoSol Research")](https://geosolresearch.com)

# 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

```bash
pip install gri-terrain
```

Or for development:

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

## Quick Start

```python
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

```python
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

```python
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:

```python
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:

```python
# 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:

```python
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:

```python
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):

```python
# 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).

```python
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

- [x] Source abstraction (TerrainSource ABC)
- [x] Terrain class with fallback chain
- [x] Ray-terrain intersection
- [x] 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](https://gitlab.com/geosol-foss/python/gri-terrain/-/blob/main/.docs_other_projects.md) we are building and maintaining.

## License

MIT License. See [LICENSE](https://gitlab.com/geosol-foss/python/gri-terrain/-/blob/main/LICENSE) for details.
