Add tile-aware harvest orchestration and CLI examples

disturbance_config.py

@@ -3,6 +3,7 @@
 import os
 import glob
 import logging
+from typing import Iterable, List, Sequence
 
 # --------------------------------------------------------------------
 # LOGGING CONFIG
@@ -22,7 +23,7 @@
 INSECT_OUTPUT_DIR = os.path.join(INSECT_GDB_DIR, "Processed")
 INSECT_FINAL_DIR = os.path.join(INSECT_GDB_DIR, "Final")
 
-# Hansen
+# Hansen (Hansen Global Forest Change harvest proxy)
 HANSEN_INPUT_DIR = os.path.join(BASE_DIR, "Hansen")
 HANSEN_OUTPUT_DIR = os.path.join(HANSEN_INPUT_DIR, "Processed")
 
@@ -71,6 +72,63 @@ def _resolve_tcc_path(year: int) -> str:
 # Thresholds (upper bounds) for assigning canopy-loss severity classes 1-4
 NLCD_TCC_SEVERITY_BREAKS = [25, 50, 75, 100]
 
+# --------------------------------------------------------------------
+# HARVEST WORKFLOW SELECTION
+# --------------------------------------------------------------------
+# Two harvest/other processing options are supported:
+#   * "hansen"              => harvest_other.py
+#   * "nlcd_tcc_severity"   => harvest_other_severity.py
+# Update HARVEST_WORKFLOW to toggle between the legacy Hansen workflow
+# and the newer NLCD Tree Canopy Cover severity workflow.
+HARVEST_WORKFLOW = "hansen"
+
+HARVEST_PRODUCTS = {
+    "hansen": {
+        "module": "harvest_other",
+        "description": "Hansen Global Forest Change based harvest/other",
+        "raster_directory": HANSEN_OUTPUT_DIR,
+        "raster_template": "hansen_{period}.tif",
+        "combined_output_subdir": "hansen",
+    },
+    "nlcd_tcc_severity": {
+        "module": "harvest_other_severity",
+        "description": "NLCD Tree Canopy Cover change severity",
+        "raster_directory": NLCD_TCC_OUTPUT_DIR,
+        "raster_template": "nlcd_tcc_severity_{period}.tif",
+        "combined_output_subdir": "nlcd_tcc_severity",
+    },
+}
+
+
+def harvest_product_config(workflow: str = None):
+    """Return the harvest configuration for the selected workflow."""
+
+    wf = workflow or HARVEST_WORKFLOW
+    if wf not in HARVEST_PRODUCTS:
+        raise ValueError(
+            f"Unknown harvest workflow '{wf}'. Valid options: {sorted(HARVEST_PRODUCTS)}"
+        )
+    return dict(HARVEST_PRODUCTS[wf])
+
+
+def harvest_raster_path(period_name: str, workflow: str = None) -> str:
+    """Return the expected harvest raster path for the given period."""
+
+    config = harvest_product_config(workflow)
+    directory = config["raster_directory"]
+    template = config["raster_template"]
+    return os.path.join(directory, template.format(period=period_name))
+
+
+def final_combined_dir(workflow: str = None) -> str:
+    """Directory for final disturbance rasters for the selected harvest workflow."""
+
+    config = harvest_product_config(workflow)
+    combined_subdir = config.get("combined_output_subdir") or (workflow or HARVEST_WORKFLOW)
+    path = os.path.join(FINAL_COMBINED_ROOT_DIR, combined_subdir)
+    os.makedirs(path, exist_ok=True)
+    return path
+
 # --------------------------------------------------------------------
 # NLCD LAND COVER (reference raster for env)
 # --------------------------------------------------------------------
@@ -98,20 +156,29 @@ def nlcd_lc_path(year: int) -> str:
 # FINAL COMBINATION OUTPUTS
 # --------------------------------------------------------------------
 INTERMEDIATE_COMBINED_DIR = os.path.join(BASE_DIR, "Intermediate")
-FINAL_COMBINED_DIR = os.path.join(BASE_DIR, "FinalCombined")
+FINAL_COMBINED_ROOT_DIR = os.path.join(BASE_DIR, "FinalCombined")
+# Backwards compatibility: legacy callers may still reference FINAL_COMBINED_DIR
+FINAL_COMBINED_DIR = FINAL_COMBINED_ROOT_DIR
 
 # Create directories if they don’t exist:
-for _dir in [
+_BASE_DIRS = [
     INSECT_OUTPUT_DIR,
     INSECT_FINAL_DIR,
     HANSEN_OUTPUT_DIR,
     NLCD_TCC_OUTPUT_DIR,
     FIRE_OUTPUT_DIR,
     INTERMEDIATE_COMBINED_DIR,
-    FINAL_COMBINED_DIR
-]:
+    FINAL_COMBINED_ROOT_DIR,
+]
+
+for _dir in _BASE_DIRS:
     os.makedirs(_dir, exist_ok=True)
 
+for _cfg in HARVEST_PRODUCTS.values():
+    combined_subdir = _cfg.get("combined_output_subdir") or ""
+    if combined_subdir:
+        os.makedirs(os.path.join(FINAL_COMBINED_ROOT_DIR, combined_subdir), exist_ok=True)
+
 # --------------------------------------------------------------------
 # REGIONS & TIME PERIODS
 # --------------------------------------------------------------------
@@ -154,3 +221,77 @@ def nlcd_lc_path(year: int) -> str:
     r"C:\GIS\Data\LEARN\Disturbances\Hansen\GFW2023_60N_120W.tif",
     r"C:\GIS\Data\LEARN\Disturbances\Hansen\GFW2023_60N_130W.tif",
 ]
+
+
+def _hansen_tile_id(tile_path: str) -> str:
+    """Return the bare tile identifier (filename without extension)."""
+
+    base = os.path.basename(tile_path)
+    stem, _ = os.path.splitext(base)
+    return stem
+
+
+_HANSEN_TILE_LOOKUP = {
+    _hansen_tile_id(path).upper(): path
+    for path in HANSEN_TILES
+}
+
+
+def normalize_tile_ids(tile_ids: Sequence[str]) -> List[str]:
+    """Normalize user-provided Hansen tile identifiers."""
+
+    normed: List[str] = []
+    for raw in tile_ids:
+        if raw is None:
+            continue
+        tile = raw.strip().upper()
+        if not tile:
+            continue
+        tile = tile.replace(".TIF", "")
+        if tile in _HANSEN_TILE_LOOKUP:
+            normed.append(tile)
+            continue
+        prefixed = tile if tile.startswith("GFW") else f"GFW2023_{tile.lstrip('_')}"
+        normed.append(prefixed)
+    return normed
+
+
+def hansen_tile_paths(tile_ids: Iterable[str] | None = None) -> List[str]:
+    """Return Hansen tile paths filtered by ``tile_ids`` if provided."""
+
+    if not tile_ids:
+        return list(HANSEN_TILES)
+
+    selected: List[str] = []
+    missing: List[str] = []
+    for tile in normalize_tile_ids(tile_ids):
+        path = _HANSEN_TILE_LOOKUP.get(tile)
+        if path:
+            selected.append(path)
+        else:
+            missing.append(tile)
+
+    if missing:
+        logging.warning(
+            "Ignored %d Hansen tile id(s) with no match: %s",
+            len(missing),
+            ", ".join(sorted(set(missing))),
+        )
+
+    # Deduplicate while preserving order of input sequence
+    seen = set()
+    ordered: List[str] = []
+    for path in selected:
+        if path not in seen:
+            ordered.append(path)
+            seen.add(path)
+
+    if not ordered:
+        logging.error("No Hansen tiles matched the provided ids. Using full tile list instead.")
+        return list(HANSEN_TILES)
+
+    logging.info("Selected %d Hansen tile(s).", len(ordered))
+    for path in ordered:
+        logging.info("  %s", path)
+
+    return ordered

final_disturbance.py

@@ -1,7 +1,16 @@
 # final_disturbance.py
 """
 Combines fire, insect, and harvest rasters via CellStatistics (MAX)
-and optionally masks out low-severity fire.
+and optionally masks out low-severity fire. The harvest raster source is
+controlled by disturbance_config.HARVEST_WORKFLOW so either the legacy
+Hansen workflow or the newer NLCD Tree Canopy Cover severity workflow can
+be used without code changes.
+
+Example
+-------
+```
+python final_disturbance.py
+```
 """
 
 import arcpy
@@ -28,33 +37,49 @@ def main():
     arcpy.env.outputCoordinateSystem = cfg.NLCD_RASTER
     arcpy.env.overwriteOutput = True
 
+    harvest_cfg = cfg.harvest_product_config()
+    logging.info(
+        "Using harvest workflow '%s' => %s",
+        cfg.HARVEST_WORKFLOW,
+        harvest_cfg.get("description", ""),
+    )
+    logging.info("Harvest rasters directory => %s", harvest_cfg.get("raster_directory"))
+
+    final_output_dir = cfg.final_combined_dir()
+    logging.info("Final disturbance outputs will be written to => %s", final_output_dir)
+
     for period in cfg.TIME_PERIODS.keys():
         # Build paths
         fire_path = os.path.join(cfg.FIRE_OUTPUT_DIR, f"fire_{period}.tif")
         insect_path = os.path.join(cfg.INSECT_FINAL_DIR, f"insect_damage_{period}.tif")
-        hansen_path = os.path.join(cfg.HANSEN_OUTPUT_DIR, f"hansen_{period}.tif")
+        harvest_path = cfg.harvest_raster_path(period)
 
         if not (arcpy.Exists(fire_path) and
                 arcpy.Exists(insect_path) and
-                arcpy.Exists(hansen_path)):
+                arcpy.Exists(harvest_path)):
             logging.warning(f"Missing one or more rasters for period={period}. Skipping.")
             continue
 
-        logging.info(f"Combining Fire={fire_path}, Insect={insect_path}, Harvest={hansen_path}")
+        logging.info(
+            "Combining Fire=%s, Insect=%s, Harvest=%s",
+            fire_path,
+            insect_path,
+            harvest_path,
+        )
 
         # Max combination
         fire_ras = Raster(fire_path)
         insect_ras = Raster(insect_path)
-        hansen_ras = Raster(hansen_path)
+        harvest_ras = Raster(harvest_path)
 
-        combined_max = CellStatistics([fire_ras, insect_ras, hansen_ras], "MAXIMUM", "DATA")
+        combined_max = CellStatistics([fire_ras, insect_ras, harvest_ras], "MAXIMUM", "DATA")
 
         # Example: mask out code=3 => 0 (low severity fire).
         # If you do NOT want to mask out code=3, skip this step.
         final_dist = Con(combined_max == 3, 0, combined_max)
 
         # Save
-        out_raster_path = os.path.join(cfg.FINAL_COMBINED_DIR, f"disturb_{period}.tif")
+        out_raster_path = os.path.join(final_output_dir, f"disturb_{period}.tif")
         final_dist.save(out_raster_path)
         logging.info(f"Final combined disturbance => {out_raster_path}")
 

fire.py

@@ -11,6 +11,12 @@
    - For each time period in cfg.TIME_PERIODS, gather the reclassified rasters
      for each year in that period, perform CellStatistics (MAX), and save
      the final as fire_{period}.tif.
+
+Example
+-------
+```
+python fire.py
+```
 """
 
 import arcpy