Skip to content

API Reference

Module for getting DEM from USGS's 3D Elevation Program (3DEP).

GeometryLike #

Bases: Protocol

A geometry exposing a shapely-style (minx, miny, maxx, maxy) envelope.

Get3DEPErrors #

Get3DEPErrors(errors, vrt_url)

Bases: Exception

Raised when one or more tiles fail to download from 3DEP.

Carries the per-tile failures in :attr:errors so callers can inspect each failure and decide whether to retry. Already-downloaded tiles are left on disk, so re-running :func:get_map will resume from the failures.

Attributes:

  • errors (list[BaseException]) –

    The exceptions raised by individual tile downloads.

  • vrt_url (str) –

    The source VRT URL the failed tiles were being read from.

decompose_bbox #

decompose_bbox(bbox, res, pixel_max, buff_npixels=0)

Divide a Bbox into equal-area sub-bboxes based on pixel count.

Parameters:

  • bbox (tuple) –

    Bounding box coordinates in decimal degrees like so: (west, south, east, north).

  • res (int) –

    Resolution of the domain in meters.

  • pixel_max (int) –

    Maximum number of pixels allowed in each sub-bbox. If None, the bbox is not decomposed.

  • buff_npixels (int, default: 0 ) –

    Number of pixels to buffer each sub-bbox by, defaults to 0. Only applied when the bbox is decomposed into multiple sub-bboxes; a bbox that fits within pixel_max is returned as a single un-buffered box (the buffer exists to give per-tile workflows a halo at internal seams, which is moot when there is only one tile).

Returns:

  • boxes ( list of tuple ) –

    List of sub-bboxes in the form (west, south, east, north).

  • sub_width ( int ) –

    Pixel width of each sub-bbox, including buffer pixels on both sides when buff_npixels > 0.

  • sub_height ( int ) –

    Pixel height of each sub-bbox, including buffer pixels on both sides when buff_npixels > 0.

get_dem #

get_dem(
    bbox,
    save_dir,
    res=10,
    pixel_max=MAX_PIXELS,
    buff_npixels=0,
    max_workers=None,
)

Get DEM from 3DEP at 10, 30, or 60 meters resolutions.

Notes

If you need a different resolution, use the :func:get_map function with map_type="DEM".

Parameters:

  • bbox (tuple) –

    Bounding box coordinates in decimal degrees: (west, south, east, north).

  • save_dir (str or Path) –

    Path to save the GeoTiff files.

  • res ((10, 30, 60), default: 10 ) –

    Target resolution of the DEM in meters, by default 10. Must be one of 10, 30, or 60.

  • pixel_max (int, default: MAX_PIXELS ) –

    Maximum number of pixels allowed in each sub-bbox for decomposing the bbox into equal-area sub-bboxes, defaults to 8 million. If None, the bbox is not decomposed and is downloaded as a single file. Values more than 8 million are not allowed.

  • buff_npixels (int, default: 0 ) –

    Number of pixels to buffer each sub-bbox by when the bbox is decomposed into multiple tiles, defaults to 0 (no buffer).

  • max_workers (int, default: None ) –

    Maximum number of concurrent tile downloads. Defaults to None which uses min(8, os.cpu_count(), n_tiles). The bottleneck is HTTPS round-trip latency to the source COG-VRT, so values up to ~16 are useful for very large requests; setting this to 1 forces serial downloading.

Returns:

  • list of pathlib.Path –

    list of GeoTiff files containing the DEM clipped to the bounding box.

Raises:

  • Get3DEPErrors –

    If one or more tiles fail to download. Already-completed tiles are kept on disk so re-running :func:get_map resumes from the failures.

get_image_server #

get_image_server(
    url,
    bbox,
    save_dir,
    name,
    image_crs,
    res,
    *,
    pixel_max=MAX_PIXELS,
    buff_npixels=0,
    rendering_rule=None,
)

Download tiled GeoTiffs from any ArcGIS ImageServer/exportImage endpoint.

Parameters:

  • url (str) –

    Full URL of an ArcGIS ImageServer exportImage endpoint, e.g. https://server/arcgis/rest/services/Layer/ImageServer/exportImage.

  • bbox (tuple) –

    Bounding box in decimal degrees (WGS84): (west, south, east, north).

  • save_dir (str or Path) –

    Directory to save the GeoTiff files.

  • name (str) –

    Prefix used when naming the output files.

  • image_crs (int, str, or rasterio.crs.CRS) –

    Output image CRS passed as imageSR. Accepts anything rasterio.crs.CRS.from_user_input understands (EPSG integer, authority string such as "EPSG:3857", or a CRS object). Resolved to an EPSG WKID when one exists, otherwise serialized as a WKT JSON string. Not all services support every CRS; check the service capabilities first.

  • res (int) –

    Target resolution in meters.

  • pixel_max (int, default: MAX_PIXELS ) –

    Maximum pixels per tile for bbox decomposition, by default 8 million. Pass None to download as a single file. Values above 8 million are not allowed.

  • buff_npixels (int, default: 0 ) –

    Pixel buffer around each tile when the bbox is decomposed, by default 0.

  • rendering_rule (dict, default: None ) –

    Esri rendering rule serialized as the renderingRule JSON parameter, e.g. {"rasterFunction": "Hillshade Gray"}. None returns raw pixel values (elevation/DEM).

Returns:

  • list of pathlib.Path –

    GeoTiff files covering the bounding box.

get_map #

get_map(
    map_type,
    bbox,
    save_dir,
    res=10,
    pixel_max=MAX_PIXELS,
    buff_npixels=0,
)

Get topo maps in 3857 coordinate system within US from 3DEP at any resolution.

Parameters:

  • map_type (MapTypes) –

    Type of map to get. Must be one of the following:

    • 'DEM'
    • 'Hillshade Gray'
    • 'Aspect Degrees'
    • 'Aspect Map'
    • 'GreyHillshade_elevationFill'
    • 'Hillshade Multidirectional'
    • 'Slope Map'
    • 'Slope Degrees'
    • 'Hillshade Elevation Tinted'
    • 'Height Ellipsoidal'
    • 'Contour 25'
    • 'Contour Smoothed 25'
  • bbox (tuple) –

    Bounding box coordinates in decimal degrees (WGS84): (west, south, east, north).

  • save_dir (str or Path) –

    Path to save the GeoTiff files.

  • res (int, default: 10 ) –

    Target resolution of the map in meters, by default 10.

  • pixel_max (int, default: MAX_PIXELS ) –

    Maximum number of pixels allowed in each sub-bbox for decomposing the bbox into equal-area sub-bboxes, defaults to 8 million. If None, the bbox is not decomposed and is downloaded as a single file. Values more than 8 million are not allowed.

  • buff_npixels (int, default: 0 ) –

    Number of pixels to buffer each sub-bbox by when the bbox is decomposed into multiple tiles, defaults to 0 (no buffer).

Returns:

  • list of pathlib.Path –

    list of GeoTiff files containing the DEM clipped to the bounding box.

get_global_dem #

get_global_dem(
    bbox,
    save_dir,
    *,
    res=30,
    pixel_max=MAX_PIXELS,
    buff_npixels=0,
)

Get the NOAA global DEM mosaic in 3857 coordinate system.

The NOAA global mosaic blends SRTM, GEBCO, ICESat, and other sources into a single seamless worldwide layer. Unlike :func:get_map, it covers the entire globe and includes sub-zero bathymetric values (down to roughly -9982 m). Native resolution is approximately 1 arc-second (~30 m) over most of the world; finer values will upsample from that source.

Parameters:

  • bbox (tuple) –

    Bounding box in decimal degrees (WGS84): (west, south, east, north).

  • save_dir (str or Path) –

    Path to save the GeoTiff files.

  • res (int, default: 30 ) –

    Target resolution in meters, by default 30 to match the global mosaic's native ~1 arc-second resolution.

  • pixel_max (int, default: MAX_PIXELS ) –

    Maximum pixels per tile for bbox decomposition, by default 8 million. Pass None to download as a single file. Values above 8 million are not allowed.

  • buff_npixels (int, default: 0 ) –

    Pixel buffer around each tile when the bbox is decomposed, by default 0.

Returns:

  • list of pathlib.Path –

    GeoTiff files covering the bounding box.

build_vrt #

build_vrt(
    vrt_path,
    tiff_files,
    *,
    pixel_function=None,
    resampling="nearest",
)

Create a VRT from a list of GeoTIFF tiles.

Notes

This function requires the installation of libgdal-core. The recommended approach is to use conda (or alternatives like mamba or micromamba). However, if using the system's package manager is the only option, ensure that the gdal-bin or gdal package is installed. For detailed instructions, refer to the GDAL documentation here. When seamless-3dep is installed from Conda, libgdal-core is installed as a dependency and this function works without any additional steps.

Parameters:

  • vrt_path (str or Path) –

    Path to save the output VRT file.

  • tiff_files (list of str or Path) –

    List of file paths to include in the VRT.

  • pixel_function (str, default: None ) –

    Name of a GDAL VRT pixel function for resolving pixel values where input tiles overlap, defaults to None. When None, the gdalbuildvrt default ("last input wins") is used. Useful when mosaicking per-tile post-processed outputs (e.g., TWI computed independently on buffered DEM tiles) where overlap pixels can differ between tiles. Common values: "first", "mean", "median", "min", "max", "mode", "geometric_mean", "harmonic_mean". Requires GDAL 3.12+.

  • resampling (str, default: 'nearest' ) –

    Resampling method to use when building the VRT, by default "nearest". "nearest" avoids propagating NaN from tile-boundary fringe pixels into the mosaic interior when source tiles have slightly different pixel grids (fractional DstRect offsets). Use "bilinear" or higher only when you need smooth blending across deliberately different resolutions. Must be one of "nearest", "bilinear", "cubic", "cubicspline", "lanczos", "average", or "mode". See the GDAL documentation for details: https://gdal.org/programs/gdalbuildvrt.html

tiffs_to_da #

tiffs_to_da(
    tiff_files,
    geometry,
    crs=4326,
    *,
    pixel_function=None,
    resampling="nearest",
    use_vrt=False,
)

Mosaic a list of tiff files and return a clipped xarray.DataArray.

Parameters:

  • tiff_files (list of Path) –

    List of file paths to convert to a DataArray.

  • geometry (geometry or sequence of float) –

    A geometry exposing a shapely-style bounds attribute (e.g. a shapely geometry), or a bounding box (west, south, east, north). A geometry is additionally clipped to its exact shape; a bounding box clips to the envelope only.

  • crs (int, str, or CRS, default: 4326 ) –

    Coordinate reference system of the input geometry, by default 4326.

  • pixel_function (str, default: None ) –

    Function for resolving overlap between tiles, by default None ("last input wins"). With use_vrt=True this is forwarded to :func:build_vrt and accepts any GDAL pixel function (e.g. "mean", requires GDAL 3.12+). With the default in-memory path (use_vrt=False) only None and "first"/"min"/"max" are supported.

  • resampling (str, default: 'nearest' ) –

    Resampling method forwarded to :func:build_vrt, by default "nearest". Only used when use_vrt=True; ignored by the in-memory mosaic. See :func:build_vrt.

  • use_vrt (bool, default: False ) –

    If True, mosaic via a lazily-read VRT built with the gdalbuildvrt CLI (requires the conda libgdal/gdal package, not just libgdal-core). If False (default), mosaic in memory with :func:rioxarray.merge.merge_arrays, which needs no GDAL command-line tool — only rasterio/rioxarray.

Returns:

  • DataArray –

    DataArray containing the clipped data.

elevation_bygrid #

elevation_bygrid(
    longs, lats, res=10, window=5, resampling=1
)

Sample elevation from 3DEP at a grid of lon/lat coordinates.

Notes

Reads directly from the USGS seamless DEM VRT (Cloud-Optimized GeoTIFFs, EPSG:4269). A small pixel window around each query point is read and downsampled to a single value using the chosen resampling kernel.

Parameters:

  • longs (array - like) –

    1D sequence of longitude values in decimal degrees.

  • lats (array - like) –

    1D sequence of latitude values in decimal degrees.

  • res ((10, 30, 60), default: 10 ) –

    Source DEM resolution in meters, by default 10. Choose a coarser source if you don't need the 10 m fidelity — coarser VRTs are smaller and faster to range-read.

  • window (int, default: 5 ) –

    Size of the read window for interpolation, must be odd, defaults to 5.

  • resampling (int, default: 1 ) –

    Resampling method from rasterio.enums.Resampling, defaults to 1 (bilinear). Methods applicable to DEM interpolation:

    • 0: nearest — fastest, no interpolation
    • 1: bilinear — good general-purpose default
    • 2: cubic — sharper than bilinear
    • 3: cubic_spline — smooth spline interpolation
    • 4: lanczos — high-quality windowed sinc

Returns:

  • ndarray –

    2D array of shape (len(lats), len(longs)) with elevation values in meters.