Elevation query using https://aws.amazon.com/public-datasets/terrain/

elevation_query.py

#!/usr/bin/env python

from fiona.transform import transform_geom
from rasterstats import point_query
import mercantile


def make_dem_url(lon: float, lat: float, z: int = 14) -> str:
    """Returns a URL referencing the GeoTiff Digitial Elevation Model (DEM)
    for the given point at a zoom level (default, max is 14).

    See https://aws.amazon.com/public-datasets/terrain/
    """
    if z > 14:
        raise ValueError("max zoom is 14")

    # Pad the point with a small epsilon to create a valid bounding box
    eps = 1e-4
    bbox = (lon - eps, lat - eps, lon + eps, lat + eps)

    # Find the parent tile at the specified zoom level
    x, y, _ = mercantile.parent(mercantile.bounding_tile(*bbox), zoom=z)

    return f"https://elevation-tiles-prod.s3.amazonaws.com/geotiff/{z}/{x}/{y}.tif"


def make_mercator_geom(lon: float, lat: float) -> dict:
    """Convert to GeoJSON-like geometry dict in Web Mercator projection
    """
    return transform_geom(
        "EPSG:4326", "EPSG:3857", {"type": "Point", "coordinates": [lon, lat]}
    )


def elevation_at_point(lon: float, lat: float) -> float:
    """Returns the elevation in meters for the given lon, lat

    In the future, we could download the geotiff locally
    and cache it to the local file system. Could improve
    performance vs going over HTTP for each point
    """
    raster_url = make_dem_url(lon, lat)
    geom = make_mercator_geom(lon, lat)
    elev = point_query(geom, raster_url)
    return elev[0]


def test_elevation_query():
    denver = (-104.9848, 39.7392)  # CO State Capitol
    elev = elevation_at_point(*denver)
    print(f"{elev} meters")
    assert abs(elev - 1609.3) < 2  # Within 2 meters of "mile high"


if __name__ == "__main__":
    test_elevation_query()