erin-glen · erin-glen · Jan 26, 2026 · Jan 26, 2026 · Jan 27, 2026
diff --git a/analysis_core.py b/analysis_core.py
@@ -153,8 +153,14 @@ def perform_analysis(
                 # Zero out other disturbances, if any
                 if "insect_damage_HA" in forest_age_df.columns:
                     forest_age_df.loc[recat_fire, "insect_damage_HA"] = 0
-                if "harvest_HA" in forest_age_df.columns:
-                    forest_age_df.loc[recat_fire, "harvest_HA"] = 0
+                for column in [
+                    "harvest_0_25_HA",
+                    "harvest_25_50_HA",
+                    "harvest_50_75_HA",
+                    "harvest_75_100_HA",
+                ]:
+                    if column in forest_age_df.columns:
+                        forest_age_df.loc[recat_fire, column] = 0
                 if "undisturbed_HA" in forest_age_df.columns:
                     forest_age_df.loc[recat_fire, "undisturbed_HA"] = 0
 

diff --git a/config.py b/config.py
@@ -4,6 +4,7 @@
 # Base directories
 DATA_FOLDER = r"C:\GIS\Data\LEARN\SourceData"
 OUTPUT_BASE_DIR = r"C:\GIS\Data\LEARN\Outputs"
+DISTURBANCE_FOLDER = r"C:\GIS\Data\LEARN\Disturbances\NLCD_harvest_severity\final_disturbances"
 
 # Valid years for analysis
 VALID_YEARS = ["2001", "2004", "2006", "2008", "2011", "2013", "2016", "2019", "2021","2023"]
@@ -67,13 +68,13 @@ def get_input_config(year1, year2, aoi_name=None, tree_canopy_source=None):
     input_config = {
         "nlcd_1": os.path.join(
             DATA_FOLDER,
-            "LandCover",
-            f"NLCD_{year1}_Land_Cover_l48_20210604.tif"
+            "NEW_NLCD",
+            f"Annual_NLCD_LndCov_{year1}_CU_C1V0.tif"
         ),
         "nlcd_2": os.path.join(
             DATA_FOLDER,
-            "LandCover",
-            f"NLCD_{year2}_Land_Cover_l48_20210604.tif"
+            "NEW_NLCD",
+            f"Annual_NLCD_LndCov_{year2}_CU_C1V0.tif"
         ),
         "forest_age_raster": os.path.join(
             DATA_FOLDER, "ForestType", "forest_raster_01062025.tif"
@@ -100,19 +101,15 @@ def get_input_config(year1, year2, aoi_name=None, tree_canopy_source=None):
     # If tree_canopy_source is set, choose appropriate canopy data
     if tree_canopy_source:
         if tree_canopy_source == "NLCD":
-            tc_folder = os.path.join(DATA_FOLDER, "TreeCanopy", "NLCD_Project")
-
-            # ─────────────────────────────────────────────────────────────────
-            # NEW LOGIC: If user picks 2021-2023 for land cover, override
-            # canopy to 2019/2021
-            # ─────────────────────────────────────────────────────────────────
-            if (year1 == 2021 and year2 == 2023):
-                input_config["tree_canopy_1"] = os.path.join(tc_folder, "nlcd_tcc_conus_2019_v2021-4_projected.tif")
-                input_config["tree_canopy_2"] = os.path.join(tc_folder, "nlcd_tcc_conus_2021_v2021-4_projected.tif")
-            else:
-                input_config["tree_canopy_1"] = os.path.join(tc_folder, f"nlcd_tcc_conus_{year1}_v2021-4_projected.tif")
-                input_config["tree_canopy_2"] = os.path.join(tc_folder, f"nlcd_tcc_conus_{year2}_v2021-4_projected.tif")
-            # ─────────────────────────────────────────────────────────────────
+            tc_folder = os.path.join(DATA_FOLDER, "TreeCanopy", "NLCD_v2023-5_project")
+            input_config["tree_canopy_1"] = os.path.join(
+                tc_folder,
+                f"nlcd_tcc_conus_wgs84_v2023-5_{year1}0101_{year1}1231_projected.tif",
+            )
+            input_config["tree_canopy_2"] = os.path.join(
+                tc_folder,
+                f"nlcd_tcc_conus_wgs84_v2023-5_{year2}0101_{year2}1231_projected.tif",
+            )
 
         elif tree_canopy_source == "CBW":
             tc_folder = os.path.join(DATA_FOLDER, "TreeCanopy", "CBW")
@@ -135,15 +132,15 @@ def get_input_config(year1, year2, aoi_name=None, tree_canopy_source=None):
 
     # Disturbance rasters info
     disturbance_rasters_info = [
-        {"name": "disturbance_0104.tif", "start_year": 2001, "end_year": 2004},
-        {"name": "disturbance_0406.tif", "start_year": 2004, "end_year": 2006},
-        {"name": "disturbance_0608.tif", "start_year": 2006, "end_year": 2008},
-        {"name": "disturbance_0811.tif", "start_year": 2008, "end_year": 2011},
-        {"name": "disturbance_1113.tif", "start_year": 2011, "end_year": 2013},
-        {"name": "disturbance_1316.tif", "start_year": 2013, "end_year": 2016},
-        {"name": "disturbance_1619.tif", "start_year": 2016, "end_year": 2019},
-        {"name": "disturbance_1921.tif", "start_year": 2019, "end_year": 2021},
-        {"name": "disturbance_2123.tif", "start_year": 2021, "end_year": 2023},
+        {"name": "disturb_abs_2001_2004.tif", "start_year": 2001, "end_year": 2004},
+        {"name": "disturb_abs_2004_2006.tif", "start_year": 2004, "end_year": 2006},
+        {"name": "disturb_abs_2006_2008.tif", "start_year": 2006, "end_year": 2008},
+        {"name": "disturb_abs_2008_2011.tif", "start_year": 2008, "end_year": 2011},
+        {"name": "disturb_abs_2011_2013.tif", "start_year": 2011, "end_year": 2013},
+        {"name": "disturb_abs_2013_2016.tif", "start_year": 2013, "end_year": 2016},
+        {"name": "disturb_abs_2016_2019.tif", "start_year": 2016, "end_year": 2019},
+        {"name": "disturb_abs_2019_2021.tif", "start_year": 2019, "end_year": 2021},
+        {"name": "disturb_abs_2021_2023.tif", "start_year": 2021, "end_year": 2023},
     ]
 
     # Pick disturbance rasters fully inside the analysis period
@@ -153,7 +150,7 @@ def get_input_config(year1, year2, aoi_name=None, tree_canopy_source=None):
         end = dist_info["end_year"]
         # If the disturbance window is entirely within the user’s chosen [year1, year2]
         if (start >= year1) and (end <= year2):
-            dist_raster_path = os.path.join(DATA_FOLDER, "Disturbances", dist_info["name"])
+            dist_raster_path = os.path.join(DISTURBANCE_FOLDER, dist_info["name"])
             selected_disturbance_rasters.append(dist_raster_path)
 
     input_config["disturbance_rasters"] = selected_disturbance_rasters

diff --git a/funcs.py b/funcs.py
@@ -416,7 +416,14 @@ def compute_disturbance_max(
     ).reset_index()
 
     # Ensure necessary columns exist
-    required_columns = ["fire_HA", "harvest_HA", "insect_damage_HA"]
+    required_columns = [
+        "fire_HA",
+        "insect_damage_HA",
+        "harvest_0_25_HA",
+        "harvest_25_50_HA",
+        "harvest_50_75_HA",
+        "harvest_75_100_HA",
+    ]
     for col in required_columns:
         if col not in disturbance_wide.columns:
             disturbance_wide[col] = 0
@@ -595,15 +602,28 @@ def fill_na_values(forest_age_df: pd.DataFrame) -> pd.DataFrame:
         pd.DataFrame: Updated DataFrame.
     """
     # Fill NA values
-    forest_age_df["fire_HA"] = forest_age_df["fire_HA"].fillna(0)
-    forest_age_df["insect_damage_HA"] = forest_age_df["insect_damage_HA"].fillna(0)
-    forest_age_df["harvest_HA"] = forest_age_df["harvest_HA"].fillna(0)
+    disturbance_columns = [
+        "fire_HA",
+        "insect_damage_HA",
+        "harvest_0_25_HA",
+        "harvest_25_50_HA",
+        "harvest_50_75_HA",
+        "harvest_75_100_HA",
+    ]
+    for column in disturbance_columns:
+        if column not in forest_age_df.columns:
+            forest_age_df[column] = 0
+        else:
+            forest_age_df[column] = forest_age_df[column].fillna(0)
 
     # Calculate undisturbed area
     forest_age_df["undisturbed_HA"] = (
         forest_age_df["Hectares"]
         - forest_age_df["fire_HA"]
-        - forest_age_df["harvest_HA"]
+        - forest_age_df["harvest_0_25_HA"]
+        - forest_age_df["harvest_25_50_HA"]
+        - forest_age_df["harvest_50_75_HA"]
+        - forest_age_df["harvest_75_100_HA"]
         - forest_age_df["insect_damage_HA"]
     )
 
@@ -632,7 +652,10 @@ def merge_age_factors(
         "Forests Remaining Forest Removal Factor",
         "Fire Emissions Factor",
         "Insect Emissions Factor",
-        "Harvest Emissions Factor",
+        "Harvest 0-25 Emissions Factor",
+        "Harvest 25-50 Emissions Factor",
+        "Harvest 50-75 Emissions Factor",
+        "Harvest 75-100 Emissions Factor",
     ]
     forest_table = pd.read_csv(forest_lookup_csv, usecols=columns_to_use)
 
@@ -672,8 +695,29 @@ def calculate_forest_removals_and_emissions(
     forest_age_df["Annual_Emissions_Fire_CO2"] = (
         forest_age_df["fire_HA"] * forest_age_df["Fire Emissions Factor"] * (44 / 12) / years_difference
     )
-    forest_age_df["Annual_Emissions_Harvest_CO2"] = (
-        forest_age_df["harvest_HA"] * forest_age_df["Harvest Emissions Factor"] * (44 / 12) / years_difference
+    forest_age_df["Annual_Emissions_Harvest_0_25_CO2"] = (
+        forest_age_df["harvest_0_25_HA"]
+        * forest_age_df["Harvest 0-25 Emissions Factor"]
+        * (44 / 12)
+        / years_difference
+    )
+    forest_age_df["Annual_Emissions_Harvest_25_50_CO2"] = (
+        forest_age_df["harvest_25_50_HA"]
+        * forest_age_df["Harvest 25-50 Emissions Factor"]
+        * (44 / 12)
+        / years_difference
+    )
+    forest_age_df["Annual_Emissions_Harvest_50_75_CO2"] = (
+        forest_age_df["harvest_50_75_HA"]
+        * forest_age_df["Harvest 50-75 Emissions Factor"]
+        * (44 / 12)
+        / years_difference
+    )
+    forest_age_df["Annual_Emissions_Harvest_75_100_CO2"] = (
+        forest_age_df["harvest_75_100_HA"]
+        * forest_age_df["Harvest 75-100 Emissions Factor"]
+        * (44 / 12)
+        / years_difference
     )
     forest_age_df["Annual_Emissions_Insect_CO2"] = (
         forest_age_df["insect_damage_HA"] * forest_age_df["Insect Emissions Factor"] * (44 / 12) / years_difference
@@ -818,7 +862,10 @@ def calculate_area(
             "Undisturbed": "undisturbed_HA",
             "Fire": "fire_HA",
             "Insect/Disease": "insect_damage_HA",
-            "Harvest/Other": "harvest_HA",
+            "Harvest 0-25%": "harvest_0_25_HA",
+            "Harvest 25-50%": "harvest_25_50_HA",
+            "Harvest 50-75%": "harvest_50_75_HA",
+            "Harvest 75-100%": "harvest_75_100_HA",
         }
         column = columns_mapping.get(type_)
         if column:
@@ -897,7 +944,10 @@ def calculate_ghg_flux(
             "Undisturbed": "Annual_Removals_Undisturbed_CO2",
             "Fire": "Annual_Emissions_Fire_CO2",
             "Insect/Disease": "Annual_Emissions_Insect_CO2",
-            "Harvest/Other": "Annual_Emissions_Harvest_CO2",
+            "Harvest 0-25%": "Annual_Emissions_Harvest_0_25_CO2",
+            "Harvest 25-50%": "Annual_Emissions_Harvest_25_50_CO2",
+            "Harvest 50-75%": "Annual_Emissions_Harvest_50_75_CO2",
+            "Harvest 75-100%": "Annual_Emissions_Harvest_75_100_CO2",
         }
         column = columns_mapping.get(type_)
         if column:
@@ -966,7 +1016,10 @@ def summarize_ghg(
         ("Forest Remaining Forest", "Undisturbed", "Removals"),
         ("Forest Remaining Forest", "Fire", "Emissions"),
         ("Forest Remaining Forest", "Insect/Disease", "Emissions"),
-        ("Forest Remaining Forest", "Harvest/Other", "Emissions"),
+        ("Forest Remaining Forest", "Harvest 0-25%", "Emissions"),
+        ("Forest Remaining Forest", "Harvest 25-50%", "Emissions"),
+        ("Forest Remaining Forest", "Harvest 50-75%", "Emissions"),
+        ("Forest Remaining Forest", "Harvest 75-100%", "Emissions"),
     ]
 
     if include_trees_outside_forest:
@@ -1041,14 +1094,19 @@ def save_results(
     landuse_numeric_cols = [
         "Hectares", "CellCount",
         "carbon_ag_bg_us", "carbon_sd_dd_lt", "carbon_so",
-        "fire_HA", "harvest_HA", "insect_damage_HA",
+        "fire_HA", "insect_damage_HA",
+        "harvest_0_25_HA", "harvest_25_50_HA", "harvest_50_75_HA", "harvest_75_100_HA",
         "TreeCanopy_HA", "TreeCanopyLoss_HA",
         "Annual Emissions Forest to Non Forest CO2",
     ]
     forest_type_numeric_cols = [
-        "Hectares", "fire_HA", "harvest_HA", "insect_damage_HA", "undisturbed_HA",
+        "Hectares", "fire_HA", "insect_damage_HA", "undisturbed_HA",
+        "harvest_0_25_HA", "harvest_25_50_HA", "harvest_50_75_HA", "harvest_75_100_HA",
         "Annual_Removals_Undisturbed_CO2", "Annual_Removals_N_to_F_CO2",
-        "Annual_Emissions_Fire_CO2", "Annual_Emissions_Harvest_CO2", "Annual_Emissions_Insect_CO2"
+        "Annual_Emissions_Fire_CO2",
+        "Annual_Emissions_Harvest_0_25_CO2", "Annual_Emissions_Harvest_25_50_CO2",
+        "Annual_Emissions_Harvest_50_75_CO2", "Annual_Emissions_Harvest_75_100_CO2",
+        "Annual_Emissions_Insect_CO2"
     ]
 
     # 2) Round columns in landuse_result
@@ -1083,4 +1141,3 @@ def save_results(
     processing_time = dt.now() - start_time
     with open(os.path.join(output_path, f"processing_time{geography_suffix}.txt"), "w") as f:
         f.write(f"Total processing time: {processing_time}\n")
-
diff --git a/lookups.py b/lookups.py
@@ -50,7 +50,10 @@
 
 disturbanceLookup = {
 	#0: 'no_disturbance_HA',
-	1: 'harvest_HA',
+	1: 'harvest_0_25_HA',
+	2: 'harvest_25_50_HA',
+	3: 'harvest_50_75_HA',
+	4: 'harvest_75_100_HA',
 	5: 'insect_damage_HA',
 	10: 'fire_HA'
 }
@@ -94,4 +97,4 @@
     2: "GAP2",
     3: "GAP3",
     4: "GAP4"
-}
+}