py3dep 0.1.3

Access USGS 3DEP database and get elevation in the US

Features

Py3DEP is a part of Hydrodata software stack and provides access to 3DEP which is a part the National Map services. The 3DEP service has multi-resolution sources and depending on the user provided resolution, the data is resampled on server-side based on all the available data sources. Py3DEP returns All the requestes as xarray.DataArray or xarray.Dataset that allows the user to process the data using all the xarray’s functionalities. The following layers are available:

• DEM

• Hillshade Gray

• Aspect Degrees

• Aspect Map

• GreyHillshade Elevation Fill

• Hillshade Multidirectional

• Slope Map

• Slope Degrees

• Hillshade Elevation Tinted

• Height Ellipsoidal

• Contour 25

• Contour Smoothed 25

Moreover, Py3DEP offers thress additonal function:

• elevation_bygrid: For getting elevations of all the grid points in a 2D array of x- and y-coordinates.

• elevation_byloc: For getting elevation of a single point based on the National Map’s Elevation Point Query Service service.

• deg2mpm: For converting slop from degree to meter per meter.

Moreover, requests for additional databases or functionalities can be submitted via issue tracker.

Installation

You can install py3dep using pip:

$ pip install py3dep

Quickstart

Py3DEP accepts Shapely’s Polygon or a bounding box (a tuple of length four) as an input geometry. We can use Hydrodata to get a watershed’s geometry, then use it to get DEM and slope data in meters/meters from Py3DEP.

The get_map function has a resolution argument that sets the target resolution in meters. Note that the highest available resolution throughout the CONUS is about 10 m, though higher resolutions are available in limited parts of the US. Note that the input geometry can be in any valid spatial reference (geo_crs). The crs argument, however, is limited to CRS:84, EPSG:4326, and EPSG:3857 since 3DEP only supports these spatial references.

import py3dep
from hydrodata import NLDI

geom = NLDI().getfeature_byid("nwissite", "USGS-01031500", basin=True).geometry[0]
dem = py3dep.get_map("DEM", geom, resolution=30, geo_crs="epsg:4326", crs="epsg:3857")
slope = py3dep.get_map("Slope Degrees", geom, resolution=30)
slope = py3dep.utils.deg2mpm(slope)

We can get the elevation for a single point within the US:

elev = py3dep.elevation_byloc((-7766049.665, 5691929.739), "epsg:3857")

Additionally, we can get the elevations of set of x- and y- coordinates of a grid. For example, let’s get the minimum temperature data within the watershed from Daymet using Hydrodata then add the elevation as a new variable to the dataset:

import hydrodata.datasets as hds
import xarray as xr
import numpy as np

clm = hds.daymet_bygeom(geom, dates=("2005-01-01", "2005-01-31"), variables="tmin")
gridxy = (clm.x.values, clm.y.values)
res = clm.res[0] * 1000
elev = py3dep.elevation_bygrid(gridxy, clm.crs, res)
clm = xr.merge([clm, elev], combine_attrs="override")
clm["elevation"] = clm.elevation.where(~np.isnan(clm.isel(time=0).tmin), drop=True)

Contributing

If you are interested in contributing to the project please get in touch. You can find information about contributing to py3dep at our Contributing page.