diff --git a/analysis_core.py b/analysis_core.py index 0a5a0d3..4ce64bc 100644 --- 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 index 3387fe5..f5b4f0a 100644 --- 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 index 90af015..4032a7a 100644 --- 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 index fe44f7b..f186e83 100644 --- 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" -} \ No newline at end of file +}