Add 3DEP_1M data source for native 1-meter resolution DEMs

The existing 3DEP source uses the ImageServer which resamples to ~30m.
3DEP_1M queries the TNM Access API for 1m tiles and fetches them as
COGs from S3 via /vsicurl/. Uses a two-step GDAL warp: first reproject
from source UTM to EPSG:4326, then convert NAVD88 to WGS84 ellipsoidal.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: scottstanie

README.md

@@ -181,6 +181,21 @@ No authentication is required.
 sardem --bbox -104 30 -103 31 --data-source 3DEP
 ```
 
+## USGS 3DEP 1-Meter Lidar DEM
+
+The `--data-source 3DEP_1M` option provides access to 1-meter resolution lidar-derived DEMs from the USGS 3DEP program. Tiles are discovered via the [TNM Access API](https://tnmaccess.nationalmap.gov/api/v1/products) and fetched as Cloud Optimized GeoTIFFs (COGs) from S3.
+
+- **US coverage only** — not all areas have 1m lidar data. If no tiles are found, the command will raise an error.
+- **No authentication required.**
+- Output resolution matches the native ~1m resolution of the source tiles.
+- `--xrate`/`--yrate` upsampling is not supported (data is already high-resolution).
+
+### Usage
+
+```bash
+sardem --bbox -118.4 33.7 -118.3 33.8 --data-source 3DEP_1M
+```
+
 ## Citations and Acknowledgments
 
 ### Copernicus DEM

sardem/cli.py

@@ -41,6 +41,7 @@ def positive_small_int(argstring):
         sardem --bbox -156 18.8 -154.7 20.3 --data-source NASA_WATER -o my_watermask.wbd # Water mask
         sardem --bbox -156 18.8 -154.7 20.3 --data COP -isce  # Generate .isce XML files as well
         sardem --bbox -104 30 -103 31 --data-source 3DEP  # USGS 3DEP lidar DEM (US only)
+        sardem --bbox -118.4 33.7 -118.3 33.8 --data-source 3DEP_1M  # USGS 3DEP 1m lidar (US only)
         sardem --bbox -104 30 -103 31 --data-source NISAR  # NISAR DEM (requires Earthdata login)
 
 

sardem/dem.py

@@ -385,6 +385,25 @@ def main(
             utils.gdal2isce_xml(output_name, keep_egm=keep_egm)
         return
 
+    # For USGS 3DEP 1m, download tiles from S3 COGs via TNM API
+    if data_source == "3DEP_1M":
+        utils._gdal_installed_correctly()
+        from sardem import usgs_3dep_1m
+
+        usgs_3dep_1m.download_and_stitch(
+            output_name,
+            bbox,
+            keep_egm=keep_egm,
+            xrate=xrate,
+            yrate=yrate,
+            output_format=output_format,
+            output_type=output_type,
+        )
+        if make_isce_xml:
+            logger.info("Creating ISCE2 XML file")
+            utils.gdal2isce_xml(output_name, keep_egm=keep_egm)
+        return
+
     if data_source == "NISAR":
         if keep_egm:
             logger.info(

sardem/download.py

@@ -199,7 +199,8 @@ class Downloader:
             "https://elevation.nationalmap.gov/arcgis/rest/services"
             "/3DEPElevation/ImageServer/exportImage"
         ),
-        "NISAR": "https://nisar.asf.earthdatacloud.nasa.gov/NISAR/DEM/v1.2/EPSG4326/EPSG4326.vrt"
+        "NISAR": "https://nisar.asf.earthdatacloud.nasa.gov/NISAR/DEM/v1.2/EPSG4326/EPSG4326.vrt",
+        "3DEP_1M": "https://tnmaccess.nationalmap.gov/api/v1/products",
     }
     VALID_SOURCES = DATA_URLS.keys()
     TILE_ENDINGS = {

sardem/tests/test_usgs_3dep_1m.py

@@ -0,0 +1,63 @@
+import json
+
+import pytest
+import responses
+
+from sardem.usgs_3dep_1m import TNM_API_URL, _find_tile_urls
+
+
+def _make_tnm_response(items):
+    """Build a fake TNM API JSON response body."""
+    return json.dumps({"items": items})
+
+
+def _make_item(url, date_created="2020-01-01"):
+    """Build a minimal TNM API item dict."""
+    return {"downloadURL": url, "dateCreated": date_created}
+
+
+@responses.activate
+def test_find_tile_urls_success():
+    """Mock TNM API response, verify URL extraction and sorting."""
+    items = [
+        _make_item("https://prd-tnm.s3.amazonaws.com/tile_newer.tif", "2022-06-15"),
+        _make_item("https://prd-tnm.s3.amazonaws.com/tile_older.tif", "2019-03-01"),
+    ]
+    responses.add(
+        responses.GET,
+        TNM_API_URL,
+        body=_make_tnm_response(items),
+        status=200,
+        content_type="application/json",
+    )
+
+    urls = _find_tile_urls((-118.4, 33.7, -118.3, 33.8))
+
+    assert len(urls) == 2
+    # Oldest first so newest overwrites in gdal.Warp
+    assert "tile_older" in urls[0]
+    assert "tile_newer" in urls[1]
+
+
+@responses.activate
+def test_find_tile_urls_empty():
+    """Mock empty response, verify RuntimeError."""
+    responses.add(
+        responses.GET,
+        TNM_API_URL,
+        body=_make_tnm_response([]),
+        status=200,
+        content_type="application/json",
+    )
+
+    with pytest.raises(RuntimeError, match="No 3DEP 1m DEM tiles found"):
+        _find_tile_urls((-118.4, 33.7, -118.3, 33.8))
+
+
+@responses.activate
+def test_find_tile_urls_http_error():
+    """Mock 500, verify exception."""
+    responses.add(responses.GET, TNM_API_URL, status=500)
+
+    with pytest.raises(Exception):
+        _find_tile_urls((-118.4, 33.7, -118.3, 33.8))

sardem/usgs_3dep_1m.py

@@ -0,0 +1,188 @@
+"""Download USGS 3DEP 1-meter DEM tiles from S3 via the TNM Access API.
+
+Uses the USGS TNM (The National Map) Access API to discover which 1-meter
+DEM tiles cover the requested bounding box, then fetches them as Cloud
+Optimized GeoTIFFs (COGs) from S3 using GDAL's /vsicurl/ virtual filesystem.
+
+Source tiles are in NAD83/UTM (varying zones) with NAVD88 heights, but without
+compound vertical CRS metadata. Horizontal reprojection to EPSG:4326 is always
+performed. When keep_egm=False (the default), a second GDAL warp converts
+NAVD88 heights to WGS84 ellipsoidal heights.
+
+Coverage: US only -- not all areas have 1m lidar data available.
+
+References:
+    https://www.usgs.gov/3d-elevation-program
+    https://tnmaccess.nationalmap.gov/api/v1/products
+"""
+
+import logging
+import os
+import tempfile
+
+import requests
+
+from sardem import utils
+
+logger = logging.getLogger("sardem")
+utils.set_logger_handler(logger)
+
+TNM_API_URL = "https://tnmaccess.nationalmap.gov/api/v1/products"
+TNM_DATASET = "Digital Elevation Model (DEM) 1 meter"
+TNM_MAX_ITEMS = 200
+
+
+def download_and_stitch(
+    output_name,
+    bbox,
+    keep_egm=False,
+    xrate=1,
+    yrate=1,
+    output_format="GTiff",
+    output_type="float32",
+):
+    """Download USGS 3DEP 1m DEM tiles and mosaic them with gdal.Warp.
+
+    Parameters
+    ----------
+    output_name : str
+        Path for the output DEM file.
+    bbox : tuple
+        (left, bottom, right, top) in decimal degrees.
+    keep_egm : bool
+        If True, keep NAVD88 geoid heights. If False (default), convert to
+        WGS84 ellipsoidal heights.
+    xrate : int
+        Upsample factor in x direction (ignored for 1m data).
+    yrate : int
+        Upsample factor in y direction (ignored for 1m data).
+    output_format : str
+        GDAL output format (default GTiff).
+    output_type : str
+        Output pixel type (default float32).
+    """
+    from osgeo import gdal
+
+    gdal.UseExceptions()
+
+    if xrate > 1 or yrate > 1:
+        logger.warning(
+            "xrate/yrate upsampling is ignored for 3DEP_1M (data is already"
+            " ~1m resolution). xrate=%d, yrate=%d.",
+            xrate,
+            yrate,
+        )
+
+    tile_urls = _find_tile_urls(bbox)
+    vsicurl_paths = ["/vsicurl/" + url for url in tile_urls]
+    logger.info("Found %d tile(s) covering the bounding box", len(vsicurl_paths))
+
+    out_type = gdal.GetDataTypeByName(output_type.title())
+
+    # Source tiles are in NAD83/UTM (various zones). Do NOT override srcSRS --
+    # GDAL must read the CRS from each file's metadata for correct reprojection.
+    reproject_opts = dict(
+        format=output_format,
+        outputBounds=list(bbox),
+        dstSRS="EPSG:4326",
+        outputType=out_type,
+        resampleAlg="nearest",
+        srcNodata=-999999,
+        multithread=True,
+        warpMemoryLimit=5000,
+        warpOptions=["NUM_THREADS=4"],
+    )
+
+    if keep_egm:
+        logger.info("Creating %s (keeping NAVD88 heights)", output_name)
+        reproject_opts["callback"] = gdal.TermProgress
+        gdal.Warp(
+            output_name, vsicurl_paths, options=gdal.WarpOptions(**reproject_opts)
+        )
+    else:
+        # Two-step warp:
+        # 1) Reproject from UTM to EPSG:4326 (horizontal only, preserves NAVD88 Z)
+        # 2) Convert NAVD88 heights to WGS84 ellipsoidal using compound CRS.
+        #    After step 1 the file is in geographic coords, so
+        #    srcSRS="EPSG:4269+5703" (NAD83 + NAVD88) is correct.
+        fd, tmp_path = tempfile.mkstemp(suffix=".tif")
+        os.close(fd)
+        try:
+            logger.info("Reprojecting tiles to EPSG:4326...")
+            reproject_opts["callback"] = gdal.TermProgress
+            gdal.Warp(
+                tmp_path, vsicurl_paths, options=gdal.WarpOptions(**reproject_opts)
+            )
+
+            logger.info("Converting NAVD88 heights to WGS84 ellipsoidal...")
+            vert_opts = dict(
+                format=output_format,
+                srcSRS="EPSG:4269+5703",
+                dstSRS="EPSG:4326",
+                outputType=out_type,
+                multithread=True,
+                callback=gdal.TermProgress,
+            )
+            gdal.Warp(
+                output_name, tmp_path, options=gdal.WarpOptions(**vert_opts)
+            )
+        finally:
+            if os.path.exists(tmp_path):
+                os.remove(tmp_path)
+
+
+def _find_tile_urls(bbox):
+    """Query TNM Access API for 1m DEM tiles covering bbox.
+
+    Returns list of S3 download URLs, sorted oldest-first so
+    newest tiles take priority in gdal.Warp overlap resolution.
+
+    Parameters
+    ----------
+    bbox : tuple
+        (left, bottom, right, top) in decimal degrees.
+
+    Returns
+    -------
+    list[str]
+        S3 URLs to COG tiles.
+
+    Raises
+    ------
+    RuntimeError
+        If no tiles are found for the requested area.
+    """
+    left, bottom, right, top = bbox
+    params = {
+        "datasets": TNM_DATASET,
+        "bbox": "{},{},{},{}".format(left, bottom, right, top),
+        "max": TNM_MAX_ITEMS,
+        "outputFormat": "JSON",
+    }
+
+    logger.info("Querying TNM API for 1m DEM tiles...")
+    response = requests.get(TNM_API_URL, params=params, timeout=60)
+    response.raise_for_status()
+
+    data = response.json()
+    items = data.get("items", [])
+    if not items:
+        raise RuntimeError(
+            "No 3DEP 1m DEM tiles found for bbox {}. "
+            "Not all US areas have 1m coverage.".format(bbox)
+        )
+
+    if len(items) >= TNM_MAX_ITEMS:
+        logger.warning(
+            "TNM API returned the maximum %d items. The bounding box may be"
+            " too large to fetch all tiles in one request.",
+            TNM_MAX_ITEMS,
+        )
+
+    # Sort by dateCreated ascending (oldest first) so that when gdal.Warp
+    # processes them in order, newer tiles overwrite older ones in overlaps
+    items.sort(key=lambda item: item.get("dateCreated", ""))
+
+    urls = [item["downloadURL"] for item in items]
+    logger.info("Found %d tiles from TNM API", len(urls))
+    return urls