USGS 3DEP Elevation Data

Overview

The USGS 3D Elevation Program (3DEP) provides high-resolution terrain elevation data for the United States (CONUS, Alaska, Hawaii, and US territories). The default resolution is 1/3 arc-second (~10m).

Data is downloaded on-demand from the USGS National Map ImageServer REST API as GeoTIFF tiles and read with TiffImages.jl (no GDAL dependency required). This is a static (time-invariant) dataset – the elevation values do not change over the simulation time span.

EarthSciData.USGS3DEP — Function

USGS3DEP(domaininfo; name, resolution, stream)

Create a ModelingToolkit System that provides terrain elevation and slope data from the USGS 3D Elevation Program (3DEP).

The system exposes three variables interpolated to the simulation coordinates:

elevation (m): terrain elevation above sea level
dzdx (dimensionless): terrain slope in the x (east) direction (rise/run)
dzdy (dimensionless): terrain slope in the y (north) direction (rise/run)

Slopes are computed from the elevation grid using central finite differences with coordinate conversion from lon/lat to metric distances.

The domain may use any coordinate reference system (lon-lat, UTM, Lambert Conformal Conic, etc.). The domain bounding box is automatically reprojected to WGS84 for the USGS API request, and coordinate transforms between the domain CRS and the data's native lon-lat grid are handled by the interpolation infrastructure.

Note that 3DEP coverage is limited to the United States (CONUS, Alaska, Hawaii, and US territories).

Arguments

domaininfo: A DomainInfo specifying the spatial and temporal domain.
name: System name (default :USGS3DEP).
resolution: Target resolution in arc-seconds (default 1/3 ≈ 10m).
stream: Whether to stream data lazily (default true).

Example

using EarthSciData, EarthSciMLBase, ModelingToolkit, Dates
domain = DomainInfo(
    DateTime(2018, 11, 8), DateTime(2018, 11, 9);
    lonrange = deg2rad(-121.7):deg2rad(0.01):deg2rad(-121.5),
    latrange = deg2rad(39.7):deg2rad(0.01):deg2rad(39.8),
    levrange = 1:1,
)
elev = USGS3DEP(domain)

source

Usage

using EarthSciData, EarthSciMLBase
using ModelingToolkit, DataFrames
using ModelingToolkit: t
using Dates
using DynamicQuantities
using DynamicQuantities: dimension

# Define a domain near Paradise, CA
domain = DomainInfo(
    DateTime(2018, 11, 8), DateTime(2018, 11, 9);
    lonrange = deg2rad(-121.7):deg2rad(0.01):deg2rad(-121.5),
    latrange = deg2rad(39.7):deg2rad(0.01):deg2rad(39.8),
    levrange = 1:1,
)

elev = USGS3DEP(domain; resolution=10.0)

\[ \begin{align} \mathtt{elevation}\left( t \right) &= \mathtt{elevation\_itp}\left( t + \mathtt{t\_ref}, \mathtt{lon}, \mathtt{lat} \right) \\ \mathtt{dzdx}\left( t \right) &= \mathtt{dzdx\_itp}\left( t + \mathtt{t\_ref}, \mathtt{lon}, \mathtt{lat} \right) \\ \mathtt{dzdy}\left( t \right) &= \mathtt{dzdy\_itp}\left( t + \mathtt{t\_ref}, \mathtt{lon}, \mathtt{lat} \right) \end{align} \]

Variables

vars = unknowns(elev)
DataFrame(
    :Name => [string(Symbolics.tosymbol(v, escape = false)) for v in vars],
    :Units => [dimension(ModelingToolkit.get_unit(v)) for v in vars],
    :Description => [ModelingToolkit.getdescription(v) for v in vars]
)

3×3 DataFrame

Row	Name	Units	Description
	String	Dimensio…	String
1	elevation	m	Terrain elevation above sea level
2	dzdx		Terrain slope in x (east) direction (dimensionless, rise/run)
3	dzdy		Terrain slope in y (north) direction (dimensionless, rise/run)

Parameters

function table(vars)
    DataFrame(
        :Name => [string(Symbolics.tosymbol(v, escape = false)) for v in vars],
        :Units => [dimension(ModelingToolkit.get_unit(v)) for v in vars],
        :Description => [ModelingToolkit.getdescription(v) for v in vars]
    )
end
table(parameters(elev))

4×3 DataFrame

Row	Name	Units	Description
	String	Dimensio…	String
1	t_ref	s	Reference time
2	elevation_itp	m	Interpolated elevation
3	dzdx_itp		Interpolated dzdx
4	dzdy_itp		Interpolated dzdy

Equations

eqs = equations(elev)

How It Works

Bounding box: The spatial extent of the domain is converted to a WGS84 (EPSG:4326) longitude/latitude bounding box with a one-pixel buffer.
Download: A single GeoTIFF tile covering the bounding box is requested from the USGS ImageServer exportImage endpoint. Pixel dimensions are computed from the requested resolution, capped at 4000x4000 to limit download size. The tile is cached locally.
Coordinate mapping: Pixel-centre coordinates are computed analytically from the bounding box and pixel dimensions (no GeoTIFF metadata parsing needed). Coordinates are stored in radians for consistency with the EarthSciML coordinate system.
Interpolation: The elevation field is provided as a DataSetInterpolator that maps simulation coordinates to elevation values via bilinear interpolation.

Coverage

US coverage only

3DEP provides elevation data for the United States only. If the domain falls outside US coverage, a warning is issued and the returned data may contain nodata values. Approximate coverage bounds: longitude -180 to -60, latitude 17 to 72.