Seamless3DEP: Streamlined Access to USGS 3DEP Topographic Data

Seamless3DEP is a lightweight Python package that simplifies access to topographic data from the USGS 3D Elevation Program (3DEP). Whether you need elevation data or its derivatives, Seamless3DEP provides an efficient interface to both static and dynamic 3DEP products. For global coverage or bathymetry-inclusive workflows, Seamless3DEP also exposes the NOAA global DEM mosaic and a generic helper for any ArcGIS ImageServer/exportImage endpoint.

Seamless3DEP utilizes connection pooling across threads (safely) to optimize service calls and minimize redundant connections. This reduces both the service load and the time required to retrieve data, making it an ideal tool for handling large-scale topographic data requests.

📚 Full documentation is available here.

Available Products

Static DEMs

• 1/3 arc-second (10 meters)
• 1 arc-second (30 meters)
• 2 arc-second (60 meters)

Dynamic Products

• Digital Elevation Model (DEM)
• Hillshade Derivatives:
  • Gray Hillshade
  • Multidirectional Hillshade
  • GreyHillshade with Elevation Fill
  • Hillshade with Elevation Tint
• Terrain Analysis:
  • Aspect (Degrees and Map)
  • Slope (Degrees and Map)
  • Height (Ellipsoidal)
• Contours:
  • Contour 25 (Dynamically generates 25 contours for the area of interest)
  • Contour Smoothed 25 (Smoothed version of the 25 contours)

Core Functions

Seamless3DEP offers eight main functions designed for efficient data retrieval and processing:

• get_dem: Retrieves static DEMs within a specified bounding box. The function automatically splits large areas into manageable tiles, downloads data as GeoTIFF files in EPSG:4326, and supports resolutions of 10m, 30m, or 60m. A buff_npixels option produces overlapping tiles to avoid NoData strips along seams when mosaicking with tiffs_to_da.
• get_map: Fetches any 3DEP product (including DEMs and terrain derivatives) for areas within the US. Works with all available product types, allows custom resolution settings, and downloads in EPSG:3857. Also accepts buff_npixels for overlapping tiles.
• get_global_dem: Fetches the NOAA global DEM mosaic (SRTM + GEBCO + ICESat) for any bounding box on Earth. Returns GeoTIFFs in EPSG:3857. Unlike 3DEP, it includes sub-zero bathymetric values (down to roughly -9982 m), making it the right choice outside the US or for coastal and marine workflows. Native resolution is approximately 30 m (1 arc-second).
• get_image_server: General-purpose downloader for any ArcGIS ImageServer/exportImage endpoint. Accepts any output CRS, an optional Esri rendering rule, and handles URL validation, bbox tiling, and resume-on-failure automatically. get_map and get_global_dem are thin wrappers over this function.
• elevation_bygrid: Samples elevation values from the 10 m seamless DEM at a grid of longitude/latitude coordinates. Reads directly from the USGS Cloud-Optimized GeoTIFFs (EPSG:4269) and supports configurable resampling methods including nearest, bilinear, cubic, cubic spline, and Lanczos.
• decompose_bbox: Handles large area requests by breaking down extensive bounding boxes into optimal sizes based on resolution and maximum pixel count, ensuring efficient data retrieval.
• build_vrt: Creates a virtual raster dataset by combining multiple GeoTIFF files via the gdalbuildvrt CLI. Supports a resampling parameter (default "nearest") and automatically sets NaN as the nodata value for float rasters so mosaic gaps appear as NaN rather than finite sentinels. Requires the libgdal/gdal conda-forge package (which ships gdalbuildvrt); note that libgdal-core provides the library but not the CLI utilities.
• tiffs_to_da: Mosaics a list of GeoTIFF files and returns a clipped xarray.DataArray. By default (use_vrt=False) tiles are merged in memory via rioxarray.merge.merge_arrays — only rioxarray is required, no gdalbuildvrt CLI or shapely. Pass use_vrt=True to mosaic through a VRT file instead, which enables GDAL pixel functions such as "mean" and "median" but requires gdalbuildvrt. The geometry parameter accepts any object with a .bounds attribute (e.g., any shapely geometry) or a plain (west, south, east, north) bounding box tuple.

Important Notes

• Bounding box coordinates should be in decimal degrees (WGS84) format: (west, south, east, north)
• Default projection for requesting maps is EPSG:3857
• get_map is restricted to US coverage; use get_global_dem or get_image_server for data outside the US

Installation

Choose your preferred installation method:

Using pip

pip install seamless-3dep

Using micromamba (recommended)

micromamba install -c conda-forge seamless-3dep

Alternatively, you can use conda or mamba.

Quick Start Guide

We can retrieve topographic data using Seamless3DEP in just a few lines of code. Then, we can visualize or even reproject the data using rioxarray.

Retrieving a DEM

from pathlib import Path
import seamless_3dep as s3dep
import rioxarray as rxr

# Define area of interest (west, south, east, north)
bbox = (-105.7006276, 39.8472777, -104.869054, 40.298293)
data_dir = Path("data")

# Download DEM
tiff_files = s3dep.get_dem(bbox, data_dir)

# Convert to xarray.DataArray
dem = s3dep.tiffs_to_da(tiff_files, bbox, crs=4326)

Retrieving a Slope Map

slope_files = s3dep.get_map("Slope Degrees", bbox, data_dir)
dem = s3dep.tiffs_to_da(slope_files, bbox)

Contributing

We welcome contributions! Please see the contributing section for guidelines and instructions.