From eaacf899e095d74f5c3a51994ee15cd26dd4ad90 Mon Sep 17 00:00:00 2001 From: "Michael Kavulich, Jr" Date: Tue, 14 Jul 2026 15:36:33 -0600 Subject: [PATCH 1/3] Read subgrid ratios from WPS file attributes, making namelist.wps optional MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit geogrid writes the per-domain subgrid_ratio_x/y to geo_em files as the sr_x/sr_y global attributes (copied to met_em by metgrid), so a namelist.wps is only needed to override them — e.g. when adding a fire grid to a domain that was originally run without one. An explicit namelist value takes precedence, with a printed note on mismatch. Also ignore __pycache__ directories. Co-Authored-By: Claude Fable 5 --- .gitignore | 1 + README.md | 2 +- fire_preprocess/cli.py | 44 ++++++++++++++++++++++++++++------- fire_preprocess/namelist.py | 12 ++++++---- fire_preprocess/wrf_domain.py | 14 +++++++++++ 5 files changed, 59 insertions(+), 14 deletions(-) diff --git a/.gitignore b/.gitignore index 9c2fa9a..80a0a29 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ *.grb2 *.grib2 FILE:* +__pycache__/ diff --git a/README.md b/README.md index 4eddeca..665d950 100644 --- a/README.md +++ b/README.md @@ -81,7 +81,7 @@ All arguments can also be supplied via a YAML config file (see below). | `--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` | +| `--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 | diff --git a/fire_preprocess/cli.py b/fire_preprocess/cli.py index 9490199..a9e38cb 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", "zsf", "fuel") # Defaults applied after CLI + config are merged, so neither source # can be mistaken for an explicit user value. @@ -110,7 +110,10 @@ def build_parser() -> argparse.ArgumentParser: ) 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( "--fuel-table", default=None, metavar="NAME|PATH", dest="fuel_table", @@ -145,7 +148,6 @@ def main(argv=None): f"The following arguments are required: {', '.join(missing)}\n" "Supply them on the command line or via --config." ) - # ── Find WPS files ───────────────────────────────────────────────────── print(f"Searching for WPS file(s): {args.wps_files}") files = sorted(glob.glob(args.wps_files)) @@ -154,12 +156,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]} ...") 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/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. From bc817e687b2e2fe1e338363130b1ac1c8167cce6 Mon Sep 17 00:00:00 2001 From: "Michael Kavulich, Jr" Date: Tue, 14 Jul 2026 15:37:09 -0600 Subject: [PATCH 2/3] Add automatic download of fuel and terrain data from LANDFIRE and USGS When --fuel / --zsf are omitted, the source rasters are fetched for the domain's bounding box (fire-grid perimeter + 2 km buffer): - NFUEL_CAT via the LANDFIRE Product Service (lfps.usgs.gov) REST API, choosing the FBFM13/FBFM40 layer to match --fuel-table and the newest full-coverage LANDFIRE version (--landfire-version pins one). LFPS requires an email address (--email). - ZSF from the USGS National Map 3DEP 1/3 arc-second DEM by default (tiles are downloaded, deduplicated to the latest revision, and merged/clipped), or from LANDFIRE 30 m elevation with --zsf-source landfire for a much smaller download. Downloads are cached in --download-dir (default ./downloads) and reused on re-runs. API queries retry transient 5xx errors, which the TNM service returns routinely. Adds the requests dependency. Co-Authored-By: Claude Fable 5 --- .gitignore | 1 + README.md | 71 ++++++-- environment.yaml | 1 + fire_preprocess/cli.py | 98 ++++++++++- fire_preprocess/download.py | 324 ++++++++++++++++++++++++++++++++++++ requirements.txt | 1 + 6 files changed, 477 insertions(+), 19 deletions(-) create mode 100644 fire_preprocess/download.py diff --git a/.gitignore b/.gitignore index 80a0a29..42da6f2 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ *.grib2 FILE:* __pycache__/ +downloads/ diff --git a/README.md b/README.md index 665d950..6dee529 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**: - +**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)**: - +**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,8 +115,12 @@ 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) | +| `--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 | | `--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) | @@ -93,12 +133,15 @@ 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 ``` 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 a9e38cb..b3924bf 100644 --- a/fire_preprocess/cli.py +++ b/fire_preprocess/cli.py @@ -13,7 +13,7 @@ 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") +_REQUIRED = ("wps_files",) # Defaults applied after CLI + config are merged, so neither source # can be mistaken for an explicit user value. @@ -22,6 +22,15 @@ "domain": 1, "namelist": "namelist.wps", "overwrite": False, + "zsf_source": "nationalmap", + "download_dir": "downloads", +} + +# 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 +75,15 @@ 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" ) parser = argparse.ArgumentParser( prog="fire_preprocess", @@ -102,11 +114,17 @@ 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", @@ -115,6 +133,28 @@ def build_parser() -> argparse.ArgumentParser: "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( "--fuel-table", default=None, metavar="NAME|PATH", dest="fuel_table", help=( @@ -148,6 +188,9 @@ 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}") files = sorted(glob.glob(args.wps_files)) @@ -205,6 +248,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) 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/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 From 884ecbdec1c6319ffbb09ddb765486d3cb5e0d71 Mon Sep 17 00:00:00 2001 From: "Michael Kavulich, Jr" Date: Tue, 21 Jul 2026 07:54:55 -0600 Subject: [PATCH 3/3] Make the ZSF missing-data fill value configurable, defaulting to 0 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Previously, fire-grid pixels with no source elevation data (outside the DEM's extent, or nodata voids within it) surfaced the DEM source raster's own nodata sentinel (e.g. -99999) in the output ZSF field. Add --zsf-fill (default 0) to control this. Also works around a GDAL/rasterio quirk where dst_nodata=0.0 is silently ignored (falsy in the underlying binding) and falls back to the source's nodata value regardless — reprojection now always fills missing pixels with NaN internally and substitutes the requested fill value afterwards, so a fill of exactly 0 works correctly. NFUEL_CAT reprojection behavior is unchanged. Co-Authored-By: Claude Sonnet 5 --- README.md | 2 ++ fire_preprocess/cli.py | 11 +++++++++- fire_preprocess/raster.py | 44 +++++++++++++++++++++++++++------------ 3 files changed, 43 insertions(+), 14 deletions(-) diff --git a/README.md b/README.md index 6dee529..1d7a7f8 100644 --- a/README.md +++ b/README.md @@ -121,6 +121,7 @@ All arguments can also be supplied via a YAML config file (see below). | `--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) | @@ -142,6 +143,7 @@ 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/fire_preprocess/cli.py b/fire_preprocess/cli.py index b3924bf..9bc69e9 100644 --- a/fire_preprocess/cli.py +++ b/fire_preprocess/cli.py @@ -24,6 +24,7 @@ "overwrite": False, "zsf_source": "nationalmap", "download_dir": "downloads", + "zsf_fill": 0.0, } # LANDFIRE product to download when --fuel is omitted, keyed by fuel table. @@ -84,6 +85,7 @@ def build_parser() -> argparse.ArgumentParser: " 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", @@ -155,6 +157,13 @@ def build_parser() -> argparse.ArgumentParser: "--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", help=( @@ -304,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/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", )