#include "BufferSpatialCalc.h"
#include "../Inputs/TerrainProcessor.h"

//References:
//Endreny, T., & Wood, E. F. (2003). Watershed Weighting of Export Coefficients to Map Critical Phosphorus Loading Areas. Journal of the American Water Resources Association, 39(1), 165-181. 
//Stephan, E. A., & Endreny, T. A. (2016). Weighting Nitrogen and Phosphorus Pixel Pollutant Loads to Represent Runoff and Buffering Likelihoods. JAWRA Journal of the American Water Resources Association, 52(2), 336-349. doi:doi:10.1111/1752-1688.12390

//Note: waters01 map assigns receiving waters 0, non-waters 1; waters0 map assigns receiving waters 0, non-waters cellsize_m
double BufferSpatialCalc::NODATA_code, BufferSpatialCalc::cellsize_m;
vector <double> BufferSpatialCalc::waters01, BufferSpatialCalc::waters0;
vector<double> BufferSpatialCalc::RI_surf, BufferSpatialCalc::RI_sub, BufferSpatialCalc::BI_surf, BufferSpatialCalc::BI_sub;
vector<double> BufferSpatialCalc::Manningn_vec, BufferSpatialCalc::Hydraulic_R_vec;
vector <double> BufferSpatialCalc::buffer_rise2, BufferSpatialCalc::buffer_rise1, BufferSpatialCalc::buffer_rise, BufferSpatialCalc::buffer_run1, BufferSpatialCalc::buffer_run, BufferSpatialCalc::buffer_slope;
vector <double> BufferSpatialCalc::FA_weighted_m;
int BufferSpatialCalc::nRows, BufferSpatialCalc::nCols;
//vector <double> waters0;
vector <double> BufferSpatialCalc::Time_surf1, BufferSpatialCalc::Time_surf2, BufferSpatialCalc::Time_surf;
vector <double> BufferSpatialCalc::Time_sub1, BufferSpatialCalc::Time_sub2, BufferSpatialCalc::Time_sub;
vector <double> BufferSpatialCalc::s_wtable_est;
vector<vector<double>>BufferSpatialCalc::slopeorganizer_temp;
vector <double> BufferSpatialCalc::p_release_vec, BufferSpatialCalc::p_trapping_vec;
vector <double> BufferSpatialCalc::filter1, BufferSpatialCalc::filter2, BufferSpatialCalc::filter, BufferSpatialCalc::filter_to_slope1, BufferSpatialCalc::filter_to_slope;
vector<long double>BufferSpatialCalc::ec_tn_vec, BufferSpatialCalc::ec_tp_vec, BufferSpatialCalc::ec_tss_vec;
vector<long double>BufferSpatialCalc::ec_tn_orig, BufferSpatialCalc::ec_tp_orig, BufferSpatialCalc::ec_tss_orig;
double BufferSpatialCalc::BI_surf_avg, BufferSpatialCalc::BI_sub_avg, BufferSpatialCalc::RI_surf_avg, BufferSpatialCalc::RI_sub_avg;
double BufferSpatialCalc::BI_surf_median, BufferSpatialCalc::BI_sub_median, BufferSpatialCalc::RI_surf_median, BufferSpatialCalc::RI_sub_median;
vector<double>BufferSpatialCalc::BI_surf_ratio, BufferSpatialCalc::BI_sub_ratio, BufferSpatialCalc::RI_surf_ratio, BufferSpatialCalc::RI_sub_ratio;
vector<double> BufferSpatialCalc::total_ec_tp_wtd, BufferSpatialCalc::total_ec_tn_wtd, BufferSpatialCalc::total_ec_tss_wtd;
vector<double> BufferSpatialCalc::total_emc_tp_wtd, BufferSpatialCalc::total_emc_tn_wtd, BufferSpatialCalc::total_emc_tss_wtd;
vector<double> BufferSpatialCalc::surf_ec_tp_wtd, BufferSpatialCalc::surf_ec_tn_wtd, BufferSpatialCalc::surf_ec_tss_wtd, BufferSpatialCalc::sub_ec_tp_wtd, BufferSpatialCalc::sub_ec_tn_wtd, BufferSpatialCalc::sub_ec_tss_wtd;
vector<int> BufferSpatialCalc::hotspot_id_temp;
vector<double> BufferSpatialCalc::hotspot_map_temp;
vector<tuple<int, int, int, double, double, double, double, double, double, double, double, double, double, double, int>> BufferSpatialCalc::hotspot_report_temp;
vector<int> BufferSpatialCalc::ContributingAreaPixels_map_writeout, BufferSpatialCalc::DispersalAreaPixels_map_writeout;
vector<double> BufferSpatialCalc::ca_concentration_tp_mg_L, BufferSpatialCalc::ca_concentration_tn_mg_L, BufferSpatialCalc::ca_concentration_tss_mg_L;
bool BufferSpatialCalc::flag_reduce_EC_EMC_after_GI = false;
long double  BufferSpatialCalc::sum_unweighted_ec_tn, BufferSpatialCalc::sum_unweighted_ec_tp, BufferSpatialCalc::sum_unweighted_ec_tss,
BufferSpatialCalc::sum_weighted_ec_tn, BufferSpatialCalc::sum_weighted_ec_tp, BufferSpatialCalc::sum_weighted_ec_tss;
vector<double> BufferSpatialCalc::slope_dem, BufferSpatialCalc::slope_wem, BufferSpatialCalc::slope_negdem, BufferSpatialCalc::slope_negwem;
vector<double> BufferSpatialCalc::groundwatertable_elevation1D, BufferSpatialCalc::neg_groundwatertable_elevation1D;
bool BufferSpatialCalc::isSubsurfaceCalculationEnabled;

//BufferSpatialCalc::ReceivingWater_ZeroValue_CreateMap_from_DEM function will assign 0 values to all receiving water pixels 
//Note: 2 versions of function. This one creates the waters01 map based on DEM analysis of flow accumulation exceeding a threshold
//Note: waters01 map assigns 0 to water and 1 to non-water; waters0 map assigns 0 to water and pixel cell size to non-water
void BufferSpatialCalc::ReceivingWater_ZeroValue_CreateMap_from_DEM() {
	nRows = Inputs::nRows;
	nCols = Inputs::nCols;
	NODATA_code = Inputs::NODATA_code;
	cellsize_m = Inputs::Length_Pixel_Side_m;
	// Ensure vectors are resized to match the dimensions of the map
	waters01.resize(nRows * nCols, NODATA_code);
	waters0.resize(nRows * nCols, NODATA_code); 

	//RiverThresh_FA = *max_element(TerrainProcessor::flowAccum1D.begin() , TerrainProcessor::flowAccum1D.end()) * 0.01; defined at 1% of FAccum
	double RiverThresh_FA = *max_element(TerrainProcessor::flowAccum1D.begin() , TerrainProcessor::flowAccum1D.end()) * 0.01;
	cout << "River networks defined at 1% of flow accumulation: " << RiverThresh_FA << endl;
	// Loop through all cells
	for (int i = 0; i < (nRows * nCols); i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (TerrainProcessor::flowAccum1D[i] > RiverThresh_FA) {
				// Assign river cells
				waters01[i] = 0;
				waters0[i] = 0;
			}
			else {
				// Assign non-river cells
				waters01[i] = 1;
				waters0[i] = cellsize_m;
			}
		}
	}
}

//BufferSpatialCalc::ReceivingWater_ZeroValue_CreateMap function will assign 0 values to all receiving water pixels 
//Note: 2 versions of function. This one reads in input->Waters0 which is a map based on National Hydrography Dataset (NHD) flowline and waterbodies 
//Note: waters01 map assigns 0 to water and 1 to non-water; waters0 map assigns 0 to water and pixel cell size to non-water
void BufferSpatialCalc::ReceivingWater_ZeroValue_CreateMap_from_NHD() {
	//auto it1 = input->Waters10.begin();
	auto it2 = Inputs::Waters0.begin();
	for (auto it = Inputs::Waters01.begin(); it != Inputs::Waters01.end(); it++) {
		if (*it != NODATA_code) {
			waters01.push_back(*it);
			waters0.push_back(*it2);
		}
		else {
			waters01.push_back(NODATA_code);
			waters0.push_back(NODATA_code);
		}
		//it1++;
		it2++;
	}
}


//BufferSpatialCalc::ManningRoughness_HydraulicRadius_CreateMap function reads Manning Roughness n and Hydraulic Radius data from inputs for each NLCD Class map input
//Note: Manning Roughness n values from Table 8.1 in Wurbs and James, Water Resources Engineering, ISBN-13: 978-0130812933, 2001, Pearson, NJ
//Note: Hydraulic Radius R values (ft) from Table 8.1 in Wurbs and James, Water Resources Engineering, ISBN-13: 978-0130812933, 2001, Pearson, NJ
//Reference: Wurbs, R. A., & James, W. P. (2001). Water Resources Engineering, 1st Edition, ISBN-13: 978-0130812933, Pages 840, Pearson, NJ
void BufferSpatialCalc::ManningRoughness_HydraulicRadius_CreateMap(Inputs* input) {
	for (double i : input->LandCover_NLCD_Class) {
		if (i != NODATA_code) {
			Manningn_vec.push_back(input->getManningn(i));
			Hydraulic_R_vec.push_back(input->getHydraulicradius(i));
		}
		else {
			Manningn_vec.push_back(NODATA_code);
			Hydraulic_R_vec.push_back(NODATA_code);
		}
	}
}

//BufferSpatialCalc::NutrientTrapping_CreateMap function creates maps of nutrient trapping based on mean (i.e., average) total phosphorus (P) release percentage (%) from the NLCD Class
//Note: p_release_vec (%) is the total phosphorus (P) the NLCD will release for a given amount that enters
//Note: p_trapping_vec (%) is the total phosphorus (P) the NLCD will trap, i.e., retain, for a given amount that enters
void BufferSpatialCalc::NutrientTrapping_CreateMap(Inputs* input) {
	double temp_val = 0;
	for (double i : input->LandCover_NLCD_Class) {
		if (i != NODATA_code) {
			temp_val = input->get_P_release_mean(i);
			p_release_vec.push_back(temp_val);
			p_trapping_vec.push_back(100- temp_val);
		}
		else {
			p_release_vec.push_back(NODATA_code);
			p_trapping_vec.push_back(NODATA_code);
		}
	}
}


//BufferSpatialCalc::Slope_DispersalArea_calc function will compute disperal area slopes for use in Buffer Index calculations of CADA theory
//Reference: Stephan, E. A., & Endreny, T. A. (2016). Weighting Nitrogen and Phosphorus Pixel Pollutant Loads to Represent Runoff and Buffering Likelihoods. JAWRA Journal of the American Water Resources Association, 52(2), 336-349. doi:doi:10.1111/1752-1688.12390
//Note: Calculation uses rise over run from local pixel to receiving water
//Note: CADA theory uses contributing area dispersal area flow analysis to identify zones with elevated runoff likelihood and buffer likelihood 
void BufferSpatialCalc::Slope_DispersalArea_calc(vector <double> slope) {
	buffer_rise1.clear();
	buffer_rise2.clear();
	buffer_rise.clear();
	buffer_run.clear();
	buffer_slope.clear();

	//start calculating buffer_rise
	double buffer_rise1_temp;
	for (int i = 0; i < (nRows * nCols); i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			//if (slope[i] != NODATA_code && waters01[i] != 0) where waters01 = 1 is land, not water
			if (slope[i] >= 0 && waters01[i] != 0) {
				buffer_rise1_temp = slope[i] * cellsize_m;
				buffer_rise1.push_back(buffer_rise1_temp);
			}
			else {
				buffer_rise1.push_back(0);
			}
		}
		else {
			buffer_rise1.push_back(NODATA_code);
		}
	}
	
	//apply the water01 map as a weight to flow accumulation to avoid collecting water pixel as -9999 during calculation
	TerrainProcessor::FlowAccum(waters01, buffer_rise1);
	//saving the 2D FA result to 1D buffer_rise vector
	for (auto row = TerrainProcessor::FlowAccumulation_Map.begin(); row != TerrainProcessor::FlowAccumulation_Map.end(); row++) {
		for (auto col = (*row).begin(); col != (*row).end(); col++) {
			buffer_rise2.push_back(*col);
		}
	}

	//exclude the nodata and water cells and save the FA result to buffer_rise
	//FA will be accumulate based on the original DEM map, rather than the elevation_m_1D (mask)
	for (int i = 0; i < (nRows * nCols); i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			//if (slope[i] != NODATA_code && waters01[i] != 0) where waters01 = 1 is land, not water
			if (buffer_rise2[i] >= 0 && waters01[i] != 0) {
				buffer_rise.push_back(buffer_rise2[i]);
			}
			else {
				buffer_rise.push_back(0);
			}
		}
		else {
			buffer_rise.push_back(NODATA_code);
		}
	}
	//finish calculating buffer rise
	
	//start calculating buffer_run
	//Determine the sum of the distance from each ridge cell (valleys in the negDem) to river cells (ridges in negDem)
	TerrainProcessor::FlowAccum(waters0);
	//saving the 2D FA result to1D vector buffer_run1
	for (auto row = TerrainProcessor::FlowAccumulation_Map.begin(); row != TerrainProcessor::FlowAccumulation_Map.end(); row++) {
		for (auto col = (*row).begin(); col != (*row).end(); col++) {
			buffer_run1.push_back(*col);
		}
	}
	//exclude the nodata and water cells and save the FA result to buffer_run
	//20231110 LI for those pixels missing NHDfdr data (-9999), FA result will be nodata
	for (int i = 0; i < (nRows * nCols); i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			//if (slope[i] != NODATA_code && waters01[i] != 0) where waters01 = 1 is land, not water
			if (buffer_run1[i] >= 0 && waters01[i] != 0) {
				buffer_run.push_back(buffer_run1[i]);
			}
			else {
				buffer_run.push_back(0);
			}
		}
		else {
			buffer_run.push_back(NODATA_code);
		}
	}

	//Determine the final slope for each pixel area
	for (int i = 0; i < (nRows * nCols); i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			//if (slope[i] != NODATA_code && waters01[i] != 0) where waters01 = 1 is land, not water
			if (buffer_run[i] > 0 && waters01[i] != 0) {
				//buffer_slope.push_back(buffer_rise[i] / buffer_run[i]) implements rise over run calculation of slope
				buffer_slope.push_back(buffer_rise[i] / buffer_run[i]);
			}
			else {
				buffer_slope.push_back(0);
			}
		}
		else {
			buffer_slope.push_back(NODATA_code);
		}
	}
}

//BufferSpatialCalc::filter_to_slope_calc function will compute disperal area slopes for use in Buffer Index calculations of CADA theory
//Reference: Endreny, T., & Wood, E. F. (2003). Watershed Weighting of Export Coefficients to Map Critical Phosphorus Loading Areas. Journal of the American Water Resources Association, 39(1), 165-181. 
//Note: CADA theory uses contributing area dispersal area flow analysis to identify zones with elevated runoff likelihood and buffer likelihood 
void BufferSpatialCalc::filter_to_slope_calc() {
	//double BI_surf_temp;
	filter_to_slope.clear();
	//Set P retention values to 0 in all cells that are classified as rivers
	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (waters01[i] != 0) {
				//note p_trapping_vec will have many 0 values for agricultural land 
				filter1.push_back(p_trapping_vec[i]);
			}
			else {
				filter1.push_back(0);
			}
		}
		else {
			filter1.push_back(NODATA_code);
		}
	}
	//Use the flow accumulation routine to determine the nutrient trapping efficiency for the downslope area
	//Note: the flow accumulation is using the flow direction with elevation_m_negative_1D
	//apply the water01 map as a weight to flow accumulation to avoid collecting water pixel as -9999 during calculation
	TerrainProcessor::FlowAccum(waters01, filter1);
	//saving the FA result to filter2
	filter2.clear();
	for (auto row = TerrainProcessor::FlowAccumulation_Map.begin(); row != TerrainProcessor::FlowAccumulation_Map.end(); row++) {
		for (auto col = (*row).begin(); col != (*row).end(); col++) {
			filter2.push_back(*col);
		}
	}

	//Use the con command to set a minimum trapping value of greater than or equal to 1 for all non river cells (add 1)
	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (filter2[i] >= 0 && waters01[i] != 0) {
				filter.push_back(filter2[i]);
			}
			else {
				filter.push_back(0);
			}
		}
		else {
			filter.push_back(NODATA_code);
		}
	}

	//Calculate slope vector (filtered/hillslope)
	double filter_to_slope_temp;
	filter_to_slope1.clear();
	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (waters01[i] != 0 && buffer_slope[i] > 0) {
				filter_to_slope_temp = filter[i] / buffer_slope[i];
				filter_to_slope1.push_back(filter_to_slope_temp);
			}
			else {
				filter_to_slope1.push_back(0);
			}
		}
		else {
			filter_to_slope1.push_back(NODATA_code);
		}
	}
	//find the median positive filter_to_slope 
	//replace any value < minimun positive filter_to_slope with minimun positive filter_to_slope to avoid filter_to_slope = 0
	//filter_to_slope = 0  -> BI=log(0+1)=0 -> BI_ratio = BI_avg/BI =BI_avg/0
	double median_filter_to_slope = calculateMedian(filter_to_slope1);
	double min_filter_to_slope = calculateMedian(filter_to_slope1);
	filter_to_slope.clear();
	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (filter_to_slope1[i] > 0 && waters01[i] != 0) {
				if (filter_to_slope1[i] < min_filter_to_slope) {
					filter_to_slope.push_back(median_filter_to_slope);
				}
				else {
					filter_to_slope.push_back(filter_to_slope1[i]);
				}
			}
			else {
				filter_to_slope.push_back(0);
			}
		}
		else {
			filter_to_slope.push_back(NODATA_code);
		}
	}
}

//BufferSpatialCalc::Depth_GroundwaterTable_calc function computes depth to the groundwater table based on Eq 13 of Stephan and Endreny (2016)
//Reference: Stephan, E. A., & Endreny, T. A. (2016). Weighting Nitrogen and Phosphorus Pixel Pollutant Loads to Represent Runoff and Buffering Likelihoods. JAWRA Journal of the American Water Resources Association, 52(2), 336-349. doi:doi:10.1111/1752-1688.12390
//Note: SSURGO data provided S_WTable_avg(cm) and f_param, DEM data provided RI_sub_iand RI_sub_avg, aka topographic index
void BufferSpatialCalc::Depth_GroundwaterTable_calc(Inputs* input) {
	double s_wtable_est_temp;
	s_wtable_est.clear();

	//compute STI first using Endreny 1999 Eq 1
	vector <double> STI;
	double STI_temp, exp_temp, STI_avg;
	double sum = 0;
	int count = 0;
	
	for (int i = 0; i < nRows * nCols; i++) {
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			//if (waters01[i] != 0 && (slope_dem[i] > 0 && input->s_trans[i] > 0) && (input->s_trans_avg > 0 &&  FA_weighted_m[i] > 0)) 
			//Note: need to be > 0 to avoid negative result after log
			if (waters01[i] != 0 && (slope_dem[i] > 0 && input->s_trans[i] > 0) && (input->s_trans_avg > 0 &&  FA_weighted_m[i] > 0)) {
				//exp_temp = (input->s_trans_avg * FA_weighted_m[i]) / (slope_dem[i] * input->s_trans[i]) + 1
				//Note: log() function in C++ returns the natural logarithm (base-e logarithm); + 1 to avoid  0 < FA/Slope_Pixel_mpm < 1
				exp_temp = (input->s_trans_avg * FA_weighted_m[i]) / (slope_dem[i] * input->s_trans[i]) + 1;
				STI_temp = log(exp_temp);
				STI.push_back(STI_temp);

				if (STI_temp >= 0) {
					sum += STI_temp;
					count++;
				}
			}
			else {
				STI.push_back(0);
			}
		}
		else {
			//STI.push_back(NODATA_code); 
			STI.push_back(NODATA_code);
		}
	}
	STI_avg = sum / count;

	//start compute watertable depth
	for (int i = 0; i < nRows * nCols; i++) {
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (waters01[i] != 0 && input->BufferSpatialWeightingNumericalParams["f_param"] > 0) {
				s_wtable_est_temp = input->BufferSpatialWeightingNumericalParams["S_WTable_avg"] - (1 / input->BufferSpatialWeightingNumericalParams["f_param"]) * (STI[i] - STI_avg);
				//Dividing by 100 to convert from cm to m
				s_wtable_est.push_back(s_wtable_est_temp / 100);
			}
			else {
				s_wtable_est.push_back(0);
			}
		}
		else {
			s_wtable_est.push_back(NODATA_code);
		}
	}
}

//BufferSpatialCalc::TravelTime_Runoff_Surface_calc function will compute runoff velocity subsurface based on Eq 9 of Stephan and Endreny (2016) with Manning velocity to estimate likelihood for nutrient removal 
//Reference: Stephan, E. A., & Endreny, T. A. (2016). Weighting Nitrogen and Phosphorus Pixel Pollutant Loads to Represent Runoff and Buffering Likelihoods. JAWRA Journal of the American Water Resources Association, 52(2), 336-349. doi:doi:10.1111/1752-1688.12390
void BufferSpatialCalc::TravelTime_Runoff_Surface_calc(vector <double> slope) {
	double R_23, S_12, Vel_surf, t_surf;
	// Clear and resize Time_surf1 to match the map dimensions
	Time_surf1.clear();
	Time_surf1.resize(nRows * nCols);
	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			//Manningn_vec[i] can't be 0 and slope can't be negative value
			if (waters01[i] != 0 && Manningn_vec[i] > 0 && slope[i] > 0) {
				S_12 = pow(slope[i], 0.5);
				R_23 = pow(Hydraulic_R_vec[i], two_three);
				Vel_surf = (1.0 / Manningn_vec[i]) * R_23 * S_12;
				if (Vel_surf > 0) {
					t_surf = cellsize_m / Vel_surf;
					Time_surf1[i] = t_surf;
				}
				else {
					Time_surf1[i] = 0;
				}
			}
			else {
				Time_surf1[i] = 0;
			}
		}
		else {
			Time_surf1[i] = NODATA_code;
		}
	}
	//Compute Time_surf_t as total accumulated travel time across the dispersal area, add in local travel time to pixel
	//apply the water01 map as a weight to flow accumulation to avoid collecting water pixel as -9999 during calculation
	TerrainProcessor::FlowAccum(waters01, Time_surf1);
	Time_surf2.clear();
	//saving the FA result to Time_surf2
	for (auto row = TerrainProcessor::FlowAccumulation_Map.begin(); row != TerrainProcessor::FlowAccumulation_Map.end(); row++) {
		for (auto col = (*row).begin(); col != (*row).end(); col++) {
			Time_surf2.push_back(*col);
		}
	}

	for (int i = 0; i < (nRows * nCols); i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (Time_surf2[i] > 0 && waters01[i] != 0) {
				//Eq 10 in Stephan and Endreny (2016)
				Time_surf.push_back(waters01[i] * Time_surf2[i]);
			}
			else {
				Time_surf.push_back(0);
			}
		}
		else {
			Time_surf.push_back(NODATA_code);
		}
	}
}

//BufferSpatialCalc::TravelTime_Runoff_Subsurface_calc function will compute runoff velocity subsurface based on Eq 12 of Stephan and Endreny (2016) with Darcy flux to estimate likelihood for nutrient removal 
//Reference: Stephan, E. A., & Endreny, T. A. (2016). Weighting Nitrogen and Phosphorus Pixel Pollutant Loads to Represent Runoff and Buffering Likelihoods. JAWRA Journal of the American Water Resources Association, 52(2), 336-349. doi:doi:10.1111/1752-1688.12390
void BufferSpatialCalc::TravelTime_Runoff_Subsurface_calc(Inputs* input, vector <double> slope) {
	double Vel_sub, t_sub;
	Time_sub1.clear();
	Time_sub1.resize(nRows * nCols);

	//auto it_slope = slope.begin();
	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (waters01[i] != 0) {
				//Unit conversion for hydraulic conductivity (um/s) multiplied by 1E-6 m/um
				Vel_sub = slope[i] * input->s_Ksat[i] * 1e-6;
				//Setting Vel_sub to a minimum 1e-5m/s or 0.8m/day to avoid creating null data values in Time_sub
				if (Vel_sub < 1e-5) {
					Vel_sub = 1e-5;
				}
				if (Vel_sub > 0) {
					//Compute t_sub (s) subsurface time for runoff as RasterLength_m / subsurface velocity
					t_sub = cellsize_m / Vel_sub;
					Time_sub1[i] = t_sub;
				}
				else {
					Time_sub1[i] = 0;
				}
			}
			else {
				Time_sub1[i] = 0;
			}
		}
		else {
			Time_sub1[i] = NODATA_code;
		}
	}

	//Compute Time_surf_t as total accumulated travel time across the dispersal area, add in local travel time to pixel
	//apply the water01 map as a weight to flow accumulation to avoid collecting water pixel as -9999 during calculation
	TerrainProcessor::FlowAccum(waters01,Time_sub1);

	Time_sub2.clear();
	//saving the FA result to Time_surf2
	for (auto row = TerrainProcessor::FlowAccumulation_Map.begin(); row != TerrainProcessor::FlowAccumulation_Map.end(); row++) {
		for (auto col = (*row).begin(); col != (*row).end(); col++) {
			Time_sub2.push_back(*col);
		}
	}

	for (int i = 0; i < (nRows * nCols); i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (Time_sub2[i] > 0 && waters01[i] != 0) {
				//Eq 10 in Stephan and Endreny (2016)
				Time_sub.push_back(Time_sub2[i]);
			}
			else {
				Time_sub.push_back(0);
			}
		}
		else {
			Time_sub.push_back(NODATA_code);
		}
	}
}

// Function to calculate the median of a vector with validation
double BufferSpatialCalc::calculateMedian(vector<double> map) {
	// Remove NODATA_code values
	map.erase(remove_if(map.begin(), map.end(), [](double v) { return v == NODATA_code; }), map.end());
	int size = map.size();
	if (size == 0) return 0.0; // Handle empty vector case

	// Sort the values
	sort(map.begin(), map.end());

	// Calculate the median
	if (size % 2 == 0) {
		return (map[size / 2 - 1] + map[size / 2]) / 2.0;
	}
	else {
		return map[size / 2];
	}
}
//BufferSpatialCalc::RunoffIndex_Surface_calc function will create Runoff Index surface values by representing upslope features of CADA theory
//Note: CADA theory uses contributing area dispersal area flow analysis to identify zones with elevated runoff likelihood and buffer likelihood 
//Note: Compute runoff index, RI_surf_i in Eq 6 of Stephan and Endreny (2016), and from impervious or Urban Topographic Index, Eq 8 of Valeo, C. et al. (2000)
//Reference: Valeo, C., & Moin, S. M. A. (2000). Variable source area modelling in urbanizing watersheds. Journal of Hydrology, 228(1-2), 68-81. 
//Reference: Stephan, E. A., & Endreny, T. A. (2016). Weighting Nitrogen and Phosphorus Pixel Pollutant Loads to Represent Runoff and Buffering Likelihoods. JAWRA Journal of the American Water Resources Association, 52(2), 336-349. doi:doi:10.1111/1752-1688.12390
void BufferSpatialCalc::RunoffIndex_Surface_calc(vector<double> FA, vector<double> Slope_Pixel_mpm) {
	double RI_surf_temp, exp_temp;
	double sum = 0;
	int count = 0;

	RI_surf.clear();

	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (FA[i] >= 0 && waters01[i] != 0 && Slope_Pixel_mpm[i] > 0) {
				//log() function in C++ returns the natural logarithm (base-e logarithm) 
				exp_temp = FA[i] / Slope_Pixel_mpm[i] + 1;
				RI_surf_temp = log(exp_temp);
				//Rescale the RI_surf_int distribution to have a minimum of 1
				//which will remove negatives and normalize indexes for cross-comparison 
				//if (RI_surf_temp < 0.0001)  RI_surf_temp = 0.0001;
				RI_surf.push_back(RI_surf_temp);
				sum += RI_surf_temp;
				count++;
			}
			else {
				//Con(IsNull('RI_surf_int'),0,'RI_surf_int')
				RI_surf.push_back(0);
			}
		}
		else {
			//Con(IsNull('RI_surf_int'),0,'RI_surf_int')
			RI_surf.push_back(NODATA_code);
		}
	}
	RI_surf_avg = sum / count;
	RI_surf_median = calculateMedian(RI_surf);
	cout << "Surface Runoff Index Mean: " << RI_surf_avg << "; Median: " << RI_surf_median <<  endl;
}

//BufferSpatialCalc::RunoffIndex_Subsurface_calc function will create Runoff Index subsurface values by representing upslope features of CADA theory
//Note: CADA theory uses contributing area dispersal area flow analysis to identify zones with elevated runoff likelihood and buffer likelihood 
//Note: Compute RI_sub in Eq 7 of Stephan and Endreny (2016), from Soil Topographic Index Eq 10 in Sivapalan et al. (1987) and Eq 1 in Endreny et al. (2000) 
//Note: Eq 7 uses average and local transmissivity, S_Trans_avg and S_Trans, and flow accumulation, FA_per, and DEM_slope, where FA_per = a and DEM_slope = tan(B)
//Note: Eq 7 separated into 2 Eqs, RI_sub = ln[(S_Trans_avg * FA_per) / (S_Trans * DEM_slope)] = ln(FA_per / DEM_slope) + ln(S_Trans_avg / S_Trans)
//Reference: Sivapalan, M., Beven, K., & Wood, E. F. (1987). On hydrologic similarity: 2. A scaled model of storm runoff production. Water Resources Research, 23(12), 2266-2278. doi:10.1029/WR023i012p02266
//Reference: Endreny, T. A., Wood, E. F., & Lettenmaier, D. P. (2000). Satellite‐derived digital elevation model accuracy: hydrogeomorphological analysis requirements. Hydrological Processes, 14(1), 1-20. doi:doi:10.1002/(SICI)1099-1085(200001)14:1<1::AID-HYP918>3.0.CO;2-#
//Reference: Stephan, E. A., & Endreny, T. A. (2016). Weighting Nitrogen and Phosphorus Pixel Pollutant Loads to Represent Runoff and Buffering Likelihoods. JAWRA Journal of the American Water Resources Association, 52(2), 336-349. doi:doi:10.1111/1752-1688.12390
void BufferSpatialCalc::RunoffIndex_Subsurface_calc(Inputs* input, vector<double> FA, vector<double> Slope_Pixel_mpm) {
	double RI_sub_temp, exp_temp;
	double sum = 0;
	int count = 0;
	RI_sub.clear();

	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (waters01[i] != 0 && Slope_Pixel_mpm[i] > 0 && input->s_trans[i] > 0 && FA[i] >= 0) {
				//when FA -> 0, log will be -> -inf
				//log() function in C++ returns the natural logarithm (base-e logarithm) 
				//+ 1 to avoid  0 < FA/Slope_Pixel_mpm < 1
				exp_temp = (input->s_trans_avg * FA[i]) / (Slope_Pixel_mpm[i] * input->s_trans[i]) + 1;

				RI_sub_temp = log(exp_temp);
				RI_sub.push_back(RI_sub_temp);
				sum += RI_sub_temp;
				count++;
			}
			else {
				//Con(IsNull('RI_sub_int'), 0, 'RI_sub_int')
				RI_sub.push_back(0);
			}
		}
		else {
			RI_sub.push_back(NODATA_code);
		}
	}
	RI_sub_avg = sum / count;
	RI_sub_median = calculateMedian(RI_sub);
	cout << "Subsurface Runoff Index Mean: " << RI_sub_avg << "; Median: " << RI_sub_median << endl;
}

//BufferSpatialCalc::BufferIndex_Subsurface_calc function will create Buffer Index subsurface values by representing downslope features of CADA theory
//Note: CADA theory uses contributing area dispersal area flow analysis to identify zones with elevated runoff likelihood and buffer likelihood 
void BufferSpatialCalc::BufferIndex_Subsurface_calc(vector<double>exp) {
	double BI_sub_temp, exp_temp;
	double sum = 0;
	int count = 0;
	BI_sub.clear();

	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (waters01[i] != 0) {
				//log() function in C++ returns the natural logarithm (base-e logarithm) 
				//+ 1 to avoid  0 < FA/slope < 1
				exp_temp = exp[i] + 1;
				BI_sub_temp = log(exp_temp);
				BI_sub.push_back(BI_sub_temp);
				sum += BI_sub_temp;
				count++;
			}
			else {
				BI_sub.push_back(0);
			}
		}
		else {
			BI_sub.push_back(NODATA_code);
		}
	}

	BI_sub_avg = sum / count;
	BI_sub_median = calculateMedian(BI_sub);
	cout << "Subsurface Runoff Index Mean: " << BI_sub_avg << "; Median: " << BI_sub_median << endl;
}

//BufferSpatialCalc::BufferIndex_Surface_calc function will create Buffer Index suurface values by representing downslope features of CADA theory
//Note: CADA theory uses contributing area dispersal area flow analysis to identify zones with elevated runoff likelihood and buffer likelihood 
void BufferSpatialCalc::BufferIndex_Surface_calc(vector<double>exp) {
	double BI_surf_temp, exp_temp;
	double sum = 0;
	int count = 0;
	BI_surf.clear();

	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (waters01[i] != 0) {
				//log() function in C++ returns the natural logarithm (base-e logarithm) 
				//+ 1 to avoid  0 < FA/slope < 1
				exp_temp = exp[i] + 1;
				BI_surf_temp = log(exp_temp);
				BI_surf.push_back(BI_surf_temp);
				sum += BI_surf_temp;
				count++;
			}
			else {
				BI_surf.push_back(0);
			}
		}
		else {
			BI_surf.push_back(NODATA_code);
		}
	}
	BI_surf_avg = sum / count;
	BI_surf_median = calculateMedian(BI_surf);
	cout << "Surface Buffer Index Mean: " << BI_surf_avg << "; Median: " << BI_surf_median << endl;
}




//BufferSpatialCalc::ExportCoefficient_Retrieve_50thPercentile function will retrieve 50th percentile (median) Export Coefficient values from USDA White et al. (2015) database
double BufferSpatialCalc::ExportCoefficient_Retrieve_50thPercentile(Inputs* input, string pollutant_type, double lc) {
	double ec_50th;
	if (lc == 21 || lc == 22 || lc == 23 || lc == 24) { //21 - 24 Urban
		if (pollutant_type == "TN") {
			ec_50th = input->EC_vec_map["Urban_TN"][3];//50th, kg/ha
		}
		if (pollutant_type == "TP") {
			ec_50th = input->EC_vec_map["Urban_TP"][3];//50th, kg/ha
		}
		if (pollutant_type == "TSS") {
			ec_50th = input->EC_vec_map["Urban_TSS"][3];//50th, kg/ha
		}
	}
	else if (lc == 41 || lc == 42 || lc == 43 || lc == 51 || lc == 52) {// LC 52 Shrub/Scrub ->LC 41 Deciduous Forest; LC 51 Dwarf scrub ->LC 41 Deciduous Forest
		if (pollutant_type == "TN") {
			ec_50th = input->EC_vec_map["Forest_TN"][3];//50th, kg/ha
		}
		if (pollutant_type == "TP") {
			ec_50th = input->EC_vec_map["Forest_TP"][3];//50th, kg/ha
		}
		if (pollutant_type == "TSS") {
			ec_50th = input->EC_vec_map["Forest_TSS"][3];//50th
		}
	}
	else if (lc == 71 || lc == 72 || lc == 73 || lc == 74 || lc == 31) {// LC 31 Barren Land -> LC71 Grassland; LC 72 Sedge/Herbaceous;LC 73 Lichens; LC 74 Moss
		if (pollutant_type == "TN") {
			ec_50th = input->EC_vec_map["Rangeland_TN"][3];//50th
		}
		if (pollutant_type == "TP") {
			ec_50th = input->EC_vec_map["Rangeland_TP"][3];//50th
		}
		if (pollutant_type == "TSS") {
			ec_50th = input->EC_vec_map["Rangeland_TSS"][3];//50th
		}
	}
	else if (lc == 81) {
		if (pollutant_type == "TN") {
			ec_50th = input->EC_vec_map["Pasture_hay_TN"][3];//50th
		}
		if (pollutant_type == "TP") {
			ec_50th = input->EC_vec_map["Pasture_hay_TP"][3];//50th
		}
		if (pollutant_type == "TSS") {
			ec_50th = input->EC_vec_map["Pasture_hay_TSS"][3];//50th
		}
	}
	else if (lc == 82) {
		if (pollutant_type == "TN") {
			ec_50th = input->EC_vec_map["Cropland_TN"][3];//50th
		}
		if (pollutant_type == "TP") {
			ec_50th = input->EC_vec_map["Cropland_TP"][3];//50th
		}
		if (pollutant_type == "TSS") {
			ec_50th = input->EC_vec_map["Cropland_TSS"][3];//50th
		}
	}
	else if (lc == 11 || lc == 95 || lc == 90) { // waters
		if (pollutant_type == "TN") {
			ec_50th = input->EC_vec_map["Water_Wetlands_TN"][3];//50th
		}
		if (pollutant_type == "TP") {
			ec_50th = input->EC_vec_map["Water_Wetlands_TP"][3];//50th
		}
		if (pollutant_type == "TSS") {
			ec_50th = input->EC_vec_map["Water_Wetlands_TSS"][3];//50th
		}
	}
	return ec_50th;
}

//BufferSpatialCalc::WaterVolume_Retrieve_50thPercentile function will retrieve 50th percentile (median) Water Volume values from USDA White et al. (2015) database
double BufferSpatialCalc::WaterVolume_Retrieve_50thPercentile(Inputs* input, double lc) {
	double wv_50th;
	//check the water volume type for EMC calculation
	string WaterVolume_Type = input->BufferNutrientsStringParams["WaterVolume_Type"];
	//vector<string> LC_name = { "Urban","Forest", "Rangeland", "Pasture_hay", "Cropland", "Water_Wetlands" };
	/*if (WaterVolume_Type == "1") {//option 1 is Water Yield
		flowmapkey = LC + "_" + "WaterYield";
	}
	if (WaterVolume_Type == "2") { //option 2 is surface runoff
		flowmapkey = LC + "_" + "SurfaceRunoff";
	}*/

	if (lc == 21 || lc == 22 || lc == 23 || lc == 24) { //21 - 24 Urban
		if (WaterVolume_Type == "1") {
			wv_50th = input->water_volume_vec_map["Urban_WaterYield"][3] * Ratio_mmHa_to_L;//50th; L
		}
		if (WaterVolume_Type == "2") {
			wv_50th = input->water_volume_vec_map["Urban_SurfaceRunoff"][3] * Ratio_mmHa_to_L;//50th; L
		}
	}
	else if (lc == 41 || lc == 42 || lc == 43 || lc == 51 || lc == 52) {// LC 52 Shrub/Scrub ->LC 41 Deciduous Forest; LC 51 Dwarf scrub ->LC 41 Deciduous Forest
		if (WaterVolume_Type == "1") {
			wv_50th = input->water_volume_vec_map["Forest_WaterYield"][3] * Ratio_mmHa_to_L;//50th; L
		}
		if (WaterVolume_Type == "2") {
			wv_50th = input->water_volume_vec_map["Forest_SurfaceRunoff"][3] * Ratio_mmHa_to_L;//50th; L
		}
	}
	else if (lc == 71 || lc == 72 || lc == 73 || lc == 74 || lc == 31) {// LC 31 Barren Land -> LC71 Grassland; LC 72 Sedge/Herbaceous;LC 73 Lichens; LC 74 Moss
		if (WaterVolume_Type == "1") {
			wv_50th = input->water_volume_vec_map["Rangeland_WaterYield"][3] * Ratio_mmHa_to_L;//50th; L
		}
		if (WaterVolume_Type == "2") {
			wv_50th = input->water_volume_vec_map["Rangeland_SurfaceRunoff"][3] * Ratio_mmHa_to_L;//50th; L
		}
	}
	else if (lc == 81) {
		if (WaterVolume_Type == "1") {
			wv_50th = input->water_volume_vec_map["Pasture_hay_WaterYield"][3] * Ratio_mmHa_to_L;//50th; L
		}
		if (WaterVolume_Type == "2") {
			wv_50th = input->water_volume_vec_map["Pasture_hay_SurfaceRunoff"][3] * Ratio_mmHa_to_L;//50th; L
		}
	}
	else if (lc == 82) {
		if (WaterVolume_Type == "1") {
			wv_50th = input->water_volume_vec_map["Cropland_WaterYield"][3] * Ratio_mmHa_to_L;//50th; L
		}
		if (WaterVolume_Type == "2") {
			wv_50th = input->water_volume_vec_map["Cropland_SurfaceRunoff"][3] * Ratio_mmHa_to_L;//50th; L
		}
	}
	else if (lc == 11 || lc == 95 || lc == 90) { // waters lc doesn't has volume 
		wv_50th = numeric_limits<double>::infinity();
	}
	return wv_50th;
}


//BufferSpatialCalc::RunoffIndex_BufferIndex_Surface_Ratio_calc function will compute the runoff index and buffer index surface ratios using CADA theory
//Note: CADA theory uses contributing area dispersal area flow analysis to identify zones with elevated runoff likelihood and buffer likelihood 
void BufferSpatialCalc::RunoffIndex_BufferIndex_Surface_Ratio_calc() {
	BI_surf_ratio.clear();
	RI_surf_ratio.clear();

	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (BI_surf[i] > 0 && RI_surf[i] > 0 && RI_surf_median > 0 && waters01[i] != 0) {
				//LI 11232024 for a AOI with most of the land cover type is agriculture (e.g. Sandusky basin has around 80% agriculture land),
				// which has trapping efficiency 0, the median value will be 0
				// to avoid the BI_ratio to be 0, use average
				// Calculate BI_surf_ratio and RI_surf_ratio using the avg
				double BI_ratio = BI_surf_median / BI_surf[i];//CADA2003 (Endrent & Wood) use median, CADA2016 (Stephen) use average
				double RI_ratio = RI_surf[i] / RI_surf_median;//CADA2003 (Endrent & Wood) use median, CADA2016 (Stephen) use average
				BI_surf_ratio.push_back(BI_ratio);
				RI_surf_ratio.push_back(RI_ratio);
			}
			else {
				BI_surf_ratio.push_back(0);
				RI_surf_ratio.push_back(0);
			}
		}
		else {
			BI_surf_ratio.push_back(NODATA_code);
			RI_surf_ratio.push_back(NODATA_code);
		}
	}
}


//BufferSpatialCalc::RunoffIndex_BufferIndex_Subsurface_Ratio_calc function will compute the runoff index and buffer index subsurface ratios using CADA theory
//Note: CADA theory uses contributing area dispersal area flow analysis to identify zones with elevated runoff likelihood and buffer likelihood 
void BufferSpatialCalc::RunoffIndex_BufferIndex_Subsurface_Ratio_calc() {
	BI_sub_ratio.clear();
	RI_sub_ratio.clear();
	//sum_of_RIBI_products_sub = 0;

	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (BI_sub[i] > 0 && RI_sub[i] > 0 && RI_sub_median > 0 && waters01[i] != 0) {
				//LI 11232024 for a AOI with most of the land cover type is agriculture (e.g. Sandusky basin has around 80% agriculture land),
				// which has trapping efficiency 0, the median value will be 0
				// to avoid the BI_ratio to be 0, use average
				// Calculate BI_surf_ratio and RI_surf_ratio using the avg
				double BI_ratio = BI_sub_median / BI_sub[i];//CADA2003 (Endrent & Wood) use median, CADA2016 (Stephen) use average
				double RI_ratio = RI_sub[i] / RI_sub_median;//CADA2003 (Endrent & Wood) use median, CADA2016 (Stephen) use average
				BI_sub_ratio.push_back(BI_ratio);
				RI_sub_ratio.push_back(RI_ratio);
			}
			else {
				BI_sub_ratio.push_back(0);
				RI_sub_ratio.push_back(0);
			}
		}
		else {
			BI_sub_ratio.push_back(NODATA_code);
			RI_sub_ratio.push_back(NODATA_code);
		}
	}
}


//BufferSpatialCalc::ExportCoefficient_unweighted_CreateMap function will create Export Coefficient maps that are unweighted, not representing CADA theory
//Note: CADA theory uses contributing area dispersal area flow analysis to identify zones with elevated runoff likelihood and buffer likelihood 
void BufferSpatialCalc::ExportCoefficient_unweighted_CreateMap(Inputs* input) {
	ec_tn_orig.clear();
	ec_tp_orig.clear();
	ec_tss_orig.clear();
	ec_tn_orig.resize(nRows * nCols, NODATA_code);
	ec_tp_orig.resize(nRows * nCols, NODATA_code);
	ec_tss_orig.resize(nRows * nCols, NODATA_code);

	// Loop through all cells
	for (int i = 0; i < (nRows * nCols); i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			ec_tn_orig[i] = ExportCoefficient_Retrieve_50thPercentile(input, "TN", input->LandCover_NLCD_Class[i]);
			ec_tp_orig[i] = ExportCoefficient_Retrieve_50thPercentile(input, "TP", input->LandCover_NLCD_Class[i]);
			ec_tss_orig[i] = ExportCoefficient_Retrieve_50thPercentile(input, "TSS", input->LandCover_NLCD_Class[i]);
		}
		else {
			ec_tn_orig[i] = NODATA_code;
			ec_tp_orig[i] = NODATA_code;
			ec_tss_orig[i] = NODATA_code;
		}
	}
}

// function to calculate surface EC sums
void BufferSpatialCalc::ExportCoefficient_Unweighted_sum_calc() {
	// Initialize temporary variables
	sum_unweighted_ec_tn = 0, sum_unweighted_ec_tp = 0, sum_unweighted_ec_tss = 0;

	// Compute original and weighted sums
	for (int i = 0; i < nRows * nCols; i++) {
		if (BI_surf_ratio[i] != NODATA_code && RI_surf_ratio[i] != NODATA_code && waters01[i] != 0) {
			sum_unweighted_ec_tn += ec_tn_orig[i];
			sum_unweighted_ec_tp += ec_tp_orig[i];
			sum_unweighted_ec_tss += ec_tss_orig[i];
		}
	}
}

// function to calculate surface EC sums
void BufferSpatialCalc::ExportCoefficient_Weighted_sum_calc(Inputs* input, bool isSubsurfaceCalculationEnabled) {
	// Initialize temporary variables
	sum_weighted_ec_tn = 0, sum_weighted_ec_tp = 0, sum_weighted_ec_tss = 0;

	// Set surface and subsurface fractions
	double surface_frac_tn = 1.0, surface_frac_tp = 1.0, surface_frac_tss = 1.0;
	if (isSubsurfaceCalculationEnabled) {
		surface_frac_tn = input->BufferSpatialWeightingNumericalParams["ECTN_surface_contribution_frac"];
		surface_frac_tp = input->BufferSpatialWeightingNumericalParams["ECTP_surface_contribution_frac"];
		surface_frac_tss = input->BufferSpatialWeightingNumericalParams["ECTSS_surface_contribution_frac"];
	}

	// Compute original and weighted sums
	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (waters01[i] != 0) {
				long double surf_ratio_product_temp = BI_surf_ratio[i] * RI_surf_ratio[i];
				long double sub_ratio_product_temp = (isSubsurfaceCalculationEnabled && BI_sub_ratio[i] != NODATA_code && RI_sub_ratio[i] != NODATA_code)
					? BI_sub_ratio[i] * RI_sub_ratio[i] : 0;
				// Apply surface and subsurface fractions from user input
				// // Calculate weighted terms with surface and subsurface fractions
				long double weighted_ec_tn = (surf_ratio_product_temp * surface_frac_tn + sub_ratio_product_temp * (1 - surface_frac_tn)) * ec_tn_orig[i];
				long double weighted_ec_tp = (surf_ratio_product_temp * surface_frac_tp + sub_ratio_product_temp * (1 - surface_frac_tp)) * ec_tp_orig[i];
				long double weighted_ec_tss = (surf_ratio_product_temp * surface_frac_tss + sub_ratio_product_temp * (1 - surface_frac_tss)) * ec_tss_orig[i];
				// Check for NaN and inf values, skipping if detected
				if (!isnan(weighted_ec_tn) && !isinf(weighted_ec_tn)) sum_weighted_ec_tn += weighted_ec_tn;
				if (!isnan(weighted_ec_tp) && !isinf(weighted_ec_tp)) sum_weighted_ec_tp += weighted_ec_tp;
				if (!isnan(weighted_ec_tss) && !isinf(weighted_ec_tss)) sum_weighted_ec_tss += weighted_ec_tss;
			}
		}
	}
}


//BufferSpatialCalc::ExportCoefficient_Weighted_Subsurface_calc function will write the CADA weighted Export Coefficient subsurface values
//Note: CADA theory uses contributing area dispersal area flow analysis to identify zones with elevated runoff likelihood and buffer likelihood 



void BufferSpatialCalc::Normalize_EC_EMC_maps(Inputs* input, bool isSubsurfaceCalculationEnabled) {

	// Clear existing vectors before recalculating
	total_ec_tn_wtd.clear();
	total_ec_tp_wtd.clear();
	total_ec_tss_wtd.clear();
	total_emc_tn_wtd.clear();
	total_emc_tp_wtd.clear();
	total_emc_tss_wtd.clear();

	surf_ec_tn_wtd.clear();
	surf_ec_tp_wtd.clear();
	surf_ec_tss_wtd.clear();
	sub_ec_tn_wtd.clear();
	sub_ec_tp_wtd.clear();
	sub_ec_tss_wtd.clear();

	// Step 1: Calculate unweighted and weighted sums
	ExportCoefficient_Unweighted_sum_calc();//sum_unweighted_ec_tn
	ExportCoefficient_Weighted_sum_calc(input, isSubsurfaceCalculationEnabled); //sum_weighted_ec_tn

	// Step 2: Calculate combined normalization factors
	double  Observed_Load_TP_kg_year = input->BufferSpatialWeightingNumericalParams["Observed_Load_TP_kg_year"];
	double  Observed_Load_TN_kg_year = input->BufferSpatialWeightingNumericalParams["Observed_Load_TN_kg_year"];
	double  Observed_Load_TSS_kg_year = input->BufferSpatialWeightingNumericalParams["Observed_Load_TSS_kg_year"];

	long double normalization_factor_for_weightedmap_tn = 1.0, normalization_factor_for_weightedmap_tp = 1.0, normalization_factor_for_weightedmap_tss = 1.0;
	long double normalization_factor_for_unweightedmap_tn = 1.0, normalization_factor_for_unweightedmap_tp = 1.0, normalization_factor_for_unweightedmap_tss = 1.0;
	long double observed_total_ec_tn = 0, observed_total_ec_tp = 0, observed_total_ec_tss = 0;

	// Step 3
	// Calculate TN normalization factor
	if (Observed_Load_TN_kg_year > 0) {
		observed_total_ec_tn = Observed_Load_TN_kg_year / (cellsize_m * cellsize_m / 10000.0);//EC1 (kg/ha) * 900/10000 + EC2 (kg/ha) * 900/10000 + ...ECn (kg/ha) * 900/10000  = sum(EC)*900/10000 = total load (kg)
		normalization_factor_for_weightedmap_tn = (sum_weighted_ec_tn != 0) ? observed_total_ec_tn / sum_weighted_ec_tn : 1.0;
		normalization_factor_for_unweightedmap_tn = (sum_unweighted_ec_tn != 0) ? observed_total_ec_tn / sum_unweighted_ec_tn : 1.0;
	}
	else {
		normalization_factor_for_weightedmap_tn = (sum_weighted_ec_tn != 0) ? sum_unweighted_ec_tn / sum_weighted_ec_tn : 1.0;
		normalization_factor_for_unweightedmap_tn = 1;
	}
	// Calculate TP normalization factor
	if (Observed_Load_TP_kg_year > 0) {
		observed_total_ec_tp = Observed_Load_TP_kg_year / (cellsize_m * cellsize_m / 10000.0);//EC1 (kg/ha) * 900/10000 + EC2 (kg/ha) * 900/10000 + ...ECn (kg/ha) * 900/10000  = sum(EC)*900/10000 = total load (kg)
		normalization_factor_for_weightedmap_tp = (sum_weighted_ec_tp != 0) ? observed_total_ec_tp / sum_weighted_ec_tp : 1.0;
		normalization_factor_for_unweightedmap_tp = (sum_unweighted_ec_tp != 0) ? observed_total_ec_tp / sum_unweighted_ec_tp : 1.0;
	}
	else {
		normalization_factor_for_weightedmap_tp = (sum_weighted_ec_tp != 0) ? sum_unweighted_ec_tp / sum_weighted_ec_tp : 1.0;
		normalization_factor_for_unweightedmap_tp = 1;
	}
	// Calculate TSS normalization factor
	if (Observed_Load_TSS_kg_year > 0) {
		observed_total_ec_tss = Observed_Load_TSS_kg_year / (cellsize_m * cellsize_m / 10000.0);//EC1 (kg/ha) * 900/10000 + EC2 (kg/ha) * 900/10000 + ...ECn (kg/ha) * 900/10000  = sum(EC)*900/10000 = total load (kg)
		normalization_factor_for_weightedmap_tss = (sum_weighted_ec_tss != 0) ? observed_total_ec_tss / sum_weighted_ec_tss : 1.0;
		normalization_factor_for_unweightedmap_tss = (sum_unweighted_ec_tss != 0) ? observed_total_ec_tss / sum_unweighted_ec_tss : 1.0;
	}
	else {
		normalization_factor_for_weightedmap_tss = (sum_weighted_ec_tss != 0) ? sum_unweighted_ec_tss / sum_weighted_ec_tss : 1.0;
		normalization_factor_for_unweightedmap_tss = 1;
	}

	cout << "Unweighted sum of Export Coefficient TP is " << sum_unweighted_ec_tp << " kg" << endl;
	cout << "Normalization factor for unweighted Export Coefficient TP is " << normalization_factor_for_unweightedmap_tp << endl;
	cout << "Weighted sum of Export Coefficient TP is " << sum_weighted_ec_tp << " kg" << endl;
	cout << "Normalization factor for weighted Export Coefficient TP is " << normalization_factor_for_weightedmap_tp << endl;
	if (Observed_Load_TP_kg_year > 0) {
		cout << "Observed sum of Export Coefficient TP is " << observed_total_ec_tp << " kg" << endl;
	}

	// Step 4: Set surface and subsurface fractions
	double surface_frac_tn = 1.0, surface_frac_tp = 1.0, surface_frac_tss = 1.0;
	if (isSubsurfaceCalculationEnabled) {
		surface_frac_tn = input->BufferSpatialWeightingNumericalParams["ECTN_surface_contribution_frac"];
		surface_frac_tp = input->BufferSpatialWeightingNumericalParams["ECTP_surface_contribution_frac"];
		surface_frac_tss = input->BufferSpatialWeightingNumericalParams["ECTSS_surface_contribution_frac"];
	}
	// Step 5: Precompute conversion factors for EMC calculation
	double conversion_factor = Ratio_kg_to_mg / Ratio_mmHa_to_L;
	
	// Resize intermediate vectors and initialize with NODATA_code
	surf_ec_tn_wtd.resize(nRows * nCols, NODATA_code);
	surf_ec_tp_wtd.resize(nRows * nCols, NODATA_code);
	surf_ec_tss_wtd.resize(nRows * nCols, NODATA_code);
	sub_ec_tn_wtd.resize(nRows * nCols, NODATA_code);
	sub_ec_tp_wtd.resize(nRows * nCols, NODATA_code);
	sub_ec_tss_wtd.resize(nRows * nCols, NODATA_code);
	ec_tn_vec.resize(nRows * nCols, NODATA_code);
	ec_tp_vec.resize(nRows * nCols, NODATA_code);
	ec_tss_vec.resize(nRows * nCols, NODATA_code);

	// Step 6: Apply normalization and populate total weighted EC and EMC maps
	for (int i = 0; i < nRows * nCols; i++) {
		// Check if the current cell is valid
		if (TerrainProcessor::elevation_m_1D[i] != NODATA_code) {
			if (waters01[i] != 0) {
				// Calculate surface and subsurface ratio products
				long double surf_ratio_product_temp = BI_surf_ratio[i] * RI_surf_ratio[i];
				long double sub_ratio_product_temp = (isSubsurfaceCalculationEnabled && BI_sub_ratio[i] != NODATA_code && RI_sub_ratio[i] != NODATA_code)
					? BI_sub_ratio[i] * RI_sub_ratio[i] : 0;

				// Store normalized weighted surface EC values
				surf_ec_tn_wtd[i] = surf_ratio_product_temp * surface_frac_tn * ec_tn_orig[i] * normalization_factor_for_weightedmap_tn;
				surf_ec_tp_wtd[i] = surf_ratio_product_temp * surface_frac_tp * ec_tp_orig[i] * normalization_factor_for_weightedmap_tp;
				surf_ec_tss_wtd[i] = surf_ratio_product_temp * surface_frac_tss * ec_tss_orig[i] * normalization_factor_for_weightedmap_tss;

				// Store normalized subsurface EC values
				sub_ec_tn_wtd[i] = sub_ratio_product_temp * (1 - surface_frac_tn) * ec_tn_orig[i] * normalization_factor_for_weightedmap_tn;
				sub_ec_tp_wtd[i] = sub_ratio_product_temp * (1 - surface_frac_tp) * ec_tp_orig[i] * normalization_factor_for_weightedmap_tp;
				sub_ec_tss_wtd[i] = sub_ratio_product_temp * (1 - surface_frac_tss) * ec_tss_orig[i] * normalization_factor_for_weightedmap_tss;

				// Calculate combined total EC values for TN, TP, and TSS
				total_ec_tn_wtd.push_back(surf_ec_tn_wtd[i] + sub_ec_tn_wtd[i]);
				total_ec_tp_wtd.push_back(surf_ec_tp_wtd[i] + sub_ec_tp_wtd[i]);
				total_ec_tss_wtd.push_back(surf_ec_tss_wtd[i] + sub_ec_tss_wtd[i]);

				//normalize unweighted EC map
				ec_tn_vec[i] = ec_tn_orig[i] * normalization_factor_for_unweightedmap_tn;
				ec_tp_vec[i] = ec_tp_orig[i] * normalization_factor_for_unweightedmap_tp;
				ec_tss_vec[i] = ec_tss_orig[i] * normalization_factor_for_unweightedmap_tss;

				// Calculate EMC values based on water volume if applicable
				double nlcd = input->LandCover_NLCD_Class[i];
				double watervolume_mmha = WaterVolume_Retrieve_50thPercentile(input, nlcd);

				if (nlcd == 11 || nlcd == 95 || nlcd == 90) {
					total_emc_tn_wtd.push_back(NODATA_code);
					total_emc_tp_wtd.push_back(NODATA_code);
					total_emc_tss_wtd.push_back(NODATA_code);
				}
				else {
					// EMC = weighted EC * conversion factor / water volume
					total_emc_tn_wtd.push_back((total_ec_tn_wtd.back() * conversion_factor) / watervolume_mmha);
					total_emc_tp_wtd.push_back((total_ec_tp_wtd.back() * conversion_factor) / watervolume_mmha);
					total_emc_tss_wtd.push_back((total_ec_tss_wtd.back() * conversion_factor) / watervolume_mmha);
				}
			}
			else {
				// Set NODATA values if conditions aren't met
				total_ec_tn_wtd.push_back(0);
				total_ec_tp_wtd.push_back(0);
				total_ec_tss_wtd.push_back(0);
				total_emc_tn_wtd.push_back(0);
				total_emc_tp_wtd.push_back(0);
				total_emc_tss_wtd.push_back(0);
			}
		}
		else {
			// Set NODATA values if conditions aren't met
			total_ec_tn_wtd.push_back(NODATA_code);
			total_ec_tp_wtd.push_back(NODATA_code);
			total_ec_tss_wtd.push_back(NODATA_code);
			total_emc_tn_wtd.push_back(NODATA_code);
			total_emc_tp_wtd.push_back(NODATA_code);
			total_emc_tss_wtd.push_back(NODATA_code);
		}
	}
}

//ReplaceElevationWithVector function to replace elevation values based on a specified vector (e.g., elevation_m_negative_1D)
void BufferSpatialCalc::ReplaceElevationWithVector(vector<double>& source_vector) {
	for (int MapPixel_ID = 0; MapPixel_ID < nRows * nCols; MapPixel_ID++) {
		//Note: Conversion between MapPixel_ID and row & column pair, with MapPixel_ID starting at 0, row & col starting at 0
		//Note: Eqs: row=MapPixel_ID / Inputs::nCols; col=MapPixel_ID % Inputs::nCols; MapPixel_ID = (row * nCols) + col
		int row = MapPixel_ID / nCols;
		//col = MapPixel_ID % nCols; indexing starts at 0, not the ESRI approach of starting row a 1
		int col = MapPixel_ID % nCols;

		//If elevation_m_2D[row][col] != NODATA_code then replace with source_vector[MapPixel_ID]
		if (TerrainProcessor::elevation_m_2D[row][col] != NODATA_code) {
			TerrainProcessor::elevation_m_2D[row][col] = source_vector[MapPixel_ID];
		}
		//Else
		else {
			//elevation_m_2D[row][col] = NODATA_code; ignored
			TerrainProcessor::elevation_m_2D[row][col] = NODATA_code;
		}
	}
}

void BufferSpatialCalc::ComputeWEM(Inputs* input){
	// Resize the WEM vector to match the 1D grid size
	groundwatertable_elevation1D.clear();
	groundwatertable_elevation1D.resize(nRows * nCols, NODATA_code);
	neg_groundwatertable_elevation1D.clear();
	neg_groundwatertable_elevation1D.resize(nRows * nCols, NODATA_code);

	//if user provide f and watertable avg, then estimateWEM using Endreny 1999 Eq 2 method
	//else read in watertable map s_wtable.asc
	if (input->BufferSpatialWeightingNumericalParams["f_param"] > 0 && input->BufferSpatialWeightingNumericalParams["S_WTable_avg"] > 0) {
		Depth_GroundwaterTable_calc(input);
		for (int index = 0; index < nRows * nCols; ++index) {
			// Check if both elevation and waterTableDepth have valid values
			if (TerrainProcessor::elevation_m_1D[index] != NODATA_code) {
				// Calculate WEM as elevation - water table depth (converted from cm to m)
				groundwatertable_elevation1D[index] = TerrainProcessor::elevation_m_1D[index] - (s_wtable_est[index] / 100.0);
				neg_groundwatertable_elevation1D[index] = -groundwatertable_elevation1D[index];
			}
		}
	}else {
		// Otherwise, use the provided water table map (s_wtable.asc)
		for (int index = 0; index < nRows * nCols; ++index) {
			// Check if both elevation and waterTableDepth have valid values
			if (TerrainProcessor::elevation_m_1D[index] != NODATA_code) {
				// Calculate WEM as elevation - water table depth (converted from cm to m)
				groundwatertable_elevation1D[index] = TerrainProcessor::elevation_m_1D[index] - (input->s_wtable[index] / 100.0);
				neg_groundwatertable_elevation1D[index] = -groundwatertable_elevation1D[index];
			}
		}
	}
}

void BufferSpatialCalc::recomputeFlowDirection( Inputs* input) {
	//note that slope calculated during the flow direction computation
	//redo the flow direction to get the new slope, which will be the watertable slope
	//re-calculate the flow direction, slopeorganizer_riserun1D has been updated accordingly
	if (input->SimulationStringParams["Algorithm_FlowDirection"] == "DInfinity")
	{
		TerrainProcessor::calcFlowDirection_DInf();
	}
	//calcFlowDirection_D8 computes the direction DEM cells will flow, using D8 algorithm (O'Callaghan and Mark, 1984)
	//Note: Slope is explicitly computed within the flow direction algorithm to find steepest single path
	else if (input->SimulationStringParams["Algorithm_FlowDirection"] == "D8")
	{
		TerrainProcessor::calcFlowDirection_D8();
	}
	//calcFlowDirection_MFD computes the direction DEM cells will flow, using multiple flow direction algorithm (Quinn et al. 1991)
	//Note: Slope is implicitly computed within the flow direction algorithm to find all downslope paths
	else if (input->SimulationStringParams["Algorithm_FlowDirection"] == "MFD")
	{
		TerrainProcessor::calcFlowDirection_MFD();
	}
}

//BufferSpatialCalc::BufferModel_IntermediateOutputs_writer function will write intermediate products from Buffer model calculations
void BufferSpatialCalc::BufferModel_IntermediateOutputs_writer(string file_name) {
	int case_number;
	string fName;
	fName = file_name + ".asc";

	if (file_name == "dem" || file_name == "neg_dem"||file_name == "slope_dem"|| file_name == "dem_filled") case_number = 1;
	if (file_name == "demfdr" || file_name == "demfdrr") case_number = 2;
	if (file_name == "nhdfdr" || file_name == "nhdfdrr") case_number = 3;
	if (file_name == "waters01" || file_name == "waters0") case_number = 4;
	if (file_name == "p_release" || file_name == "p_trapping") case_number = 5;
	if (file_name == "FA_weighted_m") case_number = 6;
	if (file_name == "RI_surf" || file_name == "RI_sub") case_number = 7;
	if (file_name == "surf_buffer_slope" || file_name == "surf_buffer_rise2" || file_name == "surf_buffer_rise1" || file_name == "surf_buffer_rise"|| file_name == "surf_buffer_run1"|| file_name == "surf_buffer_run"
		||file_name == "sub_buffer_slope" || file_name == "sub_buffer_rise2" || file_name == "sub_buffer_rise1" || file_name == "sub_buffer_rise" || file_name == "sub_buffer_run1" || file_name == "sub_buffer_run") case_number = 8;
	if (file_name == "Manningn"|| file_name == "Hydraulic_R") case_number = 9;
	if (file_name == "Time_surf" || file_name == "Time_surf1" || file_name == "Time_surf2") case_number = 10;
	if (file_name == "BI_surf") case_number = 11;
	if (file_name == "ContributingAreaPixels") {
		case_number = 12;
		fName = file_name + ".csv";
	}
	if (file_name == "s_wtable_est") case_number = 13;
	if (file_name == "slope_wem || slope_negwem") case_number = 14;
	if (file_name == "Time_sub" || file_name == "Time_sub1" || file_name == "Time_sub2") case_number = 15;
	if (file_name == "BI_sub") case_number = 16;
	if (file_name == "s_trans"|| file_name == "s_depth") case_number = 17;
	if (file_name == "surf_filter" || file_name == "surf_filter1" || file_name == "surf_filter2" || file_name == "surf_filter_to_slope"
		|| file_name == "sub_filter" || file_name == "sub_filter1" || file_name == "sub_filter2" || file_name == "sub_filter_to_slope") case_number = 18;
	if (file_name == "RI_surf_ratio" || file_name == "BI_surf_ratio" || file_name == "RI_sub_ratio" || file_name == "BI_sub_ratio") case_number = 19;
	if (file_name == "ec_tn_vec" || file_name == "ec_tp_vec" || file_name == "ec_tss_vec") case_number = 20;
	if (file_name == "surf_ec_tn_wtd" || file_name == "surf_ec_tp_wtd" || file_name == "surf_ec_tss_wtd" || file_name == "sub_ec_tn_wtd" || file_name == "sub_ec_tp_wtd" || file_name == "sub_ec_tss_wtd") case_number = 21;
	if (file_name.find("total_ec_tn_wtd") != string::npos || file_name.find("total_ec_tp_wtd") != string::npos || file_name.find ("total_ec_tss_wtd") != string::npos 
		|| file_name.find("total_emc_tn_wtd") != string::npos || file_name.find("total_emc_tp_wtd") != string::npos || file_name.find("total_emc_tss_wtd") != string::npos) case_number = 22;
	if (file_name == "hotspot_tn" || file_name == "hotspot_tp" || file_name == "hotspot_tss") case_number = 23;
	if (file_name.find("ContributingAreaPixels_map") != string::npos || file_name.find("DispersalAreaPixels_map") != string::npos) case_number = 24;

	//erase the "_surf" from the file name 
	//for avoiding confusion since Buffer 2003 doesn't have subsurface layer calculation
	if (!isSubsurfaceCalculationEnabled) {
		// Remove "_surf" from anywhere in the file name to avoid confusion
		size_t found = fName.find("_surf");
		while (found != string::npos) {
			fName.erase(found, 5); // Remove "_surf"
			found = fName.find("_surf", found); // Continue searching from the current position
		}

		// Remove "surf_" if it appears at the start of the file name
		if (fName.find("surf_") == 0) {
			fName.erase(0, 5); // Remove "surf_" from the beginning
		}
	}

	//start creating the output files
	ofstream outfile(Inputs::SimulationStringParams["OutputFolder_Path"] + fName);
	
	outfile << left << setw(14) << "ncols" << setw(15) << Inputs::nCols << endl;
	outfile << left << setw(14) << "nrows" << setw(15) << Inputs::nRows << endl;
	outfile << left << setw(14) << "xllcorner" << setw(15) << setprecision(15) << Inputs::xllcorner_m << endl;
	outfile << left << setw(14) << "yllcorner" << setw(15) << setprecision(15) << Inputs::yllcorner_m << endl;
	outfile << left << setw(14) << "cellsize" << setw(15) << setprecision(15) << Inputs::Length_Pixel_Side_m << endl;
	outfile << left << setw(14) << "NODATA_value" << setw(15) << Inputs::NODATA_code << endl;
	
	if (!outfile.good()) {
		cout << "Cannot write file: " + fName << endl;
		return;
	}

	//DataDrawers_Size defined as product of nRows * nCols
	int DataDrawers_Size =Inputs::nRows * Inputs::nCols;
	//For Loop through all map rows and columnes, defiend as DataDrawers_Size, where nRows * nCols = number of map pixels
	//Note: For map based models, DataDrawer_ID has maximum = maxRows * maxCols, and DataFolder_ID typically has maximum = 1
	for (int DataDrawer_ID = 0; DataDrawer_ID < DataDrawers_Size; ++DataDrawer_ID) {

		//Note: Below the code is converting between map vector DataDrawer_ID and pixel row, col pair
		//Note: Conversion uses 0-based indexing: row = MapPixel_ID / nCols; col = MapPixel_ID % nCols;
		//MapPixel_ID = (row * nCols) + col; where MapPixel_ID ranges from 0 to (nRows * nCols - 1)
		int col = DataDrawer_ID % Inputs::nCols;
		//newLine = ((col + 1) == Inputs::nCols); True when col+1 equals nCols
		bool newLine = ((col + 1) == Inputs::nCols);
		
		//If LandCover_NLCD_Class[DataDrawer_ID] equals NODATA_code value from HydroPlusConfig.xml then
		if (int(TerrainProcessor::elevation_m_1D[DataDrawer_ID]) == int(Inputs::NODATA_code)) {
			//write NODATA_code to map
			outfile << int(Inputs::NODATA_code);
		}
		else {
			//If needed, to set precision on output, use a line of code such as: outfile.precision(2)
			//write outputVarName variable value to map
			switch (case_number) {
			case 1:
				if (file_name == "dem"|| file_name == "dem_filled") outfile << setprecision(15) << TerrainProcessor::elevation_m_1D[DataDrawer_ID];
				if (file_name == "neg_dem") outfile << setprecision(15) << TerrainProcessor::elevation_m_negative_1D[DataDrawer_ID];
				if (file_name == "slope_dem") outfile << slope_dem[DataDrawer_ID];
				break;
			case 2:
				if (file_name == "demfdr") {
					if (Inputs::SimulationStringParams["Algorithm_FlowDirection"] == "D8") {
						outfile << Inputs::flowdircode_whitebox_to_esri(Inputs::DEMfdr[DataDrawer_ID]);
					}
					else {
						outfile << Inputs::DEMfdr[DataDrawer_ID];
					}
				}
				if (file_name == "demfdrr") {
					if (Inputs::SimulationStringParams["Algorithm_FlowDirection"] == "D8") {
						outfile << Inputs::flowdircode_whitebox_to_esri(Inputs::negDEMfdr[DataDrawer_ID]);
					}
					else {
						outfile << Inputs::negDEMfdr[DataDrawer_ID];
					}
				}
				break;
			case 3:
				if (file_name == "nhdfdr") outfile << Inputs::flowdircode_whitebox_to_esri(Inputs::NHDfdr[DataDrawer_ID]);
				if (file_name == "nhdfdrr") outfile << Inputs::flowdircode_whitebox_to_esri(Inputs::NHDfdr_reverse[DataDrawer_ID]);
				break;
			case 4:
				if (file_name == "waters01") outfile << waters01[DataDrawer_ID];
				if (file_name == "waters0") outfile << waters0[DataDrawer_ID];
				break;
			case 5:
				if (file_name == "p_release") outfile << p_release_vec[DataDrawer_ID];
				if (file_name == "p_trapping") outfile << p_trapping_vec[DataDrawer_ID];
				break;
			case 6:
				if (file_name == "FA_weighted_m") outfile << FA_weighted_m[DataDrawer_ID];
				break;
			case 7:
				if (file_name == "RI_surf") outfile << RI_surf[DataDrawer_ID];
				if (file_name == "RI_sub") outfile << RI_sub[DataDrawer_ID];
				break;
			case 8:
				if (file_name == "surf_buffer_rise1" || file_name == "sub_buffer_rise1") outfile << buffer_rise1[DataDrawer_ID];
				if (file_name == "surf_buffer_rise2"|| file_name == "sub_buffer_rise2") outfile << buffer_rise2[DataDrawer_ID];
				if (file_name == "surf_buffer_rise"|| file_name == "sub_buffer_rise") outfile << buffer_rise[DataDrawer_ID];
				if (file_name == "surf_buffer_run1"|| file_name == "sub_buffer_run1") outfile << buffer_run1[DataDrawer_ID];
				if (file_name == "surf_buffer_run"|| file_name == "sub_buffer_run") outfile << buffer_run[DataDrawer_ID];
				if (file_name == "surf_buffer_slope"|| file_name == "sub_buffer_slope") outfile << buffer_slope[DataDrawer_ID];
				break;
			case 9:
				if (file_name == "Manningn") outfile << Manningn_vec[DataDrawer_ID];
				if (file_name == "Hydraulic_R") outfile << Hydraulic_R_vec[DataDrawer_ID];
				break;
			case 10:
				if (file_name == "Time_surf1") outfile << Time_surf1[DataDrawer_ID];
				if (file_name == "Time_surf2") outfile << Time_surf2[DataDrawer_ID];
				if (file_name == "Time_surf") outfile << Time_surf[DataDrawer_ID];
				break;
			case 11:
				if (file_name == "BI_surf") outfile << BI_surf[DataDrawer_ID];
				break;
			case 12:
				if (file_name == "ContributingAreaPixels") {
					for (pair<int, int> p : TerrainProcessor::ContributingAreaPixels_FlowIntoAllPixels_vec[DataDrawer_ID]) {
						//DataDrawer_ID = (row * col_max) + col
						outfile << p.first * nCols + p.second << " ";
					}
				}
				break;
			case 13:
				if (file_name == "S_WTable") outfile << s_wtable_est[DataDrawer_ID];
				break;
			case 14:
				if (file_name == "slope_wem") outfile << slope_wem[DataDrawer_ID];
				if (file_name == "slope_negwem") outfile << slope_negwem[DataDrawer_ID];
				break;
			case 15:
				if (file_name == "Time_sub") outfile << Time_sub[DataDrawer_ID];
				if (file_name == "Time_sub1") outfile << Time_sub1[DataDrawer_ID];
				if (file_name == "Time_sub2") outfile << Time_sub2[DataDrawer_ID];
				break;
			case 16:
				if (file_name == "BI_sub") outfile << BI_sub[DataDrawer_ID];
				break;
			case 17:
				if (file_name == "s_trans") outfile << Inputs::s_trans[DataDrawer_ID];
				if (file_name == "s_depth") outfile << Inputs::s_depth[DataDrawer_ID];
				break;
			case 18:
				if (file_name == "surf_filter"|| file_name == "sub_filter") outfile << filter[DataDrawer_ID];
				if (file_name == "surf_filter1"|| file_name == "sub_filter1") outfile << filter1[DataDrawer_ID];
				if (file_name == "surf_filter2"|| file_name == "sub_filter2") outfile << filter2[DataDrawer_ID];
				if (file_name == "surf_filter_to_slope"|| file_name == "sub_filter_to_slope") outfile << filter_to_slope[DataDrawer_ID];
				break;
			case 19:
				if (file_name == "RI_surf_ratio") outfile << RI_surf_ratio[DataDrawer_ID];
				if (file_name == "BI_surf_ratio") outfile << BI_surf_ratio[DataDrawer_ID];
				if (file_name == "RI_sub_ratio") outfile << RI_sub_ratio[DataDrawer_ID];
				if (file_name == "BI_sub_ratio") outfile << BI_sub_ratio[DataDrawer_ID];
				break;
			case 20:
				if (file_name == "ec_tn_vec") outfile << ec_tn_vec[DataDrawer_ID];
				if (file_name == "ec_tp_vec") outfile << ec_tp_vec[DataDrawer_ID];
				if (file_name == "ec_tss_vec") outfile << ec_tss_vec[DataDrawer_ID];
				break;
			case 21:
				if (file_name == "surf_ec_tn_wtd") outfile << surf_ec_tn_wtd[DataDrawer_ID];
				if (file_name == "surf_ec_tp_wtd") outfile << surf_ec_tp_wtd[DataDrawer_ID];
				if (file_name == "surf_ec_tss_wtd") outfile << surf_ec_tss_wtd[DataDrawer_ID];
				if (file_name == "sub_ec_tn_wtd") outfile << sub_ec_tn_wtd[DataDrawer_ID];
				if (file_name == "sub_ec_tp_wtd") outfile << sub_ec_tp_wtd[DataDrawer_ID];
				if (file_name == "sub_ec_tss_wtd") outfile << sub_ec_tss_wtd[DataDrawer_ID];
				break;
			case 22:
				if (file_name.find("total_ec_tn_wtd") != string::npos) {
					outfile << total_ec_tn_wtd[DataDrawer_ID];
				}
				if (file_name.find("total_ec_tp_wtd") != string::npos) {
					outfile << total_ec_tp_wtd[DataDrawer_ID];
				}
				if (file_name.find("total_ec_tss_wtd") != string::npos) {
					outfile << total_ec_tss_wtd[DataDrawer_ID];
				}
				if (file_name.find("total_emc_tn_wtd") != string::npos) {
					outfile << total_emc_tn_wtd[DataDrawer_ID];
				}
				if (file_name.find("total_emc_tp_wtd") != string::npos) {
					outfile << total_emc_tp_wtd[DataDrawer_ID];
				}
				if (file_name.find("total_emc_tss_wtd") != string::npos) {
					outfile << total_emc_tss_wtd[DataDrawer_ID];
				}
				break;
			case 23:
				if (file_name == "hotspot_tn") outfile << hotspot_map_temp[DataDrawer_ID];
				if (file_name == "hotspot_tp") outfile << hotspot_map_temp[DataDrawer_ID];
				if (file_name == "hotspot_tss") outfile << hotspot_map_temp[DataDrawer_ID];
				break;
			case 24:
				if (file_name.find("ContributingAreaPixels_map") != string::npos)  outfile << ContributingAreaPixels_map_writeout[DataDrawer_ID];
				if (file_name.find("DispersalAreaPixels_map") != string::npos) outfile << DispersalAreaPixels_map_writeout[DataDrawer_ID];
				break;
			}
		}
		//If newLine boolean is true then
		if (newLine) {
			//outfile receives endl function for end of line return
			outfile << endl;
		}
		//Else If then newLine boolean is false
		else {
			if (file_name == "ContributingAreaPixels") outfile << ", ";//this is a csv file
			//outfile receives space between variables
			else outfile << " ";
		}
	}
	//outfile uses close function
	outfile.close(); 
}



//BufferSpatialCalc::Find_MinimumPositiveValue_1D_Vector_calc function finds minimum value in a 1D vector while excluding specific values like -9999 (nodata), -nan, inf, and -inf
double BufferSpatialCalc::Find_MinimumPositiveValue_1D_Vector_calc(vector<double> vec1D) {
	double min = DBL_MAX;
	for (double i : vec1D) {
		double min_temp = i;
		if (min_temp != NODATA_code && !isnan(min_temp) && isfinite(min_temp) && min_temp > 0) {
			if (min_temp < min) {
				min = min_temp;
			}
		}
	}
	return min;
}

// BufferSpatialCalc::Find_MaximumPositiveValue_1D_Vector_calc function finds the maximum value in a 1D vector 
// while excluding specific values like -9999 (NoData), NaN, inf, and -inf.
double BufferSpatialCalc::Find_MaximumPositiveValue_1D_Vector_calc(vector<double> vec1D) {
	double max = -DBL_MAX; // Initialize to the lowest possible double value
	for (double i : vec1D) {
		double max_temp = i;
		if (max_temp != NODATA_code && !isnan(max_temp) && isfinite(max_temp) && max_temp > 0) {
			if (max_temp > max) {
				max = max_temp;
			}
		}
	}
	return max;
}

//sorting and keeping track of indexes
template <typename T>
vector<size_t> sort_indexes(const vector<T>& Vector) {
	// initialize original index locations
	vector<size_t> idx(Vector.size());
	iota(idx.begin(), idx.end(), 0);
	// sort indexes based on comparing values in DistanceVector
	// using stable_sort instead of sort to avoid unnecessary index re-orderings when v contains elements of equal values 
	//stable_sort is used to sort the indices in idx based on the values in Vector. 
	//It compares the values in Vector at indices i1 and i2, sorting them in increasing order (Vector[i1] > Vector[i2]).
	stable_sort(idx.begin(), idx.end(), [&Vector](size_t i1, size_t i2) {return Vector[i1] < Vector[i2]; });
	return idx;
}

//sort the vector to find the >99th percentile pixels
void BufferSpatialCalc::Find_PercentileValue_1D_Vector_calc(Inputs* input, vector<double> weightedECvec1D, vector<long double> unweightedECvec1D, bool isSubsurfaceCalculationEnabled) {
	// Define the hotspot percentile
	int hotspot_percentile = 99;
	
	// Create a list of valid (non-NODATA) indices for ECvec1D
	vector<size_t> validIndices;
	validIndices.reserve(weightedECvec1D.size());
	for (size_t i = 0; i < weightedECvec1D.size(); ++i) {
		if (weightedECvec1D[i] != NODATA_code) {
			validIndices.push_back(i);
		}
	}
	// Calculate the number of valid data points and find the index for the 99th percentile
	int data_elements = validIndices.size();
	int percentile_index = static_cast<int>(data_elements * hotspot_percentile / 100);

	// Use nth_element to get the 99th percentile threshold value
	nth_element(validIndices.begin(), validIndices.begin() + percentile_index, validIndices.end(),
		[&weightedECvec1D](size_t i1, size_t i2) { 
			return weightedECvec1D[i1] < weightedECvec1D[i2]; 
		});
	double percentile_value = weightedECvec1D[validIndices[percentile_index]];


	// Resize hotspot_value_temp to have the same size as vec1D, filled with NODATA_code
	hotspot_map_temp.resize(weightedECvec1D.size(), NODATA_code);
	//clean up the hotspot_report before pushing back the tuples
	hotspot_report_temp.clear();
	//get the column numbers from the input
	nCols = Inputs::nCols;
	
		// Process only pixels at or above the 99th percentile
	for (size_t idx : validIndices) {
		double pixelECValue = weightedECvec1D[idx];

		if (pixelECValue >= percentile_value) {
			int pixelID = static_cast<int>(idx);
			//Note: Conversion between MapPixel_ID and row & column pair, with MapPixel_ID starting at 0, row & col starting at 0
			//Note: Eqs: row=MapPixel_ID / Inputs::nCols; col=MapPixel_ID % Inputs::nCols; MapPixel_ID = (row * nCols) + col
			int row = pixelID / nCols;
			//col = MapPixel_ID % nCols; indexing starts at 0, not the ESRI approach of starting row a 1
			int col = pixelID % nCols;
			//landcover = input->LandCover_NLCD_Class[pixelID]; initialized to NLCD_Class
			double landcover = input->LandCover_NLCD_Class[pixelID];
			double unweighted_ec = unweightedECvec1D[pixelID];
			//double pixelEMCValue;


			// Determine contributing pixel count
			int contributing_pixel_count;
			contributing_pixel_count = FA_weighted_m[pixelID]/cellsize_m;

			// Surface and subsurface runoff/buffer index values
			double runoffindex_surface = RI_surf[pixelID];
			double runoffindex_surface_median = RI_surf_median;
			double bufferindex_surface = BI_surf[pixelID];
			double bufferindex_surface_median = BI_surf_median;
			double runoffindex_subsurface = isSubsurfaceCalculationEnabled ? RI_sub[pixelID] : NODATA_code;
			double runoffindex_subsurface_median = isSubsurfaceCalculationEnabled ? RI_sub_median : NODATA_code;
			double bufferindex_subsurface = isSubsurfaceCalculationEnabled ? BI_sub[pixelID] : NODATA_code;
			double bufferindex_subsurface_median = isSubsurfaceCalculationEnabled ? BI_sub_median : NODATA_code;

			// Update the hotspot map and report
			hotspot_map_temp[pixelID] = pixelECValue;
			hotspot_report_temp.push_back({ pixelID, row, col, landcover, pixelECValue, unweighted_ec,
										   runoffindex_surface, runoffindex_surface_median, bufferindex_surface,bufferindex_surface_median,
										   runoffindex_subsurface, runoffindex_subsurface_median, bufferindex_subsurface,bufferindex_subsurface_median,
										   contributing_pixel_count });
		}
	}
}

void BufferSpatialCalc::writeHotspotReport(string filename, vector<tuple<int, int, int, double, double, double, double, double, double, double, double, double, double, double, int>>hotspot_tuple) {
	//Sort by EC descending; tie-breaker: contributing_pixel_count descending
	sort(hotspot_tuple.begin(), hotspot_tuple.end(), [](const auto& a, const auto& b) {
		// Primary: EC(element 4)
		if (get<4>(a) != get<4>(b))
			return get<4>(a) < get<4>(b);
		//Secondary: contributing_pixel_count (element 14)
		if (get<14>(a) != get<14>(b))
			return get<14>(a) < get<14>(b);
		//Tertiary: pixelID (element 0)
		return get<0>(a) > get<0>(b);
		});

	ofstream outfile(Inputs::SimulationStringParams["OutputFolder_Path"] + filename);

	if (!outfile.good()) {
		cout << "Warning: Output folder does not exist the file is open and cannot be overwritten." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Correction: Confirm the hotspotreport.txt is closed." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}
	else {
		cout << "Writing " << filename <<" to output folder. " << endl;
	}

	// Write header
	outfile << "Pixel_ID,Row,Column,Land_Cover,EC_kg_p_ha,unweighted_EC,runoffindex_surface,runoffindex_surface_median,bufferindex_surface,bufferindex_surface_median,runoffindex_subsurface,runoffindex_subsurface_median,bufferindex_subsurface,bufferindex_subsurface_median,contributing_cells_count" << endl;

	// Write hotspot data (limit to first 500 pixels)
	//int pixelCount = 0;
	for (auto it = hotspot_tuple.rbegin(); it != hotspot_tuple.rend(); ++it) {

		int pixelID = get<0>(*it);
		
		int row = get<1>(*it) + 1;//ArcGIS row col start from 1
		int col = get<2>(*it) + 1;//ArcGIS row col start from 1

		double lc = get<3>(*it);
		double ec = get<4>(*it);
		//double emc = get<5>(*it);
		double ec_unweighted = get<5>(*it);

		double ri_sur = get<6>(*it);
		double ri_sur_med = get<7>(*it);

		double bi_sur = get<8>(*it);
		double bi_sur_med = get<9>(*it);

		double ri_sub = get<10>(*it);
		double ri_sub_med = get<11>(*it);

		double bi_sub = get<12>(*it);
		double bi_sub_med = get<13>(*it);

		int count = get<14>(*it);

		outfile << pixelID << "," << row << "," << col << "," << lc << "," << ec << "," << ec_unweighted << ","
			<< ri_sur << "," << ri_sur_med << "," << bi_sur << "," << bi_sur_med << "," << ri_sub << "," << ri_sub_med << "," << bi_sub << "," << bi_sub_med << "," << count << endl;
		//++pixelCount;  // Increment counter after each pixel is processed
	}

	outfile.close();
}

//BufferSpatialCalc::ReplaceFlowDirectionOrganizerWithVector function replaces DEM flow direction values with NHD flow direction values
//Note: NHD are National Hydrography Dataset Plus that NHD corrected to ensure flow paths align with rivers as represented in planimetric maps
//Note: https://www.epa.gov/waterdata/nhdplus-national-hydrography-dataset-plus
void BufferSpatialCalc::ReplaceFlowDirectionOrganizerWithVector(vector<double>& source_vector) {
	nRows = Inputs::nRows;
	nCols = Inputs::nCols;
	NODATA_code = Inputs::NODATA_code;
	//For Loop row through nRows
	for (int row = 0; row < nRows; row++) {
		//For Loop col through nCols
		for (int col = 0; col < nCols; col++) {
			//MapPixel_ID defined as from row, col and nCols inputs
			//Note: Conversion between MapPixel_ID and row & column pair, with MapPixel_ID starting at 0, row & col starting at 0
			//Note: Eqs: row=MapPixel_ID / Inputs::nCols; col=MapPixel_ID % Inputs::nCols; MapPixel_ID = (row * nCols) + col
			int MapPixel_ID = (row * nCols) + col;

			//If TerrainProcessor::elevation_m_2D[row][col] is not a NODATA value then can replace DEM FDR with NHD FDR
			if (TerrainProcessor::elevation_m_2D[row][col] != NODATA_code) {
				//If row=0 or col=0 or row=nRows-1 or col=nCols-1 then along the boundary or edge of the map input
				//Note: Caution when along the edge of map input, as NHD flow direction may exit the map area and cause failure in subsequent function
				if (row == 0 || row == nRows - 1 || col == 0 || col == nCols - 1) {
					//If TerrainProcessor::FlowDirection_Map vector using DEM flow direction value is not equal to NHDfdr for the same MapPixel_ID then
					//Note: Comparing 2D FlowDirection_Map with 1D source_vector; they have coordinated indexing
					if (TerrainProcessor::FlowDirection_Map[row][col] != source_vector[MapPixel_ID]) {
						//Command Prompt Notification
						cout << "Note: DEM replaced NHD flow direction at map boundary pixel [row " << row << ", col " << col << "]." << endl;
					}
				}
				//If Else not along the boundary of the map input, and no caution needed about NHD flow direction 
				else {
					//TerrainProcessor::FlowDirection_Map[row][col] is set equal to NHDfdr[MapPixel_ID]; comparing 2D with 1D vectors
					TerrainProcessor::FlowDirection_Map[row][col] = source_vector[MapPixel_ID];
				}
			}
			else {
				//TerrainProcessor::FlowDirection_Map[row][col] is set equal to NODATA
				TerrainProcessor::FlowDirection_Map[row][col] = NODATA_code;
			}
		}
	}
}

void BufferSpatialCalc::Contributing_dispersal_CreateMap(int row, int col) {
	nCols = Inputs::nCols;
	int pixel_id = (row * nCols) + col;

	vector<pair<int, int>> contributing_vector_temp, dispersal_vector_temp;
	//contributing_vector_temp = BufferSpatialCalc::ContributingAreaPixels_FlowIntoAllPixels_vec[pixel_id];
	contributing_vector_temp.insert(
		contributing_vector_temp.end(),
		TerrainProcessor::ContributingAreaPixels_FlowIntoAllPixels_vec[pixel_id].begin(),
		TerrainProcessor::ContributingAreaPixels_FlowIntoAllPixels_vec[pixel_id].end()
	);
	dispersal_vector_temp.insert(
		dispersal_vector_temp.end(),
		TerrainProcessor::DispersalAreaPixels_FlowFromAllPixels_vec[pixel_id].begin(),
		TerrainProcessor::DispersalAreaPixels_FlowFromAllPixels_vec[pixel_id].end()
	);
	/*cout << pixel_id << " = (" << row << " * " << nCols << ")+" << col << endl;
	cout << contributing_vector_temp.size() << endl;
	cout << dispersal_vector_temp.size() << endl;*/
	
	for (pair<int, int> p : contributing_vector_temp) {
		//ContributingAreaPixels_map_writeout[p.first * nCols + p.second] += 1;
		ContributingAreaPixels_map_writeout[p.first * nCols + p.second] = 1;
	}
	for (pair<int, int> p : dispersal_vector_temp) {
		//DispersalAreaPixels_map_writeout[p.first * nCols + p.second] += 1;
		DispersalAreaPixels_map_writeout[p.first * nCols + p.second] = 1;
	}
}

void BufferSpatialCalc::updateLandCover(Inputs* input){
	// Check if the vectors are of the same length
	if (input->GI_list.size() != input->GI_ulcorner_rows.size() ||
		input->GI_list.size() != input->GI_ulcorner_cols.size() ||
		input->GI_list.size() != input->GI_lrcorner_rows.size() ||
		input->GI_list.size() != input->GI_lrcorner_cols.size()) {
		cout << "Warning: GI vectors must be of the same length." << endl;
		cout << "Aborting: This warning triggers the simulation to abort." << endl;
		cout << "Explanation: The vectors specifying GI types and their respective coordinates must have the same length to ensure each GI type has corresponding coordinate data." << endl;
		cout << "Correction: Check the input vectors for GI types in DataFolder and coordinates in BufferSpatialWeightingStringParams. Ensure they all have the same number of elements." << endl;
	}

	// Iterate through each GI area and update land cover
	for (size_t gi = 0; gi < input->GI_list.size(); ++gi) {
		string GI_type = input->GI_list[gi];
		int GI_ulcorner_row = input->GI_ulcorner_rows[gi];
		int GI_ulcorner_col = input->GI_ulcorner_cols[gi];
		int GI_lrcorner_row = input->GI_lrcorner_rows[gi];
		int GI_lrcorner_col = input->GI_lrcorner_cols[gi];

		// Update the land cover type in the specified area
		for (int row = GI_ulcorner_row; row <= GI_lrcorner_row; ++row) {
			for (int col = GI_ulcorner_col; col <= GI_lrcorner_col; ++col) {
				int index = row * nRows + col;

				// Determine the new land cover type based on GI type
				if (GI_type == "BioRetention" || GI_type == "RainGarden") {
					input->LandCover_NLCD_Class[index] = 51;
				}
				else if (GI_type == "GreenRoof" || GI_type == "InfiltrationTrench" || GI_type == "Swale") {
					input->LandCover_NLCD_Class[index] = 71;
				}
			}
		}
	}
}

//sum the contributing EC for a GI block
//need to run after getting weighted EC
double BufferSpatialCalc::AverageContributingAreaLoad_calc(int GI_ulcorner_row, int GI_ulcorner_col, int GI_lrcorner_row, int GI_lrcorner_col, string type) {
	double ca_weighted_ec_sum_kg_ha = 0.0;
	int count = 0;
	//loop through the ContributingAreaPixels_map_writeout to avoid adding one cells muiltiple times
	for (int pixelIndex = 0; pixelIndex < ContributingAreaPixels_map_writeout.size(); pixelIndex++) {
		double contributingPixels = ContributingAreaPixels_map_writeout[pixelIndex];
		if (contributingPixels >= 1) {
			if (type == "TP") {
				if (total_ec_tp_wtd[pixelIndex] > 0) {
					ca_weighted_ec_sum_kg_ha += total_ec_tp_wtd[pixelIndex];
				}
			}
			if (type == "TN") {
				if (total_ec_tn_wtd[pixelIndex] > 0) {
					ca_weighted_ec_sum_kg_ha += total_ec_tn_wtd[pixelIndex];
				}
			}
			if (type == "TSS") {
				if (total_ec_tss_wtd[pixelIndex] > 0) {
					ca_weighted_ec_sum_kg_ha += total_ec_tss_wtd[pixelIndex];
				}
			}
			count++;
		}
	}
	return ca_weighted_ec_sum_kg_ha/ count;
}

void BufferSpatialCalc:: ContributingConcentrationValues_calc(Inputs* input) {
	ca_concentration_tp_mg_L.clear();
	ca_concentration_tn_mg_L.clear();
	ca_concentration_tss_mg_L.clear();
	double Q_m3, M_kg_ha_yr, A_m2, t_s, M_kg, M_mg;
	double concentration_tp_mg_L, concentration_tn_mg_L, concentration_tss_mg_L;

	// Iterate through each GI area and update land cover
	for (size_t GI_Index = 0; GI_Index < input->GI_list.size(); ++GI_Index) {
		string GI_type = input->GI_list[GI_Index];
		int GI_ulcorner_row = input->GI_ulcorner_rows[GI_Index];
		int GI_ulcorner_col = input->GI_ulcorner_cols[GI_Index];
		int GI_lrcorner_row = input->GI_lrcorner_rows[GI_Index];
		int GI_lrcorner_col = input->GI_lrcorner_cols[GI_Index];

		//C = M/Q = EC (kg/ha/year) * area * simulation time period
		Q_m3 = input->RepoDict["GI_Runon_from_BulkArea_Surface_Sum_m3"];
		A_m2 = input-> GI_Contributing_Area_m2[GI_Index];
		//t_s is input->SimulationTimePeriod_seconds 
		t_s = input->SimulationTimePeriod_seconds;

		if (A_m2 <= cellsize_m * cellsize_m) {
			cout << "Warning: The selected GI location has <= 1 contributing cell, which may result in no inflow from the bulk area and NA values in the output." << endl;
			cout << "Correction: Please review the hotspot report and update the GI location in HydroPlusConfig.xml to a spot with > 1 contributing cell." << endl;
		}
		// TP concentration
		M_kg_ha_yr = AverageContributingAreaLoad_calc(GI_ulcorner_row, GI_ulcorner_col, GI_lrcorner_row, GI_lrcorner_col, "TP");
		M_kg = M_kg_ha_yr * (A_m2 / ha_to_m2) * (t_s / year_to_seconds);
		M_mg = M_kg * Ratio_kg_to_mg;
		concentration_tp_mg_L = M_mg / (Q_m3 * m3_to_L); // converting m³ to L

		// TN concentration
		M_kg_ha_yr = AverageContributingAreaLoad_calc(GI_ulcorner_row, GI_ulcorner_col, GI_lrcorner_row, GI_lrcorner_col, "TN");
		M_kg = M_kg_ha_yr * (A_m2 / ha_to_m2) * (t_s / year_to_seconds);
		M_mg = M_kg * Ratio_kg_to_mg;
		concentration_tn_mg_L = M_mg / (Q_m3 * m3_to_L); // converting m³ to L

		// TSS concentration
		M_kg_ha_yr = AverageContributingAreaLoad_calc(GI_ulcorner_row, GI_ulcorner_col, GI_lrcorner_row, GI_lrcorner_col, "TSS");
		M_kg = M_kg_ha_yr * (A_m2 / ha_to_m2) * (t_s / year_to_seconds);
		M_mg = M_kg * Ratio_kg_to_mg;
		concentration_tss_mg_L = M_mg / (Q_m3 * m3_to_L); // converting m³ to L

		// Store the concentrations
		ca_concentration_tp_mg_L.push_back(concentration_tp_mg_L);
		//cout << "push back ca_concentration_tp_mg_L" << concentration_tp_mg_L << endl;
		ca_concentration_tn_mg_L.push_back(concentration_tn_mg_L);
		//cout << "push back ca_concentration_tn_mg_L" << concentration_tn_mg_L << endl;
		ca_concentration_tss_mg_L.push_back(concentration_tss_mg_L);
		//cout << "push back ca_concentration_tss_mg_L" << concentration_tss_mg_L << endl;
	}
}

void BufferSpatialCalc::replacePollutantConcentrations(Inputs* input) {
	// Iterate through each GI device
	for (size_t GI_Index = 0; GI_Index < input->GI_list.size(); ++GI_Index) {
		// Iterate over Pollutant_Name to find and replace TP, TN (TKN), and TSS concentrations
		for (size_t i = 0; i < input->Pollutant_Name.size(); ++i) {
			if (input->Pollutant_Name[i] == "TP") {
				cout << "Replacing TP concentration from " << input->Pollutant_ConcentrationMean_mg_p_L[i] << " mg/L to " << ca_concentration_tp_mg_L[GI_Index] << " mg/L." << endl;
				input->Pollutant_ConcentrationMean_mg_p_L[i] = ca_concentration_tp_mg_L[GI_Index];
			}
			else if (input->Pollutant_Name[i] == "TKN") {  // Assuming TKN is equivalent to TN
				cout << "Replacing TN concentration from " << input->Pollutant_ConcentrationMean_mg_p_L[i] << " mg/L to " << ca_concentration_tn_mg_L[GI_Index] << " mg/L." << endl;
				input->Pollutant_ConcentrationMean_mg_p_L[i] = ca_concentration_tn_mg_L[GI_Index];
			}
			else if (input->Pollutant_Name[i] == "TSS") {
				cout << "Replacing TSS concentration from " << input->Pollutant_ConcentrationMean_mg_p_L[i] << " mg/L to " << ca_concentration_tss_mg_L[GI_Index] << " mg/L." << endl;
				input->Pollutant_ConcentrationMean_mg_p_L[i] = ca_concentration_tss_mg_L[GI_Index];
			}
		}
	}
}


bool BufferSpatialCalc::isPixelValid(int i, Inputs* input) {
	if (TerrainProcessor::elevation_m_1D[i] == Inputs::NODATA_code) return false;
	if (input->LandCover_NLCD_Class[i] == Inputs::NODATA_code) return false;
	if (input->BufferSpatialWeightingStringParams["Flag_FlowDirection_Source"] == "NHD") {
		if (input->NHDfdr[i] == Inputs::NODATA_code || input->Waters01[i] == Inputs::NODATA_code) return false;
	}
	if (isSubsurfaceCalculationEnabled) {
		if (input->s_Ksat[i] == Inputs::NODATA_code || input->s_depth[i] == Inputs::NODATA_code) return false;
	}
	if (input->BufferSpatialWeightingStringParams["Flag_ElevationMethod_Subsurface"] == "WEM" && !(input->BufferSpatialWeightingNumericalParams["f_param"] > 0 && input->BufferSpatialWeightingNumericalParams["S_WTable_avg"] > 0)) {
		if (input->s_wtable[i] == Inputs::NODATA_code) return false;
	}
	return true;
}


//BufferSpatialCalc::BufferModel_HotSpotMaps_DiffusePollution_main function implements the hot spot mapping of diffuse or non-point source runoff
//Reference: Endreny, T., & Wood, E. F. (2003). Watershed Weighting of Export Coefficients to Map Critical Phosphorus Loading Areas. Journal of the American Water Resources Association, 39(1), 165-181. 
//Reference: Stephan, E. A., & Endreny, T. A. (2016). Weighting Nitrogen and Phosphorus Pixel Pollutant Loads to Represent Runoff and Buffering Likelihoods. JAWRA Journal of the American Water Resources Association, 52(2), 336-349. doi:doi:10.1111/1752-1688.12390
void BufferSpatialCalc::BufferModel_HotSpotMaps_DiffusePollution_main(Inputs* input) {
	//* Preparation Step 1 * Read in parameters and Checking the flags 
	//reading parameters and flags from input, check the settings
	nRows = Inputs::nRows;
	nCols = Inputs::nCols;
	NODATA_code = Inputs::NODATA_code;
	cellsize_m = Inputs::Length_Pixel_Side_m;;
	//read the flags of the surface contribution, to see if subsurface simulation is needed
	isSubsurfaceCalculationEnabled = (
		input->BufferSpatialWeightingNumericalParams["ECTP_surface_contribution_frac"] < 1 ||
		input->BufferSpatialWeightingNumericalParams["ECTN_surface_contribution_frac"] < 1 ||
		input->BufferSpatialWeightingNumericalParams["ECTSS_surface_contribution_frac"] < 1
		);
	if (isSubsurfaceCalculationEnabled) {
		// Ensure flow direction source is compatible with water table slope method
		if (input -> BufferSpatialWeightingStringParams["Flag_ElevationMethod_Subsurface"] == "WEM" && 
			input-> BufferSpatialWeightingStringParams["Flag_FlowDirection_Source"] != "DEM") {
			cout << "Warning: Flag_ElevationMethod_Subsurface is set to WEM, requiring DEM as the flow direction source." << endl;
			cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
			cout << "Explanation: Flag_FlowDirection_Source is now being forced to DEM.SSURGO water table elevation (WEM) slope requires flow direction source to be DEM." << endl;
			cout << "Correction: Set Flag_FlowDirection_Source to DEM." << endl;
			abort();
		}
		//check if Flag_BufferMethod_Subsurface and Flag_ElevationMethod_Subsurface provided
		if (input->BufferSpatialWeightingStringParams["Flag_BufferMethod_Subsurface"] != "LC" && input->BufferSpatialWeightingStringParams["Flag_BufferMethod_Subsurface"] != "t") {
			cout << "Warning: Subsurface calculations are enabled, but Flag_BufferMethod_Subsurface is missing." << endl;
			cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
			cout << "Correction: Set Flag_BufferMethod_Subsurface to 'LC' or 't' to specify subsurface calculation method." << endl;
			abort();
		}
		if (input->BufferSpatialWeightingStringParams["Flag_ElevationMethod_Subsurface"]!= "DEM" && input->BufferSpatialWeightingStringParams["Flag_ElevationMethod_Subsurface"] != "WEM") {
			cout << "Warning: Subsurface calculations are enabled, but Flag_ElevationMethod_Subsurface is missing." << endl;
			cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
			cout << "Correction: Set Flag_ElevationMethod_Subsurface to 'WEM' or 'DEM' to specify elevation source for subsurface." << endl;
			abort();
		}
	}
	
	//* Preparation Step 2 *  Handle DEM slope cleanup and align NODATA across maps
	//Replaces unrealistically small or invalid slope values with the median of positive slope values, storing the cleaned slope in slope_dem 
	slope_dem.clear();
	slope_dem.resize(nRows * nCols, Inputs::NODATA_code);
	double slope_min = Find_MinimumPositiveValue_1D_Vector_calc(TerrainProcessor::slopeorganizer_riserun1D);
	double slope_median = calculateMedian(TerrainProcessor::slopeorganizer_riserun1D);

	// Iterate through each cell in the map
	for (int i = 0; i < nRows * nCols; i++) {
		//replace the negative or inf dem slope with median value
		if (TerrainProcessor::elevation_m_1D[i] != Inputs::NODATA_code) {
			slope_dem[i] = (TerrainProcessor::slopeorganizer_riserun1D[i] < slope_min)
				? slope_median
				: TerrainProcessor::slopeorganizer_riserun1D[i];
		}
		else {
			slope_dem[i] = Inputs::NODATA_code;
		}

		//Ensures spatial alignment across all input maps by marking a pixel as NODATA in all relevant datasets if it is NODATA in any one of them.
		//This includes DEM, land cover (NLCD), flow direction (NHD), soil data (Ksat, depth), and subsurface water table (if used).
		//Logic Overview:
		// - If DEM pixel is NODATA, all dependent maps are also set to NODATA.
		// - If any layer (e.g., NLCD, NHDfdr, Ksat, s_wtable) has NODATA at a pixel,that pixel is excluded from processing by propagating NODATA across others.
		// Note: This step is essential for ensuring consistent per-pixel computation across maps with different NODATA footprints.
		if (!isPixelValid(i, input)) {
			TerrainProcessor::elevation_m_1D[i] = Inputs::NODATA_code;
			input->LandCover_NLCD_Class[i] = Inputs::NODATA_code;
			if (input->BufferSpatialWeightingStringParams["Flag_FlowDirection_Source"] == "NHD") {
				input->NHDfdr[i] = Inputs::NODATA_code;
				input->Waters01[i] = Inputs::NODATA_code;
			}
			if (isSubsurfaceCalculationEnabled) {
				input->s_Ksat[i] = Inputs::NODATA_code;
				input->s_depth[i] = Inputs::NODATA_code;
			}
			if (input->BufferSpatialWeightingStringParams["Flag_ElevationMethod_Subsurface"] == "WEM" && !(input->BufferSpatialWeightingNumericalParams["f_param"] > 0 && input->BufferSpatialWeightingNumericalParams["S_WTable_avg"] > 0)) {
				input->s_wtable[i] = Inputs::NODATA_code;
			}
		}
	}

	//replace the 2D DEM vector with the latest DEM map. If not the flow accumulation process will still accumulate those nodata cells
	ReplaceElevationWithVector(TerrainProcessor::elevation_m_1D);

	//saving the original dem flow direction to Inputs::DEMfdr, in case been replaced when redo the flow direction
	Inputs::DEMfdr.clear();
	//saving flow direction results to a 1D vector
	for (auto it = TerrainProcessor::FlowDirection_Map.begin(); it != TerrainProcessor::FlowDirection_Map.end(); it++) {
		for (auto it2 = (*it).begin(); it2 != (*it).end(); it2++) {
			Inputs::DEMfdr.push_back(*it2);
		}
	}

	//* Preparation Step 3  For GI case only
	//update landcover type if GI installed
	if (input->SimulationStringParams["Model_Selection"] == "SpatialBufferwGI") {
		updateLandCover(input);
		if (Inputs::SimulationStringParams["Algorithm_FlowDirection"] == "D8") {
			//write out Contributing Area, list of contributing cell ID list for each pixel
			BufferModel_IntermediateOutputs_writer("ContributingAreaPixels");
		}
	}

	// ******************************************************************************************************
	// Compute waters, reverse fdr and compute impervious weighted grid
	// ******************************************************************************************************
	//when Flag_FlowDirection_Source = NHD, user need to provide NHD flow direction, flowline and waterbody (NHD dem is optional)
	//when Flag_FlowDirection_Source = DEM, results generated from DEM

	//LI 12042024 moving this part forward to index calculator
	//if Flag_FlowDirection_Source = DEM, compute river, else read fdr via readWaterFeaturesMap() in Inputs
	if (input->BufferSpatialWeightingStringParams["Flag_FlowDirection_Source"] == "DEM") {
		//Command Prompt Notification
		cout << "Generating river layer from DEM flow accumulation. " << endl;
		//Call ReceivingWater_ZeroValue_CreateMap_from_DEM function sending FlowAccum
		ReceivingWater_ZeroValue_CreateMap_from_DEM();
		// Re-check conditions for saving contributing/dispersal area pixels
		TerrainProcessor::setContributingDispersalFlags();
	}
	//If Else Flag_FlowDirection_Source NHD, then water layer and flow direction are derived from NHD input files
	if (input->BufferSpatialWeightingStringParams["Flag_FlowDirection_Source"] == "NHD") {
		//Command Prompt Notification
		cout << "Generating water layer from nhdwaters10.asc. " << endl;
		//Call ReceivingWater_ZeroValue_CreateMap_from_NHD function sending input
		ReceivingWater_ZeroValue_CreateMap_from_NHD();
		cout << "Updating flow direction from nhdfdr.asc. " << endl;
		//Algorithm to replace DEM derived flow direction (FDR) with NHD FDR values, whem appropriate
		ReplaceFlowDirectionOrganizerWithVector(input->NHDfdr);
		// Re-check conditions for saving contributing/dispersal area pixels
		TerrainProcessor::setContributingDispersalFlags();
	}

	//accumulate the non-water area
	TerrainProcessor::FlowAccum(waters01, TerrainProcessor::DX_vec);
	//mark ContributingArea been saved and avoid running it again in the future accumulation
	TerrainProcessor::saveContributingDispersalArea = false;

	//saving the weighted FA result to FA_weighted_m
	for (auto row = TerrainProcessor::FlowAccumulation_Map.begin(); row != TerrainProcessor::FlowAccumulation_Map.end(); row++) {
		for (auto col = (*row).begin(); col != (*row).end(); col++) {
			FA_weighted_m.push_back(*col);
		}
	}

	if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1){
		//Write some intermediate outputs
		//BufferModel_IntermediateOutputs_writer("dem");
		BufferModel_IntermediateOutputs_writer("slope_dem");
		if (Inputs::BufferSpatialWeightingStringParams["Flag_FlowDirection_Source"] == "DEM") {
			BufferModel_IntermediateOutputs_writer("demfdr");
		}
		else {
			BufferModel_IntermediateOutputs_writer("nhdfdr");
			BufferModel_IntermediateOutputs_writer("nhdfdrr");
		}

		BufferModel_IntermediateOutputs_writer("waters01");
		//BufferModel_IntermediateOutputs_writer("waters0");
	
		BufferModel_IntermediateOutputs_writer("FA_weighted_m");
	}

	// ******************************************************************************************************
	// Compute RI_surf_i in Eq 6 of Stephan and Endreny (2016), from impervious Topographic Index, Eq 8 of Valeo, C. et al. (2000). Journal of Hydrology 228(1-2): 68-81.
	// ******************************************************************************************************
	//Eq 6 is RI_surf_i = ln(FlowAccumulation(impervious surfaces) / DEM_slope) = ln(FA_imp / DEM_slope
	//calculate surface runoff index
	RunoffIndex_Surface_calc(FA_weighted_m, slope_dem);	
	//Write intermediate outputs
	BufferModel_IntermediateOutputs_writer("RI_surf");

	// ******************************************************************************************************
	//Compute RI_sub in Eq 7 of Stephan and Endreny (2016), from Soil Topographic Index, Eq 10 in Sivapalan, M., et al. (1987). Water Resources Research 23(12): 2266-2278.
	//and Eq 1 in Endreny, T.A., et al. (2000).Hydrological Processes 14(1) : 1 - 20.
	// ******************************************************************************************************
	if (isSubsurfaceCalculationEnabled) {
		//calculate subsurface runoff index using dem slope, but Eq 7 is using local subsurface watertable
		//at this point WT_slope can't be calculated without estimating watertable during RI_sub calculation
		//Endreny and W0od 1999. Eq 1. S should be the pixel slope rather thanwatertable slope. Note that in BUffer.py it's using dem slope.
		if (input->BufferSpatialWeightingStringParams["Flag_ElevationMethod_Subsurface"] == "DEM") {
			RunoffIndex_Subsurface_calc(input, FA_weighted_m, slope_dem);
		}
		if (input->BufferSpatialWeightingStringParams["Flag_ElevationMethod_Subsurface"] == "WEM") {
			//compute the WEM to get groundwatertable_elevation1D
			ComputeWEM(input);
			//update the dem and re-compute the slope and flow direction
			//replace elevation[row][col] with groundwatertable_elevation1D  
			ReplaceElevationWithVector(groundwatertable_elevation1D);
			recomputeFlowDirection(input);
			//saving the slope result 
			slope_wem.clear();
			//saving flow direction results to a 1D vector to slope_wem
			for (auto it = TerrainProcessor::Slope_Ground_rise_run_mpm_Map.begin(); it != TerrainProcessor::Slope_Ground_rise_run_mpm_Map.end(); it++) {
				for (auto it2 = (*it).begin(); it2 != (*it).end(); it2++) {
					slope_wem.push_back(*it2);
				}
			}

			//redo the flow accumulation use WEM fdr
			//accumulate the non-water area
			TerrainProcessor::FlowAccum(waters01, TerrainProcessor::DX_vec);
			//saving the weighted FA result to FA_WEMfdr_temp
			vector<double> FA_WEMfdr_weighted_temp;
			for (auto row = TerrainProcessor::FlowAccumulation_Map.begin(); row != TerrainProcessor::FlowAccumulation_Map.end(); row++) {
				for (auto col = (*row).begin(); col != (*row).end(); col++) {
					FA_WEMfdr_weighted_temp.push_back(*col);
				}
			}
			RunoffIndex_Subsurface_calc(input, FA_WEMfdr_weighted_temp, slope_wem);
		}
		
		BufferModel_IntermediateOutputs_writer("RI_sub");

		//write intermediate outputs if Flag_ExtendedOutputs turned on
		if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1){
			//Write intermediate outputs
			BufferModel_IntermediateOutputs_writer("s_trans");
		}
	}
	// ******************************************************************************************************
	// Compute BI_surf using CADA 2003 and 2016
	// ******************************************************************************************************
	//get the negative DEM fdr and slope
	//replace elevation[row][col] with elevation_m_negative_1D  
	ReplaceElevationWithVector(TerrainProcessor::elevation_m_negative_1D);
	//update the dem and re-compute the flow direciton and slope
	recomputeFlowDirection(input);
	//saving the flow direction to negDEMfdr
	Inputs::negDEMfdr.clear();
	//saving flow direction results to a 1D vector
	for (auto it = TerrainProcessor::FlowDirection_Map.begin(); it != TerrainProcessor::FlowDirection_Map.end(); it++) {
		for (auto it2 = (*it).begin(); it2 != (*it).end(); it2++) {
			Inputs::negDEMfdr.push_back(*it2);
		}
	}
	//saving the slope to slope_negdem
	slope_negdem.clear();
	//saving slope to a 1D vector
	for (auto it = TerrainProcessor::Slope_Ground_rise_run_mpm_Map.begin(); it != TerrainProcessor::Slope_Ground_rise_run_mpm_Map.end(); it++) {
		for (auto it2 = (*it).begin(); it2 != (*it).end(); it2++) {
			slope_negdem.push_back(*it2);
		}
	}

	if (input->BufferSpatialWeightingStringParams["Flag_FlowDirection_Source"] == "DEM") {
		if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {
			BufferModel_IntermediateOutputs_writer("neg_dem");
			BufferModel_IntermediateOutputs_writer("demfdrr");
		}
	}
	
	//replace the fdr with the reverse NHD fdr
	if (input->BufferSpatialWeightingStringParams["Flag_FlowDirection_Source"] == "NHD") {
		//replace the flow direction organizer with reveresed NHD fdr
		ReplaceFlowDirectionOrganizerWithVector(input->NHDfdr_reverse);
	}
	
	//done with updating the reverse flow direction for the following flow accumulation

	//Start the BI_surf calc
	if (input->BufferSpatialWeightingStringParams["Flag_BufferMethod_Surface"] == "LC" ||
		input->BufferSpatialWeightingStringParams["Flag_BufferMethod_Subsurface"] == "LC") {
		//creating nutrient trapping map
		NutrientTrapping_CreateMap(input);

		if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {
			//Write intermediate outputs
			BufferModel_IntermediateOutputs_writer("p_release");
			BufferModel_IntermediateOutputs_writer("p_trapping");
		}
	}

	//if use land cover trapping method e.g. Compute BI_surf_i  T.A. and Wood, E.F., 2003
	if (input->BufferSpatialWeightingStringParams["Flag_BufferMethod_Surface"] == "LC") {
		//Determining the average slope for each of the pixel's dispersal area Based on the basic slope = rise/run principle
		Slope_DispersalArea_calc(slope_negdem);

		if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1){
			//Write intermediate outputs
			BufferModel_IntermediateOutputs_writer("surf_buffer_rise1");
			BufferModel_IntermediateOutputs_writer("surf_buffer_rise2");
			BufferModel_IntermediateOutputs_writer("surf_buffer_rise");
			BufferModel_IntermediateOutputs_writer("surf_buffer_run1");
			BufferModel_IntermediateOutputs_writer("surf_buffer_run");
			BufferModel_IntermediateOutputs_writer("surf_buffer_slope");
		}
		//Li 20231106 use the clear() method to remove all elements from the vectors and release the memory for BI_surf 2003 partI
		buffer_rise1.clear();
		buffer_rise2.clear();
		buffer_rise.clear();
		buffer_run1.clear();
		buffer_run.clear();
		
		//calculate filtered/hillslope
		filter_to_slope_calc();
		
		if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1){
			//Write intermediate outputs
			BufferModel_IntermediateOutputs_writer("surf_filter1");
			BufferModel_IntermediateOutputs_writer("surf_filter2");
			BufferModel_IntermediateOutputs_writer("surf_filter");
			BufferModel_IntermediateOutputs_writer("surf_filter_to_slope");
		}
		//Li 20231106 use the clear() method to remove all elements from the vectors and release the memory for BI_surf 2003 partIII
		filter1.clear();
		filter2.clear();
		filter.clear();
		//Calculate the final BI, where BI = ln(filtered/hillslope)
		BufferIndex_Surface_calc(filter_to_slope);
		BufferModel_IntermediateOutputs_writer("BI_surf");

		//Li 20231106 use the clear() method to remove all elements from the vectors and release the memory for BI_surf 2003 partIV
		buffer_slope.clear();
		filter_to_slope.clear();
	}

	// if use travel time method, e.g. Compute BI_surf in Eq 10 of Stephan and Endreny (2016)
	if (input->BufferSpatialWeightingStringParams["Flag_BufferMethod_Surface"] == "t") {
		
		ManningRoughness_HydraulicRadius_CreateMap(input);
		if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1){
			//Write intermediate outputs
			BufferModel_IntermediateOutputs_writer("Manningn");
			BufferModel_IntermediateOutputs_writer("Hydraulic_R");
		}

		TravelTime_Runoff_Surface_calc(slope_dem);
		if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1){
			//Write intermediate outputs
			BufferModel_IntermediateOutputs_writer("Time_surf1");
			BufferModel_IntermediateOutputs_writer("Time_surf2");
			BufferModel_IntermediateOutputs_writer("Time_surf");
		}

		BufferIndex_Surface_calc(Time_surf);
		BufferModel_IntermediateOutputs_writer("BI_surf");

		//Li 20231106 use the clear() method to remove all elements from the vectors and release the memory for BI_surf 2016
		Time_surf1.clear();
		Time_surf2.clear();
		Time_surf.clear();
	}
	// ******************************************************************************************************
	// Compute BI_sub_i Stephan and Endreny (2016)
	// ******************************************************************************************************
	if (isSubsurfaceCalculationEnabled) {
		//Compute s_wtable (m) as groundwater table depth using Eq 13 of Stephan and Endreny (2016)
		//generate s_wtable_est estimate based on user provided f and s_wtable_avg
		//Depth_GroundwaterTable_calc(input);
		//Write intermediate outputs
		//if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) BufferModel_IntermediateOutputs_writer("s_wtable_est");
	
		//saving the 2D slope result to 1D vector
		/*WTable_slope.clear();
		for (auto it = slopeorganizer_temp.begin(); it != slopeorganizer_temp.end(); it++) {
			for (auto it2 = (*it).begin(); it2 != (*it).end(); it2++) {
				WTable_slope.push_back(*it2);
			}
		}

		//Write intermediate outputs
		if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) BufferModel_IntermediateOutputs_writer("WTable_slope");
		*/


		//at this point, the flowdirectionorganizer[][] is negDEMdsr or Inputs::NHDfdr_reverse
		//to make it easy, doing DEM case first then do the WEM case
		if (input->BufferSpatialWeightingStringParams["Flag_ElevationMethod_Subsurface"] == "DEM") {
			//if use land cover trapping method e.g. Compute BI_surf_i  T.A. and Wood, E.F., 2003
			if (input->BufferSpatialWeightingStringParams["Flag_BufferMethod_Surface"] == "LC") {
				//Determining the average slope for each of the pixel's dispersal area Based on the basic slope = rise/run principle
				Slope_DispersalArea_calc(slope_negdem);

				if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {
					//Write intermediate outputs
					BufferModel_IntermediateOutputs_writer("sub_buffer_rise1");
					BufferModel_IntermediateOutputs_writer("sub_buffer_rise2");
					BufferModel_IntermediateOutputs_writer("sub_buffer_rise");
					BufferModel_IntermediateOutputs_writer("sub_buffer_run1");
					BufferModel_IntermediateOutputs_writer("sub_buffer_run");
					BufferModel_IntermediateOutputs_writer("sub_buffer_slope");
				}
				//Li 20231106 use the clear() method to remove all elements from the vectors and release the memory for BI_surf 2003 partI
				buffer_rise1.clear();
				buffer_rise2.clear();
				buffer_rise.clear();
				buffer_run1.clear();
				buffer_run.clear();

				//calculate filtered/hillslope
				filter_to_slope_calc();

				if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {
					//Write intermediate outputs
					BufferModel_IntermediateOutputs_writer("sub_filter1");
					BufferModel_IntermediateOutputs_writer("sub_filter2");
					BufferModel_IntermediateOutputs_writer("sub_filter");
					BufferModel_IntermediateOutputs_writer("sub_filter_to_slope");
				}
				//Li 20231106 use the clear() method to remove all elements from the vectors and release the memory for BI_surf 2003 partIII
				filter1.clear();
				filter2.clear();
				filter.clear();
				//Calculate the final BI, where BI = ln(filtered/hillslope)
				BufferIndex_Subsurface_calc(filter_to_slope);
				BufferModel_IntermediateOutputs_writer("BI_sub");

				//Li 20231106 use the clear() method to remove all elements from the vectors and release the memory for BI_surf 2003 partIV
				buffer_slope.clear();
				filter_to_slope.clear();
			}

			// if use travel time method, e.g. Compute BI_surf in Eq 10 of Stephan and Endreny (2016)
			if (input->BufferSpatialWeightingStringParams["Flag_BufferMethod_Surface"] == "t") {
				//Compute Time_sub (s) subsurface time for runoff
				TravelTime_Runoff_Subsurface_calc(input, slope_negdem);
				if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {
					//Write intermediate outputs
					BufferModel_IntermediateOutputs_writer("Time_sub1");
					BufferModel_IntermediateOutputs_writer("Time_sub2");
					BufferModel_IntermediateOutputs_writer("Time_sub");
				}
				//Calculate local Buffer index as natural logarithm of time 
				BufferIndex_Subsurface_calc(Time_sub);
				//Write intermediate outputs
				BufferModel_IntermediateOutputs_writer("BI_sub");

				//Li 20231106 use the clear() method to remove all elements from the vectors and release the memory for BI_surf 2016
				Time_sub1.clear();
				Time_sub2.clear();
				Time_sub.clear();
			}
		}


		//if use WEM, we need to get the slope neg WEM through the flow direction process and update the flow direction to negWEMfdr
		if (input->BufferSpatialWeightingStringParams["Flag_ElevationMethod_Subsurface"] == "WEM") {
			//replace elevation[row][col] with neg_groundwatertable_elevation1D  
			ReplaceElevationWithVector(neg_groundwatertable_elevation1D);
			//update the dem and re-compute the flow direciton and slope
			recomputeFlowDirection(input);
			//saving the slope to slope_negwem
			slope_negwem.clear();
			//saving slope to a 1D vector
			for (auto it = TerrainProcessor::Slope_Ground_rise_run_mpm_Map.begin(); it != TerrainProcessor::Slope_Ground_rise_run_mpm_Map.end(); it++) {
				for (auto it2 = (*it).begin(); it2 != (*it).end(); it2++) {
					slope_negwem.push_back(*it2);
				}
			}

			//if use land cover trapping method e.g. Compute BI_surf_i  T.A. and Wood, E.F., 2003
			if (input->BufferSpatialWeightingStringParams["Flag_BufferMethod_Surface"] == "LC") {
				//Determining the average slope for each of the pixel's dispersal area Based on the basic slope = rise/run principle
				Slope_DispersalArea_calc(slope_negwem);

				if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {
					//Write intermediate outputs
					BufferModel_IntermediateOutputs_writer("sub_buffer_rise1");
					BufferModel_IntermediateOutputs_writer("sub_buffer_rise2");
					BufferModel_IntermediateOutputs_writer("sub_buffer_rise");
					BufferModel_IntermediateOutputs_writer("sub_buffer_run1");
					BufferModel_IntermediateOutputs_writer("sub_buffer_run");
					BufferModel_IntermediateOutputs_writer("sub_buffer_slope");
				}
				//Li 20231106 use the clear() method to remove all elements from the vectors and release the memory for BI_surf 2003 partI
				buffer_rise1.clear();
				buffer_rise2.clear();
				buffer_rise.clear();
				buffer_run1.clear();
				buffer_run.clear();

				//calculate filtered/hillslope
				filter_to_slope_calc();

				if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {
					//Write intermediate outputs
					BufferModel_IntermediateOutputs_writer("sub_filter1");
					BufferModel_IntermediateOutputs_writer("sub_filter2");
					BufferModel_IntermediateOutputs_writer("sub_filter");
					BufferModel_IntermediateOutputs_writer("sub_filter_to_slope");
				}
				//Li 20231106 use the clear() method to remove all elements from the vectors and release the memory for BI_surf 2003 partIII
				filter1.clear();
				filter2.clear();
				filter.clear();
				//Calculate the final BI, where BI = ln(filtered/hillslope)
				BufferIndex_Subsurface_calc(filter_to_slope);
				BufferModel_IntermediateOutputs_writer("BI_sub");

				//Li 20231106 use the clear() method to remove all elements from the vectors and release the memory for BI_surf 2003 partIV
				buffer_slope.clear();
				filter_to_slope.clear();
			}
				

			// if use travel time method, e.g. Compute BI_surf in Eq 10 of Stephan and Endreny (2016)
			if (input->BufferSpatialWeightingStringParams["Flag_BufferMethod_Surface"] == "t") {
				//Compute Time_sub (s) subsurface time for runoff
				TravelTime_Runoff_Subsurface_calc(input, slope_negwem);
				if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {
					//Write intermediate outputs
					BufferModel_IntermediateOutputs_writer("Time_sub1");
					BufferModel_IntermediateOutputs_writer("Time_sub2");
					BufferModel_IntermediateOutputs_writer("Time_sub");
				}
				//Calculate local Buffer index as natural logarithm of time 
				BufferIndex_Subsurface_calc(Time_sub);
				//Write intermediate outputs
				BufferModel_IntermediateOutputs_writer("BI_sub");

				//Li 20231106 use the clear() method to remove all elements from the vectors and release the memory for BI_surf 2016
				Time_sub1.clear();
				Time_sub2.clear();
				Time_sub.clear();
			}
		}
	}

	// ******************************************************************************************************
	// Calculations for buffer index ratios
	// ******************************************************************************************************
	
	//Buffering Index (BI) is ratio of BI_average to BI_local, large ratios lead to larger weights for Export Coefficient
	RunoffIndex_BufferIndex_Surface_Ratio_calc();
	if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1){
		BufferModel_IntermediateOutputs_writer("RI_surf_ratio");
		BufferModel_IntermediateOutputs_writer("BI_surf_ratio");
	}
	if (isSubsurfaceCalculationEnabled) {
		RunoffIndex_BufferIndex_Subsurface_Ratio_calc();
		if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1){
			BufferModel_IntermediateOutputs_writer("RI_sub_ratio");
			BufferModel_IntermediateOutputs_writer("BI_sub_ratio");
		}
	}

	// ******************************************************************************************************
	// Calculations for Export Coefficient values weighted
	// ******************************************************************************************************
	
	//Export Coefficient values for total nitrogen (TN) and total phosphorus (TP), unweighted, from White et al 
	ExportCoefficient_unweighted_CreateMap(input);
	//compute the weighted map
	Normalize_EC_EMC_maps(input, isSubsurfaceCalculationEnabled);
	//write out the normalized weighted EC map
	if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {
		if (isSubsurfaceCalculationEnabled && input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {
			BufferModel_IntermediateOutputs_writer("surf_ec_tn_wtd");
			BufferModel_IntermediateOutputs_writer("surf_ec_tp_wtd");
			BufferModel_IntermediateOutputs_writer("surf_ec_tss_wtd");
			BufferModel_IntermediateOutputs_writer("sub_ec_tn_wtd");
			BufferModel_IntermediateOutputs_writer("sub_ec_tp_wtd");
			BufferModel_IntermediateOutputs_writer("sub_ec_tss_wtd");
		}
	//write out the normalized unweighted EC map
		BufferModel_IntermediateOutputs_writer("ec_tn_vec");
		BufferModel_IntermediateOutputs_writer("ec_tp_vec");
		BufferModel_IntermediateOutputs_writer("ec_tss_vec");
	}
	
	BufferModel_IntermediateOutputs_writer("total_ec_tn_wtd");
	BufferModel_IntermediateOutputs_writer("total_ec_tp_wtd");
	BufferModel_IntermediateOutputs_writer("total_ec_tss_wtd");
	BufferModel_IntermediateOutputs_writer("total_emc_tn_wtd");
	BufferModel_IntermediateOutputs_writer("total_emc_tp_wtd");
	BufferModel_IntermediateOutputs_writer("total_emc_tss_wtd");


	//only find and report the hotspot when SimulationNumericalParams["Flag_ExtendedOutputs"] == 1
	if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {
		//find the hotspot of the map, find pixels >99th percentile of the whole map
		//TP
		Find_PercentileValue_1D_Vector_calc(input, total_ec_tp_wtd, ec_tp_vec, isSubsurfaceCalculationEnabled);
		//generate the hotspot map
		BufferModel_IntermediateOutputs_writer("hotspot_tp");
		//genrate a txt resport of the hotspots, including the pixel ID, row, col, and the value. This is a temp vector, that will be cleared after been wrriten out.
		writeHotspotReport("TP_hotspotreport.csv", hotspot_report_temp);

		//TN
		Find_PercentileValue_1D_Vector_calc(input, total_ec_tn_wtd, ec_tn_vec,  isSubsurfaceCalculationEnabled);
		//generate the hotspot map
		BufferModel_IntermediateOutputs_writer("hotspot_tn");
		//genrate a txt resport of the hotspots, including the pixel ID, row, col, and the value. This is a temp vector, that will be cleared after been wrriten out.
		writeHotspotReport("TN_hotspotreport.csv", hotspot_report_temp);

		//TSS
		Find_PercentileValue_1D_Vector_calc(input, total_ec_tss_wtd, ec_tss_vec, isSubsurfaceCalculationEnabled);
		//generate the hotspot map
		BufferModel_IntermediateOutputs_writer("hotspot_tss");
		//genrate a txt resport of the hotspots, including the pixel ID, row, col, and the value. This is a temp vector, that will be cleared after been wrriten out.
		writeHotspotReport("TSS_hotspotreport.csv", hotspot_report_temp);
	}

	//tracking contributing/dispersal area only works for D8 
	// Vectors are temporarily stored and cleaned up after writing output for individual pixels.
	if (input->SimulationStringParams["Algorithm_FlowDirection"] == "D8") {
		//write out contributing and dispersal area
		for (auto i : input->BufferSpatialWeightingStringParams) {
			if (i.first.find("PrintRowCol") != string::npos) {
				//Note: user provided/ArcGIS row col start from 1
				int writeout_row = atoi(i.second.substr(0, i.second.find(",")).c_str()) - 1;
				//Note: user provided/ArcGIS row col start from 1
				int writeout_col = atoi(i.second.substr(i.second.find(",") + 1, i.second.size() - i.second.find(",")).c_str()) - 1;
				//check if PrintRowCol has values
				if (writeout_row >= 0 && writeout_col >= 0 ) {
					// Check if the specified row and column are within the map boundaries
					if (writeout_row < nRows && writeout_col < nCols) {
						string C_writeout_filename = "ContributingAreaPixels_map_" + to_string(writeout_row + 1) + "_" + to_string(writeout_col + 1);
						string D_writeout_filename = "DispersalAreaPixels_map_" + to_string(writeout_row + 1) + "_" + to_string(writeout_col + 1);

						// Temporary storage of contributing and dispersal area vectors for individual GI pixels.
						// These vectors are saved as temporary files and automatically cleaned up and resized at the beginning of the new loop
						ContributingAreaPixels_map_writeout.clear();
						DispersalAreaPixels_map_writeout.clear();
						ContributingAreaPixels_map_writeout.resize(nRows * nCols, NODATA_code);
						DispersalAreaPixels_map_writeout.resize(nRows * nCols, NODATA_code);

						//create the maps for the pixels privided by user
						Contributing_dispersal_CreateMap(writeout_row, writeout_col);
						//write out the map
						BufferModel_IntermediateOutputs_writer(C_writeout_filename);
						BufferModel_IntermediateOutputs_writer(D_writeout_filename);
					}
					else {
						// The specified row and/or column is outside the map boundaries, do not write the maps
						cout << "Warning: Specified row or column is outside the map boundaries. Contributing and dispersal area maps will not be written." << endl;
					}
				}
			}
		}

		if (Inputs::SimulationStringParams["Model_Selection"] == "SpatialBufferwGI") {

			// Accumulation of contributing and dispersal area vectors for the entire GI block.
			// For multi-pixel GI blocks, the contributing and dispersal areas for each pixel are
			// combined into a single map to reflect the cumulative effect across all pixels.
			// Keep this clean up section outside of the loop to ensure it's accumulating correctly
			ContributingAreaPixels_map_writeout.clear();
			DispersalAreaPixels_map_writeout.clear();
			ContributingAreaPixels_map_writeout.resize(nRows * nCols, NODATA_code);
			DispersalAreaPixels_map_writeout.resize(nRows * nCols, NODATA_code);

			//write out contributing and dispersal area for GI pixels
			// Iterate through each GI area and update land cover
			for (size_t gi = 0; gi < input->GI_list.size(); ++gi) {
				string GI_type = input->GI_list[gi];
				int GI_ulcorner_row = input->GI_ulcorner_rows[gi];
				int GI_ulcorner_col = input->GI_ulcorner_cols[gi];
				int GI_lrcorner_row = input->GI_lrcorner_rows[gi];
				int GI_lrcorner_col = input->GI_lrcorner_cols[gi];

				//create the name of the output maps
				string C_writeout_filename = "ContributingAreaPixels_map_GI_" + to_string(GI_ulcorner_row + 1) + "_" + to_string(GI_ulcorner_col + 1) + "_" + to_string(GI_lrcorner_row + 1) + "_" + to_string(GI_lrcorner_col + 1);
				string D_writeout_filename = "DispersalAreaPixels_map_GI_" + to_string(GI_ulcorner_row + 1) + "_" + to_string(GI_ulcorner_col + 1) + "_" + to_string(GI_lrcorner_row + 1) + "_" + to_string(GI_lrcorner_col + 1);

				// Iterate over each pixel in the block
				for (int gi_row = GI_ulcorner_row; gi_row <= GI_lrcorner_row; ++gi_row) {
					for (int gi_col = GI_ulcorner_col; gi_col <= GI_lrcorner_col; ++gi_col) {

						// Check if the specified row and column are within the map boundaries
						if (gi_row >= 0 && gi_row < nRows && gi_col >= 0 && gi_col < nCols) {
							//create the maps, for  multiple gi pixels, it's cumulative
							Contributing_dispersal_CreateMap(gi_row, gi_col);
						}
						else {
							// The specified row and/or column is outside the map boundaries, do not write the maps
							cout << "Warning: GI is outside the map boundaries. Contributing and dispersal area maps will not be written." << endl;
						}
					}
				}
				//write out the map for each gi block
				BufferModel_IntermediateOutputs_writer(C_writeout_filename);
				BufferModel_IntermediateOutputs_writer(D_writeout_filename);
			}
		}
	}
	//replace the TP, TN(TKN), and TSS values with the calculated concentrations, and then show a program message for each replacement
	//for the following statistical Hydro GI modeling
	if (input->SimulationStringParams["Model_Selection"] == "SpatialBufferwGI" && input->SimulationStringParams["Algorithm_FlowDirection"] == "D8") {
		ContributingConcentrationValues_calc(input);
		replacePollutantConcentrations(input);
		//turn on the flag to allow statistical Hydro to run the process of reduce_EC_EMC_after_GI
		flag_reduce_EC_EMC_after_GI = true;
	}
}

