Terrain¶

Download DEM tiles, merge them into local rasters, and compute where a sensor’s field of view intersects the ground surface. The module docstring below documents assumptions, limitations, and the public API.

Terrain analysis using Copernicus DEM data.

Downloads, caches, and queries 30-meter Copernicus GLO-30 DEM tiles from AWS. Provides bulk elevation lookup, DEM tile merging via rasterio, and a vectorized ray-terrain intersection algorithm for computing off-nadir ground intersection points.

Assumptions and limitations¶

DEM is treated as static.: Tiles are downloaded once and cached indefinitely. There is no versioning or temporal awareness — the DEM is assumed to represent a fixed surface for the duration of a planning session.
Nearest-neighbor elevation sampling.: get_elevations rounds query coordinates to the nearest pixel center and returns that pixel’s value. There is no bilinear or bicubic interpolation. At 30 m resolution this introduces up to ~15 m of horizontal positioning error, which is negligible for flight planning but would matter for sub-pixel terrain modeling.
No nodata masking.: Pixels with nodata values (e.g. ocean voids) are returned as-is. Callers that operate over coastlines or data gaps should check for nodata sentinels (available via DEMGrid.nodata) themselves.
Out-of-bounds queries clamp to edge pixels.: Query points outside the DEM extent return the elevation of the nearest edge pixel rather than NaN. A warning is logged.
Ray-terrain intersection uses fixed-step marching.: ray_terrain_intersection steps along the slant range at a fixed interval (default 10 m) and detects the first range step where the DEM surface rises above the ray altitude. This is not a root-finding algorithm — the intersection point is located to within one step of the true crossing. Smaller precision values improve accuracy at the cost of a larger search grid.
Surface normals use 3x3 central differences.: surface_normal_at estimates slope from the two adjacent pixels in each direction, converted to metric gradients using a latitude-dependent meters-per-degree approximation. There is no sub-pixel terrain modeling or higher-order surface fitting.
Ellipsoidal altitude only.: All altitudes are heights above the WGS84 ellipsoid. The module does not apply geoid corrections (EGM96/EGM2008). For most flight-planning purposes the ~20 m geoid-ellipsoid separation is small relative to typical AGL flight altitudes (hundreds to thousands of meters).

Terrain references¶

Data source: Copernicus DEM GLO-30, European Space Agency, distributed via AWS Open Data (s3://copernicus-dem-30m).

class DEMGrid[source]¶

Bases: object

Backend-independent DEM raster.

All terrain math consumes DEMGrid instead of live rasterio datasets. The geotransform uses the GDAL convention: (origin_x, pixel_width, x_skew, origin_y, y_skew, pixel_height) where pixel_height is typically negative (north-up).

Variables:

array – 2-D elevation array (rows, cols) in meters.
geotransform – 6-element affine transform tuple.
bounds – (west, south, east, north) in degrees.
nodata – No-data sentinel value, or None.

array: ndarray¶

geotransform: Tuple[float, float, float, float, float, float]¶

bounds: Tuple[float, float, float, float]¶

nodata: float | None = None¶

property raster_min: float¶

property raster_max: float¶

property shape: Tuple[int, int]¶

__init__(array, geotransform, bounds, nodata=None)¶

Parameters:

array (ndarray)
geotransform (Tuple[float, float, float, float, float, float])
bounds (Tuple[float, float, float, float])
nodata (float | None)

Return type:

None

get_cache_root(custom_path=None)[source]¶

Get the root directory for caching files.

Return type:: str
Parameters:: custom_path (str | None)

clear_cache()[source]¶

Clears the entire cache directory after confirming it is safe to do so.

Return type:: None

clear_localdem_cache(confirm=True)[source]¶

Clears the local DEM cache directory.

This removes all files in the ‘localdem’ subdirectory of the cache root, which stores downloaded DEM tiles.

Parameters:: confirm (bool) – If True, prompt the user for confirmation before clearing the cache.
Return type:: None

build_tile_index(tile_list_file)[source]¶

Build an R-tree spatial index for DEM tiles from a tile list file.

Each line in the file is a tile name encoding lat/lon in the filename. Tiles are parsed into 1x1 degree bounding boxes and indexed spatially.

Parameters:

tile_list_file (str) – Path to the text file listing available DEM tiles.

Returns:

A tuple of (rtree_index,: tile_bboxes) where tile_bboxes is a list of (tile_name, bounding_box) pairs.

Return type:

Tuple[index.Index, List[Tuple[str, box]]]

download_dem_files(lon_min, lat_min, lon_max, lat_max, aws_dir)[source]¶

Download DEM tile files covering a geographic bounding box.

Tiles are downloaded from the specified AWS directory and cached locally. Already-downloaded tiles are reused from the cache.

Parameters:

lon_min (float) – Western longitude bound (degrees).
lat_min (float) – Southern latitude bound (degrees).
lon_max (float) – Eastern longitude bound (degrees).
lat_max (float) – Northern latitude bound (degrees).
aws_dir (str) – Base URL of the AWS-hosted DEM tile directory.

Returns:

List of local file paths to the downloaded DEM tiles.

Return type:

List[str]

merge_tiles(output_filename, tile_file_list)[source]¶

Merge multiple DEM tile files into a single GeoTIFF using rasterio.

Parameters:

output_filename (str) – Path for the merged output GeoTIFF file.
tile_file_list (List[str]) – List of file paths to DEM tiles to merge.

Raises:

ValueError – If tile_file_list is empty or contains invalid paths.
RuntimeError – If merge operation fails.

Return type:

None

generate_demfile(latitude, longitude, aws_dir='https://copernicus-dem-30m.s3.amazonaws.com/')[source]¶

Generate a DEM file covering the specified latitude and longitude extents.

Return type:

str

Parameters:

latitude (ndarray)
longitude (ndarray)
aws_dir (str)

load_dem(dem_file)[source]¶

Load a DEM file and return a cached DEMGrid.

Handles mtime lookup so callers don’t have to.

Return type:: DEMGrid
Parameters:: dem_file (str)

get_elevations(lats, lons, dem_file)[source]¶

Extract elevation values for given latitudes and longitudes from a DEM file.

Reads the entire raster band once and indexes it in bulk, rather than querying pixel-by-pixel. The raster is cached per (path, mtime) so repeated calls against the same DEM (the common case in flight planning) pay the read cost only once.

Return type:

ndarray

Parameters:

lats (ndarray)
lons (ndarray)
dem_file (str)

get_elevations_from_grid(lats, lons, dem)[source]¶

Extract elevation values from a DEMGrid (no file I/O).

Return type:

ndarray

Parameters:

lats (ndarray)
lons (ndarray)
dem (DEMGrid)

get_min_max_elevations(dem_file)[source]¶

Get the minimum and maximum elevation values from a DEM file.

Parameters:: dem_file (str) – Path to the DEM file.
Returns:: (min_elevation, max_elevation) in the DEM file.
Return type:: Tuple[float, float]

ray_terrain_intersection(lat0, lon0, h0, az, tilt, precision=10.0, dem_file=None)[source]¶

Batch computation of ray-terrain intersections using a DEM for multiple observer positions. Vectorized to handle multiple observers efficiently.

Parameters:

lat0 (np.ndarray) – Array of observer latitudes (degrees).
lon0 (np.ndarray) – Array of observer longitudes (degrees).
h0 (float) – Altitude of the observer above the ellipsoid (meters).
az (np.ndarray) – Azimuth angle of the ray (degrees).
tilt (np.ndarray) – Off-nadir tilt angle from vertical (degrees, 0-90), where 0 is nadir and 90 is horizontal.
precision (float) – Precision of the slant range sampling (meters).
dem_file (str) – Path to the DEM file. If None, one will be generated.

Returns:

(intersection_lats, intersection_lons,: intersection_alts). Observers with no terrain intersection have NaN in all three arrays.

Return type:

Tuple[np.ndarray, np.ndarray, np.ndarray]

Raises:

ValueError – If tilt or azimuth values are out of range, or if tilt is too close to horizontal (90 deg) where slant-range geometry is undefined.

terrain_elevation_along_track(flight_line, dem_file, precision=100.0)[source]¶

Min, mean, and max terrain elevation (m MSL) along a flight line’s nadir track.

Samples the DEM at evenly-spaced points along the flight line and returns summary statistics useful for Mode 3 per-line altitude planning.

Parameters:

flight_line – A FlightLine object with a track(precision) method.
dem_file (str) – Path to a DEM GeoTIFF covering the flight line.
precision (float) – Along-track sampling interval in meters. Default 100 m.

Return type:

dict

Returns:

Dict with keys "min", "mean", and "max" (all in meters MSL).

terrain_aspect_azimuth(polygon, dem_file=None)[source]¶

Dominant terrain gradient direction (degrees from north) for a polygon.

Computes the dominant downslope azimuth from the DEM gradient over the polygon area. To minimise altitude variation along each flight line (as recommended by Zhao et al. 2021 for Mode 3), orient flight lines perpendicular to the returned azimuth:

flight_azimuth = (terrain_aspect_azimuth(polygon) + 90) % 360

Parameters:

polygon – Shapely Polygon defining the survey area (WGS84 lon/lat).
dem_file (str | None) – Path to a DEM GeoTIFF. If None, one is downloaded and cached from the Copernicus GLO-30 archive.

Return type:

float

Returns:

Dominant downslope azimuth in degrees clockwise from true north (range 0-360).

surface_normal_at(lats, lons, dem_file)[source]¶

Compute outward surface normal unit vectors in the ENU frame.

Uses central-difference gradients from the 3x3 DEM neighbourhood around each query point, converted to metric slope.

Parameters:

lats (ndarray) – Latitude(s) in decimal degrees.
lons (ndarray) – Longitude(s) in decimal degrees.
dem_file (str) – Path to a GeoTIFF DEM file.

Return type:

ndarray

Returns:

(N, 3) array of [east, north, up] unit normal components. On flat terrain the normal is [0, 0, 1].