Allow terrain coverage smaller than domain extent - near fire area only

.gitignore

@@ -2,3 +2,5 @@
 *.grb2
 *.grib2
 FILE:*
+__pycache__/
+*.pyc

README.md

@@ -166,3 +166,18 @@ reprojects to the WRF fire subgrid using the projection parameters and
 `subgrid_ratio_x/y` values from `namelist.wps`, and writes the fire fields
 into the existing WPS netCDF files. The rest of the workflow
 (`real.exe` → WRF) is unchanged.
+
+### Clipped source rasters
+
+The fuel and high-resolution DEM GeoTIFFs may cover only the intended fire
+area. Outside their joint valid coverage, `NFUEL_CAT` is set to the selected
+fuel table's no-fuel category. `ZSF` remains finite over the full fire array by
+bilinearly interpolating the atmospheric `HGT_M` field to the fire grid, then
+overlaying the high-resolution DEM where both rasters are valid.
+
+`DZDXF` and `DZDYF` are calculated from the merged `ZSF` using centered
+differences in the interior, one-sided differences at array edges, and WRF
+projection map factors. Both gradients are set to zero wherever their
+two-dimensional finite-difference stencil crosses the clipped-data boundary.
+This prevents a terrain discontinuity at the boundary from entering the fire
+spread calculation.

fire_preprocess/cli.py

@@ -7,6 +7,11 @@
 import numpy as np
 import yaml
 
+from .coverage import (
+    interpolate_hgt_to_fire_grid,
+    merge_fire_coverage,
+    read_atmospheric_terrain,
+)
 from .namelist import get_fire_subgrid_ratios, get_domain_params
 from .wrf_domain import read_domain_from_file
 from .fire_grid import build_fire_grid
@@ -78,7 +83,8 @@ def build_parser() -> argparse.ArgumentParser:
     parser = argparse.ArgumentParser(
         prog="fire_preprocess",
         description=(
-            "Add CFBM fire fields (NFUEL_CAT, ZSF) directly to WPS netCDF files,\n"
+            "Add CFBM static fire fields (NFUEL_CAT, ZSF, DZDXF, DZDYF) directly\n"
+            "to WPS netCDF files,\n"
             "bypassing GEOGRID.TBL editing and geogrid binary format conversion."
         ),
         formatter_class=argparse.RawDescriptionHelpFormatter,
@@ -183,30 +189,42 @@ def main(argv=None):
 
     # ── Reproject fuel categories ─────────────────────────────────────────────
     print(f"Reprojecting fuel data: {args.fuel}")
-    nfuel_raw = reproject_fuel(args.fuel, fire_grid)
-    nfuel = fuel_table.apply(nfuel_raw)
+    nfuel_raw, fuel_valid = reproject_fuel(
+        args.fuel, fire_grid, return_mask=True
+    )
     print(
-        f"  NFUEL_CAT shape {nfuel.shape}  "
-        f"range [{int(nfuel.min())}, {int(nfuel.max())}]"
+        f"  Valid fuel coverage: {np.count_nonzero(fuel_valid) / fuel_valid.size:.2%}"
     )
 
     # ── Reproject DEM ─────────────────────────────────────────────────────────
     print(f"Reprojecting terrain DEM: {args.zsf}")
-    zsf = reproject_dem(args.zsf, fire_grid)
-    print(
-        f"  ZSF shape {zsf.shape}  "
-        f"range [{zsf.min():.1f}, {zsf.max():.1f}] m"
+    zsf_raw, dem_valid = reproject_dem(args.zsf, fire_grid, return_mask=True)
+    terrain, terrain_name = read_atmospheric_terrain(files[0])
+    terrain_background = interpolate_hgt_to_fire_grid(
+        terrain,
+        sr_x,
+        sr_y,
+        (fire_grid.ny_fire, fire_grid.nx_fire),
     )
-
-    # ── Terrain gradients ─────────────────────────────────────────────────────
-    dzdxf, dzdyf = compute_slope(zsf, fire_grid)
-    grad_max = float(np.hypot(dzdxf, dzdyf).max())
+    nfuel, zsf, high_resolution_valid = merge_fire_coverage(
+        nfuel_raw,
+        zsf_raw,
+        fuel_valid,
+        dem_valid,
+        terrain_background,
+        fuel_table,
+    )
+    dzdxf, dzdyf = compute_slope(zsf, fire_grid, high_resolution_valid)
+    coverage_fraction = np.count_nonzero(high_resolution_valid) / high_resolution_valid.size
     print(
-        f"Terrain gradients from ZSF: max |grad| = {grad_max:.3f} "
-        f"({np.degrees(np.arctan(grad_max)):.1f}° slope)"
+        f"  Joint high-resolution coverage: {coverage_fraction:.2%}\n"
+        f"  ZSF background: interpolated {terrain_name}\n"
+        f"  ZSF shape {zsf.shape} range [{zsf.min():.1f}, {zsf.max():.1f}] m\n"
+        f"  Maximum |DZDXF|={np.max(np.abs(dzdxf)):.4f}, "
+        f"|DZDYF|={np.max(np.abs(dzdyf)):.4f}"
     )
 
-    # ── Write to WPS files ─────────────────────────────────────────────────
+    # ── Write to WPS files ────────────────────────────────────────────
     skipped = 0
     for file_path in files:
         label = os.path.basename(file_path)

fire_preprocess/coverage.py

@@ -0,0 +1,120 @@
+"""Construct full-domain fire fields from clipped high-resolution rasters."""
+
+from __future__ import annotations
+
+import netCDF4 as nc
+import numpy as np
+
+
+#--------------------------------------------------------------------------------
+# Atmospheric terrain background
+#--------------------------------------------------------------------------------
+
+def read_atmospheric_terrain(path: str) -> tuple[np.ndarray, str]:
+    """Read finite mass-grid terrain from a WPS output file."""
+    with nc.Dataset(path) as dataset:
+        terrain_name = "HGT_M"
+        if terrain_name not in dataset.variables:
+            raise KeyError(f"{path} does not contain HGT_M")
+
+        terrain_var = dataset.variables[terrain_name]
+        terrain = np.ma.filled(terrain_var[0], np.nan).astype(np.float32)
+
+    if terrain.ndim != 2 or not np.isfinite(terrain).all():
+        raise ValueError(
+            f"{terrain_name} in {path} must be a finite two-dimensional field"
+        )
+    return terrain, terrain_name
+
+
+def _linear_coordinates(source_size: int, ratio: int) -> np.ndarray:
+    """Return fire-cell center coordinates in atmospheric mass-grid units."""
+    active_size = source_size * ratio
+    return (np.arange(active_size) - 0.5 * (ratio - 1)) / ratio
+
+
+def _interpolate_axis(values: np.ndarray, coordinates: np.ndarray, axis: int) -> np.ndarray:
+    """Interpolate along one axis and linearly continue through boundary cells."""
+    source = np.moveaxis(np.asarray(values, dtype=np.float32), axis, -1)
+    source_size = source.shape[-1]
+    if source_size < 2:
+        raise ValueError("atmospheric terrain needs at least two points per axis")
+
+    lower = np.floor(coordinates).astype(np.int64)
+    lower_clipped = np.clip(lower, 0, source_size - 2)
+    upper_clipped = lower_clipped + 1
+    fraction = coordinates - lower_clipped
+
+    result = (
+        np.take(source, lower_clipped, axis=-1) * (1.0 - fraction)
+        + np.take(source, upper_clipped, axis=-1) * fraction
+    )
+    return np.moveaxis(result.astype(np.float32), -1, axis)
+
+
+def interpolate_hgt_to_fire_grid(
+    terrain: np.ndarray,
+    sr_x: int,
+    sr_y: int,
+    output_shape: tuple[int, int],
+) -> np.ndarray:
+    """Bilinearly interpolate mass-grid HGT onto the complete fire array."""
+    terrain_values = np.asarray(terrain, dtype=np.float32)
+    if terrain_values.ndim != 2 or not np.isfinite(terrain_values).all():
+        raise ValueError("atmospheric terrain must be a finite two-dimensional field")
+
+    y_coordinates = _linear_coordinates(terrain_values.shape[0], sr_y)
+    x_coordinates = _linear_coordinates(terrain_values.shape[1], sr_x)
+    interpolated_x = _interpolate_axis(terrain_values, x_coordinates, axis=1)
+    active = _interpolate_axis(interpolated_x, y_coordinates, axis=0)
+
+    expected_shape = (
+        (terrain_values.shape[0] + 1) * sr_y,
+        (terrain_values.shape[1] + 1) * sr_x,
+    )
+    if output_shape != expected_shape:
+        raise ValueError(
+            f"fire grid shape {output_shape} does not match expected {expected_shape}"
+        )
+
+    background = np.empty(output_shape, dtype=np.float32)
+    active_ny, active_nx = active.shape
+    background[:active_ny, :active_nx] = active
+    background[:active_ny, active_nx:] = active[:, -1, None]
+    background[active_ny:, :] = background[active_ny - 1, :]
+    return background
+
+
+#--------------------------------------------------------------------------------
+# High-resolution fire-area overlay
+#--------------------------------------------------------------------------------
+
+def merge_fire_coverage(
+    nfuel_raw: np.ndarray,
+    zsf_raw: np.ndarray,
+    fuel_valid: np.ndarray,
+    dem_valid: np.ndarray,
+    terrain_background: np.ndarray,
+    fuel_table,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Apply high-resolution data only within joint fuel and terrain coverage."""
+    shape = terrain_background.shape
+    inputs = (nfuel_raw, zsf_raw, fuel_valid, dem_valid)
+    if any(np.asarray(item).shape != shape for item in inputs):
+        raise ValueError("fuel, DEM, masks, and terrain background must share a shape")
+    if not np.isfinite(terrain_background).all():
+        raise ValueError("interpolated HGT_M background must be finite")
+
+    joint_valid = (
+        np.asarray(fuel_valid, dtype=bool)
+        & np.asarray(dem_valid, dtype=bool)
+        & np.isfinite(nfuel_raw)
+        & np.isfinite(zsf_raw)
+    )
+
+    nfuel = fuel_table.apply(nfuel_raw)
+    nfuel[~joint_valid] = float(fuel_table.nodata_out)
+
+    zsf = np.asarray(terrain_background, dtype=np.float32).copy()
+    zsf[joint_valid] = np.asarray(zsf_raw, dtype=np.float32)[joint_valid]
+    return nfuel.astype(np.float32), zsf, joint_valid

fire_preprocess/raster.py

@@ -4,29 +4,35 @@
   reproject_fuel(src_path, fire_grid)  →  float32 array (NFUEL_CAT)
   reproject_dem (src_path, fire_grid)  →  float32 array (ZSF)
 
+Both can optionally return a Boolean source-coverage mask with the field.
 Both map onto the same fire grid (south_north_subgrid × west_east_subgrid).
 rasterio.warp.reproject handles the CRS transformation from the source
 raster's native projection to the WRF map projection automatically.
 """
+
+from __future__ import annotations
+
 import warnings
+
 import numpy as np
 import rasterio
-from rasterio.warp import reproject, Resampling
 from rasterio.crs import CRS as RasterioCRS
+from rasterio.warp import Resampling, reproject
 
 from .fire_grid import FireGrid
 
 
-def _warn_coverage(arr, nodata_val, 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
+def _warn_coverage(
+    valid_mask: np.ndarray,
+    field_name: str,
+    threshold: float = 0.05,
+) -> None:
+    """Report substantial gaps without treating clipped rasters as an error."""
+    frac = 1.0 - np.count_nonzero(valid_mask) / valid_mask.size
     if frac > threshold:
         warnings.warn(
-            f"{field_name}: {frac:.1%} of fire-grid pixels are nodata — "
-            "check that the source raster covers the full WRF domain.",
+            f"{field_name}: {frac:.1%} of fire-grid pixels are outside the "
+            "source coverage; the workflow will use its configured background.",
             stacklevel=3,
         )
 
@@ -36,9 +42,10 @@ def reproject_to_grid(
     fire_grid: FireGrid,
     resampling: Resampling,
     band: int = 1,
-    src_nodata: float = None,
+    src_nodata: float | None = None,
     field_name: str = "field",
-) -> np.ndarray:
+    return_mask: bool = False,
+) -> np.ndarray | tuple[np.ndarray, np.ndarray]:
     """Core reprojection routine.
 
     Returns a 2-D float32 array in WRF row order (south_north, west_east):
@@ -47,9 +54,9 @@ def reproject_to_grid(
     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)
+    # NaN ensures pixels rasterio does not write, including domain areas
+    # outside the source extent, are represented explicitly in the mask.
+    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
@@ -61,19 +68,29 @@ def reproject_to_grid(
             src_transform=src.transform,
             src_crs=src.crs,
             src_nodata=nodata,
+            dst_nodata=np.nan,
             dst_transform=fire_grid.transform,
             dst_crs=dst_crs,
             resampling=resampling,
+            init_dest_nodata=True,
         )
 
     # rasterio returns row 0 = north; WRF netCDF expects row 0 = south.
     dst_data = np.flipud(dst_data)
+    valid_mask = np.isfinite(dst_data)
 
-    _warn_coverage(dst_data, nodata, field_name)
+    _warn_coverage(valid_mask, field_name)
+    if return_mask:
+        return dst_data, valid_mask
     return dst_data
 
 
-def reproject_fuel(src_path: str, fire_grid: FireGrid, band: int = 1) -> np.ndarray:
+def reproject_fuel(
+    src_path: str,
+    fire_grid: FireGrid,
+    band: int = 1,
+    return_mask: bool = False,
+) -> np.ndarray | tuple[np.ndarray, np.ndarray]:
     """Reproject a fuel-category GeoTIFF to the fire grid.
 
     Uses nearest-neighbour resampling to preserve integer category values.
@@ -83,11 +100,16 @@ def reproject_fuel(src_path: str, fire_grid: FireGrid, band: int = 1) -> np.ndar
         src_path, fire_grid,
         resampling=Resampling.nearest,
         band=band,
-        field_name="NFUEL_CAT",
+        field_name="NFUEL_CAT", return_mask=return_mask,
     )
 
 
-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,
+    return_mask: bool = False,
+) -> np.ndarray | tuple[np.ndarray, np.ndarray]:
     """Reproject a terrain DEM GeoTIFF to the fire grid.
 
     Uses bilinear resampling, matching the WPS geogrid four_pt interpolation
@@ -97,5 +119,5 @@ def reproject_dem(src_path: str, fire_grid: FireGrid, band: int = 1) -> np.ndarr
         src_path, fire_grid,
         resampling=Resampling.bilinear,
         band=band,
-        field_name="ZSF",
+        field_name="ZSF", return_mask=return_mask,
     )