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,
    spatial_interp
)

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).
spatial_interp = :linear (default) does full multilinear interpolation; :nearest does spatial nearest-neighbour + time-only linear interpolation for ~8x speedup when queries are always at grid points.

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\_unit} ~ interp\_unsafe\left( \mathtt{elevation\_data}, 1 + \frac{ - \mathtt{elevation\_tstart} + t + \mathtt{t\_ref}}{\mathtt{elevation\_tstep}}, 1 + \frac{ - \mathtt{elevation\_s1start} + \mathtt{lon}}{\mathtt{elevation\_s1step}}, 1 + \frac{ - \mathtt{elevation\_s2start} + \mathtt{lat}}{\mathtt{elevation\_s2step}}, \mathtt{elevation\_extrap} \right) \\ \mathtt{dzdx}\left( t \right) &= \mathtt{dzdx\_unit} ~ interp\_unsafe\left( \mathtt{dzdx\_data}, 1 + \frac{ - \mathtt{dzdx\_tstart} + t + \mathtt{t\_ref}}{\mathtt{dzdx\_tstep}}, 1 + \frac{ - \mathtt{dzdx\_s1start} + \mathtt{lon}}{\mathtt{dzdx\_s1step}}, 1 + \frac{ - \mathtt{dzdx\_s2start} + \mathtt{lat}}{\mathtt{dzdx\_s2step}}, \mathtt{dzdx\_extrap} \right) \\ \mathtt{dzdy}\left( t \right) &= \mathtt{dzdy\_unit} ~ interp\_unsafe\left( \mathtt{dzdy\_data}, 1 + \frac{ - \mathtt{dzdy\_tstart} + t + \mathtt{t\_ref}}{\mathtt{dzdy\_tstep}}, 1 + \frac{ - \mathtt{dzdy\_s1start} + \mathtt{lon}}{\mathtt{dzdy\_s1step}}, 1 + \frac{ - \mathtt{dzdy\_s2start} + \mathtt{lat}}{\mathtt{dzdy\_s2step}}, \mathtt{dzdy\_extrap} \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))

28×3 DataFrame

Row	Name	Units	Description
	String	Dimensio…	String
1	t_ref	s	Reference time
2	elevation_s1start		Spatial grid start dim 1 for elevation
3	elevation_s1step		Spatial grid step dim 1 for elevation
4	elevation_s2start		Spatial grid start dim 2 for elevation
5	elevation_s2step		Spatial grid step dim 2 for elevation
6	elevation_extrap		Extrapolation type for elevation
7	elevation_unit	m	Unit scale for elevation
8	dzdx_s1start		Spatial grid start dim 1 for dzdx
9	dzdx_s1step		Spatial grid step dim 1 for dzdx
10	dzdx_s2start		Spatial grid start dim 2 for dzdx
11	dzdx_s2step		Spatial grid step dim 2 for dzdx
12	dzdx_extrap		Extrapolation type for dzdx
13	dzdx_unit		Unit scale for dzdx
14	dzdy_s1start		Spatial grid start dim 1 for dzdy
15	dzdy_s1step		Spatial grid step dim 1 for dzdy
16	dzdy_s2start		Spatial grid start dim 2 for dzdy
17	dzdy_s2step		Spatial grid step dim 2 for dzdy
18	dzdy_extrap		Extrapolation type for dzdy
19	dzdy_unit		Unit scale for dzdy
20	elevation_data		Interpolation data for elevation
21	elevation_tstart	s	Time grid start for elevation
22	elevation_tstep	s	Time grid step for elevation
23	dzdx_data		Interpolation data for dzdx
24	dzdx_tstart	s	Time grid start for dzdx
25	dzdx_tstep	s	Time grid step for dzdx
26	dzdy_data		Interpolation data for dzdy
27	dzdy_tstart	s	Time grid start for dzdy
28	dzdy_tstep	s	Time grid step for 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 1000x1000 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.