diff --git a/.gitignore b/.gitignore
index 9c2fa9a..42da6f2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,5 @@
 *.grb2
 *.grib2
 FILE:*
+__pycache__/
+downloads/
diff --git a/README.md b/README.md
index 4eddeca..1d7a7f8 100644
--- a/README.md
+++ b/README.md
@@ -24,15 +24,44 @@ pip install -r requirements.txt
 
 ---
 
-## 2. Retrieving input data
+## 2. Input data
+
+Two source datasets are needed: LANDFIRE fuel categories (NFUEL_CAT) and a
+high-resolution terrain DEM (ZSF). By default both are **downloaded
+automatically** for the exact WRF domain; supplying local GeoTIFF files with
+`--fuel` / `--zsf` is still supported.
+
+### Automatic download (default)
+
+If `--fuel` or `--zsf` is not given, the tool computes the domain's bounding
+box (plus a 2 km buffer) and fetches the data itself:
+
+- **NFUEL_CAT** — requested from the [LANDFIRE Product Service
+  (LFPS)](https://lfps.usgs.gov). The product matches `--fuel-table`
+  (`fbfm13` → FBFM13, `fbfm40`/`fbfm40_to_anderson13` → FBFM40), using the
+  newest full-coverage LANDFIRE version (pin one with `--landfire-version`,
+  e.g. `LF2023`). *LFPS requires an email address* (`--email`), used by
+  LANDFIRE only for usage reporting.
+- **ZSF** — from the USGS National Map (3DEP **1/3 arc-second**, ~10 m) by
+  default. The required 1°×1° tiles (~350 MB each) are downloaded and
+  merged/clipped to the domain. Alternatively `--zsf-source landfire` fetches
+  LANDFIRE's 30 m elevation through LFPS instead — a much smaller download at
+  the cost of resolution.
+
+Everything lands in `--download-dir` (default `./downloads`) and is cached:
+re-running the tool for the same domain reuses the existing files, and DEM
+tiles are shared between overlapping domains.
 
-Two datasets are required. Both are available as GeoTIFF and can be downloaded
-for any region of the contiguous US.
+```bash
+python fire_preprocess.py --wps-files geo_em.d03.nc --email you@example.org
+```
+
+### Manual download
 
-### Fuel categories (NFUEL_CAT)
+Both datasets are also available interactively as GeoTIFF for any region of
+the contiguous US.
 
-Download from the **LANDFIRE data viewer**:
-<https://landfire.gov/viewer/>
+**Fuel categories** — from the [LANDFIRE data viewer](https://landfire.gov/viewer/):
 
 1. Draw your area of interest on the map.
 2. Under *Fire Behavior Fuel Models*, select either:
@@ -44,10 +73,8 @@ Download from the **LANDFIRE data viewer**:
    rio merge LF2025_FBFM13_*.tif --output fuel.tif
    ```
 
-### High-resolution terrain (ZSF)
-
-Download from the **USGS National Map / 3D Elevation Program (3DEP)**:
-<https://apps.nationalmap.gov/downloader/>
+**High-resolution terrain** — from the [USGS National Map / 3D Elevation
+Program (3DEP)](https://apps.nationalmap.gov/downloader/):
 
 1. Under *Elevation Products (3DEP)*, select **1/3 Arc-Second DEM** (~10 m).
    1 Arc-Second (~30 m) is also available if coarser resolution is acceptable.
@@ -63,6 +90,15 @@ Download from the **USGS National Map / 3D Elevation Program (3DEP)**:
 
 ### Command-line usage
 
+Minimal — settings are read from the WPS file itself and both rasters are
+downloaded automatically:
+
+```bash
+python fire_preprocess.py --wps-files geo_em.d02.nc --email you@example.org
+```
+
+With local input files:
+
 ```bash
 python fire_preprocess.py \
     --wps-files geo_em.d02.nc \
@@ -79,9 +115,14 @@ All arguments can also be supplied via a YAML config file (see below).
 | Argument | Default | Description |
 |---|---|---|
 | `--wps-files` | — | WPS output file (`geo_em`/`met_em`) or glob pattern (e.g. `'met_em.d02.*.nc'`) |
-| `--fuel` | — | LANDFIRE fuel-category GeoTIFF (NFUEL_CAT source) |
-| `--zsf` | — | High-resolution terrain DEM GeoTIFF (ZSF source) |
-| `--namelist` | `namelist.wps` | Path to `namelist.wps` |
+| `--fuel` | *download* | LANDFIRE fuel-category GeoTIFF (NFUEL_CAT source); downloaded automatically if omitted |
+| `--zsf` | *download* | High-resolution terrain DEM GeoTIFF (ZSF source); downloaded automatically if omitted |
+| `--zsf-source` | `nationalmap` | Source for automatic ZSF download: `nationalmap` (USGS 3DEP 1/3 arc-sec, ~10 m) or `landfire` (30 m, much smaller download) |
+| `--email` | — | Email address; required by the LANDFIRE Product Service when downloading |
+| `--landfire-version` | newest full-coverage | Pin the LANDFIRE version of downloaded fuel data (e.g. `LF2023`) |
+| `--download-dir` | `downloads` | Cache directory for downloaded rasters |
+| `--zsf-fill` | `0` | Value assigned to ZSF pixels with no source elevation data (outside the DEM's extent, or nodata voids within it) |
+| `--namelist` | `namelist.wps` | Path to `namelist.wps`. Optional: if absent, `subgrid_ratio_x/y` are read from the WPS file's `sr_x`/`sr_y` global attributes (a namelist value overrides the file) |
 | `--domain` | `1` | Domain number, used to read the correct `subgrid_ratio_x/y` from the namelist |
 | `--fuel-table` | `fbfm13` | Fuel remapping table (see below) |
 | `--overwrite` | `false` | Overwrite existing fire fields without prompting |
@@ -93,12 +134,16 @@ Any argument can be set in a YAML config file. CLI flags override config values.
 
 ```yaml
 wps_files:  'met_em.d02.2024-09-08_*.nc'
-zsf:        /path/to/dem.tif
-fuel:       /path/to/fuel.tif
+zsf:        /path/to/dem.tif      # omit to download automatically
+fuel:       /path/to/fuel.tif     # omit to download automatically
 namelist:   namelist.wps
 fuel_table: fbfm13
 domain:     2
 overwrite:  false
+zsf_source: nationalmap
+email:      you@example.org       # required for LANDFIRE downloads
+download_dir: downloads
+zsf_fill:   0                     # value for missing elevation pixels
 ```
 
 Run with a config file:
diff --git a/environment.yaml b/environment.yaml
index 648ab42..b2fd1ea 100644
--- a/environment.yaml
+++ b/environment.yaml
@@ -10,3 +10,4 @@ dependencies:
   - pyproj>=3.7.2
   - f90nml>=1.5
   - pyyaml>=6.0
+  - requests>=2.28
diff --git a/fire_preprocess/cli.py b/fire_preprocess/cli.py
index 9490199..9bc69e9 100644
--- a/fire_preprocess/cli.py
+++ b/fire_preprocess/cli.py
@@ -7,13 +7,13 @@
 import yaml
 
 from .namelist import get_fire_subgrid_ratios, get_domain_params
-from .wrf_domain import read_domain_from_file
+from .wrf_domain import read_domain_from_file, read_subgrid_ratios_from_file
 from .fire_grid import build_fire_grid
 from .raster import reproject_fuel, reproject_dem
 from .fuel_tables import get_fuel_table, list_fuel_tables
 from .wps_io import check_existing_fire_vars, prompt_overwrite, write_fire_vars
 
-_REQUIRED = ("wps_files", "zsf", "fuel", "namelist")
+_REQUIRED = ("wps_files",)
 
 # Defaults applied after CLI + config are merged, so neither source
 # can be mistaken for an explicit user value.
@@ -22,6 +22,16 @@
     "domain": 1,
     "namelist": "namelist.wps",
     "overwrite": False,
+    "zsf_source": "nationalmap",
+    "download_dir": "downloads",
+    "zsf_fill": 0.0,
+}
+
+# LANDFIRE product to download when --fuel is omitted, keyed by fuel table.
+_FUEL_TABLE_PRODUCT = {
+    "fbfm13": "FBFM13",
+    "fbfm40": "FBFM40",
+    "fbfm40_to_anderson13": "FBFM40",
 }
 
 
@@ -66,12 +76,16 @@ def _merge(args: argparse.Namespace, cfg: dict) -> argparse.Namespace:
 def build_parser() -> argparse.ArgumentParser:
     config_example = (
         "  wps_files:  'met_em.d01.*.nc'\n"
-        "  zsf:        /path/to/highres_dem.tif\n"
-        "  fuel:       /path/to/landfire.tif\n"
+        "  zsf:        /path/to/highres_dem.tif   # omit to download automatically\n"
+        "  fuel:       /path/to/landfire.tif      # omit to download automatically\n"
         "  namelist:   namelist.wps\n"
         "  fuel_table: fbfm13\n"
         "  domain:     1\n"
         "  overwrite:  false\n"
+        "  zsf_source: nationalmap\n"
+        "  email:      you@example.org           # required for LANDFIRE downloads\n"
+        "  download_dir: downloads\n"
+        "  zsf_fill:   0                         # value for missing elevation pixels\n"
     )
     parser = argparse.ArgumentParser(
         prog="fire_preprocess",
@@ -102,15 +116,53 @@ def build_parser() -> argparse.ArgumentParser:
     )
     parser.add_argument(
         "--zsf", default=None, metavar="GEOTIFF",
-        help="High-resolution terrain DEM GeoTIFF (ZSF source; ≥1/3 arc-sec recommended)",
+        help=(
+            "High-resolution terrain DEM GeoTIFF (ZSF source; ≥1/3 arc-sec recommended). "
+            "If omitted, the DEM is downloaded automatically (see --zsf-source)."
+        ),
     )
     parser.add_argument(
         "--fuel", default=None, metavar="GEOTIFF",
-        help="LANDFIRE fuel-category GeoTIFF (NFUEL_CAT source)",
+        help=(
+            "LANDFIRE fuel-category GeoTIFF (NFUEL_CAT source). If omitted, the fuel "
+            "layer matching --fuel-table is downloaded automatically from LANDFIRE."
+        ),
     )
     parser.add_argument(
         "--namelist", default=None, metavar="FILE",
-        help="Path to namelist.wps (provides subgrid_ratio_x/y and projection parameters)",
+        help=(
+            "Path to namelist.wps. Optional: subgrid_ratio_x/y are read from the WPS "
+            "file's sr_x/sr_y global attributes when the namelist is absent."
+        ),
+    )
+    parser.add_argument(
+        "--zsf-source", default=None, choices=("nationalmap", "landfire"), dest="zsf_source",
+        help=(
+            "Source for automatic ZSF download: 'nationalmap' = USGS 3DEP 1/3 arc-sec "
+            "(~10 m, default), 'landfire' = LANDFIRE elevation (30 m, smaller download)"
+        ),
+    )
+    parser.add_argument(
+        "--email", default=None, metavar="ADDRESS",
+        help="Email address, required by the LANDFIRE Product Service for downloads",
+    )
+    parser.add_argument(
+        "--landfire-version", default=None, metavar="LFxxxx", dest="landfire_version",
+        help=(
+            "Pin the LANDFIRE version for downloaded fuel data (e.g. LF2023). "
+            "Default: newest version with full geographic coverage"
+        ),
+    )
+    parser.add_argument(
+        "--download-dir", default=None, metavar="DIR", dest="download_dir",
+        help="Directory where downloaded rasters are cached. Default: ./downloads",
+    )
+    parser.add_argument(
+        "--zsf-fill", type=float, default=None, metavar="VALUE", dest="zsf_fill",
+        help=(
+            "Value assigned to ZSF pixels with no source elevation data "
+            "(outside the DEM's extent, or nodata voids within it). Default: 0"
+        ),
     )
     parser.add_argument(
         "--fuel-table", default=None, metavar="NAME|PATH", dest="fuel_table",
@@ -145,6 +197,8 @@ def main(argv=None):
             f"The following arguments are required: {', '.join(missing)}\n"
             "Supply them on the command line or via --config."
         )
+    if args.zsf_source not in ("nationalmap", "landfire"):
+        parser.error(f"zsf_source must be 'nationalmap' or 'landfire', got '{args.zsf_source}'")
 
     # ── Find WPS files ─────────────────────────────────────────────────────
     print(f"Searching for WPS file(s): {args.wps_files}")
@@ -154,12 +208,36 @@ def main(argv=None):
     basenames = [os.path.basename(f) for f in files]
     print(f"  Found {len(files)} file(s): {basenames}")
 
-    # ── Read namelist ─────────────────────────────────────────────────────────
+    # ── Determine fire subgrid ratios and domain parameters ──────────────────
+    # An explicit namelist value wins (e.g. adding a fire grid to a geo_em file
+    # generated without one); otherwise the WPS file's sr_x/sr_y attributes are used.
     domain_index = args.domain - 1  # namelist arrays are 0-indexed, domain numbers are 1-indexed
-    print(f"Reading namelist: {args.namelist}")
-    sr_x, sr_y = get_fire_subgrid_ratios(args.namelist, domain_index=domain_index)
-    print(f"  subgrid_ratio_x = {sr_x},  subgrid_ratio_y = {sr_y}")
-    nml_params = get_domain_params(args.namelist, domain_index=domain_index)
+    file_sr = read_subgrid_ratios_from_file(files[0])
+    nml_sr = (None, None)
+    nml_params = None
+    if os.path.isfile(args.namelist):
+        print(f"Reading namelist: {args.namelist}")
+        nml_sr = get_fire_subgrid_ratios(args.namelist, domain_index=domain_index)
+        nml_params = get_domain_params(args.namelist, domain_index=domain_index)
+
+    if nml_sr[0] is not None:
+        sr_x, sr_y = nml_sr
+        sr_source = args.namelist
+        if file_sr[0] is not None and file_sr != nml_sr:
+            print(
+                f"  Note: {basenames[0]} has sr_x/sr_y = {file_sr[0]}/{file_sr[1]}; "
+                f"using the namelist values {sr_x}/{sr_y} instead."
+            )
+    elif file_sr[0] is not None:
+        sr_x, sr_y = file_sr
+        sr_source = f"{basenames[0]} global attributes"
+    else:
+        parser.error(
+            f"Fire subgrid ratios not found: '{basenames[0]}' has no sr_x/sr_y global "
+            f"attributes and no namelist with subgrid_ratio_x/y was found at "
+            f"'{args.namelist}'. Supply --namelist."
+        )
+    print(f"  subgrid_ratio_x = {sr_x},  subgrid_ratio_y = {sr_y}  (from {sr_source})")
 
     # ── Build WRF domain geometry ─────────────────────────────────────────────
     print(f"Reading domain geometry from {basenames[0]} ...")
@@ -179,6 +257,51 @@ def main(argv=None):
     fuel_table = get_fuel_table(args.fuel_table)
     print(f"Fuel table: {fuel_table.name}  —  {fuel_table.description}")
 
+    # ── Download any missing source rasters ───────────────────────────────────
+    if args.fuel is None or args.zsf is None:
+        # Imported lazily so environments without `requests` can still run
+        # with locally supplied rasters.
+        from .download import (
+            fire_grid_bbox_wgs84, bbox_tag, resolve_landfire_layer,
+            download_landfire, download_usgs_dem,
+        )
+        bbox = fire_grid_bbox_wgs84(fire_grid)
+        tag = bbox_tag(bbox)
+        os.makedirs(args.download_dir, exist_ok=True)
+        print(
+            f"Download area (WGS84): W={bbox[0]:.4f}  S={bbox[1]:.4f}  "
+            f"E={bbox[2]:.4f}  N={bbox[3]:.4f}"
+        )
+
+        if args.fuel is None:
+            product = _FUEL_TABLE_PRODUCT.get(fuel_table.name)
+            if product is None:
+                parser.error(
+                    "--fuel is required when using a custom fuel table "
+                    "(cannot infer which LANDFIRE product to download)."
+                )
+            print(f"No fuel raster given; downloading {product} from LANDFIRE ...")
+            layer = resolve_landfire_layer(product, args.landfire_version)
+            args.fuel = download_landfire(
+                layer, bbox, os.path.join(args.download_dir, f"{layer}_{tag}.tif"),
+                args.email,
+            )
+
+        if args.zsf is None:
+            if args.zsf_source == "landfire":
+                print("No terrain DEM given; downloading elevation from LANDFIRE ...")
+                layer = resolve_landfire_layer("Elev")
+                args.zsf = download_landfire(
+                    layer, bbox, os.path.join(args.download_dir, f"{layer}_{tag}.tif"),
+                    args.email,
+                )
+            else:
+                print("No terrain DEM given; downloading 1/3 arc-second DEM from the USGS National Map ...")
+                args.zsf = download_usgs_dem(
+                    bbox, args.download_dir,
+                    os.path.join(args.download_dir, f"USGS_3DEP_13as_{tag}.tif"),
+                )
+
     # ── Reproject fuel categories ─────────────────────────────────────────────
     print(f"Reprojecting fuel data: {args.fuel}")
     nfuel_raw = reproject_fuel(args.fuel, fire_grid)
@@ -190,7 +313,7 @@ def main(argv=None):
 
     # ── Reproject DEM ─────────────────────────────────────────────────────────
     print(f"Reprojecting terrain DEM: {args.zsf}")
-    zsf = reproject_dem(args.zsf, fire_grid)
+    zsf = reproject_dem(args.zsf, fire_grid, fill_value=args.zsf_fill)
     print(
         f"  ZSF shape {zsf.shape}  "
         f"range [{zsf.min():.1f}, {zsf.max():.1f}] m"
diff --git a/fire_preprocess/download.py b/fire_preprocess/download.py
new file mode 100644
index 0000000..51ed048
--- /dev/null
+++ b/fire_preprocess/download.py
@@ -0,0 +1,324 @@
+"""Download source rasters for NFUEL_CAT and ZSF from public web services.
+
+Fuel categories (NFUEL_CAT) come from the LANDFIRE Product Service (LFPS):
+  https://lfps.usgs.gov  — REST API; asynchronous jobs return a zipped GeoTIFF
+  clipped to the requested bounding box. LFPS requires an email address with
+  every request (used by LANDFIRE for usage reporting only).
+
+Terrain (ZSF) comes from either:
+  - USGS National Map / 3DEP via the TNM Access API (default; 1/3 arc-second,
+    ~10 m): https://tnmaccess.nationalmap.gov/api/v1/products
+    Tiles are 1°x1° GeoTIFFs (~350 MB each); they are downloaded, cached, and
+    merged/clipped to the domain.
+  - LANDFIRE elevation (LF2020_Elev, 30 m) through the same LFPS API — much
+    smaller download, but coarser than 3DEP.
+
+All downloads land in a cache directory; a file that already exists there
+(same layer + bounding box) is reused instead of re-downloaded.
+"""
+import json
+import math
+import os
+import re
+import time
+import warnings
+import zipfile
+
+import numpy as np
+import requests
+import rasterio
+from rasterio.merge import merge as rio_merge
+from pyproj import CRS, Transformer
+
+from .fire_grid import FireGrid
+from .wrf_domain import WRF_SPHERE_RADIUS
+
+LFPS_API = "https://lfps.usgs.gov/api"
+TNM_API = "https://tnmaccess.nationalmap.gov/api/v1/products"
+TNM_DEM_DATASET = "National Elevation Dataset (NED) 1/3 arc-second"
+
+# Used only if the live LFPS product listing cannot be queried.
+FALLBACK_LAYERS = {
+    "FBFM13": "LF2024_FBFM13",
+    "FBFM40": "LF2024_FBFM40",
+    "Elev": "LF2020_Elev",
+}
+
+_CHUNK = 1 << 20  # 1 MiB
+
+
+# ── Bounding box ──────────────────────────────────────────────────────────────
+
+def fire_grid_bbox_wgs84(fire_grid: FireGrid, buffer_m: float = 2000.0):
+    """Return (west, south, east, north) lat/lon bounds covering the fire grid.
+
+    Samples points along the grid perimeter (projected edges are curved in
+    lat/lon) and pads by *buffer_m* so edge pixels are fully covered by the
+    downloaded data. WRF's spherical lat/lon values are used directly as
+    WGS84, following the standard WRF convention of ignoring the
+    sphere-vs-ellipsoid datum shift.
+    """
+    t = fire_grid.transform
+    x0, y_north = t.c, t.f
+    x1 = x0 + t.a * fire_grid.nx_fire
+    y_south = y_north + t.e * fire_grid.ny_fire
+
+    n = 50
+    xs = np.linspace(x0, x1, n)
+    ys = np.linspace(y_south, y_north, n)
+    edge_x = np.concatenate([xs, xs, np.full(n, x0), np.full(n, x1)])
+    edge_y = np.concatenate([np.full(n, y_south), np.full(n, y_north), ys, ys])
+
+    geo = CRS.from_proj4(
+        f"+proj=longlat +a={WRF_SPHERE_RADIUS} +b={WRF_SPHERE_RADIUS} +no_defs"
+    )
+    lon, lat = Transformer.from_crs(fire_grid.crs, geo, always_xy=True).transform(edge_x, edge_y)
+
+    blat = buffer_m / 111320.0
+    blon = buffer_m / (111320.0 * math.cos(math.radians(float(np.abs(lat).max()))))
+    return (
+        float(lon.min()) - blon,
+        float(lat.min()) - blat,
+        float(lon.max()) + blon,
+        float(lat.max()) + blat,
+    )
+
+
+def bbox_tag(bbox) -> str:
+    """Filename-friendly tag identifying a bounding box (used for cache names)."""
+    return "_".join(f"{v:.4f}" for v in bbox)
+
+
+def _get_json(url: str, params: dict = None, timeout: int = 60, retries: int = 4):
+    """GET a JSON endpoint, retrying transient failures (5xx, timeouts).
+
+    Both LFPS and the TNM Access API intermittently return gateway errors;
+    a few retries with increasing back-off rides those out.
+    """
+    delay = 5
+    for attempt in range(retries):
+        try:
+            r = requests.get(url, params=params, timeout=timeout)
+            if r.status_code >= 500:
+                raise requests.HTTPError(f"{r.status_code} for {r.url}", response=r)
+            r.raise_for_status()
+            return r.json()
+        except (requests.ConnectionError, requests.Timeout, requests.HTTPError) as exc:
+            status = getattr(getattr(exc, "response", None), "status_code", None)
+            if status is not None and status < 500:
+                raise  # 4xx: our request is wrong; retrying won't help
+            if attempt == retries - 1:
+                raise
+            print(f"    Transient error from {url} ({exc}); retrying in {delay} s ...")
+            time.sleep(delay)
+            delay *= 2
+
+
+def _stream_download(url: str, dest: str, label: str = None, size: int = None):
+    """Download *url* to *dest*, atomically via a .part file."""
+    label = label or os.path.basename(dest)
+    size_txt = f" ({size / 1e6:.0f} MB)" if size else ""
+    print(f"  Downloading {label}{size_txt} ...", flush=True)
+    tmp = dest + ".part"
+    with requests.get(url, stream=True, timeout=120) as r:
+        r.raise_for_status()
+        with open(tmp, "wb") as fh:
+            for chunk in r.iter_content(_CHUNK):
+                fh.write(chunk)
+    os.replace(tmp, dest)
+
+
+# ── LANDFIRE Product Service ──────────────────────────────────────────────────
+
+def resolve_landfire_layer(product: str, version: str = None) -> str:
+    """Return the LFPS layer name (e.g. 'LF2024_FBFM13') for a product acronym.
+
+    Queries the live LFPS product listing. *version* (e.g. 'LF2023') pins a
+    specific LANDFIRE version; otherwise the newest version with full ('All')
+    geographic coverage is preferred, because the newest release sometimes
+    covers only part of CONUS (e.g. LF2025 covers only SW/NW GeoAreas).
+    """
+    try:
+        products = _get_json(f"{LFPS_API}/products", timeout=30)["products"]
+    except Exception as exc:
+        layer = FALLBACK_LAYERS[product]
+        warnings.warn(f"Could not query the LFPS product list ({exc}); falling back to '{layer}'.")
+        return layer
+
+    # Match on theme too: excludes 'Seasonal Fuels' variants like LF2025_FBFM40_SP26.
+    cands = [
+        p for p in products
+        if p.get("acronym") == product and p.get("theme") in ("Fuels", "Topographic")
+    ]
+    if not cands:
+        raise ValueError(f"No LFPS product found with acronym '{product}'.")
+
+    if version:
+        for p in cands:
+            if p["version"].lower() == version.lower():
+                return p["layerName"]
+        avail = sorted({p["version"] for p in cands})
+        raise ValueError(
+            f"LFPS has no {product} product for version '{version}'. "
+            f"Available versions: {', '.join(avail)}"
+        )
+
+    full_coverage = [p for p in cands if p.get("geoAreas") == "All"]
+    best = max(full_coverage or cands, key=lambda p: p["version"])
+    return best["layerName"]
+
+
+def download_landfire(layer: str, bbox, out_path: str, email: str,
+                      poll_seconds: int = 10, timeout_seconds: int = 1800) -> str:
+    """Fetch *layer* clipped to *bbox* via an asynchronous LFPS job.
+
+    Returns *out_path* (reused as-is if it already exists).
+    """
+    if os.path.isfile(out_path):
+        print(f"  Using cached download: {out_path}")
+        return out_path
+    if not email:
+        raise ValueError(
+            "LANDFIRE downloads require an email address (LFPS 'Email' parameter). "
+            "Supply one with --email or the 'email' config key."
+        )
+
+    west, south, east, north = bbox
+    params = {
+        "Email": email,
+        "Layer_List": layer,
+        "Area_of_Interest": f"{west} {south} {east} {north}",
+    }
+    job = _get_json(f"{LFPS_API}/job/submit", params=params)
+    job_id = job.get("jobId")
+    if not job_id:
+        raise RuntimeError(f"LFPS job submission failed: {job}")
+    print(f"  LFPS job submitted: {job_id}  (layer {layer})")
+
+    deadline = time.time() + timeout_seconds
+    while True:
+        status_json = _get_json(f"{LFPS_API}/job/status", params={"JobId": job_id})
+        status = status_json.get("status", "")
+        if status == "Succeeded":
+            break
+        if status in ("Failed", "Canceled"):
+            msgs = status_json.get("messages", [])
+            raise RuntimeError(f"LFPS job {job_id} {status.lower()}: {msgs}")
+        if time.time() > deadline:
+            raise TimeoutError(
+                f"LFPS job {job_id} did not finish within {timeout_seconds} s "
+                f"(last status: {status})."
+            )
+        pos = status_json.get("queuePosition")
+        queue_txt = f", queue position {pos}" if pos not in (None, -1) else ""
+        print(f"    LFPS job status: {status}{queue_txt} — retrying in {poll_seconds} s")
+        time.sleep(poll_seconds)
+
+    m = re.search(r'https?://[^"\s\\]+?\.zip', json.dumps(status_json))
+    if not m:
+        raise RuntimeError(
+            f"LFPS job {job_id} succeeded but no download URL was found "
+            f"in the status response: {status_json}"
+        )
+
+    zip_path = out_path + ".zip"
+    _stream_download(m.group(0), zip_path, label=f"LFPS bundle {job_id}.zip")
+    with zipfile.ZipFile(zip_path) as zf:
+        tifs = [n for n in zf.namelist() if n.lower().endswith(".tif")]
+        if not tifs:
+            raise RuntimeError(f"No GeoTIFF found in the LFPS bundle {zip_path}")
+        with zf.open(tifs[0]) as src, open(out_path, "wb") as dst:
+            while True:
+                chunk = src.read(_CHUNK)
+                if not chunk:
+                    break
+                dst.write(chunk)
+    os.remove(zip_path)
+    print(f"  Saved {out_path}")
+    return out_path
+
+
+# ── USGS National Map (3DEP) ──────────────────────────────────────────────────
+
+def _tnm_query(bbox, dataset: str) -> list:
+    """Return all TNM product items for *dataset* intersecting *bbox*."""
+    west, south, east, north = bbox
+    items, offset = [], 0
+    while True:
+        params = {
+            "datasets": dataset,
+            "bbox": f"{west},{south},{east},{north}",
+            "prodFormats": "GeoTIFF",
+            "outputFormat": "JSON",
+            "max": 100,
+            "offset": offset,
+        }
+        payload = _get_json(TNM_API, params=params, timeout=120)
+        batch = payload.get("items", [])
+        items.extend(batch)
+        offset += len(batch)
+        if not batch or offset >= int(payload.get("total", 0)):
+            return items
+
+
+def _latest_per_tile(items: list) -> list:
+    """Keep only the most recent product per 1°x1° tile (drops superseded versions)."""
+    tiles = {}
+    for it in items:
+        url = it.get("downloadURL") or ""
+        if not url.lower().endswith(".tif"):
+            continue
+        m = re.search(r"[ns]\d{2,3}[ew]\d{3}", url)
+        key = m.group(0) if m else url
+        prev = tiles.get(key)
+        if prev is None or (it.get("publicationDate") or "") > (prev.get("publicationDate") or ""):
+            tiles[key] = it
+    return list(tiles.values())
+
+
+def download_usgs_dem(bbox, download_dir: str, out_path: str,
+                      dataset: str = TNM_DEM_DATASET) -> str:
+    """Download the 3DEP DEM tiles covering *bbox* and merge/clip into *out_path*.
+
+    Returns *out_path* (reused as-is if it already exists). Individual tiles
+    are cached in *download_dir* under their original names.
+    """
+    if os.path.isfile(out_path):
+        print(f"  Using cached download: {out_path}")
+        return out_path
+
+    print(f"  Querying the USGS National Map: {dataset}")
+    tiles = _latest_per_tile(_tnm_query(bbox, dataset))
+    if not tiles:
+        raise RuntimeError(
+            f"The USGS National Map returned no '{dataset}' products for bbox {bbox}. "
+            "Check that the domain is inside 3DEP coverage, or use --zsf-source landfire."
+        )
+    print(f"  {len(tiles)} DEM tile(s) cover the domain")
+
+    tile_paths = []
+    for it in tiles:
+        url = it["downloadURL"]
+        dest = os.path.join(download_dir, os.path.basename(url))
+        if os.path.isfile(dest):
+            print(f"  Using cached tile: {os.path.basename(dest)}")
+        else:
+            _stream_download(url, dest, size=it.get("sizeInBytes"))
+        tile_paths.append(dest)
+
+    print(f"  Merging {len(tile_paths)} tile(s), clipped to the domain")
+    sources = [rasterio.open(p) for p in tile_paths]
+    try:
+        data, transform = rio_merge(sources, bounds=bbox)
+        profile = sources[0].profile.copy()
+        profile.update(
+            height=data.shape[1], width=data.shape[2],
+            transform=transform, driver="GTiff", compress="deflate",
+        )
+        with rasterio.open(out_path, "w", **profile) as dst:
+            dst.write(data)
+    finally:
+        for s in sources:
+            s.close()
+    print(f"  Saved {out_path}")
+    return out_path
diff --git a/fire_preprocess/namelist.py b/fire_preprocess/namelist.py
index 07ebe73..423be55 100644
--- a/fire_preprocess/namelist.py
+++ b/fire_preprocess/namelist.py
@@ -26,12 +26,16 @@ def get_fire_subgrid_ratios(path, domain_index=0):
         domain_index: 0-based domain index (default 0 = first domain)
 
     Returns:
-        (sr_x, sr_y) as integers; defaults to (1, 1) if not present
+        (sr_x, sr_y) as integers, or (None, None) if the namelist does not
+        set subgrid_ratio_x/y (the caller may then fall back to the sr_x/sr_y
+        global attributes of the WPS file).
     """
     params = read_namelist_wps(path)
-    sr_x = int(_scalar(params.get("subgrid_ratio_x", 1), domain_index))
-    sr_y = int(_scalar(params.get("subgrid_ratio_y", 1), domain_index))
-    return sr_x, sr_y
+    sr_x = params.get("subgrid_ratio_x")
+    sr_y = params.get("subgrid_ratio_y")
+    if sr_x is None or sr_y is None:
+        return None, None
+    return int(_scalar(sr_x, domain_index)), int(_scalar(sr_y, domain_index))
 
 
 def get_domain_params(path, domain_index=0):
diff --git a/fire_preprocess/raster.py b/fire_preprocess/raster.py
index ef636d6..b87b363 100644
--- a/fire_preprocess/raster.py
+++ b/fire_preprocess/raster.py
@@ -17,12 +17,9 @@
 from .fire_grid import FireGrid
 
 
-def _warn_coverage(arr, nodata_val, field_name, threshold=0.05):
+def _warn_coverage(missing_mask, field_name, threshold=0.05):
     """Emit a warning if more than *threshold* fraction of pixels are nodata."""
-    if nodata_val is None:
-        return
-    mask = (arr == nodata_val) | np.isnan(arr)
-    frac = mask.sum() / arr.size
+    frac = missing_mask.sum() / missing_mask.size
     if frac > threshold:
         warnings.warn(
             f"{field_name}: {frac:.1%} of fire-grid pixels are nodata — "
@@ -37,19 +34,29 @@ def reproject_to_grid(
     resampling: Resampling,
     band: int = 1,
     src_nodata: float = None,
+    dst_fill: float = None,
     field_name: str = "field",
 ) -> np.ndarray:
     """Core reprojection routine.
 
+    Pixels not covered by the source raster (domain edges beyond its
+    extent) or that were themselves nodata in the source end up at
+    *dst_fill*. If *dst_fill* is None, it defaults to the source's own
+    nodata value (or 0 if the source declares none).
+
+    Reprojection itself always fills missing pixels with NaN rather than
+    *dst_fill* directly: rasterio/GDAL silently ignores a dst_nodata of
+    exactly 0.0 (falsy in the underlying binding) and falls back to the
+    source's nodata value, which would leak it through whenever a caller
+    wants a fill of 0. Substituting *dst_fill* ourselves afterwards avoids
+    that pitfall for any fill value, including 0.
+
     Returns a 2-D float32 array in WRF row order (south_north, west_east):
     row 0 is the *southernmost* row, matching WRF/WPS netCDF convention.
     """
     height, width = fire_grid.ny_fire, fire_grid.nx_fire
     dst_crs = RasterioCRS.from_user_input(fire_grid.crs.to_wkt())
-
-    # np.zeros ensures any pixels rasterio does not write (e.g. at domain
-    # edges outside the source extent) get a known value rather than garbage.
-    dst_data = np.zeros((height, width), dtype=np.float32)
+    dst_data = np.full((height, width), np.nan, dtype=np.float32)
 
     with rasterio.open(src_path) as src:
         nodata = src_nodata if src_nodata is not None else src.nodata
@@ -63,13 +70,18 @@ def reproject_to_grid(
             src_nodata=nodata,
             dst_transform=fire_grid.transform,
             dst_crs=dst_crs,
+            dst_nodata=np.nan,
             resampling=resampling,
         )
 
     # rasterio returns row 0 = north; WRF netCDF expects row 0 = south.
     dst_data = np.flipud(dst_data)
 
-    _warn_coverage(dst_data, nodata, field_name)
+    missing = np.isnan(dst_data)
+    fill = dst_fill if dst_fill is not None else (nodata if nodata is not None else 0.0)
+    dst_data = np.where(missing, np.float32(fill), dst_data)
+
+    _warn_coverage(missing, field_name)
     return dst_data
 
 
@@ -77,7 +89,9 @@ def reproject_fuel(src_path: str, fire_grid: FireGrid, band: int = 1) -> np.ndar
     """Reproject a fuel-category GeoTIFF to the fire grid.
 
     Uses nearest-neighbour resampling to preserve integer category values.
-    LANDFIRE categorical data must not be interpolated.
+    LANDFIRE categorical data must not be interpolated. Missing pixels are
+    left at the source's own nodata value, which the fuel table's
+    nodata_values remap to nodata_out (see fuel_tables/base.py).
     """
     return reproject_to_grid(
         src_path, fire_grid,
@@ -87,15 +101,19 @@ def reproject_fuel(src_path: str, fire_grid: FireGrid, band: int = 1) -> np.ndar
     )
 
 
-def reproject_dem(src_path: str, fire_grid: FireGrid, band: int = 1) -> np.ndarray:
+def reproject_dem(src_path: str, fire_grid: FireGrid, band: int = 1,
+                   fill_value: float = 0.0) -> np.ndarray:
     """Reproject a terrain DEM GeoTIFF to the fire grid.
 
     Uses bilinear resampling, matching the WPS geogrid four_pt interpolation
-    specified for ZSF in GEOGRID.TBL.FIRE.
+    specified for ZSF in GEOGRID.TBL.FIRE. Pixels with no source elevation
+    data (outside the DEM's extent, or nodata voids within it) are set to
+    *fill_value*.
     """
     return reproject_to_grid(
         src_path, fire_grid,
         resampling=Resampling.bilinear,
         band=band,
+        dst_fill=fill_value,
         field_name="ZSF",
     )
diff --git a/fire_preprocess/wrf_domain.py b/fire_preprocess/wrf_domain.py
index ded3a4f..f45e122 100644
--- a/fire_preprocess/wrf_domain.py
+++ b/fire_preprocess/wrf_domain.py
@@ -92,6 +92,20 @@ def _sw_corner_from_namelist(params, crs):
     return x_sw_mass, y_sw_mass
 
 
+def read_subgrid_ratios_from_file(file_path):
+    """Return (sr_x, sr_y) from a WPS file's global attributes, or (None, None).
+
+    geogrid writes the per-domain subgrid_ratio_x/y values to geo_em files as
+    the global attributes sr_x/sr_y (and metgrid copies them to met_em files),
+    so a namelist.wps is not needed when these attributes are present.
+    """
+    with nc.Dataset(file_path) as ds:
+        try:
+            return int(ds.getncattr("sr_x")), int(ds.getncattr("sr_y"))
+        except AttributeError:
+            return None, None
+
+
 def read_domain_from_file(file_path, sr_x, sr_y, namelist_params=None):
     """Build a WRFDomain by reading a WPS netCDF file.
 
diff --git a/requirements.txt b/requirements.txt
index 6525888..72a4d11 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -4,3 +4,4 @@ rasterio>=1.5.0
 pyproj>=3.7.2
 f90nml>=1.5
 pyyaml>=6.0
+requests>=2.28