Download DEMs for a study area
For a US study area, download DEMs through The National Map: draw or enter your area, select an elevation product, inspect the returned footprints, then download each GeoTIFF with its metadata1. The TNMAccess API provides the same kind of product information for a repeatable acquisition record. The worked example below downloads two neighbouring one-metre 3DEP tiles, records their identities and hashes, and inspects their actual raster headers.
Find the product that matches the work
Start at The National Map’s data-delivery page and open its download application. Search for your place or set the map extent, choose the Elevation Products/3DEP category, select the DEM resolution you need, and search the area. Examine each result’s footprint and metadata before adding its download to your project. A search rectangle is a selection aid; it is not necessarily the final boundary of your analysis.
Choose the product family deliberately. A one-metre DEM is a different delivery from a one-third-arc-second or one-arc-second DEM. “One metre” describes grid spacing, not a promise that every elevation is accurate to one metre. Keep the product description, collection dates and accuracy information alongside the resolution choice. For a broad overview, the elevation product comparison helps distinguish families before downloading.
The example uses two tiles from one named collection, CO_ArapahoRooseveltPikeNF_D23, near Boulder, Colorado. Both were published on 20 March 2026. Their ScienceBase product identifiers end in 33c9a and 33c98; the full identifiers are recorded by the script. Keeping identifiers matters because a similar filename or a newer publication may otherwise look like the same input.
Keep the API response, not just the download link
The request supplies the one-metre dataset name, a longitude/latitude bounding box and a maximum result count. Its response contains product IDs, publication dates, download addresses, metadata addresses and reported file sizes2. The script selects two explicit IDs from that response and stops if either is absent. It therefore does not silently replace a missing tile with the next search result.
For another study area, change the search and select its returned IDs together. Inspect the full result set rather than assuming this small max=5 example exhausts a larger area. Save the request and response before processing; they explain why those particular files entered the project. If a publisher updates an object at the same address, your stored SHA-256 distinguishes the bytes you used.
The filenames also carry useful identity. USGS_1M_13_x47y443_CO_ArapahoRooseveltPikeNF_D23.tif names the one-metre product, zone and tile within the collection. Its neighbour changes y443 to y444. Use that pattern to organise files, but obtain the coordinate reference and extent from metadata and the header rather than deriving them from the filename alone.
Download and inspect the exact inputs
Use Bash with curl, Python and GDAL available in the active environment. Set DEM_WORK to a new local directory with room for the source rasters and later outputs. The two downloads here total 36,235,737 bytes; the subsequent mosaic needs additional working space. The script saves the API response, selected records, XML, text and JSON inspections, and a download manifest in that directory.
The loop obtains both the raster and its vendorMetaUrl, then calls gdalinfo on the local file. The last section computes hashes from downloaded bytes rather than copying a catalogue value. It extracts NAVD88 and metres from the actual XML abstract, and records that location explicitly because these XML files do not supply separate altdatum and altunits elements.
# SPDX-License-Identifier: MIT
set -euo pipefail
: "${DEM_WORK:?Set DEM_WORK to a new local working directory}"
mkdir -p "$DEM_WORK"
cd "$DEM_WORK"
gdalinfo --version
api='https://tnmaccess.nationalmap.gov/api/v1/products?datasets=Digital%20Elevation%20Model%20%28DEM%29%201%20meter&bbox=-105.28,40.01,-105.27,40.02&max=5'
curl -fsSL --retry 2 "$api" -o tnm-products.json
python - <<'PY'
import json
from pathlib import Path
wanted = {'69c20ed7b66b01d6ebf33c9a', '69c20ed6b66b01d6ebf33c98'}
items = [x for x in json.loads(Path('tnm-products.json').read_text())['items'] if x['sourceId'] in wanted]
assert len(items) == 2, 'Expected both selected USGS products'
Path('selected-products.json').write_text(json.dumps(items, indent=2))
with Path('downloads.tsv').open('w') as f:
for item in items:
name = item['downloadURL'].rsplit('/', 1)[1]
f.write('\t'.join([name, item['downloadURL'], item['vendorMetaUrl']]) + '\n')
PY
while IFS=$'\t' read -r name url metadata; do
curl -fsSL --retry 2 "$url" -o "$name"
curl -fsSL --retry 2 "$metadata" -o "${name%.tif}.xml"
gdalinfo "$name" > "${name%.tif}-gdalinfo.txt"
gdalinfo -json "$name" > "${name%.tif}-gdalinfo.json"
done < downloads.tsv
python - <<'PY'
import hashlib, json, datetime, xml.etree.ElementTree as ET
from pathlib import Path
rows = []
for item in json.loads(Path('selected-products.json').read_text()):
p = Path(item['downloadURL'].rsplit('/', 1)[1])
xml = p.with_suffix('.xml'); tree = ET.parse(xml)
info = json.loads(p.with_name(p.stem + '-gdalinfo.json').read_text())
band = info['bands'][0]
abstract = tree.findtext('.//abstract') or ''
row = dict(product_id=item['sourceId'], publication_date=item['publicationDate'],
url=item['downloadURL'], metadata_url=item['vendorMetaUrl'],
bytes=p.stat().st_size, sha256=hashlib.sha256(p.read_bytes()).hexdigest(),
metadata_sha256=hashlib.sha256(xml.read_bytes()).hexdigest(),
size=info['size'], geotransform=info['geoTransform'],
band_type=band['type'], nodata=band.get('noDataValue'),
crs=info['coordinateSystem']['wkt'].split(chr(34))[1],
vertical_datum='NAVD88' if 'NAVD88' in abstract else None,
elevation_units='metres' if 'elevation values are in meters' in abstract else None,
vertical_source='XML abstract; no dedicated altdatum/altunits elements')
rows.append(row)
Path('download-manifest.json').write_text(json.dumps({'retrieved_utc':datetime.datetime.now(datetime.timezone.utc).isoformat(), 'tiles':rows}, indent=2))
print(json.dumps(rows, indent=2))
PY
Read the output before using the surface
The relevant output was:
GDAL 3.11.5 "Eganville", released 2025/10/29
South tile: 10,425,274 bytes; 10012 × 10012; Float32
North tile: 25,810,463 bytes; 10012 × 10012; Float32
CRS: NAD83 / UTM zone 13N (EPSG:26913)
South geotransform: (469993.9999821812, 1, 0, 4430005.999955563, 0, -1)
North geotransform: (469993.99998021487, 1, 0, 4440005.999968549, 0, -1)
Band NoData: -999999
XML abstract: NAVD88 elevations; metres
The two negative values have different meanings. The geotransform’s -1 is the north-up grid’s row step; -999999 is the band’s missing-value marker. Neither should be treated as a measured terrain height. The saved JSON keeps the full coordinate system and numeric transform, while the compact output above exposes the fields most likely to affect the next operation3.
The horizontal reference in the raster is NAD83/UTM zone 13N. The vertical reference comes from the XML’s description of elevations: NAVD88, in metres. Keeping both prevents a horizontal CRS label from being mistaken for a complete explanation of the height values. Record the metadata wording and its hash even when your GIS application displays only the horizontal system.
The footprints overlap slightly, but that alone does not establish continuous terrain coverage. Continue to the worked mosaic and NoData check before assuming that the shared border contains valid heights. That check reads the actual overlap and compares the materialised mosaic with both source footprints.
Leave an acquisition record someone else can follow
Keep the original filenames, full product IDs, publication dates, request URL, downloaded metadata and hashes together. Add your study boundary and the reason for selecting this resolution. Preserve the unmodified inputs separately from derived mosaics so a later correction can be reproduced without guessing which files were edited.
Use the publication date to identify the delivered product and collection information to understand the observation period. Those dates answer different questions. A 2026 publication does not by itself mean that the terrain was collected in 2026. Before analysing change over time, inspect the metadata’s collection dates and processing description rather than substituting the download date.
Not covered on this page
- Selecting an accuracy threshold for a scientific or engineering decision.
- Transforming horizontal or vertical coordinates.
- Downloading every tile in a large region or managing subsequent product revisions.
- Validating terrain against independent survey observations.
Sources
Code on this page is original and MIT licensed.
How we test: methodology. Found a mistake? Corrections.