Reading USGS 3DEP LiDAR from Public Cloud Storage

TL;DR: The USGS 3DEP LiDAR collection is published on AWS as Entwine Point Tile (EPT) datasets, one per project, in the public usgs-lidar-public bucket. Look up which datasets cover your area in the public boundaries index, then point readers.ept at https://s3-us-west-2.amazonaws.com/usgs-lidar-public/<project>/ept.json with a bounds or polygon for your area. PDAL fetches only the octree nodes that intersect it, so you can read a few square kilometres from a project of billions of points in seconds.

# Context and Motivation

This guide is part of S3 Cloud Storage I/O. Downloading 3DEP data used to mean finding the right project, downloading whole tiles through a web map, and clipping locally. The EPT copies on AWS change that: each project is an octree of LAZ files on S3, indexed by a small ept.json, and PDAL’s readers.ept can query it by area and level of detail. No account or credentials are needed for these public datasets, and nothing outside your area is downloaded.

That makes 3DEP a practical input for batch jobs — a site-suitability screen over hundreds of parcels, a flood model needing a fresh DTM for one catchment — without keeping a local copy of the national collection.

# Prerequisites and Assumptions

• PDAL 2.4 or newer with readers.ept (standard in conda-forge builds) and python-pdal.
• GeoPandas and Shapely to query the boundaries index.
• Internet access to s3-us-west-2.amazonaws.com; compute in AWS us-west-2 is fastest.

# Step-by-Step Implementation

# Step 1 — Find datasets that cover your area

Load the public resources.geojson boundaries index, maintained alongside the collection, and intersect it with your area of interest. Each feature has a name and an url pointing at its ept.json.

# Step 2 — Express the area in the dataset’s CRS

The EPT copies are stored in Web Mercator (EPSG:3857). Either reproject your polygon to 3857 or append the CRS to the bounds string, for example ([xmin, xmax], [ymin, ymax])/EPSG:4326.

# Step 3 — Read with readers.ept

Pass filename, bounds or polygon, and optionally resolution to stop at a coarser octree level for previews.

# Step 4 — Reproject and filter

Add filters.reprojection to a local projected CRS, and filters.range to drop noise classes 7 and 18.

# Step 5 — Write locally or to your own bucket

Write LAZ or COPC for reuse, or go straight to a DTM with writers.gdal.

# Complete Working Example

"""Fetch 3DEP points for a polygon from public EPT and write a local COPC."""
import json

import geopandas as gpd
import pdal
from shapely.geometry import box

INDEX = "https://raw.githubusercontent.com/hobuinc/usgs-lidar/master/boundaries/resources.geojson"
aoi = gpd.GeoDataFrame(geometry=[box(-105.29, 40.00, -105.26, 40.02)], crs="EPSG:4326")

# 1. which projects cover the area?
index = gpd.read_file(INDEX).to_crs(4326)
hits = index[index.intersects(aoi.geometry.iloc[0])]
print(hits[["name", "count"]].to_string(index=False))
ept_url = hits.iloc[0]["url"]                      # pick the most recent or densest

# 2. polygon in the dataset's CRS (EPSG:3857)
poly_3857 = aoi.to_crs(3857).geometry.iloc[0].wkt

# 3-5. read, clean, reproject, write
spec = {"pipeline": [
    {"type": "readers.ept", "filename": ept_url, "polygon": poly_3857, "threads": 8},
    {"type": "filters.range", "limits": "Classification![7:7],Classification![18:18]"},
    {"type": "filters.reprojection", "out_srs": "EPSG:26913+5703"},
    {"type": "writers.copc", "filename": "boulder_aoi.copc.laz", "forward": "all"},
]}
p = pdal.Pipeline(json.dumps(spec))
n = p.execute()
print(f"{n:,} points written")

A quick preview at reduced density uses resolution, in the dataset’s units (metres in EPSG:3857):

pdal translate https://s3-us-west-2.amazonaws.com/usgs-lidar-public/<project>/ept.json preview.laz \
  --readers.ept.bounds="([-11722000, -11718000], [4865000, 4869000])" \
  --readers.ept.resolution=5

# Choosing Between Overlapping Projects

Many places are covered by several 3DEP projects: an older county collection, a newer statewide one, sometimes a specialised survey. The boundaries index lists them all, and the choice matters. Prefer the most recent project unless you need a specific date for change analysis, and check its point count against the area to estimate density — a Quality Level 1 project has several times the density of an older QL2 one. Project names usually carry a year, and the metadata linked from the USGS project pages gives the collection dates, accuracy reports and classification scheme, which you should read before trusting a classification for anything important. Reading two projects over the same area and differencing their DTMs is also an easy first change-detection experiment.

Remember that project boundaries are approximate footprints, not exact coverage; gaps and irregular edges inside a boundary are normal, so check point counts after reading.

# Key Parameter Table

# Verification

• Point count. Compare the count with the project’s approximate density times your area; a much smaller number usually means the polygon was in the wrong CRS.
• Bounds. pdal info --summary on the output shows bounds inside your area after reprojection.
• Classification. pdal info --stats shows the classes present; some older projects lack a full ASPRS classification.

# Gotchas and Edge Cases

Wrong CRS for the query. Passing a longitude and latitude box without a CRS suffix is interpreted in Web Mercator and selects a tiny area near the null island. Always reproject or append the CRS.

Mixed vertical datums. Projects differ in vertical datum and geoid; when combining two, reproject both to the same compound CRS before comparing heights.

Service changes. The public index location and bucket layout are maintained by the community project behind the collection; if a URL moves, check the project’s repository for the current index.

# Frequently Asked Questions

Do I need AWS credentials to read 3DEP EPT data?

No. The Entwine Point Tile copies of 3DEP are in a public bucket and are read over HTTPS, so readers.ept works without an AWS account or credentials.

How do I find which 3DEP project covers my area?

Load the public boundaries GeoJSON that lists every EPT dataset with its footprint and URL, and intersect it with your area of interest. Several projects may overlap.

What CRS are the 3DEP EPT datasets in?

Web Mercator, EPSG:3857. Transform your query polygon to that CRS, or append the CRS to the bounds string, and reproject the output to a local projected CRS for analysis.

How can I get a quick low-density preview?

Set the resolution option on readers.ept to a few metres. PDAL then stops at a coarser octree level and downloads far fewer nodes.

Can I process 3DEP tiles in batch without downloading the whole project?

Yes. Split your area into tiles, and have each batch job read its own tile bounds, plus a small buffer, from the EPT dataset. Each job downloads only the octree nodes for its tile.

# Related

• S3 Cloud Storage I/O — reading and writing LiDAR in object storage
• Streaming LAZ from S3 with PDAL — reading individual LAZ tiles
• Reading COPC in Python with laspy — the single-file cloud format
• Choosing a Projected CRS for a LiDAR Project — picking out_srs