SILC v8-v10: AGM MORTCAL year-100, FVS-NE + OSM-ACD comparison, refined inputs

Key additions in calibration/output/comparisons_overstory/:

10 year empirical CFI scorecard (silc_cfi/):
- Long horizon CFI scorecard at 10 yr horizon (n=6 routine pairs)
- AGM MORTCAL cuts BA bias from +18.5% to +7.8%, sawlog BdFt R^2 0.84
- Mechanism decomposition figure (under mortality plus over growth)
- Triple-check verification: pipeline reproduces to 1e-13

100 year strata projections:
- AcadianGY MORTCAL=TRUE on 11 byStrata stands (Cardinal 11076533)
- FVS-NE default + calibrated on same stands (Cardinal 11078908)
- OSM-ACD 100 yr with refined BGI=3902 (Cardinal 11080104, 6 of 11 stands)
- CFI plot backfill for Cedar A+B and Mixedwood A+B (Cardinal 11077925)
- Year 100 BA: AGM MORTCAL 98-108, FVS-NE cal 125-225, OSM-ACD 172-221

Refined inputs from SILC CFI v3 lat/long (46.4628 N, 68.4253 W, Davistown):
- BGI: 3000 to 3902 (ME_BGI_V1.tif raster sample)
- CSI: 12.0 to 15.78 m (CSI_2030.tif raster sample)
- FVS-NE site index: 42 to 35 ft (empirical from dominant heights)

Deliverables:
- SILC_4Model_Benchmark_v10_Weiskittel.pdf (18 slides, operational deck)
- SILC_4Model_Benchmark_Report_v8_Weiskittel.docx (companion report)
- silc_strata_5x2_AGM_MORTCAL_BA.png (100 yr trajectory figure)
- silc_cfi_long_mechanism.png, silc_cfi_long_mortcal_compare.png
- All raw trajectory CSVs and aggregated cell tables

Operational recommendation: AGM with MORTCAL on anchors the long horizon
floor at 35 to 39 cords/ac; FVS-NE calibrated brackets the upper bound
at 44 to 70 cords/ac; OSM-ACD provides additional cross check.

Author: holoros

calibration/output/comparisons_overstory/lookup_bgi_csi.py

@@ -0,0 +1,45 @@
+#!/usr/bin/env python3
+"""lookup_bgi_csi.py - Extract BGI and CSI from rasters at SILC CFI lat/long
+and the 11 byStrata stand locations on Cardinal.
+Runs on Cardinal."""
+import os, sys
+try:
+    import rasterio
+    from rasterio.warp import transform as rio_transform
+except Exception:
+    print("rasterio missing; install with: pip install --user rasterio")
+    sys.exit(1)
+import pandas as pd
+from pyproj import Transformer
+
+BGI_TIF = os.path.expanduser("~/raster_layers/bgi/ME_BGI_V1.tif")
+# Find CSI raster
+import glob
+CSI_TIFS = glob.glob(os.path.expanduser("~/raster_layers/csi/**/*.tif"), recursive=True)
+print("BGI:", BGI_TIF, "exists:", os.path.exists(BGI_TIF))
+print("CSI:", CSI_TIFS)
+
+def sample_raster(path, lat, lon):
+    with rasterio.open(path) as src:
+        # transform lat/long (EPSG:4326) to raster CRS
+        xs, ys = rio_transform("EPSG:4326", src.crs, [lon], [lat])
+        try:
+            v = next(src.sample([(xs[0], ys[0])]))[0]
+            return float(v) if v is not None else None
+        except StopIteration:
+            return None
+
+# === SILC CFI: single approximate location ===
+print("\n=== SILC CFI plots (approximate lat 46.4628, lon -68.4253) ===")
+bgi_cfi = sample_raster(BGI_TIF, 46.4628, -68.4253)
+print(f"  BGI: {bgi_cfi}")
+for c in CSI_TIFS[:3]:
+    v = sample_raster(c, 46.4628, -68.4253)
+    print(f"  CSI from {os.path.basename(c)}: {v}")
+
+# === byStrata stands: from StandInit lat/long if any ===
+print("\n=== byStrata stands ===")
+si = pd.read_csv(os.path.expanduser("~/silc_strata/Acadian_Matrix_StandInit_2023.csv"))
+# Stand init has Climate Site Index already
+print(si[["STAND_ID","INV_PLOT_SIZE","ClimateSiteIndexMeters","ElevationMeters"]].to_string(index=False))
+print("\nNote: byStrata stands have CSI already in StandInit (ClimateSiteIndexMeters column)")

calibration/output/comparisons_overstory/osm_strata_StandList.csv

@@ -0,0 +1,105 @@
+#!/usr/bin/env python3
+"""run_fvs_strata_100yr.py
+=====================================================================
+Run FVS-NE (default + calibrated) for 100 years on the 11 SILC
+byStrata stands, starting from the GrownDB year-2023 snapshot
+(matches the AGM MORTCAL 100-yr runner so trajectories are comparable).
+NUMCYCLE=20 with TIMEINT=5 (i.e. 5-yr cycles x 20 = 100 yr).
+"""
+from __future__ import annotations
+import os, sys, sqlite3, subprocess, tempfile
+import pandas as pd
+import numpy as np
+
+PROJECT_ROOT = os.environ.get("FVS_PROJECT_ROOT", "/users/PUOM0008/crsfaaron/fvs-modern")
+FVS_LIB_DIR  = os.environ.get("FVS_LIB_DIR", os.path.join(PROJECT_ROOT, "lib"))
+sys.path.insert(0, "/users/PUOM0008/crsfaaron/silc_cfi")
+from run_fvs_on_cfi import (build_standinit, KEYFILE, calibrated_keywords)
+
+# AGM alpha species code to FIA SPCD (Acadian / Maine region)
+AGY_TO_SPCD = {
+    "BF":  12, "RS":  97, "PB": 375, "YB": 371, "RM": 316,
+    "SM": 318, "WA": 541, "RO": 833, "BS":  95, "WS":  91,
+    "JP": 105, "RP": 125, "WP": 129, "EH": 261, "WC": 241,
+    "BC": 531, "QA": 746, "GB": 934, "BA": 543, "AB": 531,
+    "ST": 951, "TA": 547, "OH": 998, "OS": 299, "HW": 998,
+    "PC": 117, "PR": 125, "AS": 543, "BT": 951, "EC": 241,
+    "HH": 701, "NS": 94, "RB": 372, "RN": 94, "BP": 98,
+    "SW": 998, "SB":  95, "TL": 547, "WB": 375,
+}
+
+def build_treeinit_strata(td, sid):
+    """GrownDB year-2023 snapshot already has imperial DBH (in), Ht (ft),
+    and TPA (per-acre). Map alpha species to FIA SPCD."""
+    rows = []
+    for i, r in td.iterrows():
+        spcd = AGY_TO_SPCD.get(str(r["Species"]), 998)
+        dbh = float(r["DBH"])
+        if not (np.isfinite(dbh) and dbh > 0):
+            continue
+        ht = float(r["Ht"]) if pd.notna(r["Ht"]) and r["Ht"] > 0 else 0.0
+        tpa = float(r["TPA"])
+        rows.append({
+            "stand_id": sid,
+            "plot_id": 1,
+            "tree_id": i + 1,
+            "tree_count": round(tpa, 6),
+            "species": spcd,
+            "diameter": round(dbh, 3),
+            "ht": round(ht, 1),
+            "crratio": 40,
+        })
+    return pd.DataFrame(rows)
+
+def run_one_100yr(sid, tree_df, inv_year, variant, config):
+    """Run FVS for 20 cycles of 5 yr each = 100 yr horizon."""
+    binary = os.path.join(FVS_LIB_DIR, f"FVS{variant.lower()}")
+    with tempfile.TemporaryDirectory() as d:
+        db = os.path.join(d, "FVS_Data.db")
+        con = sqlite3.connect(db)
+        build_standinit(sid, inv_year, variant).to_sql("fvs_standinit", con, if_exists="replace", index=False)
+        tree_df.to_sql("fvs_treeinit", con, if_exists="replace", index=False)
+        con.close()
+        calib = calibrated_keywords(variant) if config == "calibrated" else "** DEFAULT"
+        key = os.path.join(d, "strata.key")
+        open(key, "w").write(KEYFILE.format(sid=sid, db=db, clen=5, ncyc=20, calib=calib))
+        try:
+            subprocess.run([binary, f"--keywordfile={key}"], cwd=d, capture_output=True, timeout=600)
+        except subprocess.TimeoutExpired:
+            return None
+        try:
+            df = pd.read_sql_query("SELECT * FROM FVS_Summary2", sqlite3.connect(db))
+            df["variant"] = variant.upper()
+            df["config"]  = config
+            return df
+        except Exception:
+            return None
+
+def main():
+    wd = os.getcwd()
+    gr = pd.read_csv(os.path.join(wd, "GrownDB_byStrata_ALL.csv"))
+    gr2023 = gr[(gr.Year == 2023) & gr.DBH.notna() & gr.TPA.notna() & (gr.TPA > 0)]
+    stand_ids = sorted(gr2023["StandID"].unique())
+    print(f"Running FVS-NE 100-yr on {len(stand_ids)} byStrata stands")
+
+    rows = []
+    for sid in stand_ids:
+        td = gr2023[gr2023.StandID == sid]
+        tree_df = build_treeinit_strata(td, sid)
+        for variant in ("ne",):  # FVS-NE only
+            for config in ("default", "calibrated"):
+                df = run_one_100yr(sid, tree_df, 2023, variant, config)
+                if df is None or len(df) == 0:
+                    print(f"  {sid} {variant} {config}: FAILED")
+                    continue
+                df["StandID"] = sid
+                rows.append(df)
+                last = df.iloc[-1]
+                print(f"  {sid} {variant} {config}: yr{int(last.get('Year', 2123))} BA={last.get('BA', 0):.1f} TPA={last.get('Tpa', 0):.0f} Cords={last.get('MCuFt', 0)/79:.1f}")
+
+    out = pd.concat(rows, ignore_index=True) if rows else pd.DataFrame()
+    out.to_csv(os.path.join(wd, "silc_strata_100yr_fvsne_results.csv"), index=False)
+    print(f"\nWrote {len(out)} rows for {out['StandID'].nunique() if len(out) else 0} stands")
+
+if __name__ == "__main__":
+    main()

calibration/output/comparisons_overstory/run_fvs_strata_100yr.py

@@ -0,0 +1,103 @@
+#!/usr/bin/env python3
+"""run_fvs_strata_100yr_v2.py
+v2: refined site index. Default SI=42 (balsam fir) was too optimistic
+for these mature spruce/fir/mixed stands. Empirically observed dominant
+heights of 25-35 ft on the SILC CFI at ~30-50 yr stand age put true SI
+at ~35 ft. Lowering to SI=35 gives less optimistic growth projections.
+"""
+from __future__ import annotations
+import os, sys, sqlite3, subprocess, tempfile
+import pandas as pd, numpy as np
+
+PROJECT_ROOT = os.environ.get("FVS_PROJECT_ROOT", "/users/PUOM0008/crsfaaron/fvs-modern")
+FVS_LIB_DIR  = os.environ.get("FVS_LIB_DIR", os.path.join(PROJECT_ROOT, "lib"))
+sys.path.insert(0, "/users/PUOM0008/crsfaaron/silc_cfi")
+from run_fvs_on_cfi import KEYFILE, calibrated_keywords
+
+SITE_INDEX_REFINED = 35  # was 42; empirical estimate for mature spruce/fir
+
+# Override build_standinit to use refined SI
+def build_standinit_v2(sid, inv_year, variant):
+    return pd.DataFrame([{
+        "stand_id": sid, "variant": variant.upper(), "inv_year": int(inv_year),
+        "latitude": 46.4628, "longitude": -68.4253, "region": 9,
+        "forest": 0, "district": 0,
+        "basal_area_factor": 0.0, "inv_plot_size": 1.0, "brk_dbh": 999.0,
+        "num_plots": 1,
+        "age": 60, "aspect": 0, "slope": 5, "elevft": 1000,
+        "site_species": 12, "site_index": SITE_INDEX_REFINED,
+        "state": 23, "county": 21,
+        "forest_type": 121, "sam_wt": 1.0,
+    }])
+
+AGY_TO_SPCD = {
+    "BF":  12, "RS":  97, "PB": 375, "YB": 371, "RM": 316,
+    "SM": 318, "WA": 541, "RO": 833, "BS":  95, "WS":  91,
+    "JP": 105, "RP": 125, "WP": 129, "EH": 261, "WC": 241,
+    "BC": 531, "QA": 746, "GB": 934, "BA": 543, "AB": 531,
+    "ST": 951, "TA": 547, "OH": 998, "OS": 299, "HW": 998,
+    "PC": 117, "PR": 125, "AS": 543, "BT": 951, "EC": 241,
+    "HH": 701, "NS": 94, "RB": 372, "RN": 94, "BP": 98,
+    "SW": 998, "SB":  95, "TL": 547, "WB": 375,
+}
+
+def build_treeinit_strata(td, sid):
+    rows = []
+    for i, r in td.iterrows():
+        spcd = AGY_TO_SPCD.get(str(r["Species"]), 998)
+        dbh = float(r["DBH"])
+        if not (np.isfinite(dbh) and dbh > 0): continue
+        ht = float(r["Ht"]) if pd.notna(r["Ht"]) and r["Ht"] > 0 else 0.0
+        tpa = float(r["TPA"])
+        rows.append({"stand_id": sid, "plot_id": 1, "tree_id": i + 1,
+                     "tree_count": round(tpa, 6), "species": spcd,
+                     "diameter": round(dbh, 3), "ht": round(ht, 1), "crratio": 40})
+    return pd.DataFrame(rows)
+
+def run_one_100yr(sid, tree_df, inv_year, variant, config):
+    binary = os.path.join(FVS_LIB_DIR, f"FVS{variant.lower()}")
+    with tempfile.TemporaryDirectory() as d:
+        db = os.path.join(d, "FVS_Data.db")
+        con = sqlite3.connect(db)
+        build_standinit_v2(sid, inv_year, variant).to_sql("fvs_standinit", con, if_exists="replace", index=False)
+        tree_df.to_sql("fvs_treeinit", con, if_exists="replace", index=False)
+        con.close()
+        calib = calibrated_keywords(variant) if config == "calibrated" else "** DEFAULT"
+        key = os.path.join(d, "strata.key")
+        open(key, "w").write(KEYFILE.format(sid=sid, db=db, clen=5, ncyc=20, calib=calib))
+        try:
+            subprocess.run([binary, f"--keywordfile={key}"], cwd=d, capture_output=True, timeout=600)
+        except subprocess.TimeoutExpired:
+            return None
+        try:
+            df = pd.read_sql_query("SELECT * FROM FVS_Summary2", sqlite3.connect(db))
+            df["variant"] = variant.upper(); df["config"]  = config
+            return df
+        except Exception:
+            return None
+
+def main():
+    wd = os.getcwd()
+    gr = pd.read_csv(os.path.join(wd, "GrownDB_byStrata_ALL.csv"))
+    gr2023 = gr[(gr.Year == 2023) & gr.DBH.notna() & gr.TPA.notna() & (gr.TPA > 0)]
+    stand_ids = sorted(gr2023["StandID"].unique())
+    print(f"Running FVS-NE 100-yr v2 (SI={SITE_INDEX_REFINED}) on {len(stand_ids)} byStrata stands")
+    rows = []
+    for sid in stand_ids:
+        td = gr2023[gr2023.StandID == sid]
+        tree_df = build_treeinit_strata(td, sid)
+        for variant in ("ne",):
+            for config in ("default", "calibrated"):
+                df = run_one_100yr(sid, tree_df, 2023, variant, config)
+                if df is None or len(df) == 0:
+                    print(f"  {sid} {variant} {config}: FAILED"); continue
+                df["StandID"] = sid
+                rows.append(df)
+                last = df.iloc[-1]
+                print(f"  {sid} {variant} {config}: yr2123 BA={last.get('BA',0):.1f} Cords={last.get('MCuFt',0)/79:.1f}")
+    out = pd.concat(rows, ignore_index=True) if rows else pd.DataFrame()
+    out.to_csv(os.path.join(wd, "silc_strata_100yr_fvsne_si35_results.csv"), index=False)
+    print(f"\nWrote {len(out)} rows for {out['StandID'].nunique() if len(out) else 0} stands")
+
+if __name__ == "__main__":
+    main()

calibration/output/comparisons_overstory/run_fvs_strata_100yr_v2.py

@@ -0,0 +1,134 @@
+#!/usr/bin/env python3
+"""run_osm_strata_100yr.py - OSM-ACD 100-yr driver for SILC byStrata
+stands. Uses GrownDB year-2023 snapshot for the starting tree list
+(matches AGM/FVS-NE), with refined per-stand BGI=3902 from ME_BGI_V1
+raster at Davistown lat/long and CSI from StandInit."""
+import sqlite3, subprocess
+from io import StringIO
+from pathlib import Path
+import numpy as np, pandas as pd
+
+OSM_BINARY = Path.home() / "OSM/v2.26.1/OSMv2.26.1_Linux64/OSM.ConsoleApp"
+BGI_REFINED = 3902  # from ME_BGI_V1.tif sampled at lat 46.4628, lon -68.4253
+BA_PER_TREE_CM = 0.00007854
+FT2_AC_PER_M2_HA = 4.35
+ACRES_PER_HA = 2.4710538147
+N_YEARS = 100
+
+OSM_REMAP = {"PB": "WB", "TA": "TL", "PC": "PR"}
+acadgy_species = {"AB","AS","BA","BC","BF","BP","BS","BT","EC","EH","GA","GB",
+                  "HH","HW","JP","NS","OH","OS","PB","PC","PR","QA","RB","RM",
+                  "RN","RO","RP","RS","SB","SM","ST","SW","TA","WA","WC","WP",
+                  "WS","YB"}
+
+def to_osm_sp(sp):
+    s = str(sp)
+    if s in acadgy_species:
+        return OSM_REMAP.get(s, s)
+    return "OS"
+
+def build_inputs(stand_ids, gr2023, si_lookup, out_dir):
+    out_dir.mkdir(parents=True, exist_ok=True)
+    sqlite_path = out_dir / "osm_strata.sqlite"
+    osmc_path   = out_dir / "osm_strata.osmc"
+    stand_csv   = out_dir / "osm_strata_StandList.csv"
+    tree_csv    = out_dir / "osm_strata_TreeList.csv"
+    if sqlite_path.exists(): sqlite_path.unlink()
+
+    sid_map = {}
+    stand_rows = []
+    tree_rows  = []
+    for i, sid in enumerate(stand_ids):
+        n = i + 1
+        sid_map[sid] = n
+        stand_rows.append((n, 1, "ME", "None", BGI_REFINED, 0))
+        td = gr2023[gr2023["StandID"] == sid]
+        for _, r in td.iterrows():
+            dbh = float(r["DBH"]) if pd.notna(r["DBH"]) else 0
+            if dbh <= 0: continue
+            dbh_cm = dbh * 2.54
+            ht = float(r["Ht"]) if pd.notna(r["Ht"]) and r["Ht"] > 0 else 0
+            if ht <= 0:
+                ht_ft = max(6, 4.27 + 82 * (1 - np.exp(-0.04 * dbh_cm)))
+                ht_m = ht_ft * 0.3048
+            else:
+                ht_m = ht * 0.3048
+            expf_ha = float(r["TPA"]) * ACRES_PER_HA
+            sp = to_osm_sp(str(r["Species"]))
+            tree_rows.append((n, sp, dbh_cm, ht_m, expf_ha))
+
+    con = sqlite3.connect(str(sqlite_path)); cur = con.cursor()
+    cur.execute("CREATE TABLE OSM_StandList (SurveyID INT, Plots INT, Zone TEXT, "
+                "Management TEXT, BGI INT, PoorSite SMALLINT)")
+    cur.execute("CREATE TABLE OSM_TreeList (SurveyID INT, Species TEXT, DBH REAL, "
+                "HT REAL, Stems REAL)")
+    cur.execute("CREATE INDEX idx_tl_sid ON OSM_TreeList(SurveyID)")
+    cur.executemany("INSERT INTO OSM_StandList VALUES (?,?,?,?,?,?)", stand_rows)
+    cur.executemany("INSERT INTO OSM_TreeList (SurveyID, Species, DBH, HT, Stems) VALUES (?,?,?,?,?)", tree_rows)
+    con.commit(); con.close()
+
+    lines = ["SIMULATION", f" YEARS {N_YEARS}", " YPC 5", "",
+             f'INPUTS.SOURCE "{sqlite_path}"', "", "OUTPUTS", " Messages FALSE",
+             " StandSummary.ConsoleOn FALSE",
+             f' StandSummary.FilePath "{stand_csv}"',
+             f' TreeList.FilePath "{tree_csv}"', "", "SIMULATION.Scenario BASE"]
+    for sid in sorted(sid_map.values()): lines.append(f"SIMULATE {sid}")
+    osmc_path.write_text("\n".join(lines) + "\n")
+    return sqlite_path, osmc_path, stand_csv, tree_csv, sid_map
+
+def parse_stand(stand_csv):
+    raw = Path(stand_csv).read_text(errors="replace").replace("\x00", "")
+    df = pd.read_csv(StringIO(raw))
+    def col(*names):
+        for n in names:
+            if n in df.columns: return n
+    c_sid = col("SurveyID","Survey_ID","Id")
+    c_yr  = col("Year","Yr","Cycle","Age")
+    c_ba  = col("BA","BasalArea","BasalArea_m2ha","GBA")
+    c_tph = col("Trees","TPH","Stems","Stems_per_ha","Density")
+    c_qmd = col("QMD","Qmd","Qmd_cm")
+    df = df[df[c_sid].notna()].copy()
+    df["SurveyID"] = df[c_sid].astype(float).astype(int)
+    df["yr_off"]   = df.groupby("SurveyID")[c_yr].transform(
+        lambda s: (s.astype(float) - s.astype(float).min()).round().astype(int))
+    df["BA_m2ha"]  = df[c_ba].astype(float)
+    df["TPH_ha"]   = df[c_tph].astype(float)
+    df["QMD_cm"]   = (df[c_qmd].astype(float) if c_qmd
+                     else np.sqrt((df["BA_m2ha"] / df["TPH_ha"].replace(0, np.nan)) / BA_PER_TREE_CM))
+    return df[["SurveyID","yr_off","BA_m2ha","TPH_ha","QMD_cm"]]
+
+def main():
+    wd = Path.cwd()
+    gr = pd.read_csv(wd / "GrownDB_byStrata_ALL.csv")
+    gr2023 = gr[(gr["Year"] == 2023) & gr["DBH"].notna() & gr["TPA"].notna() & (gr["TPA"] > 0)]
+    si = pd.read_csv(wd / "Acadian_Matrix_StandInit_2023.csv")
+    si_lookup = dict(zip(si["STAND_ID"], si["ClimateSiteIndexMeters"]))
+    stand_ids = sorted(gr2023["StandID"].unique())
+    print(f"Running OSM-ACD 100-yr on {len(stand_ids)} byStrata stands, BGI={BGI_REFINED}")
+
+    out_dir = wd / "osm_strata_work"
+    sqlite_path, osmc_path, stand_csv, tree_csv, sid_map = build_inputs(
+        stand_ids, gr2023, si_lookup, out_dir)
+    print(f"Built inputs for {len(sid_map)} stands")
+    res = subprocess.run([str(OSM_BINARY), "Acadian", str(osmc_path)],
+                         capture_output=True, text=True, timeout=1800)
+    if res.returncode != 0:
+        print(f"OSM failed: {res.stderr[-500:]}")
+        raise SystemExit(1)
+    sp = parse_stand(stand_csv)
+    print(f"Parsed {len(sp)} stand-year rows")
+
+    # Convert to imperial and map back to StandID
+    sid_df = pd.DataFrame([(s, sid) for sid, s in sid_map.items()],
+                           columns=["StandID","SurveyID"])
+    sp = sp.merge(sid_df, on="SurveyID", how="left")
+    sp["Year"]     = 2023 + sp["yr_off"]
+    sp["BA"]       = sp["BA_m2ha"]   * FT2_AC_PER_M2_HA
+    sp["Tpa"]      = sp["TPH_ha"]    / ACRES_PER_HA
+    sp["QMD"]      = sp["QMD_cm"]    / 2.54
+    sp = sp[["StandID","Year","BA","Tpa","QMD"]]
+    sp.to_csv(wd / "silc_strata_100yr_osmacd_results.csv", index=False)
+    print(f"Wrote silc_strata_100yr_osmacd_results.csv: {len(sp)} rows")
+
+if __name__ == "__main__":
+    main()