Add TNM elevation product and source-coverage support

pyproject.toml

@@ -35,6 +35,7 @@ dependencies = [
     "pyproj",
     "shapely",
     "pyogrio",
+    "rasterio>1.4.0",
 ]
 
 keywords = ["Geospatial"]
@@ -128,6 +129,7 @@ module = [
     "filelock",
     "pyogrio.*",
     "pyogrio.raw",
+    "rasterio",
     "pyproj",
     "sentinelsat",
     "shapefile",

src/fetchez/cli/pipeline.py

@@ -81,6 +81,14 @@ def get_command(self, ctx, name):
         if bundle_yml:
             help_text = bundle_yml.get("description", "")
             mod_args = []
+            if bundle_yml.get("products"):
+                mod_args.append(
+                    [
+                        "--products",
+                        "Comma-separated products to include (default: all)",
+                        "all",
+                    ]
+                )
 
         @click.command(name=name, help=help_text, hidden=True, cls=FetchezMainCommand)
         @click.option("--weight", type=float, default=1.0)

src/fetchez/core.py

@@ -405,10 +405,11 @@ def __init__(self, url, session=None, callback=None):
         self.size = self._get_size()
 
     def _get_size(self):
-        resp = self.session.head(self.url)
-        if "Content-Length" not in resp.headers:
-            return 0
-        return int(resp.headers["Content-Length"])
+        with self.session.head(self.url, timeout=(10, 60)) as resp:
+            resp.raise_for_status()
+            if "Content-Length" not in resp.headers:
+                raise OSError("HTTP file has no Content-Length")
+            return int(resp.headers["Content-Length"])
 
     def seek(self, offset, whence=io.SEEK_SET):
         if whence == io.SEEK_SET:
@@ -436,10 +437,19 @@ def read(self, size=-1):
 
         # Fetch ONLY the specific bytes requested
         headers = {"Range": f"bytes={self.offset}-{end}"}
-        response = self.session.get(self.url, headers=headers, timeout=(10, 60))
-        response.raise_for_status()
-
-        data = response.content
+        with self.session.get(
+            self.url, headers=headers, timeout=(10, 60), stream=True
+        ) as response:
+            response.raise_for_status()
+            expected = f"bytes {self.offset}-{end}/{self.size}"
+            if (
+                response.status_code != 206
+                or response.headers.get("Content-Range") != expected
+            ):
+                raise OSError("HTTP server did not honor the requested byte range")
+            data = response.raw.read(end - self.offset + 2, decode_content=True)
+            if len(data) != end - self.offset + 1:
+                raise OSError("HTTP server returned an incomplete byte range")
 
         if self.callback:
             self.callback(len(data))

src/fetchez/hooks/stream_init.py

@@ -14,6 +14,10 @@
 import math
 import logging
 
+from pyproj import CRS
+from pyproj.crs import CompoundCRS
+from pyproj.exceptions import CRSError
+
 from fetchez.spatial import Region
 from fetchez.hooks import FetchHook
 from fetchez.registry import ReaderRegistry, ProfileRegistry
@@ -124,8 +128,20 @@ def run(self, entries):
                         base_srs = reader.get_srs() or base_srs
 
                     vert_srs = kwargs_copy.get("vert_srs")
-                    if vert_srs and "+" not in base_srs:
-                        base_srs = f"{base_srs}+{vert_srs}"
+                    if vert_srs:
+                        try:
+                            horizontal = CRS.from_user_input(base_srs)
+                            vertical = CRS.from_user_input(vert_srs)
+                        except CRSError:
+                            # Retain custom vertical references such as global:mss.
+                            if "+" not in base_srs:
+                                base_srs = f"{base_srs}+{vert_srs}"
+                        else:
+                            if len(horizontal.axis_info) == 2:
+                                base_srs = CompoundCRS(
+                                    f"{horizontal.name} + {vertical.name}",
+                                    [horizontal, vertical],
+                                ).to_wkt()
 
                     logger.debug(f"[{self.name}] Using SRS: {base_srs}")
 

src/fetchez/modules/tnm.py

@@ -20,6 +20,8 @@
 from fetchez import utils
 from fetchez import spatial
 from fetchez import cli
+from fetchez.modules import tnm_ned, tnm_raster
+from fetchez.modules.tnm_wesm import WESM
 
 logger = logging.getLogger(__name__)
 
@@ -56,12 +58,19 @@
     "US Topo Historical",
     "Land Cover - Woodland",
     "3D Hydrography Program (3DHP)",
+    "Seamless 1-m DEM (S1M)",
 ]
 DATASET_ALIASES = {
     "1m": 2,
     "1_9as": 4,
     "1_3as": 3,
     "1_as": 1,
+    "2_as": 5,
+    "5m": 6,
+    "s1m": 29,
+}
+DATASET_PRODUCTS = {
+    DATASET_CODES[index]: alias for alias, index in DATASET_ALIASES.items()
 }
 
 
@@ -76,6 +85,8 @@
     q="Free text search query",
     date_start="Start date (YYYY-MM-DD)",
     date_end="End date (YYYY-MM-DD)",
+    strict_datasets="Fail instead of broadening a rejected dataset query",
+    source_coverage="Attach authoritative USGS source coverage",
 )
 class TheNationalMap(FetchModule):
     name = "tnm"
@@ -114,6 +125,8 @@ def __init__(
         date_start: Optional[str] = None,
         date_end: Optional[str] = None,
         dedupe: bool = True,
+        strict_datasets: bool = False,
+        source_coverage: bool = False,
         **kwargs,
     ):
         super().__init__(name="tnm", **kwargs)
@@ -127,6 +140,12 @@ def __init__(
         self.dedupe = utils.str2bool(dedupe)
         if self.dedupe is None:
             self.dedupe = True
+        self.strict_datasets = utils.str2bool(strict_datasets)
+        if self.strict_datasets is None:
+            self.strict_datasets = False
+        self.source_coverage = utils.str2bool(source_coverage)
+        if self.source_coverage is None:
+            self.source_coverage = False
 
     def run(self):
         """Run the TNM fetching module."""
@@ -138,7 +157,7 @@ def run(self):
         bbox_str = f"{w},{s},{e},{n}"
 
         offset = 0
-        total = 0
+        expected_total: Optional[int] = None
 
         # Determine Datasets to query
         dataset_names = []
@@ -160,6 +179,9 @@ def run(self):
 
         if not dataset_names:
             dataset_names = ["National Elevation Dataset (NED) 1 arc-second"]
+        product = (
+            DATASET_PRODUCTS.get(dataset_names[0]) if len(dataset_names) == 1 else None
+        )
 
         best_tiles = {}
         all_tiles = []
@@ -188,10 +210,12 @@ def run(self):
             req = core.Fetch(TNM_API_PRODUCTS_URL).fetch_req(params=params)
 
             if (
-                req
+                req is not None
                 and "All dataset queries failed" in req.text
                 and "datasets" in params
             ):
+                if self.strict_datasets or self.source_coverage:
+                    raise RuntimeError("TNM API rejected the strict dataset query")
                 logger.warning(
                     "USGS rejected the strict dataset strings. Retrying with broad text search..."
                 )
@@ -207,24 +231,44 @@ def run(self):
                 req = core.Fetch(TNM_API_PRODUCTS_URL).fetch_req(params=params)
 
             if req is None or req.status_code != 200:
-                logger.error(
-                    f"TNM API Failed: {req.status_code if req else 'No Response'}"
-                )
-                break
-
-            if req.text.strip().startswith("{errorMessage"):
-                logger.error(f"TNM API Error: {req.text}")
-                break
+                status = req.status_code if req is not None else "no response"
+                raise RuntimeError(f"TNM API request failed: {status}")
 
             try:
-                data = req.json()
-                total = data.get("total", 0)
-                items = data.get("items", [])
+                try:
+                    data = req.json()
+                except Exception as exc:
+                    raise RuntimeError("TNM API returned invalid JSON") from exc
+                if not isinstance(data, dict):
+                    raise RuntimeError("TNM API returned an invalid response")
+                if data.get("errorMessage"):
+                    raise RuntimeError(f"TNM API error: {data['errorMessage']}")
+                if "total" not in data or "items" not in data:
+                    raise RuntimeError("TNM API response is missing total or items")
+                total = int(data["total"])
+                items = data["items"]
+                if total < 0 or not isinstance(items, list):
+                    raise RuntimeError("TNM API returned an invalid result page")
+                if expected_total is None:
+                    expected_total = total
+                elif total != expected_total:
+                    raise RuntimeError(
+                        "TNM API result total changed during pagination: "
+                        f"{expected_total} to {total}"
+                    )
+                if offset > total or len(items) > total - offset:
+                    raise RuntimeError("TNM API returned an invalid result page")
+                if offset < total and not items:
+                    raise RuntimeError(
+                        f"TNM API pagination stopped after {offset} of {total} products"
+                    )
 
                 for item in items:
+                    if not isinstance(item, dict):
+                        raise RuntimeError("TNM API returned an invalid product entry")
                     url = item.get("downloadURL")
-                    if not url:
-                        continue
+                    if not isinstance(url, str) or not url.strip():
+                        raise RuntimeError("TNM product is missing its download URL")
 
                     filename = url.split("/")[-1]
                     fmt = item.get("format", "Unknown")
@@ -279,6 +323,7 @@ def run(self):
                         "remote_size": item.get("sizeInBytes"),
                         "title": item.get("title"),
                         "tnm_project": project,
+                        "tnm_product": product,
                         "tnm_source_id": item.get("sourceId"),
                         "tnm_publication_date": item.get("publicationDate"),
                         "tnm_last_updated": item.get("lastUpdated"),
@@ -297,15 +342,39 @@ def run(self):
                         if date and date > existing_date:
                             best_tiles[fn_bn] = item_data
 
-            except Exception as e:
-                logger.exception(f"Error parsing TNM JSON: {e}")
-                break
+            except RuntimeError:
+                raise
+            except Exception as exc:
+                raise RuntimeError("Error parsing TNM API response") from exc
 
-            offset += 100
+            offset += len(items)
             if offset >= total:
                 break
 
-        tiles = best_tiles.values() if self.dedupe else all_tiles
+        tiles = list(best_tiles.values()) if self.dedupe else all_tiles
+        if self.source_coverage and tiles:
+            unsupported = set(dataset_names).difference(
+                {
+                    DATASET_CODES[2],
+                    DATASET_CODES[4],
+                    DATASET_CODES[6],
+                    DATASET_CODES[29],
+                }
+            )
+            if unsupported:
+                raise ValueError(
+                    "TNM source coverage supports only S1M, 1 m, 5 m and 1/9 arc-second DEMs"
+                )
+            if len(dataset_names) != 1:
+                raise ValueError("TNM source coverage requires one dataset per module")
+            if product in {"s1m", "5m"}:
+                tiles = tnm_raster.add_source_coverage(tiles, self.wgs_region)
+            elif product == "1_9as":
+                tiles = tnm_ned.add_source_coverage(tiles, self.wgs_region)
+            else:
+                tiles = WESM.add_source_coverage(
+                    tiles, self.wgs_region, require_year=True
+                )
         for tile_data in tiles:
             self.add_entry_to_results(**tile_data)