#include "WeatherProcessor.h"
#include "Inputs.h"


//Note: Flag_HottestDaySimulated = 1 searches all weather files for hottest hour, allowing differences between stations across years, months, and days
//Note: Model output file z_StationInfoMap.csv contains variable SimulationDate_Hottest_GD that reports the hottest date used for each station
//Note: Flag_MultipleStations = 1 computes distance between stations and map pixels that use the same projection, e.g., Albers_Conical_Equal_Area, Lambert_Azimuthal_Equal_Area, or UTM

//Variables defined in WeatherProcessor.h need to be declared in WeatherProcessor.cpp
int WeatherProcessor::Actual_hottest_GDH_index, WeatherProcessor::SimulationDateStart_GDH, WeatherProcessor::SimulationDateStop_GDH, WeatherProcessor::Date_YYYYMMDD, WeatherProcessor::Time_HHMMSS, WeatherProcessor::timeStepMetData, WeatherProcessor::previous_YYYYMMDD, WeatherProcessor::previous_HMS, WeatherProcessor::SimulationDate_Hottest_GD, WeatherProcessor::DH_data_Cooling_total_count, WeatherProcessor::DH_data_Heating_total_count, WeatherProcessor::Climate_Year_Start_YYYY;
double WeatherProcessor::DH_Cooling_total_sum, WeatherProcessor::DH_Heating_total_sum, WeatherProcessor::DegreeHour_Threshold_K, WeatherProcessor::Tair_mesoScale_final_weighted_K, WeatherProcessor::H_total_weighted_W_p_m2;
string WeatherProcessor::ReadStation_path, WeatherProcessor::ReadStationID;
long double  AlbersEasting, AlbersNorthing;
bool WeatherProcessor::Flag_DegreeHour_Simulated = false, WeatherProcessor::Flag_Tair_degrees_F = false, WeatherProcessor::Flag_Tdew_degrees_F = false;
vector<string> WeatherProcessor::OBJECTID_shortest, WeatherProcessor::USAF_shortest, WeatherProcessor::WBAN_shortest, WeatherProcessor::Path_folder_shortest, WeatherProcessor::TI_option_shortest, WeatherProcessor::SimulationTime_HMS_AllRecord;
vector<long double> WeatherProcessor::Easting_shortest, WeatherProcessor::Northing_shortest, WeatherProcessor::Elev_m_shortest, WeatherProcessor::Distance_shortest, WeatherProcessor::Aspect_deg_shortest, WeatherProcessor::Slope_deg_shortest, WeatherProcessor::NLCD_Class_shortest, WeatherProcessor::TC_per_shortest, WeatherProcessor::IC_per_shortest, WeatherProcessor::SVeg_in_Gap_frac_shortest, WeatherProcessor::Station_AH_Flux_Wpm2_shortest, WeatherProcessor::Distance_Pixel_to_Station_vec_km, WeatherProcessor::Weight_StationInverseDistance_vec_frac, WeatherProcessor::Distance_InverseSquared_vec;
vector<double> WeatherProcessor::RadiationNet_Wpm2_AllRecord, WeatherProcessor::Tair_C_AllRecord, WeatherProcessor::Tdew_C_AllRecord, WeatherProcessor::Rain_m_AllRecord, WeatherProcessor::Snow_m_AllRecord, WeatherProcessor::WindSpd_mps_AllRecord, WeatherProcessor::WindDir_deg_AllRecord, WeatherProcessor::AtmPres_kPa_AllRecord, WeatherProcessor::Precip_mpdt_AllRecord, WeatherProcessor::DegreeHour_Cooling_All_abs_K_vec, WeatherProcessor::DegreeHour_Heating_All_abs_K_vec, WeatherProcessor::DegreeHour_Year_Cooling_Ratio_vec, WeatherProcessor::DegreeHour_Year_Heating_Ratio_vec, WeatherProcessor::DegreeHour_Hour_Cooling_K_vec, WeatherProcessor::DegreeHour_Hour_Heating_K_vec, WeatherProcessor::DegreeHour_Month_Heating_Ratio_vec, WeatherProcessor::DegreeHour_Month_Cooling_Ratio_vec;
vector<pair<int, int>> WeatherProcessor::DaylightSaving_StartStopPairs_YYYYMMDDHH;
vector<int> WeatherProcessor::SimulationDate_GD_AllRecord, WeatherProcessor::SimulationDate_GDH_AllRecord, WeatherProcessor::SimulationDate_HH_AllRecord, WeatherProcessor::SimulationDate_HH;
string WeatherProcessor::SimulationTime_Hottest_HMS, WeatherProcessor::Time_HHMMSS_str, WeatherProcessor::DegreeHour_statisticMethod;
map<string, vector<double>> WeatherProcessor::WeatherMap;
map<int, double>  WeatherProcessor::DH_Cooling_Year_sum, WeatherProcessor::DH_Cooling_Month_sum, WeatherProcessor::DH_Cooling_Hour_sum, WeatherProcessor::DH_Heating_Year_sum, WeatherProcessor::DH_Heating_Month_sum, WeatherProcessor::DH_Heating_Hour_sum;
map<int, int>  WeatherProcessor::DH_Year_Cooling_count, WeatherProcessor::DH_Year_Heating_count, WeatherProcessor::DH_Month_Cooling_count, WeatherProcessor::DH_Month_Heating_count, WeatherProcessor::DH_Hour_Cooling_count, WeatherProcessor::DH_Hour_Heating_count;
map<int, vector<double>>  WeatherProcessor::DegreeHour_Cooling_Year_abs_K_vec, WeatherProcessor::DegreeHour_Cooling_Month_abs_K_vec, WeatherProcessor::DegreeHour_Cooling_Hour_abs_K_vec, WeatherProcessor::DegreeHour_Heating_Year_abs_K_vec, WeatherProcessor::DegreeHour_Heating_Month_abs_K_vec, WeatherProcessor::DegreeHour_Heating_Hour_abs_K_vec, WeatherProcessor::DegreeHour_Year_Cooling_Ratio_StationID_vec_map, WeatherProcessor::DegreeHour_Year_Heating_Ratio_StationID_vec_map, WeatherProcessor::DegreeHour_Month_Cooling_Ratio_StationID_vec_map, WeatherProcessor::DegreeHour_Month_Heating_Ratio_StationID_vec_map, WeatherProcessor::DegreeHour_Hour_Cooling_K_StationID_vec_map, WeatherProcessor::DegreeHour_Hour_Heating_K_StationID_vec_map, WeatherProcessor::DegreeHour_Year_Cooling_Ratio_Weighted_vec_map, WeatherProcessor::DegreeHour_Year_Heating_Ratio_Weighted_vec_map, WeatherProcessor::DegreeHour_Month_Cooling_Ratio_Weighted_vec_map, WeatherProcessor::DegreeHour_Month_Heating_Ratio_Weighted_vec_map, WeatherProcessor::DegreeHour_Hour_Cooling_K_Weighted_vec_map, WeatherProcessor::DegreeHour_Hour_Heating_K_Weighted_vec_map;
//Variables defined in WeatherProcessor.h need to be declared in WeatherProcessor.cpp; relates to static vs non-static function
map<string, vector<string> > WeatherProcessor::StationInfoMap;
map<int, map<int, vector<string> > >  WeatherProcessor::PixelNearestStation;
int  WeatherProcessor::Global_Adjusted_Timesteps = INT_MAX;

//Note: WeatherStationAttributeTable.csv Header is: OBJECTID,USAF,WBAN,Station_Easting_m,Station_Northing_m,Station_BlockGroup_ID,Station_TI_option,Station_Elev_m,Station_Slope_deg,Station_Aspect_deg,Station_NLCD_Class,Station_TC_percent,Station_IC_percent,Station_SVeg_in_Gap_frac,Station_AnthropogenicHeat_Flux_Qtot_Avg_Wpm2,Station_Path_folder
const string WeatherStationAttributeTable_header = "OBJECTID,USAF,WBAN,Station_Easting_m,Station_Northing_m,Station_BlockGroup_ID,Station_TI_option,Station_Elev_m,Station_Slope_deg,Station_Aspect_deg,Station_NLCD_Class,Station_TC_percent,Station_IC_percent,Station_SVeg_in_Gap_frac,Station_AnthropogenicHeat_Flux_Qtot_Avg_Wpm2,Station_Path_folder";

//Inputs::readMultipleWeatherStations_List_CSV function reads from WeatherStationAttributeTable.csv the weather station ID, location, and characteristics
//Note: Flag_MultipleStations = 1 computes distance between stations and map pixels that use the same projection, e.g., Albers_Conical_Equal_Area, Lambert_Azimuthal_Equal_Area, or UTM
//Note: Consider refactor to clearly name the input variables, e.g., changing Station_TC_percent to TreeCanopy_Cover_percent; not all variables (e.g., Station_SVeg_in_Gap_frac (%)) are used in function
void WeatherProcessor::readMultipleWeatherStations_List_CSV(Inputs* input)
{
	//RefWeatherLocationAttributeFile_str set to value of HydroPlusConfig.xml element RefWeatherLocationAttributeFile
	string RefWeatherLocationAttributeFile_str = input->SimulationStringParams["RefWeatherLocationAttributeFile"];
	// Read station coordinates based on RefWeatherLocationCoordinateFormat
	string coord_format_flag = "PCS";  // Default to PCS

	// Check if user provided the flag
	if (input->SimulationStringParams.find("RefWeatherLocationCoordinateFormat") != input->SimulationStringParams.end()) {
		coord_format_flag = input->SimulationStringParams["RefWeatherLocationCoordinateFormat"];
	}
	else {
		cout << "Warning: RefWeatherLocationCoordinateFormat not found in HydroPlusConfig.xml. Defaulting to PCS (Easting/Northing in meters)." << endl;
	}
	
	//If RefWeatherLocationAttributeFile_str is not populated, then HydroPlusConfig.xml is missing file, and look in local directory
	//Note: This option is for those our NUCFAC National runs
	if (RefWeatherLocationAttributeFile_str.size() < 1) {
		RefWeatherLocationAttributeFile_str = "WeatherStationAttributeTable.csv";
	}
	//fstream C++ function used to assign WeatherStation as the RefWeatherLocationAttributeFile_str file
	fstream WeatherStation(RefWeatherLocationAttributeFile_str);
	//If WeatherStation does not exist, then alert user
	if (!WeatherStation) {
		cout << endl;
		cout << "Warning: This program cannot open the required file " << RefWeatherLocationAttributeFile_str << "." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: The HydroPlusConfig.xml element RefWeatherLocationAttributeFile should point to WeatherStationAttributeTable.csv." << endl;
		cout << "Explanation: The RefWeatherLocationAttributeFile should not be open in Excel when trying to access it with HydroPlus.exe." << endl;
		cout << "Explanation: The WeatherStationAttributeTable.csv requires the following columns when Flag_MultipleStations is 1:" << endl;
		cout << WeatherStationAttributeTable_header << endl;
		cout << "Correction: Update the WeatherStationAttributeTable.csv to have the correct values." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	string Station_OBJECTID, Station_USAF, Station_WBAN, Station_Location_x_str, Station_Location_y_str, Station_BlockGroup_ID_str, Station_TI_option, Station_Elevation_m_str, Station_Slope_deg_str, Station_Aspect_deg_str, Station_NLCD_Class_str, Station_TC_percent_str, Station_IC_percent_str, Station_SVeg_in_Gap_per_str, Station_AH_Flux_Wpm2_str, Station_Path_folder;
	double Station_Easting_m, Station_Northing_m, Station_Latitude_dd, Station_Longitude_dd, Station_Elevation_m, Station_Slope_deg, Station_Aspect_deg, Station_NLCD_Class, Station_TC_percent, Station_IC_percent, Station_SVeg_in_Gap_frac, Station_AH_Flux_Wpm2;
	int Station_BlockGroup_ID;
	string header, temp_line, temp_string;

	//expected_column_count of 16 defined as columns in WeatherStationAttributeTable.csv 
	const int expected_column_count = 16;
	//line_num is intialized to 0, the header line, but increments beyond header line below
	int line_num = 0;
	getline(WeatherStation, header);
	//While the next data exist, read a line in as the temp_line
	while (getline(WeatherStation, temp_line)) {
		line_num = line_num + 1;
		//read cells within this line
		istringstream readlinedata(temp_line);

		if (!check_column_count(temp_line, expected_column_count, line_num)) {
			cout << "Warning: WeatherStationAttributeTable.csv has the wrong number of colunns." << endl;
			cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
			cout << "Explanation: The WeatherStationAttributeTable.csv requires the following columns when Flag_MultipleStations is 1:" << endl;
			cout << WeatherStationAttributeTable_header << endl;
			cout << "Correction: Update the WeatherStationAttributeTable.csv to have the correct values." << endl;
			//Call abort function, which ends the HydroPlus.exe simulation
			abort();
		}

		getline(readlinedata, Station_OBJECTID, ',');
		getline(readlinedata, Station_USAF, ',');
		getline(readlinedata, Station_WBAN, ',');
		getline(readlinedata, Station_Location_x_str, ',');
		getline(readlinedata, Station_Location_y_str, ',');
		getline(readlinedata, Station_BlockGroup_ID_str, ',');
		getline(readlinedata, Station_TI_option, ',');
		getline(readlinedata, Station_Elevation_m_str, ',');
		getline(readlinedata, Station_Slope_deg_str, ',');
		getline(readlinedata, Station_Aspect_deg_str, ',');
		getline(readlinedata, Station_NLCD_Class_str, ',');
		getline(readlinedata, Station_TC_percent_str, ',');
		getline(readlinedata, Station_IC_percent_str, ',');
		getline(readlinedata, Station_SVeg_in_Gap_per_str, ',');
		getline(readlinedata, Station_AH_Flux_Wpm2_str, ',');
		getline(readlinedata, Station_Path_folder, ',');

		//Station_Path_folder is updated by EnsureTrailingBackslash function to ensure it has no parenthesis and ends with backslash
		Station_Path_folder = EnsureTrailingBackslash(Station_Path_folder);
		//safe_stod function checks with isMissing function, then converts to double
		if (coord_format_flag == "PCS") {
			// User provided Easting/Northing in meters
			Station_Easting_m = safe_stod(Station_Location_x_str, "Station_Easting_m", line_num);
			Station_Northing_m = safe_stod(Station_Location_y_str, "Station_Northing_m", line_num);
		}
		else if (coord_format_flag == "GCS") {
			// User provided Lat/Lon in decimal degrees
			Station_Latitude_dd = safe_stod(Station_Location_x_str, "Station_Latitude_dd", line_num);
			Station_Longitude_dd = safe_stod(Station_Location_y_str, "Station_Longitude_dd", line_num);
			
			//Convert Latitude/Longitude to Easting/Northing using Albers projection parameters
			//MapProjection_WKID is read from HydroPlusConfig.xml as SimulationNumericalParams["MapProjection_WKID"]
			int MapProjection_WKID = Inputs::SimulationNumericalParams["MapProjection_WKID"];

			// Look up projection parameters from supported WKID map
			if (Inputs::AlbersParamMap.find(MapProjection_WKID) != Inputs::AlbersParamMap.end()) {
				//Create a local copy of the projection parameters for this WKID
				ProjectionParams_Albers projectionParams = Inputs::AlbersParamMap[MapProjection_WKID];
				//Call the projection function to convert Lat/Lon (dd) to Easting/Northing (meters)
				pair<double, double> projected = Inputs::Projection_LatitudeLongitude_to_Albers(Station_Latitude_dd, Station_Longitude_dd, projectionParams);
				//Assign projected coordinates
				Station_Easting_m = projected.first;
				Station_Northing_m = projected.second;
			}
			else if (Inputs::LambertParamMap.find(MapProjection_WKID) != Inputs::LambertParamMap.end()) {
				//Create a local copy of the projection parameters for this WKID
				ProjectionParams_Lambert projectionParams = Inputs::LambertParamMap[MapProjection_WKID];
				//Call the projection function to convert Lat/Lon (dd) to Easting/Northing (meters)
				pair<double, double> projected = Inputs::Projection_LatitudeLongitude_to_Lambert(Station_Latitude_dd, Station_Longitude_dd, projectionParams);
				//Assign projected coordinates
				Station_Easting_m = projected.first;
				Station_Northing_m = projected.second;
			}
			else if (((MapProjection_WKID >= 32601 && MapProjection_WKID <= 32660) || (MapProjection_WKID >= 32701 && MapProjection_WKID <= 32760))){
				//Create a local copy of the projection parameters for this WKID
				ProjectionParams_UTM projectionParams = Inputs::initializeUTMParams(MapProjection_WKID);
				//Call the projection function to convert Lat/Lon (dd) to Easting/Northing (meters)
				pair<double, double> projected = Inputs::Projection_LatitudeLongitude_to_UTM(Station_Latitude_dd, Station_Longitude_dd, projectionParams);
				//Assign projected coordinates
				Station_Easting_m = projected.first;
				Station_Northing_m = projected.second;
			
			}

			//cout << Station_USAF << " " << Station_USAF << " Station_Latitude_dd " << Station_Latitude_dd << " Station_Longitude_dd " << Station_Longitude_dd << " Station_Easting_m " << to_string(Station_Easting_m) << " Station_Northing_m " << to_string(Station_Northing_m) << endl;

		}
		else {
			cout << "Error: Invalid RefWeatherLocationCoordinateFormat value (" << coord_format_flag << "). Should be PCS or GCS." << endl;
			abort();
		}

		//safe_stoi function checks with isMissing function, then converts to int
		Station_BlockGroup_ID = safe_stoi(Station_BlockGroup_ID_str, "Station_BlockGroup_ID", line_num);
		Station_Elevation_m = safe_stod(Station_Elevation_m_str, "Station_Elevation_m", line_num);
		Station_Slope_deg = safe_stod(Station_Slope_deg_str, "Station_Slope_deg", line_num);
		Station_Aspect_deg = safe_stod(Station_Aspect_deg_str, "Station_Aspect_deg", line_num);
		Station_NLCD_Class = safe_stod(Station_NLCD_Class_str, "Station_NLCD_Class", line_num);
		Station_TC_percent = safe_stod(Station_TC_percent_str, "Station_TC_percent", line_num);
		Station_IC_percent = safe_stod(Station_IC_percent_str, "Station_IC_percent", line_num);
		Station_SVeg_in_Gap_frac = safe_stod(Station_SVeg_in_Gap_per_str, "Station_SVeg_in_Gap_frac", line_num);
		Station_AH_Flux_Wpm2 = safe_stod(Station_AH_Flux_Wpm2_str, "Station_AH_Flux_Wpm2", line_num);

		check_range(Station_BlockGroup_ID, "Station_BlockGroup_ID", line_num, -10000, 10000);
		check_range(Station_Elevation_m, "Station_Elevation_m", line_num, -1000, 50000);
		check_range(Station_Slope_deg, "Station_Slope_deg", line_num, 0, 180);
		check_range(Station_NLCD_Class, "Station_NLCD_Class", line_num, 11, 95);
		check_range(Station_Aspect_deg, "Station_Aspect_deg", line_num, 0, 360);
		check_range(Station_TC_percent, "Station_TC_percent", line_num, 0, 100);
		check_range(Station_IC_percent, "Station_IC_percent", line_num, 0, 100);
		check_range(Station_SVeg_in_Gap_frac, "Station_SVeg_in_Gap_frac", line_num, 0, 1);
		check_range(Station_AH_Flux_Wpm2, "Station_AH_Flux_Wpm2", line_num, 0, 10000); 

		input->Station_OBJECTID_vec.push_back(Station_OBJECTID);
		input->Station_USAF_vec.push_back(Station_USAF);
		input->Station_WBAN_vec.push_back(Station_WBAN);
		input->Station_Easting_m_vec.push_back(Station_Easting_m);
		input->Station_Northing_m_vec.push_back(Station_Northing_m);
		input->Station_BlockGroup_ID_vec.push_back(Station_BlockGroup_ID);
		//Station_TI_option_vec directly takes Station_TI_option and is always treated as 'avg'
		input->Station_TI_option_vec.push_back(Station_TI_option);
		input->Station_Elev_m_vec.push_back(Station_Elevation_m);
		input->Station_Slope_deg_vec.push_back(Station_Slope_deg);
		input->Station_Aspect_deg_vec.push_back(Station_Aspect_deg);
		input->Station_NLCD_Class_vec.push_back(Station_NLCD_Class);
		input->Station_TC_per_vec.push_back(Station_TC_percent);
		input->Station_IC_per_vec.push_back(Station_IC_percent);
		input->Station_SVeg_in_Gap_frac_vec.push_back(Station_SVeg_in_Gap_frac);
		input->Station_AH_Flux_Wpm2_vec.push_back(Station_AH_Flux_Wpm2);
		input->Station_Path_folder_vec.push_back(Station_Path_folder);
	}
	//Check if the MultipleStations_NumberOfStationsAveraged <= weather station provided in the attribute table
	if (input->Station_OBJECTID_vec.size() < int(input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"])) {
		cout << "Warning: The number of weather stations in the attribute table is less than the number requested for averaging." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: The simulation requires at least "
			<< int(input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"])
			<< " weather stations to perform averaging, but only "
			<< input->Station_OBJECTID_vec.size() << " were provided in the attribute table." << endl;
		cout << "Correction: Add more stations to WeatherStationAttributeTable.csv or reduce MultipleStations_NumberOfStationsAveraged in HydroPlusConfig.xml." << endl;
		abort();
	}
}

//Distance_Calc_EastingNorthing function computes distance between stations and map pixels in 2 steps, 1st distance between each point, 2nd the square root of those distances squared and summed
//Note: Flag_MultipleStations = 1 computes distance between stations and map pixels that use the same projection, e.g., Albers_Conical_Equal_Area, Lambert_Azimuthal_Equal_Area, or UTM
//Note: The Station_easting_m and Station_northing_m are read from the RefWeatherLocationAttributeFile, and PixelCenter_easting_m and PixelCenter_northing_m are read from the input maps
long double WeatherProcessor::Distance_Calc_EastingNorthing(long double PixelCenter_easting_m, long double PixelCenter_northing_m, long double Station_easting_m, long double Station_northing_m)
{
	//Distance_easting_m is distance along easting between station and pixel center
	long double Distance_easting_m = Station_easting_m - PixelCenter_easting_m;
	//Distance_northing_m is distance along northing between station and pixel center
	long double Distance_northing_m = Station_northing_m - PixelCenter_northing_m;

	//Distance_Pixel_to_Station_m (m) uses standard distance formula, taking the square root of the sum of squared distances in the northing and easting directions 
	double Distance_Pixel_to_Station_m = sqrt(Distance_easting_m * Distance_easting_m + Distance_northing_m * Distance_northing_m);
	return Distance_Pixel_to_Station_m * Ratio_km_to_m; //convert to km
}

//sorting and keeping track of indexes
//Example vector<double> data = { 4.2, 2.1, 5.9, 3.3 };
//vector<size_t> sorted_idx = sort_indexes(data);
// sorted_idx will be: {1, 3, 0, 2}
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(idx.begin(), idx.end(), [&Vector](size_t i1, size_t i2) {return Vector[i1] < Vector[i2]; });

	return idx;
}

//convertHrMinSecToInt function converts string HrMinSec (HH:MM:SS) to an integer (HHMMSS)
//Note: Function used by readWeatherFile
int convertHrMinSecToInt(const string& timeStr) {
	//timeHHMMSS_str makes a copy of the timeStr 
	string timeHHMMSS_str = timeStr;
	//timeHHMMSS_str is parsed to remove colons
	timeHHMMSS_str.erase(remove(timeHHMMSS_str.begin(), timeHHMMSS_str.end(), ':'), timeHHMMSS_str.end());
	//If timeHHMMSS_str.empty() throw error
	if (timeHHMMSS_str.empty()) {
		throw invalid_argument("Empty time string in convertHrMinSecToInt");
	}
	//timeHHMMSS_str is converted to integer
	int timeHHMMSS_int = stoi(timeHHMMSS_str);
	return timeHHMMSS_int;
}

//WeatherStationDistance_sort function will find nearest weather stations
//Note: The Station_Easting_m and Station_Northing_m are read from the RefWeatherLocationAttributeFile, and PixelCenter_easting_m and PixelCenter_northing_m are read from the input maps
//Note: The function will work when stations and map pixels use the same projection, e.g., Albers_Conical_Equal_Area, Lambert_Azimuthal_Equal_Area, or UTM
void WeatherProcessor::WeatherStationDistance_sort(Inputs* input, long double PixelCenter_easting_m, long double PixelCenter_northing_m) {
	//MultipleStations_NumberOfStationsAveraged defined as input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"], number of stations used in spatial weighting
	int MultipleStations_NumberOfStationsAveraged = int(input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]);

	//For Loop through all members in Station_Easting_m_vec.size()
	//Note: Station_Easting_m_vec is same size as Station_Northing_m_vec, and these contain the data from the reference weather stations
	for (int i = 0; i < input->Station_Easting_m_vec.size(); i++) {

		//Station_Easting_x_m (m) defined by input-> Station_Easting_m_vec (m) vector, obtained from RefWeatherLocationAttributeFile in coordinates identical to those of map
		long double Station_Easting_x_m = input->Station_Easting_m_vec[i];
		//Station_Northing_y_m (m) defined by input-> Station_Northing_m_vec (m) vector, obtained from RefWeatherLocationAttributeFile in coordinates identical to those of map
		long double Station_Northing_y_m = input->Station_Northing_m_vec[i];
		//Call Distance_Calc_EastingNorthing function, which requires the four variables passed to contain the same map projection (e.g., Albers, Lampert, UTM, etc.)
		//Distance_Pixel_to_Station_km (km) is computed with Distance_Calc_LatLong function, which sends back distance in km, taking inputs in distance m
		long double Distance_Pixel_to_Station_km = Distance_Calc_EastingNorthing(PixelCenter_easting_m, PixelCenter_northing_m, Station_Easting_x_m, Station_Northing_y_m);
		//Distance_Pixel_to_Station_vec_km has Distance_Pixel_to_Station_km added
		Distance_Pixel_to_Station_vec_km.push_back(Distance_Pixel_to_Station_km);
	}

	// Sort once outside the loop
	vector<size_t> sorted_station_indices = sort_indexes(Distance_Pixel_to_Station_vec_km);
	//For Loop through all MultipleStations_NumberOfStationsAveraged, which constains the number of nearest stations to include in weighting
	for (int i = 0; i < MultipleStations_NumberOfStationsAveraged; i++) {
		//StationID contains the weather station index of the nearest station, e.g. 12
		//Note: sort_indexes returns the location in the vector of the ith smallest distance; e.g., index locaiton of first, second, and third nearest
		int StationID = static_cast<int>(sorted_station_indices[i]);
		//Distance_WeatherStation_km contains weather station ditance to the nearest station, e.g. 16093
		long double Distance_WeatherStation_km = Distance_Pixel_to_Station_vec_km[StationID];

		//Variable vectors have members added for each of the weather stations identified to influence the pixel
		Distance_shortest.push_back(Distance_WeatherStation_km);
		OBJECTID_shortest.push_back(input->Station_OBJECTID_vec[StationID]);
		USAF_shortest.push_back(input->Station_USAF_vec[StationID]);
		WBAN_shortest.push_back(input->Station_WBAN_vec[StationID]);
		Easting_shortest.push_back(input->Station_Easting_m_vec[StationID]);
		Northing_shortest.push_back(input->Station_Northing_m_vec[StationID]);
		TI_option_shortest.push_back(input->Station_TI_option_vec[StationID]);
		Elev_m_shortest.push_back(input->Station_Elev_m_vec[StationID]);
		Slope_deg_shortest.push_back(input->Station_Slope_deg_vec[StationID]);
		Aspect_deg_shortest.push_back(input->Station_Aspect_deg_vec[StationID]);
		NLCD_Class_shortest.push_back(input->Station_NLCD_Class_vec[StationID]);
		TC_per_shortest.push_back(input->Station_TC_per_vec[StationID]);
		IC_per_shortest.push_back(input->Station_IC_per_vec[StationID]);
		SVeg_in_Gap_frac_shortest.push_back(input->Station_SVeg_in_Gap_frac_vec[StationID]);
		Station_AH_Flux_Wpm2_shortest.push_back(input->Station_AH_Flux_Wpm2_vec[StationID]);
		Path_folder_shortest.push_back(input->Station_Path_folder_vec[StationID]);
		//albersEasting_shortest.push_back(input->Station_Easting_m_vec[StationID]);
		//albersNorthing_shortest.push_back(input->Station_Northing_m_vec[StationID]);
	}
	//Distance_Pixel_to_Station_vec_km vector is cleared
	Distance_Pixel_to_Station_vec_km.clear();
}

//WeatherDataReadandProcessing function called from SimulationCoordinator to process weather within each time step of each pixel
//Note: Flag_MultipleStations = 1 computes distance between stations and map pixels that use the same projection
void WeatherProcessor::WeatherDataReadandProcessing(Inputs* input, DataOrganizer* organizer) {
	
	//If Model_Selection equals SpatialTemperatureHydro then 
	if (input->SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro") {
		//MapPixelCalcWeightsandRead will give centroid for each pixel
		MapPixelCalcWeightsandRead(input, organizer);
	}
	//Else If Model_Selection equals StatisticalHydro then 
	else if (input->SimulationStringParams["Model_Selection"] == "StatisticalHydro") {
		//MapCenterCalcWeightsandRead will give centroid for map
		MapCenterCalcWeightsandRead(input);
	}
	//MapPixelNearestWeatherStationWriter function called to Write output of weather preprocessor
	MapPixelNearestWeatherStationWriter(input);
	//WeatherStationInformationWriter function called to Write output of weather preprocessor
	WeatherStationInformationWriter(input);
	//resizeWeatherVectors function called to resize weather vectors to the size of total time step
	resizeWeatherVectors(input);

	//write DH result for each station
	if (Flag_DegreeHour_Simulated) {
		//WriteDegreeHourToCSV(bool isMultipleStation);
		WriteDegreeHourToCSV(input->SimulationNumericalParams["Flag_MultipleStations"] == 1);
	}
}

//readWeatherFile function reads in Weather.csv time series data for HydroPlus models for all models and conditions 
//Note: Header 9 variables of Weather.csv: YYYYMMDD,HH:MM:SS,Temperature_Air_C(C),Temperature_DewPoint_C(C),Radiation_Net_Wpm2(W/m^2),Wind_Speed_mps(m/s),Pressure_Atmosphere_kPa(kPa),Precipitation_Rate_mph(m/h),Snow_DepthOnGround_m(m)
//Note: i-Tree Energy Weather.csv header row has additional field after Wind_Speed_mps(m/s), Wind_Direction_deg(deg)
//Note: Radiation data aretypically computed in i-Tree WeatherPrep using NREL National Solar Radiation Database methods (1995)
//Note: NCEI Global Historical Climatology Network (GHCN) hourly data is expected to replace their integrated surface database (ISD) hourly weather data
void WeatherProcessor::readWeatherFile(Inputs* input) {
	
	//if ReadStation_path.empty() then assign to InputDirectory for calls from Input
	if (ReadStation_path.empty()) {
		ReadStation_path = input->SimulationStringParams["InputDirectory"];
	}

	//readWeather created as part of fstream, defined as the ReadStation_path and file name for Weather.csv
	fstream readWeather(ReadStation_path + "Weather.csv");

	//if readWeather does not exist, alert user
	if (!readWeather) {
		cout << "Warning: Weather.csv could not be opened." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: When the Weather.csv file is missing there is no meteorological data to drive the simulation." << endl;
		cout << "Correction: Create a Weather.csv file, either manually or perhaps using the i-Tree WeatherPrep.exe utility." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}
	//Initialize variables
	input->TotalRain_m = 0.0;
	input->TotalPrecip_m = 0.0;
	input->TotalSnow_m = 0.0;
	string lineFromCSV;
	string tempStr = "";
	string temp;
	string Date_YYYYMMDD_str;
	int date_InputDataRow_HH;
	int date_InputDataRow_GDH;
	int Counter_inputRow = 0;
	//timeStep_InputDataRow_sec initialized to HydroPlusConfig.xml TimeStep_sec
	bool Flag_SavingSimulationTimeRange = false;
	bool Flag_ProcessTimeStep = false;
	bool Flag_HottestDaySimulated = false;
	bool Flag_Update_SimulationDate_GD_GDH = true;
	Inputs::InitializeLAIVariables();

	//getline reads in header from readWeather, stores in tempStr
	//Note: getline function of fstream library takes characters from 1st item & stores them in 2nd item, until delimitation character found
	getline(readWeather, tempStr);

	//Flag_Tair_degrees_F is true if tempStr contains Temperature_Air_F or Temperature_DewPoint_F
	Flag_Tair_degrees_F = (tempStr.find("Temperature_Air_F") != string::npos);
	Flag_Tdew_degrees_F = (tempStr.find("Temperature_DewPoint_F") != string::npos);

	//numColumnsInFile within tempStr using .begin and .end, adding 1 to move from base 0 to 1 for count
	int numColumnsInFile = count(tempStr.begin(), tempStr.end(), ',') + 1;
	//Weather_Columns_count = 9 for all Model_Selection types except EnergyLoad
	int Weather_Columns_count = 9;
	//Weather_Column_names_F contain Tair_F, Tdew_F
	string Weather_Column_names_F = "YYYYMMDD,HH:MM:SS,Temperature_Air_F(F),Temperature_DewPoint_F(F),"
		"Radiation_Net_Wpm2(W/m^2),Wind_Speed_mps(m/s),Pressure_Atmosphere_kPa(kPa),"
		"Precipitation_Rate_mph(m/h),Snow_DepthOnGround_m(m)";
	//Weather_Column_names_C contain Tair_C, Tdew_C
	string Weather_Column_names_C = "YYYYMMDD,HH:MM:SS,Temperature_Air_C(C),Temperature_DewPoint_C(C),"
		"Radiation_Net_Wpm2(W/m^2),Wind_Speed_mps(m/s),Pressure_Atmosphere_kPa(kPa),"
		"Precipitation_Rate_mph(m/h),Snow_DepthOnGround_m(m)";

	//If Model_Selection is CoolBuilding set to Weather_Columns_count = 10
	if (input->SimulationStringParams["Model_Selection"] == "CoolBuilding") {
		Weather_Columns_count = 10;
		//Weather_Column_names_F also contain WindDirection_deg(deg)
		Weather_Column_names_F = "YYYYMMDD,HH:MM:SS,Temperature_Air_F(F),Temperature_DewPoint_F(F),"
			"Radiation_Net_Wpm2(W/m^2),Wind_Speed_mps(m/s),WindDirection_deg(deg),"
			"Pressure_Atmosphere_kPa(kPa),Precipitation_Rate_mph(m/h),Snow_DepthOnGround_m(m)";
		//Weather_Column_names_C also contain WindDirection_deg(deg)
		Weather_Column_names_C = "YYYYMMDD,HH:MM:SS,Temperature_Air_C(C),Temperature_DewPoint_C(C),"
			"Radiation_Net_Wpm2(W/m^2),Wind_Speed_mps(m/s),WindDirection_deg(deg),"
			"Pressure_Atmosphere_kPa(kPa),Precipitation_Rate_mph(m/h),Snow_DepthOnGround_m(m)";
	}
	//Weather_Column_names for error message
	string Weather_Column_names;
	//If (Flag_Tair_degrees_F || Flag_Tdew_degrees_F) then error-message template matches the header units
	if (Flag_Tair_degrees_F || Flag_Tdew_degrees_F) {
		Weather_Column_names = Weather_Column_names_F;
	}
	else {
		Weather_Column_names = Weather_Column_names_C;
	}

	//If numColumnsInFile is not equal to Weather_Columns_count then warn and abort
	if (numColumnsInFile != Weather_Columns_count) {
		cout << endl;
		cout << "Error: Weather.csv has an incorrect number of data columns." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: The HydroPlus model needs" << Weather_Columns_count << "columns of data to properly assign Weather.csv variables." << endl;
		cout << "Correction: Update the Weather.csv file to contain the following comma separated variables:" << endl;
		cout << Weather_Column_names << endl;
		//abort simulation due to incorrect CSV format
		abort();
	}

	//Flag_OutputTime_DST initialized to 0
	int Flag_OutputTime_DST = 0;
	try {
		//if (input->SimulationNumericalParams.count("Flag_OutputTime_DST") > 0); attempt to convert if key exists; count > 0 if element 0 or 1
		if (input->SimulationNumericalParams.count("Flag_OutputTime_DST") > 0) {
			//Flag_OutputTime_DST = input->SimulationNumericalParams["Flag_OutputTime_DST"]; provides value
			Flag_OutputTime_DST = input->SimulationNumericalParams["Flag_OutputTime_DST"];
		}
	}
	//catch (const exception& e) 
	catch (const exception& e) {
		//cerr << "Warning: Invalid Flag_OutputTime_DST value; defaulting to 0. Error: " << e.what() << endl;
		cerr << "Warning: Invalid Flag_OutputTime_DST value; defaulting to 0. Error: " << e.what() << endl;
		Flag_OutputTime_DST = 0;
	}

	//If Flag_OutputTime_DST == 1 then enter
	if (Flag_OutputTime_DST == 1) {
		//Get_DaylightSavingTime_Periods(input->SimulationDateStart_YYYYMMDDHH, input->SimulationDateStop_YYYYMMDDHH, DaylightSaving_StartStopPairs_YYYYMMDDHH)
		//Note: Call Get_DaylightSavingTime_Periods to coordinate vectors
		Get_DaylightSavingTime_Periods(input->SimulationDateStart_YYYYMMDDHH, input->SimulationDateStop_YYYYMMDDHH, DaylightSaving_StartStopPairs_YYYYMMDDHH);
	}

	//Flag_DegreeHour_Simulated initialized to false
	Flag_DegreeHour_Simulated = false;
	//If Flag_CoolAir_AnthropogenicHeat_Flux,Flag_AH_Flux_Qcr_Qncr_not_Qtot, and Flag_DegreeHour_Dependent_AH_Flux_Qcr are 1 then ensure DegreeHour defined ...
	if (input->SimulationNumericalParams["Flag_CoolAir_AnthropogenicHeat_Flux"] == 1 &&
		input->SimulationNumericalParams["Flag_AH_Flux_Qcr_Qncr_not_Qtot"] == 1 &&
		input->SimulationNumericalParams["Flag_DegreeHour_Dependent_AH_Flux_Qcr"] == 1) {

		//If DegreeHour_Threshold_K is missing then alert user
		if (input->SimulationNumericalParams.count("DegreeHour_Threshold_K") == 0) {
			cout << "Warning: HydroPlusConfig.xml element Flag_DegreeHour_Dependent_AH_Flux_Qcr is 1 but element DegreeHour_Threshold_K is missing." << endl;
			cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
			cout << "Explanation: The HydroPlus model requires DegreeHour_Threshold_K if Flag_DegreeHour_Dependent_AH_Flux_Qcr is 1." << endl;
			cout << "Explanation: DegreeHour_Threshold_K (Kelvin) is the threshold for cooling or heating degree days." << endl;
			cout << "Correction: Include element DegreeHour_Threshold_K in HydroPlusConfig.xml; 291.483 K is often used." << endl;
			abort();
		}
		//Else proceed to assign DegreeHour_Threshold_K 
		else {
			Flag_DegreeHour_Simulated = true;
			//DegreeHour_Threshold_K = input->SimulationNumericalParams["DegreeHour_Threshold_K"]
			//Note: DegreeHour_Threshold_K is typically set to the equivalent of 65 F
			DegreeHour_Threshold_K = input->SimulationNumericalParams["DegreeHour_Threshold_K"];
			//DegreeHour_statisticMethod = "median" calculation method; alternative is mean
			DegreeHour_statisticMethod = "median";
		}
	}

	//while (getline(readWeather, lineFromCSV)) to read and process each row of data as lineFromCSV
	while (getline(readWeather, lineFromCSV)) {

		//If lineFromCSV (string) is empty then continue, which reads next line of file; this skips empty lines
		if (lineFromCSV.empty()) {
			continue;
		}

		//Quick check: if any two commas appear with nothing between them → malformed row
		if (lineFromCSV.find(",,") != string::npos) {
			cout << "\nERROR: Weather.csv row " << Counter_inputRow
				<< " contains empty fields (\",,\").\n";
			cout << "Line:\n  \"" << lineFromCSV << "\"\n";
			abort();
		}

		//istringstream lineWeatherFile(lineFromCSV) creates object lineWeatherFile from string lineFromCSV
		istringstream lineWeatherFile(lineFromCSV);

		//Date_YYYYMMDD_str is obtained as 1st string before comma, using getline lineWeatherFile
		getline(lineWeatherFile, Date_YYYYMMDD_str, ',');
		//Time_HHMMSS_str is obtained as 2nd string before comma, using getline lineWeatherFile
		getline(lineWeatherFile, Time_HHMMSS_str, ',');

		//Counter_inputRow advanced
		Counter_inputRow = Counter_inputRow + 1;
		
		//Date_YYYYMMDD is obtained using stoi function with Date_YYYYMMDD_str
		Date_YYYYMMDD = stoi(Date_YYYYMMDD_str);
		//Time_HHMMSS extracted with stoi function with Time_HHMMSS_str and .substr functions, 1st number is position, 2nd is length
		Time_HHMMSS = stoi(Time_HHMMSS_str.substr(0, 2) + Time_HHMMSS_str.substr(3, 2) + Time_HHMMSS_str.substr(6, 2));
		//date_InputDataRow_GDH is created by Date_YYYYMMDD * 100 creating space for adding first two digits of Time_HHMMSS
		date_InputDataRow_GDH = (Date_YYYYMMDD * 100) + (Time_HHMMSS / 10000);
		//date_InputDataRow_HH integer from string with atoi function, extracting 2 digits of temp w substr function; c_str() points to temp
		date_InputDataRow_HH = atoi(Time_HHMMSS_str.substr(0, 2).c_str());

		//If Counter_inputRow equals 1 then check if Weather.csv first date is not greater than HydroPlusConfig.xml StartDate_YYYYMMDD
		if (Counter_inputRow == 1) {
			//timeStepMetData set to 0
			timeStepMetData = 0;
			//Check if year matches HydroPlusConfig.xml StartDate_YYYYMMDD Year (and below if Weather.csv contains full YYYYMMDD called for in HydroPlusConfig.xml)
			//HydroPlusConfig_StartDay_str_YYYYMMDD is string form of HydroPlusConfig.xml input StartDate_YYYYMMDD
			string HydroPlusConfig_StartDay_str_YYYYMMDD = to_string(input->SimulationNumericalParams["StartDate_YYYYMMDD"]);
			//WeatherInput_StartDay_YYYY is string form of Weather.csv input date column
			string WeatherInput_StartDay_str_YYYYMMDD = to_string(Date_YYYYMMDD);
			//WeatherInput_StartDay_YYYYMMDD defined as integer of Date_YYYYMMDD input
			int WeatherInput_StartDay_YYYYMMDD = Date_YYYYMMDD;
			//HydroPlusConfig_StartDay_YYYYMMDD is integer form of substring from HydroPlusConfig_StartDay_str_YYYYMMDD, taking first 8 elements
			int HydroPlusConfig_StartDay_YYYYMMDD = atoi(HydroPlusConfig_StartDay_str_YYYYMMDD.substr(0, 8).c_str());

			//if HydroPlusConfig_StartDay_YYYYMMDD > HydroPlusConfig_StartDay_YYYYMMDD then error with inputs or expectations
			if (WeatherInput_StartDay_YYYYMMDD > HydroPlusConfig_StartDay_YYYYMMDD) {
				cout << "Warning: The Weather.csv 1st column for date, YYYYMMDD, has a starting date of " << WeatherInput_StartDay_YYYYMMDD << " ..." << endl;
				cout << "... which is later than the HydroPlusConfig.xml parameter StartDate_YYYYMMDD, set to " << HydroPlusConfig_StartDay_YYYYMMDD << "." << endl;
				cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
				cout << "Explanation: The HydroPlus model needs the Weather.csv to contain the date given in HydroPlusConfig.xml." << endl;
				cout << "Correction: Either generate new Weather.csv or change the StartDate_YYYYMMDD in the HydroPlusConfig.xml file." << endl;
				//Call abort function, which ends the HydroPlus.exe simulation
				abort();
			}
		}

		//If Flag_HottestDaySimulated equals 1 then 
		//Note: Flag_HottestDaySimulated = 1 searches all weather files for hottest hour, allowing differences between stations across years, months, and days
		//Note: Model output file z_StationInfoMap.csv contains variable SimulationDate_Hottest_GD that reports the hottest date used for each station
		//Note: If left DegreeHour_Threshold_K empty, the model will use only the temperatures within the simulation date range. If a value is provided (>0), the model will read all available temperature data and computeaverage cooling and heating Degree Hours relative to this threshold.
		if (input->SimulationNumericalParams["Flag_HottestDaySimulated"] == 1 ) {
			//SimulationDate_GD_AllRecord vector contains Date_YYYYMMDD, list of Gregorian dates (YYYYMMDD) from input file
			SimulationDate_GD_AllRecord.push_back(static_cast<int>(Date_YYYYMMDD));
			//SimulationDate_GDH_AllRecord vector contains date_InputDataRow_GDH, list of Gregorian dates and hour (YYYYMMDDHH) from input file
			SimulationDate_GDH_AllRecord.push_back(static_cast<int>(date_InputDataRow_GDH));
			//SimulationDate_HH_AllRecord vector contains date_InputDataRow_HH, hour data from input file
			SimulationDate_HH_AllRecord.push_back(static_cast<int>(date_InputDataRow_HH));
			//SimulationTime_HMS_AllRecord vector contains Time_HHMMSS_str, time as HH:MM:SS string from input file
			SimulationTime_HMS_AllRecord.push_back(Time_HHMMSS_str);
			//Flag_SavingSimulationTimeRange set to false
			Flag_SavingSimulationTimeRange = false;
			//Flag_ProcessTimeStep set to false; to avoid writing time steps when not within time range
			Flag_ProcessTimeStep = false;
			//Save data to _AllRecord vectors
			Flag_HottestDaySimulated = true;
			//If date_InputDataRow_GDH <= SimulationDateStop_GDH && date_InputDataRow_GDH >= SimulationDateStart_GDH then
			//Note: date_InputDataRow_GDH is from Weather.csv & SimulationDateStop_YYYYMMDDHH is from HydroPlusConfig.xml
			if (date_InputDataRow_GDH >= input->SimulationDateStart_YYYYMMDDHH && date_InputDataRow_GDH <= input->SimulationDateStop_YYYYMMDDHH) {

				//Flag_ProcessTimeStep set to true; to allow writing time steps when within time range
				Flag_ProcessTimeStep = true;
				//Flag_Update_SimulationDate_GD_GDH is false to not update SimulationDate_YYYYMMDD and SimulationDate_GDH
				Flag_Update_SimulationDate_GD_GDH = false;
				//Call Update_SimulationTime_Vectors(input, Date_YYYYMMDD, Time_HHMMSS_str, date_InputDataRow_HH, date_InputDataRow_GDH, Flag_OutputTime_DST, DaylightSaving_StartStopPairs_YYYYMMDDHH, Flag_Update_SimulationDate_GD_GDH)
				Update_SimulationTime_Vectors(input, Date_YYYYMMDD, Time_HHMMSS_str, date_InputDataRow_HH, date_InputDataRow_GDH, Flag_OutputTime_DST, DaylightSaving_StartStopPairs_YYYYMMDDHH, Flag_Update_SimulationDate_GD_GDH);

				//SimulationTime_HH vector will contain the HH:MM:SS values
				//Note: Variable only used in this WeatherProcessor.cpp when Flag_HottestDaySimulated is true
				SimulationDate_HH.push_back(static_cast<int>(date_InputDataRow_HH));
			}
			
			//VectorsOfWeatherData function called for Flag_SavingSimulationTimeRange as false, Flag_ProcessTimeStep as true, Flag_HottestDaySimulated as true
			//Note: This is for multiple stations with flag for hottest day simulated
			//Note: Each line of Weather.csv is read for hottest hour settings, not only those within the start and stop dates
			VectorsOfWeatherData(input, lineWeatherFile, Flag_SavingSimulationTimeRange, Flag_ProcessTimeStep, Flag_HottestDaySimulated);
		}
		//Else if Flag_HottestDaySimulated is false ...
		//... If Flag_DegreeHour_Simulated is true & date_InputDataRow_GDH not within SimulationDateStart_YYYYMMDDHH and SimulationDateStop_YYYYMMDDHH
		//Note: When Flag_DegreeHour_Simulated is true it requires reading entire Weather.csv to determine degree hour statistics
		//Note: This is to read weather.csv that is outside the simulation date range 
		else if (Flag_DegreeHour_Simulated && ((date_InputDataRow_GDH < input->SimulationDateStart_YYYYMMDDHH) || (date_InputDataRow_GDH > input->SimulationDateStop_YYYYMMDDHH))) {
			//Flag_SavingSimulationTimeRange set to false
			Flag_SavingSimulationTimeRange = false;
			//Flag_ProcessTimeStep set to false
			Flag_ProcessTimeStep = false;
			//VectorsOfWeatherData function called for Flag_SavingSimulationTimeRange as false, Flag_ProcessTimeStep as false, Flag_HottestDaySimulated as false
			//Note: This is for multiple stations with flag for hottest day not simulated
			VectorsOfWeatherData(input, lineWeatherFile, Flag_SavingSimulationTimeRange, Flag_ProcessTimeStep, Flag_HottestDaySimulated);
		}
		//Else if Flag_HottestDaySimulated is false and ...
		//... date_InputDataRow_GDH is within SimulationDateStart_YYYYMMDDHH and SimulationDateStop_YYYYMMDDHH
		else if (date_InputDataRow_GDH >= input->SimulationDateStart_YYYYMMDDHH && date_InputDataRow_GDH <= input->SimulationDateStop_YYYYMMDDHH) {

			//Flag_Update_SimulationDate_GD_GDH is true to update SimulationDate_YYYYMMDD and SimulationDate_GDH
			Flag_Update_SimulationDate_GD_GDH = true;
			//Call Update_SimulationTime_Vectors(input, Date_YYYYMMDD, Time_HHMMSS_str, date_InputDataRow_HH, date_InputDataRow_GDH, Flag_OutputTime_DST, DaylightSaving_StartStopPairs_YYYYMMDDHH, Flag_Update_SimulationDate_GD_GDH)
			Update_SimulationTime_Vectors(input, Date_YYYYMMDD, Time_HHMMSS_str, date_InputDataRow_HH, date_InputDataRow_GDH, Flag_OutputTime_DST, DaylightSaving_StartStopPairs_YYYYMMDDHH, Flag_Update_SimulationDate_GD_GDH);

			//Flag_SavingSimulationTimeRange set to true
			Flag_SavingSimulationTimeRange = true;
			//Flag_ProcessTimeStep set to true
			Flag_ProcessTimeStep = true;
			//VectorsOfWeatherData function called for Flag_SavingSimulationTimeRange as true, Flag_ProcessTimeStep as true, Flag_HottestDaySimulated as false
			//Note: This is for single stations with flag for hottest day not simulated
			VectorsOfWeatherData(input, lineWeatherFile, Flag_SavingSimulationTimeRange, Flag_ProcessTimeStep, Flag_HottestDaySimulated);
		}
	} //while loop of all time ends

	//if Flag_DegreeHour_Simulated is true then call DegreeHour_CreateRatios_From_Sums_by_YearMonthHour function
	//Note: Degree Hours used in Anthropogenic Heat Flux algorithm, adjusting flux rates based on expected heating and cooling
	if (Flag_DegreeHour_Simulated) {
		if (input->SimulationNumericalParams["Flag_MultipleStations"] != 1) {
			//DegreeHour_CreateRatios_From_Sums_by_YearMonthHour function computes the degree hours using DegreeHour_statisticMethod
			//Note: If StationID is passed as -99999, it indicates Flag_MultipleStations != 1, not multiple station
			DegreeHour_CreateRatios_From_Sums_by_YearMonthHour(input, DegreeHour_statisticMethod, -99999);
			//WriteDegreeHourToCSV(bool isMultipleStation) for single station
			WriteDegreeHourToCSV(input->SimulationNumericalParams["Flag_MultipleStations"] == 1);
		}
		else {
			//DegreeHour_CreateRatios_From_Sums_by_YearMonthHour function computes the degree hours using DegreeHour_statisticMethod
			DegreeHour_CreateRatios_From_Sums_by_YearMonthHour(input, DegreeHour_statisticMethod, stoi(ReadStationID));
		}
	}
	

	//If SimulationDate_YYYYMMDD.size() <= 0 and input->SimulationNumericalParams["Flag_MultipleStations"] != 1 then error with inputs or expectations about simulation dates
	if (input->SimulationDate_YYYYMMDD.size() <= 0 && input->SimulationNumericalParams["Flag_MultipleStations"] != 1) {
		cout << "Warning: The HydroPlusConfig.xml StartDate_YYYYMMDD and StopDate_YYYYMMDD are not within Weather.csv dates." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: The HydroPlus model needs the Weather.csv to contain the dates given in HydroPlusConfig.xml." << endl;
		cout << "Correction: Change the dates in Weather.csv or in HydroPlusConfig.xml StartDate_YYYYMMDD and StopDate_YYYYMMDD." << endl;
		cout << "Correction: Change the dates in Weather.csv or in HydroPlusConfig.xml StartDate_YYYYMMDD and StopDate_YYYYMMDD." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	readWeather.close();

	//SimulationTimePeriod_seconds (sec) is total simulated time
	//Note: Uses function accumulate(input->SimulationTimeStep_Duration_sec.begin(), input->SimulationTimeStep_Duration_sec.end(),0)
	input->SimulationTimePeriod_seconds = accumulate(input->SimulationTimeStep_Duration_sec.begin(), input->SimulationTimeStep_Duration_sec.end(),0);
	//SimulationTimePeriod_timeSteps (counter of time steps) is SimulationTimeStep_Duration_sec.size()
	input->SimulationTimePeriod_timeSteps = input->SimulationTimeStep_Duration_sec.size();

	int station_Adjusted_Timesteps;
	//If Flag_HottestDaySimulated equals 1, search for the hottest date
	//Note: Flag_HottestDaySimulated = 1 searches all weather files for hottest hour, allowing differences between stations across years, months, and days
	//Note: Model output file z_StationInfoMap.csv contains variable SimulationDate_Hottest_GD that reports the hottest date used for each station
	if (input -> SimulationNumericalParams["Flag_HottestDaySimulated"] == 1 ) {
		//Tair_Max_index is sorted as the maximum value within the Tair_C_AllRecord, using sort_indexes function
		int Tair_Max_index = int(sort_indexes(Tair_C_AllRecord).back());
		//Tair_Max_C assigned Tair_C_AllRecord[Tair_Max_index], the hottest air temp value
		long double Tair_Max_C = Tair_C_AllRecord[Tair_Max_index];

		//SimulationDate_Hottest_GD (GD) is date of hottest hour
		SimulationDate_Hottest_GD = SimulationDate_GD_AllRecord[Tair_Max_index];
		//SimulationTime_Hottest_HMS (HMS) is time of hottest hour (real hottest hour and time)
		SimulationTime_Hottest_HMS = SimulationTime_HMS_AllRecord[Tair_Max_index];
		//Actual_hottest_GDH_index = Tair_Max_index
		Actual_hottest_GDH_index = Tair_Max_index;

		//total_timestep is input->SimulationTimePeriod_timeSteps
		int total_timestep = input->SimulationTimePeriod_timeSteps;
		//adjusted_timesteps initialized to total_timestep; will be reduced if index is too early
		station_Adjusted_Timesteps = total_timestep;
		
		//If the hottest hour occurs too early in the weather file (before enough previous timesteps exist)
		if (Actual_hottest_GDH_index < total_timestep - 1) {
			//Reduce number of timesteps to avoid invalid index access. For example the total_timestep = 5, the hottest hour happend in timestep 3, the station_Adjusted_Timesteps will be 3+1 = 4, shorten by 1
			//...if Actual_hottest_GDH_index = 0, station_Adjusted_Timesteps will be 1, only one timestep will be simulated, and trigger the warning message
			station_Adjusted_Timesteps = Actual_hottest_GDH_index + 1;

			cout << "Warning: Station " << ReadStationID << " hottest hour occurs too early in weather file(Index = " << Actual_hottest_GDH_index << "), unable to extract " << total_timestep << " timesteps.\n";
			cout << "Adjusting simulation to use " << station_Adjusted_Timesteps << " timesteps instead.\n";
			
			//update Global_Adjusted_Timesteps
			Global_Adjusted_Timesteps = min(Global_Adjusted_Timesteps, station_Adjusted_Timesteps);
		}
		//Loop through each timestep leading up to and including the hottest hour
		for (int timeStep = 0; timeStep < station_Adjusted_Timesteps; ++timeStep) {
			//Index within the full weather record for this timestep
			int current_index = Actual_hottest_GDH_index - (station_Adjusted_Timesteps - 1 - timeStep);
			//SimulationDate_YYYYMMDD.push_back(SimulationDate_Hottest_GD) appended for tracking simulation
			input->SimulationDate_YYYYMMDD.push_back(SimulationDate_Hottest_GD);
			//SimulationDate_GDH.push_back(SimulationDate_Hottest_GD * 100 + SimulationDate_HH[timeStep])
			input->SimulationDate_GDH.push_back(SimulationDate_Hottest_GD * 100 + SimulationDate_HH[timeStep]);

			//input->Tair_C.push_back(Tair_C_AllRecord[current_index]); grabbing hottest hour 
			//Note: Within Flag_HottestDaySimulated is true, and total_timestep should be limited to StopDate_YYYYMMDD - StartDate_YYYYMMDD
			input->Tair_C.push_back(Tair_C_AllRecord[current_index]);
			input->Tdew_C.push_back(Tdew_C_AllRecord[current_index]);
			input->RadiationNet_Wpm2.push_back(RadiationNet_Wpm2_AllRecord[current_index]);
			input->WindSpd_mps.push_back(WindSpd_mps_AllRecord[current_index]);
			//If Model_Selection is CoolBuilding then add WindDir_deg to list
			if (input->SimulationStringParams["Model_Selection"] == "CoolBuilding") {
				input->WindDir_deg.push_back(WindDir_deg_AllRecord[current_index]);
			}
			input->AtmPres_kPa.push_back(AtmPres_kPa_AllRecord[current_index]);
			input->Precip_mpdt.push_back(Precip_mpdt_AllRecord[current_index]);
			input->Rain_m.push_back(Rain_m_AllRecord[current_index]);
			input->Snow_m.push_back(Snow_m_AllRecord[current_index]);
			//total items need to accumulate
			input->TotalPrecip_m += Precip_mpdt_AllRecord[current_index];
			input->TotalRain_m += Rain_m_AllRecord[current_index];
			input->TotalSnow_m += Snow_m_AllRecord[current_index];
		}
	}

	// Create map keys and weather map only when Flag_MultipleStations is set to 1
	if (input->SimulationNumericalParams["Flag_MultipleStations"] == 1) {
		//MeteorologicalVariableName_string string contains variable names listed in header column of Weather.csv 
		vector<string> MeteorologicalVariableName_string = { "RadiationNet_Wpm2", "Tair_C", "Tdew_C", "WindSpd_mps", "AtmPres_kPa", "Precip_mpdt", "Rain_m", "Snow_m","TotalPrecip_m", "TotalRain_m", "TotalSnow_m"};

		//If Model_Selection is CoolBuilding then add WindDir_deg to list
		if (input->SimulationStringParams["Model_Selection"] == "CoolBuilding") {
			//MeteorologicalVariableName_string push_back WindDir_deg
			MeteorologicalVariableName_string.push_back("WindDir_deg");
		}
		//StationID_string string contains OBJECTID_shortest[i] item, station ID of nearest weather station
		string StationID_string = ReadStationID;
		string StationID_VariableName;
	
		//For Loop with string name through MeteorologicalVariableName_Weather_string vector list
		for (string Station_Variable_name : MeteorologicalVariableName_string) {
			//StationID_VariableName is string appended by StationID_string and Station_Variable_name, e.g., 05_RadiationNet_Wpm2
			StationID_VariableName = StationID_string + "_" + Station_Variable_name;
		
			if (Station_Variable_name == "RadiationNet_Wpm2") {
				WeatherMap[StationID_VariableName] = input->RadiationNet_Wpm2;
			}
			else if (Station_Variable_name == "Tair_C") {
				WeatherMap[StationID_VariableName] = input->Tair_C;
			}
			else if (Station_Variable_name == "Tdew_C") {
				WeatherMap[StationID_VariableName] = input->Tdew_C;
			}
			else if (Station_Variable_name == "WindSpd_mps") {
				WeatherMap[StationID_VariableName] = input->WindSpd_mps;
			}
			//If Model_Selection is CoolBuilding
			else if (Station_Variable_name == "WindDir_deg") {
				WeatherMap[StationID_VariableName] = input->WindDir_deg;
			}
			else if (Station_Variable_name == "AtmPres_kPa") {
				WeatherMap[StationID_VariableName] = input->AtmPres_kPa;
			}
			else if (Station_Variable_name == "Precip_mpdt") {
				WeatherMap[StationID_VariableName] = input->Precip_mpdt;
			}
			else if (Station_Variable_name == "Rain_m") {
				WeatherMap[StationID_VariableName] = input->Rain_m;
			}
			else if (Station_Variable_name == "Snow_m") {
				WeatherMap[StationID_VariableName] = input->Snow_m;
			}
			else if (Station_Variable_name == "TotalPrecip_m") {
				WeatherMap[StationID_VariableName] = { input->TotalPrecip_m };
			}
			else if (Station_Variable_name == "TotalRain_m") {
				WeatherMap[StationID_VariableName] = { input->TotalRain_m };
			}
			else if (Station_Variable_name == "TotalSnow_m") {
				WeatherMap[StationID_VariableName] = { input->TotalSnow_m };
			}
		}

		//If Flag_HottestDaySimulated is 1 then convert SimulationDate_GDH_AllRecord to SimulationDateStart_GDH and SimulationDateStop_GDH
		//Note this step of adjusting SimulationDateStart_GDH and SimulationDateStop_GDH is necessary, because it controls the TS read in for the radiation reading and evaporation reading
		if (input->SimulationNumericalParams["Flag_HottestDaySimulated"] == 1) {
			//SimulationDateStart_GDH = SimulationDate_GDH_AllRecord[Actual_hottest_GDH_index - (station_Adjusted_Timesteps - 1)]
			SimulationDateStart_GDH = SimulationDate_GDH_AllRecord[Actual_hottest_GDH_index - (station_Adjusted_Timesteps - 1)];
			//SimulationDateStop_GDH = SimulationDate_GDH_AllRecord[Actual_hottest_GDH_index]
			SimulationDateStop_GDH = SimulationDate_GDH_AllRecord[Actual_hottest_GDH_index];
		
			//clean up the AllRecord vectors which used to save the years long full record to find the hottest day and hour
			SimulationDate_GD_AllRecord.clear();
			SimulationDate_GDH_AllRecord.clear();
			SimulationTime_HMS_AllRecord.clear();
			SimulationDate_HH_AllRecord.clear();
			RadiationNet_Wpm2_AllRecord.clear();
			Tair_C_AllRecord.clear();
			Tdew_C_AllRecord.clear();
			Rain_m_AllRecord.clear();
			Snow_m_AllRecord.clear();
			WindSpd_mps_AllRecord.clear();
			WindDir_deg_AllRecord.clear();
			AtmPres_kPa_AllRecord.clear();
			Precip_mpdt_AllRecord.clear();
		}
	}	
}

//DegreeHour_SumHeatingCooling_by_YearMonthHour will determine magnitude of heating and cooling by year, month, and hour
//Note: It is called from within VectorsOfWeatherData, for each time step during a while loop through the Weather.csv file
void WeatherProcessor::DegreeHour_SumHeatingCooling_by_YearMonthHour(Inputs* input, double tempTair_K, int Year_YYYY, int Month_MM, int Hour_HH, string DegreeHour_statisticMethod) {
	//DegreeHour_Delta_abs_K = abs(tempTair_K - DegreeHour_Threshold_K)
	//Note: DegreeHour_Threshold_K often set to equivalent of 65 F for degree day heating and cooling calculations
	double DegreeHour_Delta_abs_K = abs(tempTair_K - DegreeHour_Threshold_K);

	//DegreeHour_Cooling_abs_K = tempTair_K > DegreeHour_Threshold_K ? DegreeHour_Delta_abs_K : 0.0 based on ternary function
	//Note: DegreeHour_Cooling_abs_K should be set to zero if not active
	double DegreeHour_Cooling_abs_K = tempTair_K > DegreeHour_Threshold_K ? DegreeHour_Delta_abs_K : 0.0;
	//DegreeHour_Heating_abs_K = tempTair_K < DegreeHour_Threshold_K ? DegreeHour_Delta_abs_K : 0.0 based on ternary function
	//Note: DegreeHour_Heating_abs_K should be set to zero if not active
	double DegreeHour_Heating_abs_K = tempTair_K < DegreeHour_Threshold_K ? DegreeHour_Delta_abs_K : 0.0;

	//If Inputs::Convert_String_to_Lower_Case(DegreeHour_statisticMethod) == "median"
	//Note: Median is the default method, and could likely remove option for mean below
	if (Inputs::Convert_String_to_Lower_Case(DegreeHour_statisticMethod) == "median") {
		//If (DegreeHour_Cooling_abs_K > 0) then building cooling is expected based on DegreeHour_Threshold_K
		if (DegreeHour_Cooling_abs_K > 0) {
			//DegreeHour_Cooling_Year_abs_K_vec[Year_YYYY].push_back(DegreeHour_Cooling_abs_K); append signal of cooling
			//Note: DegreeHour_Year_Cooling_Ratio_vec is derived from DegreeHour_Cooling_Year_abs_K_vec
			DegreeHour_Cooling_Year_abs_K_vec[Year_YYYY].push_back(DegreeHour_Cooling_abs_K);
			DegreeHour_Cooling_Month_abs_K_vec[Month_MM].push_back(DegreeHour_Cooling_abs_K);
			//DegreeHour_Cooling_Hour_abs_K_vec[Hour_HH].push_back(DegreeHour_Cooling_abs_K); append signal of cooling
			//Note: DegreeHour_Hour_Cooling_K_vec is derived from DegreeHour_Cooling_Hour_abs_K_vec
			DegreeHour_Cooling_Hour_abs_K_vec[Hour_HH].push_back(DegreeHour_Cooling_abs_K);
			DegreeHour_Cooling_All_abs_K_vec.push_back(DegreeHour_Cooling_abs_K);
		}
		//If DegreeHour_Heating_abs_K > 0 then building heating is expected based on DegreeHour_Threshold_K
		if (DegreeHour_Heating_abs_K > 0) {
			//creating the heating vector
			DegreeHour_Heating_Year_abs_K_vec[Year_YYYY].push_back(DegreeHour_Heating_abs_K);
			DegreeHour_Heating_Month_abs_K_vec[Month_MM].push_back(DegreeHour_Heating_abs_K);
			DegreeHour_Heating_Hour_abs_K_vec[Hour_HH].push_back(DegreeHour_Heating_abs_K);
			DegreeHour_Heating_All_abs_K_vec.push_back(DegreeHour_Heating_abs_K);
		}
	}

	//If if (Inputs::Convert_String_to_Lower_Case(DegreeHour_statisticMethod) == "mean")
	//Note: Method not currently used and could likely be deleted
	if (Inputs::Convert_String_to_Lower_Case(DegreeHour_statisticMethod) == "mean") {
		//If DegreeHour_Cooling_abs_K > 0 then building cooling is expected based on DegreeHour_Threshold_K
		if (DegreeHour_Cooling_abs_K > 0) {
			//add the count to the cooling count
			DH_Year_Cooling_count[Year_YYYY]++;
			DH_Month_Cooling_count[Month_MM]++;
			DH_Hour_Cooling_count[Hour_HH]++;
			DH_data_Cooling_total_count++;

			//compute the cooling sum
			DH_Cooling_Year_sum[Year_YYYY] += DegreeHour_Cooling_abs_K;
			DH_Cooling_Month_sum[Month_MM] += DegreeHour_Cooling_abs_K;
			DH_Cooling_Hour_sum[Hour_HH] += DegreeHour_Cooling_abs_K;
			DH_Cooling_total_sum += DegreeHour_Cooling_abs_K;
		}
		//If DegreeHour_Heating_abs_K > 0 then building heating is expected based on DegreeHour_Threshold_K
		if (DegreeHour_Heating_abs_K > 0) {
			//add the count to the heating count
			DH_Year_Heating_count[Year_YYYY]++;
			DH_Month_Heating_count[Month_MM]++;
			DH_Hour_Heating_count[Hour_HH]++;
			DH_data_Heating_total_count++;

			//compute the heating sum
			DH_Heating_Year_sum[Year_YYYY] += DegreeHour_Heating_abs_K;
			DH_Heating_Month_sum[Month_MM] += DegreeHour_Heating_abs_K;
			DH_Heating_Hour_sum[Hour_HH] += DegreeHour_Heating_abs_K;
			DH_Heating_total_sum += DegreeHour_Heating_abs_K;
		}
	}
}

//WriteDegreeHourToCSV function to write cooling/heating DH results to CSV
void WeatherProcessor::WriteDegreeHourToCSV(bool isMultipleStation) {

	ofstream WriteDegreeHour(Inputs::SimulationStringParams["OutputFolder_Path"] + "z_DegreeHour_Ratios.csv");
	if (!WriteDegreeHour.is_open()) {
		cout << "Failed to open " << "DegreeHour_Results.csv" << " for writing.\n";
		return;
	}

	//write out in setprecision(4)
	WriteDegreeHour << fixed << setprecision(4);

	if (isMultipleStation) {
		for (const auto& entry : DegreeHour_Year_Cooling_Ratio_StationID_vec_map) {
			int StationID = entry.first;

			WriteDegreeHour << "\n--- Station ID: " << StationID << endl;

			// Yearly Cooling Ratio
			WriteDegreeHour << "Yearly Cooling Ratio: ";
			for (const auto& val : DegreeHour_Year_Cooling_Ratio_StationID_vec_map[StationID])
				WriteDegreeHour << val << " ";
			WriteDegreeHour << endl;

			// Yearly Heating Ratio
			WriteDegreeHour << "Yearly Heating Ratio: ";
			for (const auto& val : DegreeHour_Year_Heating_Ratio_StationID_vec_map[StationID])
				WriteDegreeHour << val << " ";
			WriteDegreeHour << endl;

			// Monthly Cooling Ratio
			WriteDegreeHour << "Monthly Cooling Ratio: ";
			for (const auto& val : DegreeHour_Month_Cooling_Ratio_StationID_vec_map[StationID])
				WriteDegreeHour << val << " ";
			WriteDegreeHour << endl;

			// Monthly Heating Ratio
			WriteDegreeHour << "Monthly Heating Ratio: ";
			for (const auto& val : DegreeHour_Month_Heating_Ratio_StationID_vec_map[StationID])
				WriteDegreeHour << val << " ";
			WriteDegreeHour << endl;

			// Hourly Cooling (K)
			WriteDegreeHour << "Hourly Cooling (K): ";
			for (const auto& val : DegreeHour_Hour_Cooling_K_StationID_vec_map[StationID])
				WriteDegreeHour << val << " ";
			WriteDegreeHour << endl;

			// Hourly Heating (K)
			WriteDegreeHour << "Hourly Heating (K): ";
			for (const auto& val : DegreeHour_Hour_Heating_K_StationID_vec_map[StationID])
				WriteDegreeHour << val << " ";
			WriteDegreeHour << endl;
			WriteDegreeHour << endl;
		}
	}
	else {
		WriteDegreeHour << "Single Station Results\n";

		// Yearly Cooling Ratio
		WriteDegreeHour << "Yearly Cooling Ratio:\n";
		for (const auto& val : DegreeHour_Year_Cooling_Ratio_vec)
			WriteDegreeHour << val << ",";
		WriteDegreeHour << "\n";

		// Yearly Heating Ratio
		WriteDegreeHour << "Yearly Heating Ratio:\n";
		for (const auto& val : DegreeHour_Year_Heating_Ratio_vec)
			WriteDegreeHour << val << ",";
		WriteDegreeHour << "\n";

		// Monthly Cooling Ratio
		WriteDegreeHour << "Monthly Cooling Ratio:\n";
		for (const auto& val : DegreeHour_Month_Cooling_Ratio_vec)
			WriteDegreeHour << val << ",";
		WriteDegreeHour << "\n";

		// Monthly Heating Ratio
		WriteDegreeHour << "Monthly Heating Ratio:\n";
		for (const auto& val : DegreeHour_Month_Heating_Ratio_vec)
			WriteDegreeHour << val << ",";
		WriteDegreeHour << "\n";

		// Hourly Cooling (K)
		WriteDegreeHour << "Hourly Cooling (K):\n";
		for (const auto& val : DegreeHour_Hour_Cooling_K_vec)
			WriteDegreeHour << val << ",";
		WriteDegreeHour << "\n";

		// Hourly Heating (K)
		WriteDegreeHour << "Hourly Heating (K):\n";
		for (const auto& val : DegreeHour_Hour_Heating_K_vec)
			WriteDegreeHour << val << ",";
		WriteDegreeHour << "\n";
	}
	WriteDegreeHour.close();
}

//VectorsOfWeatherData function will assign Weather.csv lines of data to model vectors for single and multiple stations
//Note: Called from readWeatherFile function, for each time step within while loop 
//Note: Flag_SavingSimulationTimeRange when true will save simulation time range to a vector
//Note: Flag_ProcessTimeStep when true then lineWeatherFile is within the simulation range
//Note: Flag_HottestDaySimulated when true the weather variables are appended to _AllRecord variables for searching hottest hour
void WeatherProcessor::VectorsOfWeatherData(Inputs* input, istream& lineWeatherFile, bool Flag_SavingSimulationTimeRange, bool Flag_ProcessTimeStep, bool Flag_HottestDaySimulated) {

	string temp;
	//initialize date as -1, other variables to -99999
	double tempTair_C = -99999;
	double tempTdew_C = -99999;
	double tempTair_F = -99999;
	double tempTdew_F = -99999;
	double tempNetRadiation_Wpm2 = -99999;
	double tempWindSpd_mps = -99999;
	//tempWindDir_deg for EnergyLoad model
	double tempWindDir_deg = -99999;
	double tempAtmPres_kPa = -99999;
	double tempPrecip_mph = -99999;
	double tempPrecip_m = -99999;
	double tempSnowDepth_m = -99999;

	int date_InputDataRow_HH;
	int date_InputDataRow_GDH;
	int Date_JulianDay;
	int Counter_inputRow = 0;
	//timeStep_InputDataRow_sec initialized to HydroPlusConfig.xml TimeStep_sec
	int timeStep_InputDataRow_sec = input->SimulationNumericalParams["TimeStep_sec"];
	int SimulationTimeStep_Duration_sec;

	//If Flag_ProcessTimeStep is true then Time_HHMMSS_str is within the simulation range
	//Note: If Flag_HottestDaySimulated is true then Flag_ProcessTimeStep is only true when within StopDate_YYYYMMDD - StartDate_YYYYMMDD period	
	if (Flag_ProcessTimeStep == true ) {
		//Multiple Station implmentation; note vectors are cleared prior to each station, so only last station is kept
		//Time_HHMMSS returned from convertHrMinSecToInt(Time_HHMMSS_str) int function
		Time_HHMMSS = convertHrMinSecToInt(Time_HHMMSS_str);
		//timeStep_InputDataRow_sec returned from calculateTimeStep(Date_YYYYMMDD, Time_HHMMSS) function
		timeStep_InputDataRow_sec = calculateTimeStep(Date_YYYYMMDD, Time_HHMMSS, previous_YYYYMMDD, previous_HMS, input->SimulationNumericalParams["TimeStep_sec"]);
		//If timeStep_InputDataRow_sec <= 0 then throw warning
		if (timeStep_InputDataRow_sec <= 0) {
			cout << "Warning: The Weather.csv has problematic date or time chronology at row: " << Counter_inputRow + 1 << endl;
			cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
			cout << "Explanation: The Weather.csv likely has 2 or more sections with the same date: " << Date_YYYYMMDD << endl;
			cout << "Correction: Search and remove from Weather.csv any instances of problematic chronological order." << endl;
			//Call abort function, which ends the HydroPlus.exe simulation
			abort();
		}
		//SimulationTimeStep_Duration_sec vector populated with timeStep_InputDataRow_sec to record duration or time of each timeStep
		input->SimulationTimeStep_Duration_sec.push_back(timeStep_InputDataRow_sec);
		//timeStepMetData advanced when within simulation date range for use with input->SimulationTimeStep_Duration_sec[timeStep] 
		timeStepMetData = timeStepMetData + 1;
		//final area within while loop reading Weather.csv
		//previous_GDH updated to date_InputDataRow_GDH
		previous_YYYYMMDD = Date_YYYYMMDD;
		//previous_HMS updated to Time_HHMMSS
		previous_HMS = Time_HHMMSS;
	}
	//if Flag_SavingSimulationTimeRange or Flag_HottestDaySimulated are true then calculate CalcTreeLAI and CalcShortVegLAI
	//Note: If Flag_HottestDaySimulated is true then regardless of actual hottest day, output is given date from HydroPlusConfig.xml
	if (Flag_SavingSimulationTimeRange || Flag_HottestDaySimulated) {
		//Date_JulianDay = input->Date_YYYYMMDD_to_JulianDay(Date_YYYYMMDD); function obtains Julian dates
		Date_JulianDay = input->Date_YYYYMMDD_to_JulianDay(Date_YYYYMMDD);
		//CalcTreeLAI function obtains LAI for Date_JulianDay
		Inputs::CalcTreeLAI(Date_JulianDay);
		//CalcShortVegLAI function obtains LAI for Date_JulianDay
		Inputs::CalcShortVegLAI(Date_JulianDay);
	}

	//Reading air temperature input (F)
	getline(lineWeatherFile, temp, ',');
	istringstream ssTair(temp);
	//If Flag_Tair_degrees_F is true then input header contains Temperature_Air_F, read Tair_F 
	if (Flag_Tair_degrees_F) {
		ssTair >> tempTair_F;
		//Convert Temperature from Fahrenheit to Celcius and Kelvin
		tempTair_C = (tempTair_F - 32.0) * 5.0 / 9.0;
	}
	//Else default is input is Tair_C
	else {
		ssTair >> tempTair_C;
	}
	//tempTair_K converted from tempTair_C
	double tempTair_K = tempTair_C + 273.15;

	//If Flag_SavingSimulationTimeRange is true then save a single station to input->variable
	if (Flag_SavingSimulationTimeRange) {
		input->Tair_C.push_back(tempTair_C);
	}
	//If Flag_HottestDaySimulated is true, save to variable with _AllRecord to search for hottest hour
	if (Flag_HottestDaySimulated) {
		Tair_C_AllRecord.push_back(tempTair_C);
	}

	//If Flag_DegreeHour_Simulated is true then enter to compute degree Hour_HH (DegreeHour_C) cooling and heating
	//Note; if input->SimulationNumericalParams["Flag_MultipleStations"] == 1, compute DH later after reading all weather files
	if (Flag_DegreeHour_Simulated) {
		//Year_YYYY = Date_YYYYMMDD / 10000 based on division removing last 4 digits of Date_YYYYMMDD
		int Year_YYYY = Date_YYYYMMDD / 10000;
		//Month_MM = Date_YYYYMMDD / 100 based on division removing last 2 digits of Date_YYYYMMDD, mod removing 1st 4 digits
		int Month_MM = (Date_YYYYMMDD / 100) % 100 - 1;//note: 0-based index
		//Hour_HH = Time_HHMMSS / 10000 based on division removing last 4 digits of Time_HHMMSS
		int Hour_HH = Time_HHMMSS / 10000;
		//DegreeHour_SumHeatingCooling_by_YearMonthHour called to determine magnitude of heating and cooling by year, month, and hour
		DegreeHour_SumHeatingCooling_by_YearMonthHour(input, tempTair_K, Year_YYYY, Month_MM, Hour_HH, DegreeHour_statisticMethod);
		//if the Climate_Year_Start_YYYY never been defined, assign Year_YYYY to Climate_Year_Start_YYYY
		if (Climate_Year_Start_YYYY <= 0) {
			Climate_Year_Start_YYYY = Year_YYYY;
		}
	}
	//Note, when computing Degree Hour, the rest of reading is not needed;

	//Reading dew point temperature input (F)
	getline(lineWeatherFile, temp, ',');
	istringstream ssDew(temp);
	//If Flag_Tdew_degrees_F is true then input header contains Temperature_DewPoint_F, read Tdew_F 
	if (Flag_Tdew_degrees_F) {
		ssDew >> tempTdew_F;
		//Convert Temperature from Fahrenheit to Celcius and Kelvin
		tempTdew_C = (tempTdew_F - 32.0) * 5.0 / 9.0;
	}
	//Else default is input is Tdew_C
	else {
		ssDew >> tempTair_C;
	}

	//If Flag_SavingSimulationTimeRange is true then save a single station to input->variable
	if (Flag_SavingSimulationTimeRange) {
		input->Tdew_C.push_back(tempTdew_C);
	}
	//If Flag_HottestDaySimulated is true, save to variable with _AllRecord to search for hottest hour
	if (Flag_HottestDaySimulated) {
		Tdew_C_AllRecord.push_back(tempTdew_C);
	}

	//Reading net radiation input (W/m^2)
	getline(lineWeatherFile, temp, ',');
	istringstream ssRn(temp);
	ssRn >> tempNetRadiation_Wpm2;
	//If Flag_SavingSimulationTimeRange is true then save a single station to input->variable
	if (Flag_SavingSimulationTimeRange) {
		input->RadiationNet_Wpm2.push_back(tempNetRadiation_Wpm2);
	}
	//If Flag_HottestDaySimulated is true, save to variable with _AllRecord to search for hottest hour
	if (Flag_HottestDaySimulated) {
		RadiationNet_Wpm2_AllRecord.push_back(tempNetRadiation_Wpm2);
	}

	//Reading wind speed input (m/s)
	getline(lineWeatherFile, temp, ',');
	istringstream ssWind(temp);
	ssWind >> tempWindSpd_mps;
	//If Flag_SavingSimulationTimeRange is true then save a single station to input->variable
	if (Flag_SavingSimulationTimeRange) {
		input->WindSpd_mps.push_back(tempWindSpd_mps);
	}
	//If Flag_HottestDaySimulated is true, save to variable with _AllRecord to search for hottest hour
	if (Flag_HottestDaySimulated) {
		WindSpd_mps_AllRecord.push_back(tempWindSpd_mps);
	}

	//If Model_Selection is CoolBuilding then read WindDir_deg
	if (input->SimulationStringParams["Model_Selection"] == "CoolBuilding") {
		//Reading wind direction input (deg)
		getline(lineWeatherFile, temp, ',');
		istringstream ssWindDir(temp);
		ssWindDir >> tempWindDir_deg;
		//If Flag_SavingSimulationTimeRange is true then save a single station to input->variable
		if (Flag_SavingSimulationTimeRange) {
			input->WindDir_deg.push_back(tempWindDir_deg);
		}
		//If Flag_HottestDaySimulated is true, save to variable with _AllRecord to search for hottest hour
		if (Flag_HottestDaySimulated) {
			WindDir_deg_AllRecord.push_back(tempWindDir_deg);
		}
	}

	//Reading atmospheric pressure input (kPa)
	getline(lineWeatherFile, temp, ',');
	istringstream ssAtmPres(temp);
	ssAtmPres >> tempAtmPres_kPa;
	//If Flag_SavingSimulationTimeRange is true then save a single station to input->variable
	if (Flag_SavingSimulationTimeRange) {
		input->AtmPres_kPa.push_back(tempAtmPres_kPa);
	}
	//If Flag_HottestDaySimulated is true, save to variable with _AllRecord to search for hottest hour
	if (Flag_HottestDaySimulated) {
		AtmPres_kPa_AllRecord.push_back(tempAtmPres_kPa);
	}

	//Reading precipitation input (m/hr)
	getline(lineWeatherFile, temp, ',');
	istringstream ssPrecip(temp);
	ssPrecip >> tempPrecip_mph;
	//tempPrecip_m (m) = tempPrecip_mph (m/h) * Ratio_Hour_to_Second * timeStep_InputDataRow_sec
	tempPrecip_m = tempPrecip_mph * Ratio_Hour_to_Second * timeStep_InputDataRow_sec;
	//If Flag_SavingSimulationTimeRange is true then save a single station to input->variable
	if (Flag_SavingSimulationTimeRange) {
		//Note: For multiple stations, these vectors are updated later
		input->Precip_mpdt.push_back(tempPrecip_m);
		input->TotalPrecip_m += tempPrecip_m;
	}
	//If Flag_HottestDaySimulated is true, save to variable with _AllRecord to search for hottest hour
	if (Flag_HottestDaySimulated) {
		Precip_mpdt_AllRecord.push_back(tempPrecip_m);
	}

	//Threshold_Precipitation_m (m) catches faulty data read for precipitation, triggered perhaps by wrong ordered columns 
	//Note: Reading CoolBuilding Weather.csv can cause an error, mistaking atmospheric pressure for precipition
	int Threshold_Precipitation_m = 50;
	//If tempPrecip_m > Threshold_Precipitation_m, then abort program
	if (tempPrecip_m > Threshold_Precipitation_m) {
		cout << "Warning: The Weather.csv seems to ahve Precipitation_Rate_mph values > " << Threshold_Precipitation_m << "." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: The HydroPlus does not simulate the hydrology of such large precipitation rates." << endl;
		cout << "Explanation: This error may arise if Weather.csv erroneously organizes data columns." << endl;
		cout << "Explanation: The 8th data column of the Weather.csv file should contain precipitation data." << endl;
		cout << "Explanation: Weather.csv data columns for i-Tree have the following sequence:" << endl << endl;
		cout << "YYYYMMDD,HH:MM:SS,Temperature_Air_C(C),Temperature_DewPoint_C(C),Radiation_Net_Wpm2(W/m^2),Wind_Speed_mps(m/s),Pressure_Atmosphere_kPa(kPa),Precipitation_Rate_mph(m/h),Snow_DepthOnGround_m(m)" << endl << endl;
		cout << "Correction: Generate new Weather.csv with WeatherPrepConfig.xml file, setting parameter model to Hydro." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	//If tempTair_C > 0 C then store tempPrecip_m to Rain_m vector, and increment TotalRain_m
	if (tempTair_C > 0.0) {
		//If Flag_SavingSimulationTimeRange is true then save a single station to input->variable
		if (Flag_SavingSimulationTimeRange) {
			input->Rain_m.push_back(tempPrecip_m);
			input->Snow_m.push_back(0);
			input->TotalRain_m += tempPrecip_m;
		}
		//Else save as multiple station to variable with _AllRecord 
		//If Flag_HottestDaySimulated is true, save to variable with _AllRecord to search for hottest hour
		if (Flag_HottestDaySimulated) {
			Rain_m_AllRecord.push_back(tempPrecip_m);
			Snow_m_AllRecord.push_back(0);
		}
	}
	//Else If tempTair_C <= 0 C then store tempPrecip_m to Snow_m vector, and increment TotalSnow_m
	else {
		//If Flag_SavingSimulationTimeRange is true then save a single station to input->variable
		if (Flag_SavingSimulationTimeRange) {
			input->Rain_m.push_back(0);
			input->Snow_m.push_back(tempPrecip_m);
			input->TotalSnow_m += tempPrecip_m;
		}
		//Else save as multiple station to variable with _AllRecord 
		//If Flag_HottestDaySimulated is true, save to variable with _AllRecord to search for hottest hour
		if (Flag_HottestDaySimulated) {
			Rain_m_AllRecord.push_back(0);
			Snow_m_AllRecord.push_back(tempPrecip_m);
		}
	}

	//Reading Snow_DepthOnGround_m or snow depth on ground (m) but not using variable in model
	//Note: Linebreak \n is read to signify end of row
	//Note: if Tair_C <= 0 then tempPrecip_m converted to Snow_m as liquid equivalent
	getline(lineWeatherFile, temp, '\n');
	istringstream ssSnowDepth(temp);
	ssSnowDepth >> tempSnowDepth_m;
}

//Variables defined in WeatherProcessor.h need to be declared in WeatherProcessor.cpp; relates to static vs non-static function
map<string, vector<double>> WeatherProcessor::EvaporationMap;
//readEvaporationFile function will read the input files containing meteorological data
//Note: Header 7 variables of Evaporation.csv: YYYYMMDD,HH:MM:SS,PotentialEvapotranspiration_TreeCover_mph(m/h),PotentialEvaporation_WaterOnGround_mph(m/h),PotentialEvaporation_WaterOnTree_mph(m/h),PotentialSublimation_SnowOnGround_mph(m/h),PotentialSublimation_SnowOnTree_mph(m/h)
//Note: readEvaporationFile is called only when Flag_ReadEvaporationData = 1 and Model = StatisticalHydro
void WeatherProcessor::readEvaporationFile(Inputs* input)
{
	//readEvaporation created as part of fstream, defined as the ReadStation_path and file name for Evaporation.csv
	fstream readEvaporation(ReadStation_path + "Evaporation.csv");

	//If meteorlogical input data file does not exist, then alert user and abort
	if (!readEvaporation) {
		cout << "Warning: Evaporation.csv could not be opened." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: When the Evaporation.csv file is missing there is no evapotranspiration or sublimation data to drive the simulation." << endl;
		cout << "Correction: Create an Evaporation.csv file, either manually or perhaps using the i-Tree WeatherPrep.exe utility." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	string tempStr = "";
	double tempDbl = -1.0;
	string temp;
	double Date_YYYYMMDD = -1.0;
	double tempPotentialEvapotranspiration_TreeCover_m = -1.0;
	double tempPotentialEvaporation_WaterOnGround_m = -1.0;
	double tempPotentialEvaporation_WaterOnTree_m = -1.0;
	double tempPotentialSublimation_SnowOnGround_m = -1.0;
	double tempPotentialSublimation_SnowOnTree_m = -1.0;
	double date_InputDataRow_HH;
	double date_InputDataRow_GDH;
	string HrMinSec;
	int Counter_inputRow = 0;
	int timeStep = 0;

	getline(readEvaporation, tempStr);

	//numColumnsInFile within tempStr using .begin and .end, adding 1 to move from base 0 to 1 for count
	int numColumnsInFile = count(tempStr.begin(), tempStr.end(), ',') + 1;
	//Weather_Columns_count = 7 for all Model_Selection types
	int Evaporation_Columns_count = 7;
	string Evaporation_Column_names = "YYYYMMDD,HH:MM:SS,PotentialEvapotranspiration_TreeCover_mph(m/h),PotentialEvaporation_WaterOnGround_mph(m/h),PotentialEvaporation_WaterOnTree_mph(m/h),PotentialSublimation_SnowOnGround_mph(m/h),PotentialSublimation_SnowOnTree_mph(m/h)";
	//If numColumnsInFile is not equal to Evaporation_Columns_count then warn and abort
	if (numColumnsInFile != Evaporation_Columns_count) {
		cout << endl;
		cout << "Error: Evaporation.csv has an incorrect number of data columns." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: The HydroPlus model needs" << Evaporation_Columns_count << "columns of data to properly assign Evaporation.csv variables." << endl;
		cout << "Correction: Update the Evaporation.csv file to contain the following comma separated variables:" << endl;
		cout << Evaporation_Column_names << endl;
		//abort simulation due to incorrect CSV format
		abort();
	}

	while (readEvaporation.good())
	{
		getline(readEvaporation, temp, ',');

		//temp.erase will trim leading and then trailing whitespace (optional but recommended)
		temp.erase(0, temp.find_first_not_of(" \t\r\n"));
		temp.erase(temp.find_last_not_of(" \t\r\n") + 1);

		//if (temp.empty()) then continue, skip the line if it's empty or contains only whitespace
		if (temp.empty()) {
			continue;
		}

		//Counter_inputRow advances
		Counter_inputRow = Counter_inputRow + 1;
		//istringstream stores an instance of the string temp in the new ss string
		istringstream ss(temp);
		//ss string is passed into Date_YYYYMMDD as double
		ss >> Date_YYYYMMDD;

		getline(readEvaporation, temp, ',');
		istringstream ssTime(temp);
		ssTime >> HrMinSec;

		//If Counter_inputRow equals 1 then at start of file then check if Weather.csv first date is not greater than HydroPlusConfig.xml StartDate_YYYYMMDD
		if (Counter_inputRow == 1) {
			//Check if year matches HydroPlusConfig.xml StartDate_YYYYMMDD Year (and below if Weather.csv contains full YYYYMMDD called for in HydroPlusConfig.xml)
			//HydroPlusConfig_StartDay_str_YYYYMMDD is string form of HydroPlusConfig.xml input StartDate_YYYYMMDD
			string HydroPlusConfig_StartDay_str_YYYYMMDD = to_string(input->SimulationNumericalParams["StartDate_YYYYMMDD"]);
			//WeatherInput_StartDay_YYYY is string form of Weather.csv input date column
			string WeatherInput_StartDay_str_YYYYMMDD = to_string(Date_YYYYMMDD);
			//WeatherInput_StartDay_YYYYMMDD defined as integer of Date_YYYYMMDD input
			int WeatherInput_StartDay_YYYYMMDD = Date_YYYYMMDD;
			//HydroPlusConfig_StartDay_YYYYMMDD is integer form of substring from HydroPlusConfig_StartDay_str_YYYYMMDD, taking first 8 elements
			int HydroPlusConfig_StartDay_YYYYMMDD = atoi(HydroPlusConfig_StartDay_str_YYYYMMDD.substr(0, 8).c_str());

			//if HydroPlusConfig_StartDay_YYYYMMDD > HydroPlusConfig_StartDay_YYYYMMDD then error with inputs or expectations
			if (WeatherInput_StartDay_YYYYMMDD > HydroPlusConfig_StartDay_YYYYMMDD) {
				cout << "Warning: The Evaporation.csv 1st column for date, YYYYMMDD, has a starting date of " << WeatherInput_StartDay_YYYYMMDD << " ..." << endl;
				cout << "... which is later than the HydroPlusConfig.xml parameter StartDate_YYYYMMDD, set to " << HydroPlusConfig_StartDay_YYYYMMDD << "." << endl;
				cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
				cout << "Explanation: The HydroPlus model needs the Evaporation.csv to contain the date given in HydroPlusConfig.xml." << endl;
				cout << "Correction: Either generate new Evaporation.csv or change the StartDate_YYYYMMDD in the HydroPlusConfig.xml file." << endl;
				//Call abort function, which ends the HydroPlus.exe simulation
				abort();
			}
		}

		//Algorithm to extract hour integer from section of temp string, used to refine location where input rows are read
		//date_InputDataRow_HH integer from string with atoi function, extracting 2 digits of temp w substr function; c_str() points to temp
		date_InputDataRow_HH = atoi(temp.substr(0, 2).c_str());

		//date_InputDataRow_GDH contains Gregorian Date and Hour, extending 8 to 10 digits by multiplication with 100, adding hour
		date_InputDataRow_GDH = (Date_YYYYMMDD * 100) + date_InputDataRow_HH;

		//If (Flag_HottestDaySimulated equals 1 & SimulationDateStop_GDH and SimulationDateStart_GDH are for hottest day) OR ..
		//... If (Flag_HottestDaySimulated equals 0 & input->["SimulationDateStop_GDH"] .. are for HydroPlusConfig input)
		//Note: If Flag_HottestDaySimulated = 1 then Weather.csv for hottest hour, with YYYYMMDD differences between stations
		//Note: Output file z_StationInfoMap.csv contains SimulationDate_Hottest_GD reporting hottest date used for each station
		if (((input->SimulationNumericalParams["Flag_HottestDaySimulated"] == 1) && (date_InputDataRow_GDH <= SimulationDateStop_GDH && date_InputDataRow_GDH >= SimulationDateStart_GDH)) || ((input->SimulationNumericalParams["Flag_HottestDaySimulated"] != 1) && (date_InputDataRow_GDH <= input->SimulationDateStop_YYYYMMDDHH && date_InputDataRow_GDH >= input->SimulationDateStart_YYYYMMDDHH))) {

			getline(readEvaporation, temp, ',');
			istringstream ssPet(temp);
			ssPet >> tempPotentialEvapotranspiration_TreeCover_m;
			//tempPotentialEvapotranspiration_TreeCover_m (m) adjusted by TimeStep_sec and Ratio_Hour_to_Second
			tempPotentialEvapotranspiration_TreeCover_m = (tempPotentialEvapotranspiration_TreeCover_m * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep]);
			input->PotentialEvapotranspiration_TreeCover_m.push_back(tempPotentialEvapotranspiration_TreeCover_m);

			getline(readEvaporation, temp, ',');
			istringstream ssPe(temp);
			ssPe >> tempPotentialEvaporation_WaterOnGround_m;
			//tempPotentialEvaporation_WaterOnGround_m (m) adjusted by TimeStep_sec and Ratio_Hour_to_Second
			tempPotentialEvaporation_WaterOnGround_m = (tempPotentialEvaporation_WaterOnGround_m * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep]);
			input->PotentialEvaporation_WaterOnGround_m.push_back(tempPotentialEvaporation_WaterOnGround_m);

			getline(readEvaporation, temp, ',');
			istringstream ssPetree(temp);
			ssPetree >> tempPotentialEvaporation_WaterOnTree_m;
			//tempPotentialEvaporation_WaterOnTree_m (m) adjusted by TimeStep_sec and Ratio_Hour_to_Second
			tempPotentialEvaporation_WaterOnTree_m = (tempPotentialEvaporation_WaterOnTree_m * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep]);
			input->PotentialEvaporation_WaterOnTree_m.push_back(tempPotentialEvaporation_WaterOnTree_m);

			getline(readEvaporation, temp, ',');
			istringstream ssPesnow(temp);
			ssPesnow >> tempPotentialSublimation_SnowOnGround_m;
			//tempPotentialSublimation_SnowOnGround_m (m) adjusted by TimeStep_sec and Ratio_Hour_to_Second
			tempPotentialSublimation_SnowOnGround_m = (tempPotentialSublimation_SnowOnGround_m * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep]);
			input->PotentialSublimation_SnowOnGround_m.push_back(tempPotentialSublimation_SnowOnGround_m);

			//Linebreak \n is read to signify end of row
			getline(readEvaporation, temp, '\n');
			istringstream ssPotentialSublimation_SnowOnTree_m(temp);
			ssPotentialSublimation_SnowOnTree_m >> tempPotentialSublimation_SnowOnTree_m;
			//tempPotentialSublimation_SnowOnTree_m (m) adjusted by TimeStep_sec and Ratio_Hour_to_Second
			tempPotentialSublimation_SnowOnTree_m = (tempPotentialSublimation_SnowOnTree_m * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep]);
			input->PotentialSublimation_SnowOnTree_m.push_back(tempPotentialSublimation_SnowOnTree_m);
			//timeStep advanced when within simulation date range for use with input->SimulationTimeStep_Duration_sec[timeStep] 
			timeStep = timeStep + 1;
		}
		//If Flag_HottestDaySimulated != 1 or Model_Selection != SpatialTemperatureHydro && date_InputDataRow_GDH is between SimulationDateStop_GDH and SimulationDateStart_GDH then proceed to read input
		else if ((input->SimulationNumericalParams["Flag_HottestDaySimulated"] != 1 || input->SimulationStringParams["Model_Selection"] != "SpatialTemperatureHydro") && (date_InputDataRow_GDH <= input->SimulationDateStop_YYYYMMDDHH &&
			date_InputDataRow_GDH >= input->SimulationDateStart_YYYYMMDDHH)) {
			
			getline(readEvaporation, temp, ',');
			istringstream ssPet(temp);
			ssPet >> tempPotentialEvapotranspiration_TreeCover_m;
			//Multiplying tempPotentialEvapotranspiration_TreeCover_m (m) by TimeStep_sec and Ratio_Hour_to_Second
			tempPotentialEvapotranspiration_TreeCover_m = (tempPotentialEvapotranspiration_TreeCover_m * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep]);
			input->PotentialEvapotranspiration_TreeCover_m.push_back(tempPotentialEvapotranspiration_TreeCover_m);

			getline(readEvaporation, temp, ',');
			istringstream ssPe(temp);
			ssPe >> tempPotentialEvaporation_WaterOnGround_m;
			//Multiplying tempPotentialEvaporation_WaterOnGround_m (m) by TimeStep_sec and Ratio_Hour_to_Second
			tempPotentialEvaporation_WaterOnGround_m = (tempPotentialEvaporation_WaterOnGround_m * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep]);
			input->PotentialEvaporation_WaterOnGround_m.push_back(tempPotentialEvaporation_WaterOnGround_m);

			getline(readEvaporation, temp, ',');
			istringstream ssPetree(temp);
			ssPetree >> tempPotentialEvaporation_WaterOnTree_m;
			//Multiplying tempPotentialEvaporation_WaterOnTree_m (m) by TimeStep_sec and Ratio_Hour_to_Second
			tempPotentialEvaporation_WaterOnTree_m = (tempPotentialEvaporation_WaterOnTree_m * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep]);
			input->PotentialEvaporation_WaterOnTree_m.push_back(tempPotentialEvaporation_WaterOnTree_m);

			getline(readEvaporation, temp, ',');
			istringstream ssPesnow(temp);
			ssPesnow >> tempPotentialSublimation_SnowOnGround_m;
			//Multiplying tempPotentialSublimation_SnowOnGround_m (m) by TimeStep_sec and Ratio_Hour_to_Second
			tempPotentialSublimation_SnowOnGround_m = (tempPotentialSublimation_SnowOnGround_m * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep]);
			input->PotentialSublimation_SnowOnGround_m.push_back(tempPotentialSublimation_SnowOnGround_m);

			//Linebreak \n is read to signify end of row
			getline(readEvaporation, temp, '\n');
			istringstream ssPetreesnow(temp);
			ssPetreesnow >> tempPotentialSublimation_SnowOnTree_m;
			//Multiplying tempPotentialSublimation_SnowOnTree_m (m) by TimeStep_sec and Ratio_Hour_to_Second
			tempPotentialSublimation_SnowOnTree_m = (tempPotentialSublimation_SnowOnTree_m * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep]);
			input->PotentialSublimation_SnowOnTree_m.push_back(tempPotentialSublimation_SnowOnTree_m);
			//timeStep advanced when within simulation date range for use with input->SimulationTimeStep_Duration_sec[timeStep] 
			timeStep = timeStep + 1;

		}
		//Else If date_InputDataRow_GDH is not between SimulationDateStop_GDH and SimulationDateStart_GDH then advance to next date
		else {
			getline(readEvaporation, temp);
		}
	}	

	//If PotentialEvapotranspiration_TreeCover_m.size() <= 0 and input->SimulationNumericalParams["Flag_MultipleStations"] != 1 then error with inputs or expectations about simulation dates
	if (input->PotentialEvapotranspiration_TreeCover_m.size() <= 0 && input->SimulationNumericalParams["Flag_MultipleStations"] != 1) {
		cout << "Warning: The HydroPlusConfig.xml StartDate_YYYYMMDD and StopDate_YYYYMMDD are not within Evaporation.csv dates." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: The HydroPlus model needs the Evaporation.csv to contain the dates given in HydroPlusConfig.xml." << endl;
		cout << "Correction: Change the dates in Evaporation.csv or in HydroPlusConfig.xml StartDate_YYYYMMDD and StopDate_YYYYMMDD." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	readEvaporation.close();

	// Create map keys and evaporation map only when Flag_MultipleStations is set to 1
	if (input->SimulationNumericalParams["Flag_MultipleStations"] == 1) {
		//saving all data to map
		vector<string> MeteorologicalVariableName_string = { "PotentialEvapotranspiration_TreeCover_m", "PotentialEvaporation_WaterOnGround_m", "PotentialEvaporation_WaterOnTree_m", "PotentialSublimation_SnowOnGround_m", "PotentialSublimation_SnowOnTree_m" };
		string 	StationID_string = ReadStationID;
		string StationID_VariableName;

		for (string name : MeteorologicalVariableName_string) {
			//key is like 05_RadiationNet_Wpm2
			StationID_VariableName = StationID_string + "_" + name;

			if (name == "PotentialEvapotranspiration_TreeCover_m") {
				EvaporationMap[StationID_VariableName] = input->PotentialEvapotranspiration_TreeCover_m;
			}
			else if (name == "PotentialEvaporation_WaterOnGround_m") {
				EvaporationMap[StationID_VariableName] = input->PotentialEvaporation_WaterOnGround_m;
			}
			else if (name == "PotentialEvaporation_WaterOnTree_m") {
				EvaporationMap[StationID_VariableName] = input->PotentialEvaporation_WaterOnTree_m;
			}
			else if (name == "PotentialSublimation_SnowOnGround_m") {
				EvaporationMap[StationID_VariableName] = input->PotentialSublimation_SnowOnGround_m;
			}
			else if (name == "PotentialSublimation_SnowOnTree_m") {
				EvaporationMap[StationID_VariableName] = input->PotentialSublimation_SnowOnTree_m;
			}
		}
	}
}

//Variables defined in WeatherProcessor.h need to be declared in WeatherProcessor.cpp; relates to static vs non-static function
map<string, vector<double>> WeatherProcessor::RadiationMap;
map<string, double> WeatherProcessor::LatitudeMap;

//readRadiationFile function reads in time series data for HydroPlus models
//Note: Header of Radiation.csv: YYYYMMDD,HH:MM:SS,Radiation_Shortwave_Direct_Wpm2(W/m^2),Radiation_Shortwave_Diffuse_Wpm2(W/m^2),Radiation_Longwave_Downwelling_Wpm2(W/m^2),Radiation_Longwave_Upwelling_Wpm2(W/m^2),Latitude_rad(Radian),HourAngle_rad(Radian),DeclinationAngle_rad(Radian)
//Note: Direct and diffuse solar radiation is typically computed in i-Tree WeatherPrep using NREL National Solar Radiation Database methods (1995)
void WeatherProcessor::readRadiationFile(Inputs* input) {

	//Initialize variables
	int Date_YYYYMMDD = 0;
	double tempDirSW = 0;
	double tempDiffSW = 0;
	double tempDownLW = 0;
	double tempUpLW = 0;
	double tempLatitude = 0;
	double ntemphourAngle = 0.0;
	double ntempDecline = 0.0;
	double date_InputDataRow_HH;
	double date_InputDataRow_GDH;
	int Counter_inputRow = 0;

	string HrMinSec;
	string temphourAngle;
	string tempDecline;
	string tempStr = "";
	string temp;

	//readRadiation created as part of fstream, defined as the ReadStation_path and file name for Radiation.csv
	ifstream readRadiation(ReadStation_path + "Radiation.csv");

	//If meteorlogical input data file does not exist, then alert user and abort
	if (!readRadiation) {
		cout << "Warning: Radiation.csv could not be opened." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: When the Radiation.csv file is missing there is no short or longwave radiation data to drive the simulation." << endl;
		cout << "Correction: Create a Radiation.csv file, either manually or perhaps using the i-Tree WeatherPrep.exe utility." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	getline(readRadiation, tempStr);

	//numColumnsInFile within tempStr using .begin and .end, adding 1 to move from base 0 to 1 for count
	int numColumnsInFile = count(tempStr.begin(), tempStr.end(), ',') + 1;
	//Radiation_Columns_count = 9 for all Model_Selection types
	int Radiation_Columns_count = 9;
	string Radiation_Column_names = "YYYYMMDD,HH:MM:SS,Radiation_Shortwave_Direct_Wpm2(W/m^2),Radiation_Shortwave_Diffuse_Wpm2(W/m^2),Radiation_Longwave_Downwelling_Wpm2(W/m^2),Radiation_Longwave_Upwelling_Wpm2(W/m^2),Latitude_rad(Radian),HourAngle_rad(Radian),DeclinationAngle_rad(Radian)";
	//If numColumnsInFile is not equal to Evaporation_Columns_count then warn and abort
	if (numColumnsInFile != Radiation_Columns_count) {
		cout << endl;
		cout << "Error: Radiation.csv has an incorrect number of data columns." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: The HydroPlus model needs" << Radiation_Columns_count << "columns of data to properly assign Radiation.csv variables." << endl;
		cout << "Correction: Update the Radiation.csv file to contain the following comma separated variables:" << endl;
		cout << Radiation_Column_names << endl;
		//abort simulation due to incorrect CSV format
		abort();
	}

	while (readRadiation.good()) {
		getline(readRadiation, temp, ',');

		//temp.erase will trim leading and then trailing whitespace (optional but recommended)
		temp.erase(0, temp.find_first_not_of(" \t\r\n"));
		temp.erase(temp.find_last_not_of(" \t\r\n") + 1);

		//if (temp.empty()) then continue, skip the line if it's empty or contains only whitespace
		if (temp.empty()) {
			continue;
		}

		//Counter_inputRow advances
		Counter_inputRow = Counter_inputRow + 1;
		//istringstream stores an instance of the string temp in the new ss string
		istringstream ss(temp);
		//ss string is passed into Date_YYYYMMDD as double
		ss >> Date_YYYYMMDD;

		getline(readRadiation, temp, ',');
		istringstream time(temp);
		time >> HrMinSec;

		//If Counter_inputRow equals 1 then at start of file then check if Weather.csv first date is not greater than HydroPlusConfig.xml StartDate_YYYYMMDD
		if (Counter_inputRow == 1) {
			//Check if year matches HydroPlusConfig.xml StartDate_YYYYMMDD Year (and below if Weather.csv contains full YYYYMMDD called for in HydroPlusConfig.xml)
			//HydroPlusConfig_StartDay_str_YYYYMMDD is string form of HydroPlusConfig.xml input StartDate_YYYYMMDD
			string HydroPlusConfig_StartDay_str_YYYYMMDD = to_string(input->SimulationNumericalParams["StartDate_YYYYMMDD"]);
			//WeatherInput_StartDay_YYYY is string form of Weather.csv input date column
			string WeatherInput_StartDay_str_YYYYMMDD = to_string(Date_YYYYMMDD);
			//WeatherInput_StartDay_YYYYMMDD defined as integer of Date_YYYYMMDD input
			int WeatherInput_StartDay_YYYYMMDD = Date_YYYYMMDD;
			//HydroPlusConfig_StartDay_YYYYMMDD is integer form of substring from HydroPlusConfig_StartDay_str_YYYYMMDD, taking first 8 elements
			int HydroPlusConfig_StartDay_YYYYMMDD = atoi(HydroPlusConfig_StartDay_str_YYYYMMDD.substr(0, 8).c_str());

			//if HydroPlusConfig_StartDay_YYYYMMDD > HydroPlusConfig_StartDay_YYYYMMDD then error with inputs or expectations
			if (WeatherInput_StartDay_YYYYMMDD > HydroPlusConfig_StartDay_YYYYMMDD) {
				cout << "Warning: The Radiation.csv 1st column for date, YYYYMMDD, has a starting date of " << WeatherInput_StartDay_YYYYMMDD << " ..." << endl;
				cout << "... which is later than the HydroPlusConfig.xml parameter StartDate_YYYYMMDD, set to " << HydroPlusConfig_StartDay_YYYYMMDD << "." << endl;
				cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
				cout << "Explanation: The HydroPlus model needs the Radiation.csv to contain the date given in HydroPlusConfig.xml." << endl;
				cout << "Correction: Either generate new Radiation.csv or change the StartDate_YYYYMMDD in the HydroPlusConfig.xml file." << endl;
				//Call abort function, which ends the HydroPlus.exe simulation
				abort();
			}
		}

		//Algorithm to extract hour integer from section of temp string, used to refine location where input rows are read
		//date_InputDataRow_HH integer from string with atoi function, extracting 2 digits of temp w substr function; c_str() points to temp
		date_InputDataRow_HH = atoi(temp.substr(0, 2).c_str());
		//date_InputDataRow_GDH contains Gregorian Date and Hour, extending 8 to 10 digits by multiplication with 100, adding hour
		date_InputDataRow_GDH = (Date_YYYYMMDD * 100) + date_InputDataRow_HH;


		//If (Flag_HottestDaySimulated equals 1 & SimulationDateStop_GDH and SimulationDateStart_GDH are for hottest day) OR ..
		//... If (Flag_HottestDaySimulated equals 0 & input->["SimulationDateStop_GDH"] .. are for HydroPlusConfig input)
		//Note: If Flag_HottestDaySimulated = 1 then Weather.csv for hottest hour, with YYYYMMDD differences between stations
		//Note: Output file z_StationInfoMap.csv contains SimulationDate_Hottest_GD reporting hottest date used for each station
		if (((input->SimulationNumericalParams["Flag_HottestDaySimulated"] == 1) && (date_InputDataRow_GDH <= SimulationDateStop_GDH && date_InputDataRow_GDH >= SimulationDateStart_GDH)) || ((input->SimulationNumericalParams["Flag_HottestDaySimulated"] != 1) && (date_InputDataRow_GDH <= input->SimulationDateStop_YYYYMMDDHH && date_InputDataRow_GDH >= input->SimulationDateStart_YYYYMMDDHH))) {

			getline(readRadiation, temp, ',');
			istringstream ssDirsw(temp);
			ssDirsw >> tempDirSW;
			input->Radiation_Shortwave_Direct_Wpm2.push_back(tempDirSW);

			getline(readRadiation, temp, ',');
			istringstream ssDiffsw(temp);
			ssDiffsw >> tempDiffSW;
			input->Radiation_Shortwave_Diffuse_Wpm2.push_back(tempDiffSW);

			getline(readRadiation, temp, ',');
			istringstream ssDownlw(temp);
			ssDownlw >> tempDownLW;
			//Note: Radiation is positive downward, w/ Radiation_Longwave_Downwelling_Wpm2 (W/m2) added within net radiation balances
			input->Radiation_Longwave_Downwelling_Wpm2.push_back(tempDownLW);

			getline(readRadiation, temp, ',');
			istringstream ssUplw(temp);
			ssUplw >> tempUpLW;
			//Note: Radiation is positive downward, w/ Radiation_Longwave_Upwelling_Wpm2 (W/m2) subtracted within net radiation balances
			input->Radiation_Longwave_Upwelling_Wpm2.push_back(tempUpLW);

			//Latitude_rad is a constant scalar, needing one value from Radiation.csv
			//Note: This value is backup if HydroPlusConfig.xml does not provide WKID for calculation of Latitude and Longitude
			getline(readRadiation, temp, ',');
			istringstream latitude(temp);
			latitude >> tempLatitude;
			input->Latitude_rad = tempLatitude;

			getline(readRadiation, temp, ',');
			istringstream hourangle(temp);
			hourangle >> ntemphourAngle;
			input->HourAngle_Solar_rad.push_back(ntemphourAngle);

			//Linebreak \n is read to signify end of row
			getline(readRadiation, temp, '\n');
			istringstream tempdecline(temp);
			tempdecline >> ntempDecline;
			input->DeclinationAngle_Solar_rad.push_back(ntempDecline);
		}
		//If Flag_HottestDaySimulated != 1 or Model_Selection != SpatialTemperatureHydro && date_InputDataRow_GDH is between SimulationDateStop_GDH and SimulationDateStart_GDH then proceed to read input
		else if ((input->SimulationNumericalParams["Flag_HottestDaySimulated"] != 1 || input->SimulationStringParams["Model_Selection"] != "SpatialTemperatureHydro") && (date_InputDataRow_GDH <= input->SimulationDateStop_YYYYMMDDHH &&
			date_InputDataRow_GDH >= input->SimulationDateStart_YYYYMMDDHH)) {

			getline(readRadiation, temp, ',');
			istringstream ssDirsw(temp);
			ssDirsw >> tempDirSW;
			input->Radiation_Shortwave_Direct_Wpm2.push_back(tempDirSW);

			getline(readRadiation, temp, ',');
			istringstream ssDiffsw(temp);
			ssDiffsw >> tempDiffSW;
			input->Radiation_Shortwave_Diffuse_Wpm2.push_back(tempDiffSW);

			getline(readRadiation, temp, ',');
			istringstream ssDownlw(temp);
			ssDownlw >> tempDownLW;
			//Note: Radiation is positive downward, w/ Radiation_Longwave_Downwelling_Wpm2 (W/m2) added within net radiation balances
			input->Radiation_Longwave_Downwelling_Wpm2.push_back(tempDownLW);

			getline(readRadiation, temp, ',');
			istringstream ssUplw(temp);
			ssUplw >> tempUpLW;
			//Note: Radiation is positive downward, w/ Radiation_Longwave_Upwelling_Wpm2 (W/m2) subtracted within net radiation balances
			input->Radiation_Longwave_Upwelling_Wpm2.push_back(tempUpLW);

			//Note: Consider refactor to read and store only one value, as it is constant, or remove from file given we use map Easting and Northing
			getline(readRadiation, temp, ',');
			istringstream latitude(temp);
			latitude >> tempLatitude;
			input->Latitude_rad = tempLatitude;

			getline(readRadiation, temp, ',');
			istringstream hourangle(temp);
			hourangle >> ntemphourAngle;
			input->HourAngle_Solar_rad.push_back(ntemphourAngle);

			//Linebreak \n is read to signify end of row
			getline(readRadiation, temp, '\n');
			istringstream tempdecline(temp);
			tempdecline >> ntempDecline;
			input->DeclinationAngle_Solar_rad.push_back(ntempDecline);

		}
		//Else If date_InputDataRow_GDH is not between SimulationDateStop_GDH and SimulationDateStart_GDH then advance to next date
		else {
			getline(readRadiation, temp, '\n');
		}
	}

	//If Radiation_Shortwave_Direct_Wpm2.size() <= 0 and input->SimulationNumericalParams["Flag_MultipleStations"] != 1 then error with inputs or expectations about simulation dates
	if (input->Radiation_Shortwave_Direct_Wpm2.size() <= 0 && input->SimulationNumericalParams["Flag_MultipleStations"] != 1) {
		cout << "Warning: The HydroPlusConfig.xml StartDate_YYYYMMDD and StopDate_YYYYMMDD are not within Radiation.csv dates." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: The HydroPlus model needs the Radiation.csv to contain the dates given in HydroPlusConfig.xml." << endl;
		cout << "Correction: Change the dates in Radiation.csv or in HydroPlusConfig.xml StartDate_YYYYMMDD and StopDate_YYYYMMDD." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	readRadiation.close();

	// Create map keys and radiation map only when Flag_MultipleStations is set to 1
	if (input->SimulationNumericalParams["Flag_MultipleStations"] == 1) {
		//MeteorologicalVariableName_string string vector defined with variable names in header row of Radiation.csv file
		vector<string> MeteorologicalVariableName_string = { "Radiation_Shortwave_Direct_Wpm2", "Radiation_Shortwave_Diffuse_Wpm2", "Radiation_Longwave_Downwelling_Wpm2", "Radiation_Longwave_Upwelling_Wpm2", "Latitude_rad", "HourAngle_Solar_rad", "DeclinationAngle_Solar_rad" };
		//StationID_string string defined with OBJECTID_shortest[i], StationID for nearest station
		string StationID_string = ReadStationID;
		string StationID_VariableName;

		//For Loop through MeteorologicalVariableName_string
		for (string Variable_name : MeteorologicalVariableName_string) {
			//StationID_VariableName is appended StationID_string and Variable_name, e.g., 05_DirectShortWaveRadiation
			StationID_VariableName = StationID_string + "_" + Variable_name;

			//If Variable_name equals Radiation_Shortwave_Direct_Wpm2 then
			if (Variable_name == "Radiation_Shortwave_Direct_Wpm2") {
				//RadiationMap[StationID_VariableName] equals Radiation_Shortwave_Direct_Wpm2
				RadiationMap[StationID_VariableName] = input->Radiation_Shortwave_Direct_Wpm2;
			}
			else if (Variable_name == "Radiation_Shortwave_Diffuse_Wpm2") {
				RadiationMap[StationID_VariableName] = input->Radiation_Shortwave_Diffuse_Wpm2;
			}
			else if (Variable_name == "Radiation_Longwave_Downwelling_Wpm2") {
				RadiationMap[StationID_VariableName] = input->Radiation_Longwave_Downwelling_Wpm2;
			}
			else if (Variable_name == "Radiation_Longwave_Upwelling_Wpm2") {
				RadiationMap[StationID_VariableName] = input->Radiation_Longwave_Upwelling_Wpm2;
			}
			//Else if Variable_name == "Latitude_rad" then use LatitudeMap for scalar value, without timeStep variability
			else if (Variable_name == "Latitude_rad") {
				LatitudeMap[StationID_VariableName] = input->Latitude_rad;
			}
			else if (Variable_name == "HourAngle_Solar_rad") {
				RadiationMap[StationID_VariableName] = input->HourAngle_Solar_rad;
			}
			else if (Variable_name == "DeclinationAngle_Solar_rad") {
				RadiationMap[StationID_VariableName] = input->DeclinationAngle_Solar_rad;
			}
		}
	}
}

//resizeWeatherVectors function called via SimulationCoordinator, after SimulationTimePeriod_timeSteps has been determined
void WeatherProcessor::resizeWeatherVectors(Inputs* input) {

	//Variable vectors are resized based on SimulationTimePeriod_timeSteps, containing value -9999 which is NODATA_code
	//Note: .resize function only needed when vectors are filled without .pushback function
	input->RadiationNet_Wpm2.resize(input->SimulationTimePeriod_timeSteps, -9999);
	input->Tair_C.resize(input->SimulationTimePeriod_timeSteps, -9999);
	input->Tdew_C.resize(input->SimulationTimePeriod_timeSteps, -9999);
	input->WindSpd_mps.resize(input->SimulationTimePeriod_timeSteps, -9999);
	input->WindDir_deg.resize(input->SimulationTimePeriod_timeSteps, -9999);
	input->AtmPres_kPa.resize(input->SimulationTimePeriod_timeSteps, -9999);
	input->Precip_mpdt.resize(input->SimulationTimePeriod_timeSteps, -9999);
	input->Rain_m.resize(input->SimulationTimePeriod_timeSteps, -9999);
	input->Snow_m.resize(input->SimulationTimePeriod_timeSteps, -9999);

	input->Radiation_Shortwave_Direct_Wpm2.resize(input->SimulationTimePeriod_timeSteps, -9999);
	input->Radiation_Shortwave_Diffuse_Wpm2.resize(input->SimulationTimePeriod_timeSteps, -9999);
	input->Radiation_Longwave_Downwelling_Wpm2.resize(input->SimulationTimePeriod_timeSteps, -9999);
	input->Radiation_Longwave_Upwelling_Wpm2.resize(input->SimulationTimePeriod_timeSteps, -9999);
	input->HourAngle_Solar_rad.resize(input->SimulationTimePeriod_timeSteps, -9999);
	input->DeclinationAngle_Solar_rad.resize(input->SimulationTimePeriod_timeSteps, -9999);

	//If Flag_ReadEvaporationData = 1 and SimulationStringParams["Model_Selection"] != SpatialTemperatureHydro then proceed to call readEvaporationFile function
	//Note: Vectors should be zero size in order to compute evaporation data in PotentialEvaporationCalc::Method_PenmanMonteith function
	if (input->SimulationNumericalParams["Flag_ReadEvaporationData"] == 1 && input->SimulationStringParams["Model_Selection"] != "SpatialTemperatureHydro") {
		input->PotentialEvapotranspiration_TreeCover_m.resize(input->SimulationTimePeriod_timeSteps, -9999);
		input->PotentialEvaporation_WaterOnGround_m.resize(input->SimulationTimePeriod_timeSteps, -9999);
		input->PotentialEvaporation_WaterOnTree_m.resize(input->SimulationTimePeriod_timeSteps, -9999);
		input->PotentialSublimation_SnowOnGround_m.resize(input->SimulationTimePeriod_timeSteps, -9999);
		input->PotentialSublimation_SnowOnTree_m.resize(input->SimulationTimePeriod_timeSteps, -9999);
	}
}

//clearWeatherVectors function will clear the vectors between station calls
void WeatherProcessor::clearWeatherVectors(Inputs* input) {
	//Time vectors should likely be cleared before filled by any weather station in multiple station calls
	input->SimulationTime_HMS.clear();
	input->SimulationDate_YYYYMMDD.clear();
	input->SimulationDate_GDH.clear();
	input->SimulationTime_Output_HMS.clear();
	input->SimulationDate_Output_GD.clear();
	input->SimulationDate_Output_GDH.clear();
	SimulationDate_HH.clear();
	input->SimulationTimeStep_Duration_sec.clear();

	//Meterological vectors need to be cleared before filled by any weather station in multiple station calls
	input->RadiationNet_Wpm2.clear();
	input->Tair_C.clear();
	input->Tdew_C.clear();
	input->WindSpd_mps.clear();
	input->WindDir_deg.clear();
	input->AtmPres_kPa.clear();
	input->Precip_mpdt.clear();
	input->Rain_m.clear();
	input->Snow_m.clear();

	input->Radiation_Shortwave_Direct_Wpm2.clear();
	input->Radiation_Shortwave_Diffuse_Wpm2.clear();
	input->Radiation_Longwave_Downwelling_Wpm2.clear();
	input->Radiation_Longwave_Upwelling_Wpm2.clear();
	input->HourAngle_Solar_rad.clear();
	input->DeclinationAngle_Solar_rad.clear();

	input->PotentialEvapotranspiration_TreeCover_m.clear();
	input->PotentialEvaporation_WaterOnGround_m.clear();
	input->PotentialEvaporation_WaterOnTree_m.clear();
	input->PotentialSublimation_SnowOnGround_m.clear();
	input->PotentialSublimation_SnowOnTree_m.clear();
}

//CalcWeightsandRead function is called to compute the distance based weight between pixel and weather station
//Note: Flag_MultipleStations = 1 computes distance between stations and map pixels that use the same projection, e.g., Albers_Conical_Equal_Area, Lambert_Azimuthal_Equal_Area, or UTM
void WeatherProcessor::CalcWeightsandRead(long double PixelCenter_easting_m, long double PixelCenter_northing_m, Inputs* input, int MapPixel_ID) {

	//WeatherStationDistance_sort function is called get the station distance vector for the pixel center (units convert between decimal degree and radians)
	WeatherProcessor::WeatherStationDistance_sort(input, PixelCenter_easting_m, PixelCenter_northing_m);
	//MultipleStations_InverseDistancePower variable defined with SimulationNumericalParams["MultipleStations_InverseDistancePower"]
	long double MultipleStations_InverseDistancePower = input->SimulationNumericalParams["MultipleStations_InverseDistancePower"];
	long double Weight_StationInverseDistance_frac, Distance_InverseSquared;
	//found_zero_distance = false; initialize to false
	bool found_zero_distance = false;
	//index_zero_distance = -1; initialize to -1
	int index_zero_distance = -1;
	//Distance_InverseSquared_vec.clear(); clear vector as a precaution
	Distance_InverseSquared_vec.clear();

	//For Loop through SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]
	for (int i = 0; i < input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]; i++) {
		//If Distance_shortest is 0 than this station is on the pixel, likely a reference station
		if (Distance_shortest[i] == 0.0) {
			//found_zero_distance = true is activated
			found_zero_distance = true;
			index_zero_distance = i;
			//break given you found a station with 100% of the weight
			break;
		}
	}
	//Weight_StationInverseDistance_vec_frac.clear; precaution
	Weight_StationInverseDistance_vec_frac.clear();
	//If found_zero_distance is true then 1 station gets all weight, others 0 weight
	if (found_zero_distance) {
		//For loop through MultipleStations_NumberOfStationsAveraged
		for (int i = 0; i < input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]; i++) {
			//Weight_StationInverseDistance_vec_frac updated with weight of 1 or 0
			//Note: Only the zero-distance station gets weight 1.0, others get 0.0
			Weight_StationInverseDistance_vec_frac.push_back((i == index_zero_distance) ? 1.0 : 0.0);
		}
	}
	//Else no station gets all weight, and use inverse distance squared
	else {
		//Distance_InverseSquared_sum initialized to 0
		long double Distance_InverseSquared_sum = 0.0;
		//For loop through MultipleStations_NumberOfStationsAveraged
		for (int i = 0; i < input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]; i++) {
			//Distance_InverseSquared computed as 1.0 / pow(Distance_shortest[i], MultipleStations_InverseDistancePower); power squared
			long double Distance_InverseSquared = 1.0 / pow(Distance_shortest[i], MultipleStations_InverseDistancePower);
			Distance_InverseSquared_vec.push_back(Distance_InverseSquared);
			Distance_InverseSquared_sum += Distance_InverseSquared;
		}
		//For loop through MultipleStations_NumberOfStationsAveraged
		for (int i = 0; i < input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]; i++) {
			//weight is computed if Distance_InverseSquared_sum > 0.0, ...
			//Note: weight is as Distance_InverseSquared_vec[i] / Distance_InverseSquared_sum
			double Station_weight = (Distance_InverseSquared_sum > 0.0)
				? Distance_InverseSquared_vec[i] / Distance_InverseSquared_sum
				: 0.0;
			Weight_StationInverseDistance_vec_frac.push_back(Station_weight);
		}
	}

	//Distance_InverseSquared_sum is summation of all Distance_InverseSquared values in vector
	long double Distance_InverseSquared_sum = accumulate(Distance_InverseSquared_vec.begin(), Distance_InverseSquared_vec.end(), 0.0);

	//For Loop through SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]
	for (int i = 0; i < input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]; i++) {
		//Weight_StationInverseDistance_frac (fraction) is Distance_InverseSquared_vec / Distance_InverseSquared_sum
		Weight_StationInverseDistance_frac = input->SafeDivide(Distance_InverseSquared_vec[i], Distance_InverseSquared_sum);
		//If Distance_shortest is 0 then assign weight of 1 when distance is 0
		if (Distance_shortest[i] == 0.0) {
			// Special case: distance is zero -> give full weight to this station
			Weight_StationInverseDistance_frac = 1.0;
		}
		//Weight_StationInverseDistance_vec_frac vector receives Weight_StationInverseDistance_frac
		Weight_StationInverseDistance_vec_frac.push_back(Weight_StationInverseDistance_frac);
		
		//If conditional met (StationInfoMap.find(OBJECTID_shortest[i]) == StationInfoMap.end()) then OBJECTID_shortest[i] not in StationInfoMap, and new item is added
		if (StationInfoMap.find(OBJECTID_shortest[i]) == StationInfoMap.end()) {

			//ReadStation_path is set to Path_folder_shortest
			ReadStation_path = Path_folder_shortest[i];
			//ReadStationID is set to OBJECTID_shortest[i]
			ReadStationID = OBJECTID_shortest[i];

			//clearWeatherVectors called between stations to clear vectors of prior data
			clearWeatherVectors(input);

			//Call functions to read meteorological data readWeatherFile, readRadiationFile, and possibly readEvaporationFile
			readWeatherFile(input);
			readRadiationFile(input);

			//Check if Radiation and Weather files have same length
			//Length_RadiationInputs is integer length of 1st variable from Radiation.csv
			int Length_RadiationInputs = input->Radiation_Shortwave_Direct_Wpm2.size();
			//Length_WeatherInputs is integer length of 1st variable from Weather.csv
			int Length_WeatherInputs = input->Tair_C.size();
			//If Length_WeatherInputs not equal to Length_RadiationInputs then alert user and abort
			//Flag_HottestDaySimulated == 1, don't trigger this warning. The weather, radiation and evaporation file may nned to be trimmed later
			if (Length_WeatherInputs != Length_RadiationInputs && input->SimulationNumericalParams["Flag_HottestDaySimulated"] != 1) {
				cout << "Warning: The Weather.csv has " << Length_WeatherInputs << " items, but Radiation.csv has " << Length_RadiationInputs << " items." << endl;
				cout << "This program cannot verify if they are correct datasets properly overlapping in time." << endl;
				cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
				cout << "Explanation: The HydroPlus model needs Weather.csv and Radiation.csv to contain identical rows of data." << endl;
				cout << "Correction: Generate new Weather.csv or Radiation.csv data, or make those datasets have the same row count and dates." << endl;
				//Call abort function, which ends the HydroPlus.exe simulation
				abort();
			}

			//If Flag_ReadEvaporationData = 1 and Model_Selection != SpatialTemperatureHydro then proceed to call readEvaporationFile function
			if (input->SimulationNumericalParams["Flag_ReadEvaporationData"] == 1 && input->SimulationStringParams["Model_Selection"] != "SpatialTemperatureHydro") {
				//readEvaporationFile called
				//Note: Consider refactor to test if Length_EvaporationInputs != Length_RadiationInputs
				readEvaporationFile(input);
			}
			
			//If SimulationNumericalParams["Flag_HottestDaySimulated"] equals 1 then create vector with instructions to simulate hottest day
			if (input->SimulationNumericalParams["Flag_HottestDaySimulated"] == 1) {
				//StationInfoMap[OBJECTID_shortest[i]] is created, contianing data to map once the station is read
				//Note: Vector contains 2 additional variables, SimulationDate_Hottest_GD[i]), SimulationTime_Hottest_HMS[i] 
				StationInfoMap[OBJECTID_shortest[i]] = { USAF_shortest[i],WBAN_shortest[i],to_string(Easting_shortest[i]),to_string(Northing_shortest[i]), to_string(Elev_m_shortest[i]), to_string(Slope_deg_shortest[i]), to_string(Aspect_deg_shortest[i]),TI_option_shortest[i], to_string_fixed(NLCD_Class_shortest[i],0), to_string(TC_per_shortest[i]), to_string(IC_per_shortest[i]), to_string(SVeg_in_Gap_frac_shortest[i]), to_string(Station_AH_Flux_Wpm2_shortest[i]), to_string(SimulationDate_Hottest_GD), SimulationTime_Hottest_HMS};
			}
			//Else SimulationNumericalParams["Flag_HottestDaySimulated"] equals 0 then create vector without instructions to simulate hottest day
			else {
				//StationInfoMap[OBJECTID_shortest[i]] is created, contianing data to map once the station is read
				//Note: Vector does not contain 2 variables, SimulationDate_Hottest_GD[i]), SimulationTime_Hottest_HMS[i] 
				StationInfoMap[OBJECTID_shortest[i]] = { USAF_shortest[i],WBAN_shortest[i],to_string(Easting_shortest[i]),to_string(Northing_shortest[i]), to_string(Elev_m_shortest[i]), to_string(Slope_deg_shortest[i]), to_string(Aspect_deg_shortest[i]),TI_option_shortest[i], to_string_fixed(NLCD_Class_shortest[i],0), to_string(TC_per_shortest[i]), to_string(IC_per_shortest[i]), to_string(SVeg_in_Gap_frac_shortest[i]), to_string(Station_AH_Flux_Wpm2_shortest[i]), "NA(Flag_HottestDaySimulated turned off)", "NA(Flag_HottestDaySimulated turned off)"};
			}
		}
	}	

	//finish reading all maps, do the trim
	//trim only been called if Flag_HottestDaySimulated turned on and available timestep < simulation timestep
	if (input->SimulationNumericalParams["Flag_HottestDaySimulated"] == 1 && Global_Adjusted_Timesteps < input->SimulationTimePeriod_timeSteps) {
		TrimWeatherMapToAdjustedTimesteps(Global_Adjusted_Timesteps, WeatherMap);
		TrimWeatherMapToAdjustedTimesteps(Global_Adjusted_Timesteps, RadiationMap);
		TrimWeatherMapToAdjustedTimesteps(Global_Adjusted_Timesteps, EvaporationMap);
		//trim the input->SimulationDate_YYYYMMDD and input->SimulationDate_GDH.
		input->SimulationDate_YYYYMMDD.erase(input->SimulationDate_YYYYMMDD.begin(), input->SimulationDate_YYYYMMDD.end() - Global_Adjusted_Timesteps);
		input->SimulationDate_GDH.erase(input->SimulationDate_GDH.begin(), input->SimulationDate_GDH.end() - Global_Adjusted_Timesteps);
		input->SimulationTimeStep_Duration_sec.erase(input->SimulationTimeStep_Duration_sec.begin(), input->SimulationTimeStep_Duration_sec.end() - Global_Adjusted_Timesteps);
		input->SimulationTimePeriod_timeSteps = input->SimulationTimeStep_Duration_sec.size();
	}


	//temp_subPixelNearestStationMap vector created as temporary storage for weather station data
	//Note: temp_subPixelNearestStationMap[a][b], where [a] contains 1 of 3 stations, and [b] contains station ID, distance, and weight
	map<int, vector<string> > temp_subPixelNearestStationMap;
	//For Loop through MultipleStations_NumberOfStationsAveraged to populate temp_subPixelNearestStationMap with ...
	for (int i = 0; i < input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]; i++) {
		//temp_subPixelNearestStationMap receives string for OBJECTID_shortest, containing weather station ID
		temp_subPixelNearestStationMap[i].push_back(OBJECTID_shortest[i]);
		//temp_subPixelNearestStationMap receives string for Distance_shortest, containing weather station distance
		temp_subPixelNearestStationMap[i].push_back(to_string(Distance_shortest[i]));
		//temp_subPixelNearestStationMap receives string for Weight_StationInverseDistance_vec_frac, containing weather station weight
		temp_subPixelNearestStationMap[i].push_back(to_string(Weight_StationInverseDistance_vec_frac[i]));
	}
	//PixelNearestStation[MapPixel_ID] receives data from vector temp_subPixelNearestStationMap
	//Note: PixelNearestStation[MapPixel_ID][a][b], where [a] contains 1 of n stations, and [b] contains station ID, distance, and weight
	//Note: PixelNearestStation outer map key is MapPixel_ID, inner map key is reference station ranking, e.g., 1, 2, 3
	//Note: PixelNearestStation vector value holds station ID, distance (km) and weight (fraction)
	//Note: Use of loop for .first and .second to get key-value pairs
	PixelNearestStation[MapPixel_ID] = temp_subPixelNearestStationMap;
	
	/*
	//Script to view reference station properities
	if (MapPixel_ID == 149 || (MapPixel_ID > (Inputs::nCols * Inputs::nRows) - 1)) {
		cout << "PixelNearestStation contents for MapPixel_ID = " << MapPixel_ID << endl;
		for (size_t a = 0; a < PixelNearestStation[MapPixel_ID].size(); ++a) {
			cout << "  NearestStation " << a << ": ";
			for (size_t b = 0; b < PixelNearestStation[MapPixel_ID][a].size(); ++b) {
				cout << "[" << b << "]=" << PixelNearestStation[MapPixel_ID][a][b] << " ";
			}
			cout << endl;
		}
	}
	*/

	//For Loop through MultipleStations_NumberOfStationsAveraged to clear the temporary vector temp_subPixelNearestStationMap
	for (int i = 0; i < input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]; i++) {
		//temp_subPixelNearestStationMap emptied with C++ clear command
		temp_subPixelNearestStationMap[i].clear();
	}

	//Vectors will contain each map pixel, and contain data for stations with the shortest distance and best weights
	//Vectors are cleared at start
	Distance_shortest.clear();
	OBJECTID_shortest.clear();
	USAF_shortest.clear();
	WBAN_shortest.clear();
	Easting_shortest.clear();
	Northing_shortest.clear();
	Elev_m_shortest.clear();
	Slope_deg_shortest.clear();
	Aspect_deg_shortest.clear();
	TI_option_shortest.clear();
	NLCD_Class_shortest.clear();
	TC_per_shortest.clear();
	IC_per_shortest.clear();
	SVeg_in_Gap_frac_shortest.clear();
	Station_AH_Flux_Wpm2_shortest.clear();
	Path_folder_shortest.clear();
	Distance_InverseSquared_vec.clear();
	Weight_StationInverseDistance_vec_frac.clear();
}

//to_string_fixed helper function will control precision of output
string WeatherProcessor::to_string_fixed(double value_sent, int precision_desired) {
	//creates string_managed as ostringstream
	ostringstream string_managed;
	//oss takes a fixed precision the value_sent
	string_managed << fixed << setprecision(precision_desired) << value_sent;
	//return string_managed
	return string_managed.str();
}

//Variables defined in WeatherProcessor.h need to be declared in WeatherProcessor.cpp; relates to static vs non-static function
long double WeatherProcessor::PixelCenter_X_MapLowerLeft_m, WeatherProcessor::PixelCenter_Y_MapLowerLeft_m;
vector<long double> WeatherProcessor::Map_Easting_X_m_vec, WeatherProcessor::Map_Northing_Y_m_vec;

//get the (x,y) for each pixel center, then calculate weight and read weather data
void WeatherProcessor::MapPixelCalcWeightsandRead(Inputs* input, DataOrganizer* organizer) {
	
	//1. To get corner location, left lower pixel
	//PixelCenter_X_MapLowerLeft_m is xllcorner_m + 0.5 * Inputs::Length_Pixel_Side_m
	PixelCenter_X_MapLowerLeft_m = Inputs::xllcorner_m + 0.5 * Inputs::Length_Pixel_Side_m;
	PixelCenter_Y_MapLowerLeft_m = Inputs::yllcorner_m + 0.5 * Inputs::Length_Pixel_Side_m;

	//2. For each pixel, get Latitude_DecimalDegree, Longitude_DecimalDegree location, looping through SimulationNumericalParams["NumberOfLocations"]
	for (int MapPixel_ID = 0; MapPixel_ID < input->SimulationNumericalParams["NumberOfLocations"]; ++MapPixel_ID) {

		//if (Folder_count == 0) then a NODATA_code pixel with empty DataDrawer; continue to next MapPixel_ID
		if (Inputs::LandCover_NLCD_Class[MapPixel_ID] == Inputs::NODATA_code) { continue; }

		//Note: Conversion with MapPixel_ID starting at 0, row & col starting at 1 for map input & output
		//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 row = MapPixel_ID / Inputs::nCols;
		int col = MapPixel_ID % Inputs::nCols;

		long double Map_Northing_Y_m;
		long double Map_Easting_X_m;

		//If MapPixel_ID < Inputs::nCols * Inputs::nRows then it is within the map
		if (MapPixel_ID < Inputs::nCols * Inputs::nRows) {
			//Northing decreases from top to bottom, so we subtract row index from total heigh. Map_Northing_Y_m (m) is PixelCenter_Y_MapLowerLeft_m + (Inputs::nRows - row) * Inputs::Length_Pixel_Side_m
			Map_Northing_Y_m = PixelCenter_Y_MapLowerLeft_m + (Inputs::nRows - row - 1) * Inputs::Length_Pixel_Side_m;
			//Map_Easting_X_m (m) is PixelCenter_X_MapLowerLeft_m + (col - 1) * Inputs::Length_Pixel_Side_m;
			Map_Easting_X_m = PixelCenter_X_MapLowerLeft_m + col * Inputs::Length_Pixel_Side_m;
		}
		//Else it is outside the map and a reference station
		else {
			//RefStation_Index = input->getReferenceStationIndex(MapPixel_ID); obtained from search function
			int RefStation_Index = input->getReferenceStationIndex(MapPixel_ID);
			//Map_Easting_X_m = input->Station_Easting_m_vec[RefStation_Index]; Station_Easting_m_vec holds reference station location
			Map_Easting_X_m = input->Station_Easting_m_vec[RefStation_Index];
			//Map_Northing_Y_m = input->Station_Northing_m_vec[RefStation_Index]; Station_Northing_m_vec holds reference station location
			Map_Northing_Y_m = input->Station_Northing_m_vec[RefStation_Index];
		}

		//CalcWeightsandRead function called with decimal degree data to compute station spatial weights
		//CalcWeightsandRead(Latitude_DecimalDegree, Longitude_DecimalDegree, input, MapPixel_ID);
		CalcWeightsandRead(Map_Easting_X_m, Map_Northing_Y_m, input, MapPixel_ID);
		//Map_Easting_X_m_vec (m) contains vector of Map_Easting_X_m (m)
		Map_Easting_X_m_vec.push_back(Map_Easting_X_m);
		//Map_Northing_Y_m_vec (m) contains vector of Map_Northing_Y_m (m)
		Map_Northing_Y_m_vec.push_back(Map_Northing_Y_m);
	}
}

int WeatherProcessor::centerMapPixel_ID;

//get the (x,y) for map center pixel, then calculate weight and read weather data
void WeatherProcessor::MapCenterCalcWeightsandRead(Inputs* input) {
	//center_East_m (m) equals xllcorner_m + 0.5 * nCols * Inputs::Length_Pixel_Side_m
	long double center_East_m = Inputs::xllcorner_m + 0.5 * Inputs::nCols * Inputs::Length_Pixel_Side_m;
	//center_North_m (m) equals yllcorner_m + 0.5 * nRows * Inputs::Length_Pixel_Side_m
	long double center_North_m = Inputs::yllcorner_m + 0.5 * Inputs::nRows * Inputs::Length_Pixel_Side_m;
	//calculatethe center pixel ID
	centerMapPixel_ID = (Inputs::nRows / 2) * Inputs::nCols + (Inputs::nCols / 2);
	//CalcWeightsandRead function called with decimal degree data to compute station spatial weights
	CalcWeightsandRead(center_East_m, center_North_m, input, centerMapPixel_ID);
	// Store the center coordinates in the vectors
	Map_Easting_X_m_vec.push_back(center_East_m);
	Map_Northing_Y_m_vec.push_back(center_North_m);
}

//WeightedInputMap defined within WeatherProcessor::WeightedInputMap structure, after declaration in WeatherProcessor.h
map<int, map<string, vector<double> > > WeatherProcessor::WeightedInputMap;
//ComputeSpatiallyWeighted_Weather_Inputs function computes spatially weighted values for each MapPixel_ID to input->
//Note: Multiple Station simulation will write spatially weighted values of meteorological data to the input-> pointer
void WeatherProcessor::ComputeSpatiallyWeighted_Weather_Inputs(Inputs* input, int MapPixel_ID, int timeStep) {
	//MeteorologicalVariableName_Weather_string string contains variable names listed in header column of Weather.csv 
	vector<string> MeteorologicalVariableName_Weather_string = { "Tair_C", "Tdew_C", "RadiationNet_Wpm2", "WindSpd_mps", "AtmPres_kPa", "Precip_mpdt","Rain_m", "Snow_m","TotalPrecip_m", "TotalRain_m", "TotalSnow_m" };
	//MeteorologicalVariableName_Evaporation_string string contains variable names listed in header column of Evaporation.csv 
	vector<string> MeteorologicalVariableName_Evaporation_string = { "PotentialEvapotranspiration_TreeCover_m", "PotentialEvaporation_WaterOnGround_m", 
		"PotentialEvaporation_WaterOnTree_m", "PotentialSublimation_SnowOnGround_m", "PotentialSublimation_SnowOnTree_m" };
	//MeteorologicalVariableName_Radiation_string string contains variable names listed in header column of Radiation.csv 
	vector<string> MeteorologicalVariableName_Radiation_string = { "Radiation_Shortwave_Direct_Wpm2", "Radiation_Shortwave_Diffuse_Wpm2", 
		"Radiation_Longwave_Downwelling_Wpm2", "Radiation_Longwave_Upwelling_Wpm2", "Latitude_rad", "HourAngle_Solar_rad", "DeclinationAngle_Solar_rad" };
	
	//If Model_Selection is CoolBuilding then add WindDir_deg to list
	if (input->SimulationStringParams["Model_Selection"] == "CoolBuilding") {
		//MeteorologicalVariableName_Weather_string push_back WindDir_deg
		MeteorologicalVariableName_Weather_string.push_back("WindDir_deg");
	}

	//1. applying weights to weather data and saving to map WeightedInputMap
	//For Loop with string name through MeteorologicalVariableName_Weather_string vector list
	for (string Station_Variable_name : MeteorologicalVariableName_Weather_string) {
		//MeteorologicalVariable_weightedSum initialized to 0
		double MeteorologicalVariable_weightedSum = 0;
		double MeteorologicalVariable_Total_weightedSum = 0;
		//For Loop with integer i through MultipleStations_NumberOfStationsAveraged
		for (int i = 0; i < input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]; i++) {
			//StationID_string defined as PixelNearestStation[MapPixel_ID][i][0], which is station i's StationID
			string StationID_string = PixelNearestStation[MapPixel_ID][i][0];
			//StationID_VariableName defined as combination of keynamepart1_name and Station_Variable_name
			string StationID_VariableName = StationID_string + "_" + Station_Variable_name;

			//STID_TotalPrecip_m just have one value in this vector
			if (Station_Variable_name == "TotalPrecip_m" || Station_Variable_name == "TotalRain_m" || Station_Variable_name == "TotalSnow_m") {
				MeteorologicalVariable_weightedSum = MeteorologicalVariable_weightedSum + WeatherMap[StationID_VariableName][0] * stod(PixelNearestStation[MapPixel_ID][i][2]);
			}
			else {
				//MeteorologicalVariable_weightedSum is added to MeteorologicalVariable_weightedSum and product of WeatherMap[StationID_VariableName] Value at timeStep and Station Weight, contained in PixelNearestStation[MapPixel_ID][i][2]
				MeteorologicalVariable_weightedSum = MeteorologicalVariable_weightedSum + WeatherMap[StationID_VariableName][timeStep] * stod(PixelNearestStation[MapPixel_ID][i][2]);
			}		 
		}

		//If Station_Variable_name equals Tair_C then
		if (Station_Variable_name == "Tair_C") {
			//Tair_C[timeStep] equals MeteorologicalVariable_weightedSum
			input->Tair_C[timeStep] = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "Tdew_C") {
			input->Tdew_C[timeStep] = MeteorologicalVariable_weightedSum;
		}
		else if(Station_Variable_name == "RadiationNet_Wpm2") {
			input->RadiationNet_Wpm2[timeStep] = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "WindSpd_mps") {
			input->WindSpd_mps[timeStep] = MeteorologicalVariable_weightedSum;
		}
		//If Model_Selection is CoolBuilding
		else if (Station_Variable_name == "WindDir_deg") {
			input->WindDir_deg[timeStep] = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "AtmPres_kPa") {
			input->AtmPres_kPa[timeStep] = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "Precip_mpdt") {
			input->Precip_mpdt[timeStep] = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "Rain_m") {
			input->Rain_m[timeStep] = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "Snow_m") {
			input->Snow_m[timeStep] = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "TotalPrecip_m") {
			input->TotalPrecip_m = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "TotalRain_m") {
			input->TotalRain_m = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "TotalSnow_m") {
			input->TotalSnow_m = MeteorologicalVariable_weightedSum;
		}
		WeightedInputMap[MapPixel_ID][Station_Variable_name].push_back(MeteorologicalVariable_weightedSum);
	}
	//If Flag_ReadEvaporationData = 1 and SimulationStringParams["Model_Selection"] != SpatialTemperatureHydro then proceed to call readEvaporationFile function
	if (input->SimulationNumericalParams["Flag_ReadEvaporationData"] == 1 && input->SimulationStringParams["Model_Selection"] != "SpatialTemperatureHydro") {
		//2. applying weights to evaporation data and saving to map WeightedInputMap
		//For Loop with string name through MeteorologicalVariableName_Evaporation_string vector list
		for (string Station_Variable_name : MeteorologicalVariableName_Evaporation_string) {
			//MeteorologicalVariable_weightedSum initialized to 0
			double MeteorologicalVariable_weightedSum = 0;
			//For Loop with integer i through MultipleStations_NumberOfStationsAveraged
			for (int i = 0; i < input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]; i++) {
				//StationID_string defined as PixelNearestStation[MapPixel_ID][i][0], which is station i's StationID
				string StationID_string = PixelNearestStation[MapPixel_ID][i][0];
				//StationID_VariableName defined as combination of keynamepart1_name and Station_Variable_name
				string StationID_VariableName = StationID_string + "_" + Station_Variable_name;
				//MeteorologicalVariable_weightedSum is added to MeteorologicalVariable_weightedSum and product of EvaporationMap[StationID_VariableName] Value at timeStep and Station Weight, contained in PixelNearestStation[MapPixel_ID][i][2]
				MeteorologicalVariable_weightedSum = MeteorologicalVariable_weightedSum + EvaporationMap[StationID_VariableName][timeStep] * stod(PixelNearestStation[MapPixel_ID][i][2]);
			}

			//If Station_Variable_name equals PotentialEvapotranspiration_TreeCover_m then
			if (Station_Variable_name == "PotentialEvapotranspiration_TreeCover_m") {
				//PotentialEvapotranspiration_TreeCover_m[timeStep] equals MeteorologicalVariable_weightedSum
				input->PotentialEvapotranspiration_TreeCover_m[timeStep] = MeteorologicalVariable_weightedSum;
			}
			else if (Station_Variable_name == "PotentialEvaporation_WaterOnGround_m") {
				input->PotentialEvaporation_WaterOnGround_m[timeStep] = MeteorologicalVariable_weightedSum;
			}
			else if (Station_Variable_name == "PotentialEvaporation_WaterOnTree_m") {
				input->PotentialEvaporation_WaterOnTree_m[timeStep] = MeteorologicalVariable_weightedSum;
			}
			else if (Station_Variable_name == "PotentialSublimation_SnowOnGround_m") {
				input->PotentialSublimation_SnowOnGround_m[timeStep] = MeteorologicalVariable_weightedSum;
			}
			else if (Station_Variable_name == "PotentialSublimation_SnowOnTree_m") {
				input->PotentialSublimation_SnowOnTree_m[timeStep] = MeteorologicalVariable_weightedSum;
			}
			WeightedInputMap[MapPixel_ID][Station_Variable_name].push_back(MeteorologicalVariable_weightedSum);
		}
	}
	//3. applying weights to radiation data and saving to map WeightedInputMap
	//For Loop with string name through MeteorologicalVariableName_Radiation_string vector list
	for (string Station_Variable_name : MeteorologicalVariableName_Radiation_string) {
		//MeteorologicalVariable_weightedSum initialized to 0
		double MeteorologicalVariable_weightedSum = 0;
		//For Loop with integer i through MultipleStations_NumberOfStationsAveraged
		for (int i = 0; i < input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]; i++) {
			//StationID_string defined as PixelNearestStation[MapPixel_ID][i][0], which is station i's StationID
			string StationID_string = PixelNearestStation[MapPixel_ID][i][0];
			//StationID_VariableName defined as combination of keynamepart1_name and Station_Variable_name
			string StationID_VariableName = StationID_string + "_" + Station_Variable_name;
			//If Station_Variable_name not LatitudeMap then variable has vector timeStep structure
			if (Station_Variable_name != "Latitude_rad") {
				//MeteorologicalVariable_weightedSum is added to MeteorologicalVariable_weightedSum and product of RadiationMap[StationID_VariableName] Value at timeStep and Station Weight, contained in PixelNearestStation[MapPixel_ID][i][2]
				MeteorologicalVariable_weightedSum = MeteorologicalVariable_weightedSum + RadiationMap[StationID_VariableName][timeStep] * stod(PixelNearestStation[MapPixel_ID][i][2]);
			}
			//Else Station_Variable_name is LatitudeMap and has no vector timeStep structure
			else {
				//MeteorologicalVariable_weightedSum is added to MeteorologicalVariable_weightedSum and product of LatitudeMap[StationID_VariableName] Value and Station Weight, contained in PixelNearestStation[MapPixel_ID][i][2]
				MeteorologicalVariable_weightedSum = MeteorologicalVariable_weightedSum + LatitudeMap[StationID_VariableName] * stod(PixelNearestStation[MapPixel_ID][i][2]);
			}
		}
		//If Station_Variable_name equals Radiation_Shortwave_Direct_Wpm2 then
		if (Station_Variable_name == "Radiation_Shortwave_Direct_Wpm2") {
			//Radiation_Shortwave_Direct_Wpm2[timeStep] equals MeteorologicalVariable_weightedSum
			input->Radiation_Shortwave_Direct_Wpm2[timeStep] = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "Radiation_Shortwave_Diffuse_Wpm2") {
			input->Radiation_Shortwave_Diffuse_Wpm2[timeStep] = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "Radiation_Longwave_Downwelling_Wpm2") {
			input->Radiation_Longwave_Downwelling_Wpm2[timeStep] = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "Radiation_Longwave_Upwelling_Wpm2") {
			input->Radiation_Longwave_Upwelling_Wpm2[timeStep] = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "Latitude_rad") {
			input->Latitude_rad = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "HourAngle_Solar_rad") {
			input->HourAngle_Solar_rad[timeStep] = MeteorologicalVariable_weightedSum;
		}
		else if (Station_Variable_name == "DeclinationAngle_Solar_rad") {
			input->DeclinationAngle_Solar_rad[timeStep] = MeteorologicalVariable_weightedSum;
		}
		//WeightedInputMap[MapPixel_ID][Station_Variable_name] contains weighted average value of variable using distance to all weather stations
		//Note: WeightedInputMap[MapPixel_ID]["Radiation_Shortwave_Direct_Wpm2"][timeStep] is how it is written when accessed.
		WeightedInputMap[MapPixel_ID][Station_Variable_name].push_back(MeteorologicalVariable_weightedSum);
	}

	//4. applying weights to elevation
	//MeteorologicalVariable_weightedSum initialized to zero
	double MeteorologicalVariable_weightedSum = 0;
	//For Loop through MultipleStations_NumberOfStationsAveraged
	for (int i = 0; i < input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]; i++) {
		//StationID_int defined as PixelNearestStation[MapPixel_ID][i][0], which is station i's StationID
		string StationID_int = PixelNearestStation[MapPixel_ID][i][0];
		//StationWeight_frac (fraction) is influence of the station on the local pixel, computed as the inverse distance squared
		double StationWeight_frac = stod(PixelNearestStation[MapPixel_ID][i][2]);
		//StationElevation_m (m) is elevation of the station 
		double StationElevation_m = stod(StationInfoMap[StationID_int][4]);

		//MeteorologicalVariable_weightedSum is added to MeteorologicalVariable_weightedSum and product of StationInfoMap [4=Elevation] and Station Weight, contained in PixelNearestStation[MapPixel_ID][StationID_string][2]
		MeteorologicalVariable_weightedSum = MeteorologicalVariable_weightedSum + StationElevation_m * StationWeight_frac;
	}
	//WeightedInputMap[MapPixel_ID]["Elevation_WeatherStations_weighted_m"] receives value, at each time step which is not needed
	WeightedInputMap[MapPixel_ID]["Elevation_WeatherStations_weighted_m"].push_back(MeteorologicalVariable_weightedSum);
}

//MapPixelNearestWeatherStationWriter function writes output for each map pixel, its nearest weather stations distances and weights
//Note: PixelNearestStation[MapPixel_ID][a][b], where [a] contains 1 of n stations, and [b] contains station ID, distance, and weight
//Note: PixelNearestStation outer map key is MapPixel_ID, inner map key is reference station ranking, e.g., 1, 2, 3
//Note: PixelNearestStation vector value holds station ID, distance (km) and weight (fraction)
void WeatherProcessor::MapPixelNearestWeatherStationWriter( Inputs* input ) {

	//If Flag_ExtendedOutputs equals 1 then write z_WeightedInputMap.csv output, which for each pixel has time series of each weather variable
	if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {

		//Note: File z_NearestStationWeight.csv will contain weighted weather station data values, computed as spatially weighted averages of observed data, not computed by HeatFluxCal.cpp
		string File_NearestStationWeight_csv = input->SimulationStringParams["OutputFolder_Path"] + "z_NearestStationWeight.csv";
		ofstream outfileNearestStationWeight(File_NearestStationWeight_csv);
		if (!outfileNearestStationWeight.good()) {
			cout << "Warning: The file " << File_NearestStationWeight_csv << " could not be accessed for writing output." << endl;
			return;
		}

		//outfileNearestStationWeight file has header written
		outfileNearestStationWeight << "Map_pixelID,Map_Easting_X,Map_Northing_Y,Rank,WeatherStationID,Distance_km,Weight_frac" << endl;

		//For loop with auto it from PixelNearestStation.begin() to not PixelNearestStation.end()
		auto it_x = Map_Easting_X_m_vec.begin();
		auto it_y = Map_Northing_Y_m_vec.begin();


		//For loop with it through contents of map container with key and value pairs
		//Note: Use of loop for .first and .second to get key-value pairs
		for (auto it = PixelNearestStation.begin(); it != PixelNearestStation.end(); it++) {
			//For loop with auto it2 from ((*it).second).begin() to not ((*it).second).end()
			for (auto it2 = ((*it).second).begin(); it2 != ((*it).second).end(); it2++) {
				//outfileNearestStationWeight file has the following lines written ...
				outfileNearestStationWeight
					//it.first is the map pixel ID, the vector form of map row and col numbers
					<< (*it).first << ","
					//write out Albers xy location
					<< to_string(*it_x) << ","
					<< to_string(*it_y) << ","
					//it2.first is the station rank, with 1 being nearest
					<< (*it2).first << ","
					//it2.second[0] is the station ID, first column in input file WeatherStationAttributeTable.csv
					<< (*it2).second[0] << ","
					//it2.second[1] is the station distance (km)
					<< (*it2).second[1] << "," 
					//it2.second[2] is the station weight (fraction)
					<< (*it2).second[2] << endl;
			}
			it_x++;
			it_y++;
		}
		//outfileNearestStationWeight is closed at the last pixel
		outfileNearestStationWeight.close();
	}
}

//WeatherStationInformationWriter function writes optional output with ID, location, and characteristics of weather stations
void WeatherProcessor::WeatherStationInformationWriter(Inputs* input) {

	string File_StationInfoMap_csv = input->SimulationStringParams["OutputFolder_Path"] + "z_StationInfoMap.csv";
	ofstream outfileWeatherStationInfo(File_StationInfoMap_csv);
	if (!outfileWeatherStationInfo.good()) {
		cout << "Warning: The file " << File_StationInfoMap_csv << " could not be accessed for writing output." << endl;
		return;
	}

	//outfileWeatherStationInfo file has header written
	outfileWeatherStationInfo << "OBJECTID" << "," << "USAF" << "," << "WBAN" << "," << "Easting" << ", " << "Northing" << ", " << "Elev_m" << ", " << "Slope_deg" << "," << "Aspect_deg" << "," << "Station_TI_option" << "," << "NLCD_Class" << "," << "Station_TC_percent" << "," << "Station_IC_percent" << "," << "Station_SVeg_in_Gap_frac" << "," << "SimulationDate_Hottest_GD" << "," << "SimulationTime_Hottest_HMS"  << endl;
	
	//For loop with auto it from StationInfoMap.begin() to not StationInfoMap.end(); 
	//Note: StationInfoMap contains all weather stations for simulation, listed in file identified by HydroPlusConfig.xml element RefWeatherLocationAttributeFile 
	//Note: StationInfoMap is a map that has key value pairs, accessed in .first or .second position
	for (auto it = StationInfoMap.begin(); it != StationInfoMap.end(); it++) {
		//outfileWeatherStationInfo file has the following line written ...
		outfileWeatherStationInfo
			//(*it).first is weather station OBJECTID
			<< (*it).first;
			//For loop with auto it2 from (*it).second.begin() to not (*it).second.end()
			for (auto it2 = (*it).second.begin(); it2 != (*it).second.end(); it2++) {
				//outfileWeatherStationInfo file has all variables listed in header written ...
				outfileWeatherStationInfo
					<< "," << *it2;
			}
			outfileWeatherStationInfo << endl;
	}
	//outfileWeatherStationInfo file is closed 
	outfileWeatherStationInfo.close();
}

//MapPixelWeightedWeatherInputsWriter function writes optional output with row for each map pixel and its weighted meteorological inputs 
//Note: Each pixel has as many rows as there are meteorological variables, and each variable is written for each time step. 
//Note: Header is: Map_pixelID, Variable, TimeStep, etc. 
void WeatherProcessor::MapPixelWeightedWeatherInputsWriter(Inputs* input) {

	//If Flag_ExtendedOutputs equals 1 then write z_WeightedInputMap.csv output, which for each pixel has time series of each weather variable
	if (input->SimulationNumericalParams["Flag_ExtendedOutputs"] == 1) {

		//File_WeightedInputMap_csv defined as OutputDirectory and contains the file z_WeightedInputMap.csv
		string File_WeightedInputMap_csv = input->SimulationStringParams["OutputFolder_Path"] + "z_WeightedInputMap.csv";
		ofstream outfileInputMap(File_WeightedInputMap_csv);
		if (!outfileInputMap.good()) {
			cout << "Warning: The file " << File_WeightedInputMap_csv << " could not be accessed for writing output." << endl;
			return;
		}
		//outfileInputMap receiving header row with Map_pixelID and Variable and ...
		outfileInputMap << "Map_pixelID" << "," << " Variable" << ",";
		//For loop used to write to header TimeStep value for SimulationTimePeriod_timeSteps
		for (int i = 0; i < input->SimulationTimePeriod_timeSteps; i++) {
			//outfileInputMap receives TimeStep header
			//Note: TimeStep is appended to i, counter of SimulationTimePeriod_timeSteps by converting integer to string with to_string function
			outfileInputMap << "TimeStep" + to_string(i) << ",";
		}
		//outfileInputMap receives endl to close line
		outfileInputMap << endl;

		//For loop through all map pixels with auto it variable from WeightedInputMap.begin to WeightedInputMap.end
		for (auto it = WeightedInputMap.begin(); it != WeightedInputMap.end(); it++) {
			//Note: Variables .first and .second access key-value pairs in map or unordered_map containers
			//Note: (*it) and (*it2) are iterators pointing to elements of a map, and ...
			//Note: ... .first and .second are used to access the key and value of those elements.
			//Note: (*it).first: Accesses the key of the map element pointed to by the iterator it.
			//Note: (*it).second: Accesses the value associated with the key of the map element pointed to by the iterator it.
			//Note: (*it2).first: Accesses the key of the map element (likely an inner map) pointed to by the iterator it2.
			//Note: (*it2).second: Accesses the value associated with the key of the map element (likely an inner map) pointed to by the iterator it2.

			//For loop through all variables with auto it2 variable using pointer to *it.second to define begin and end of variable list
			for (auto it2 = ((*it).second).begin(); it2 != ((*it).second).end(); it2++) {
				//outfileInputMap receives content of 
				outfileInputMap
					//(*it).first is pointer to pixel ID
					<< (*it).first << ","
					//(*it2).first is pointer to variable name
					<< (*it2).first << ",";
				//For loop through all time steps to append the variable values
				for (int i = 0; i < input->SimulationTimePeriod_timeSteps; i++) {
					//outfileInputMap receives (*it2).second as pointer to variable value
					outfileInputMap << (*it2).second[i] << ",";
				}
				//outfileInputMap receives endl to close line
				outfileInputMap << endl;
			}
		}
		//outfileInputMap is closed at the last pixel
		outfileInputMap.close();
	}
}

//Tair_Mesoscale_byStationID_K and AbsHumidity_Mesoscale_byStationID_kg_p_m3 need declaration and definition
map<int, vector<double>> WeatherProcessor::Tair_Mesoscale_byStationID_K, WeatherProcessor::AbsHumidity_Mesoscale_byStationID_kg_p_m3;
//Tair_weighted_K, Tdew_weighted_K need declaration and definition
double WeatherProcessor::Tair_weighted_K, WeatherProcessor::Tdew_weighted_K, WeatherProcessor::Tair_mesoScale_weighted_K, 
WeatherProcessor::AbsHumidity_mesoScale_weighted_kg_p_m3;
//Tair_mesoScale_final_byStationID_K is 2D vector with (ID)(Tair_mesoScale_final_K)
map<int, vector<double>> WeatherProcessor::Tair_mesoScale_final_byStationID_K;
//H_total_byStationID_W_p_m2 is 2D vector with (ID)(H_total_W_p_m2)
map<int, vector<double>> WeatherProcessor::H_total_byStationID_W_p_m2;

//WeatherProcessor::ComputeSpatiallyWeighted_Tair_DEM_Variables function updates Tair variables by weight and elevation
//Note: For multiple stations, WeatherProcessor::ComputeSpatiallyWeighted_Tair_DEM_Variables uses ELR to update ...
//Note: ... Tair_weighted_K, which updates Tair_meso_K
void WeatherProcessor::ComputeSpatiallyWeighted_Tair_DEM_Variables(Inputs* input, int StationID_int, int MapPixel_ID, int timeStep) {

	//If PixelNearestStation.find(MapPixel_ID) == PixelNearestStation.end() then PixelNearestStation
	if (PixelNearestStation.find(MapPixel_ID) == PixelNearestStation.end()) {
		cerr << "Warning: PixelNearestStation does not contain MapPixel_ID " << MapPixel_ID << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: The HydroPlus model needs proper weather station data for the multiple station runs." << endl;
		cout << "Correction: Provide reasonable values in WeatherStationAttributeTable.csv." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	//Tair_mesoScale_weighted_K will contain a weighted value from multiple stations 
	Tair_mesoScale_weighted_K = 0;
	//AbsHumidity_mesoScale_weighted_kg_p_m3 will contain a weighted value from multiple stations 
	AbsHumidity_mesoScale_weighted_kg_p_m3 = 0;
	//Tair_weighted_K will contain a weighted value from multiple stations, presuming no variation in land cover between stations
	Tair_weighted_K = 0;
	//Tdew_weighted_K will contain a weighted value from multiple stations, presuming no variation in land cover between stations 
	Tdew_weighted_K = 0;
	//initilize Tair_mesoScale_final_weighted_K to 0 before each weighting process
	Tair_mesoScale_final_weighted_K = 0;
	//initilize H_total_weighted_W_p_m2 to 0 before each weighting process
	H_total_weighted_W_p_m2 = 0;
	//keyname_Tair_C and keyname_Tdew_C initialized
	string keyname_Tair_C, keyname_Tdew_C;

	//useSingleStation = (StationID_int >= 0) is true, but if StationID_int is negative then false
	bool useSingleStation = (StationID_int >= 0);

	//station_map holds reference to PixelNearestStation[MapPixel_ID], which has nearest station and its ID, weight, distance 
	auto& station_map = PixelNearestStation[MapPixel_ID];
	
	//If useSingleStation is true, then sending a single reference station, which is MapPixel_ID
	if (useSingleStation) {
		//stationData = station_map.at(0); this takes nearest station, with weight of 1
		//Note: Given reference station is MapPixel_ID, use only the top-ranked (nearest) station for this pixel
		const vector<string>& stationData = station_map.at(0);
		string StationID_string = stationData[0];
		StationID_int = stoi(StationID_string);
		//StationWeight_frac = 1.0; a weight of 1 is assigned when MapPixel_ID sits within reference station
		double StationWeight_frac = 1.0;

		//keyname_Tair_C is StationID_string + "_Tair_C"; Construct weather variable names
		string keyname_Tair_C = StationID_string + "_Tair_C";
		//keyname_Tdew_C is StationID_string + "_Tdew_C"; Construct weather variable names
		string keyname_Tdew_C = StationID_string + "_Tdew_C";

		//Tair_station_K is (WeatherMap[keyname_Tair_C][timeStep] + 273.15; Retrieve weather values
		double Tair_station_K = WeatherMap[keyname_Tair_C][timeStep] + 273.15;
		//Tdew_station_K is (WeatherMap[keyname_Tdew_C][timeStep] + 273.15; Retrieve weather values
		double Tdew_station_K = WeatherMap[keyname_Tdew_C][timeStep] + 273.15;

		//Tair_weighted_K += Tair_station_K * StationWeight_frac; Apply weighting
		Tair_weighted_K += Tair_station_K * StationWeight_frac;
		Tdew_weighted_K += Tdew_station_K * StationWeight_frac;
		Tair_mesoScale_weighted_K += Tair_Mesoscale_byStationID_K[StationID_int][timeStep] * StationWeight_frac;
		AbsHumidity_mesoScale_weighted_kg_p_m3 += AbsHumidity_Mesoscale_byStationID_kg_p_m3[StationID_int][timeStep] * StationWeight_frac;
		Tair_mesoScale_final_weighted_K += Tair_mesoScale_final_byStationID_K[StationID_int][timeStep] * StationWeight_frac;
		H_total_weighted_W_p_m2 += H_total_byStationID_W_p_m2[StationID_int][timeStep] * StationWeight_frac;
	}
	//Else useSingleStation is false, then sending a all reference station, which is MapPixel_ID
	else {
		// Loop through all nearby stations
		for (const auto& rank_entry : station_map) {
			//& stationData = rank_entry.second; 
			const vector<string>& stationData = rank_entry.second;

			//StationID_string = stationData[0]
			string StationID_string = stationData[0];
			StationID_int = stoi(StationID_string);
			//StationWeight_frac = stod(stationData[2]); location in vector of weight
			double StationWeight_frac = stod(stationData[2]);

			//keyname_Tair_C is StationID_string + "_Tair_C"; Construct weather variable names
			string keyname_Tair_C = StationID_string + "_Tair_C";
			string keyname_Tdew_C = StationID_string + "_Tdew_C";

			//Tair_station_K is WeatherMap[keyname_Tair_C][timeStep] + 273.15; Retrieve weather values
			double Tair_station_K = WeatherMap[keyname_Tair_C][timeStep] + 273.15;
			double Tdew_station_K = WeatherMap[keyname_Tdew_C][timeStep] + 273.15;

			//Tair_weighted_K += Tair_station_K * StationWeight_frac; Apply weighting
			Tair_weighted_K += Tair_station_K * StationWeight_frac;
			Tdew_weighted_K += Tdew_station_K * StationWeight_frac;
			Tair_mesoScale_weighted_K += Tair_Mesoscale_byStationID_K[StationID_int][timeStep] * StationWeight_frac;
			AbsHumidity_mesoScale_weighted_kg_p_m3 += AbsHumidity_Mesoscale_byStationID_kg_p_m3[StationID_int][timeStep] * StationWeight_frac;
			Tair_mesoScale_final_weighted_K += Tair_mesoScale_final_byStationID_K[StationID_int][timeStep] * StationWeight_frac;
			H_total_weighted_W_p_m2 += H_total_byStationID_W_p_m2[StationID_int][timeStep] * StationWeight_frac;
		}
	}
	//Taken from HeatFluxCal::AdiabaticLapseRates function:
	//Note: Even if Flag_MultipleStations=1 and WeatherProcessor lapse rate adjusted Tair_weighted_K, HeatFluxCal lapse rate adjustment is needed 
	//Note: Function calculates the effect of elevation difference with the reference station on local cell air temperature above ground
	//Note: Function modifies Stull (2000) theory of cooling by dry adiabatic lapse rate (ALR) with environmental lapse rate (ELR) of ~5.5 C/km, ...
	//Note: Assumes climbing 1 km up mountain air on mountain cools at ELR, while air lifting 1 km above ground cools at ALR.
	//Note: Stull (2000) gives dry adiabatic lapse rate Gamma_dry_C_p_km = 9.8 C/km; dew point lapse rate Gamma_dew_C_p_km = 1.8 C/km derived Eq 5.8
	//Note: Elevation_lcl_km = 1/(9.8 - 1.8) * (Tair_K - Tdew_K), Elevation_lcl_km = lifting condensation level based on Eq 5.8 Stull (2000)
	//Note: Elevation_lcl_km = 1/(5.5 - 1.8) * (Tair_K - Tdew_K), where Gamma_EnvironmentalLapseRate_C_p_km = 5.5 C/km, based on modification
	//Note: Elevation_delta_km = (Elevation_pixel - Elevation_referenceStation) and sets Tair_K = Tdew_K at Elevation_lcl_km
	//Note: If Elevation_delta_km <= Elevation_lcl_km Then Tair_update = Tair_K - Gamma_dry_C_p_km * Elevation_delta_km, based on Stull (2000)
	//Note: Else Tair_update = Tair_K - (Gamma_dry_C_p_km * Elevation_lcl_km) - [Gamma_dew_C_p_km * (Elevation_delta_km - Elevation_lcl_km)]; Stull (2000)
	//Note: If Elevation_delta_km <= Elevation_lcl_km Then Tair_update = Tair_K - Gamma_EnvironmentalLapseRate_C_p_km * Elevation_delta_km
	//Note: Else Tair_update = Tair_K - (Gamma_EnvironmentalLapseRate_C_p_km * Elevation_lcl_km) - [Gamma_dew_C_p_km * (Elevation_delta_km - Elevation_lcl_km)] based modification
	//Note: Tdew_update = Tdew_K - Gamma_dew_C_p_km * Elevation_delta_km based on Stull (2000) and modification
	//Note: For differencing temperatures, Kelvin can be used in place of Celsius, as they are linearly related
	//Stull, R. (2000). Meteorology for Scientists and Engineers, 2nd Edition. Brooks Cole. ISBN 0-534-37214-7

	//Elevation_Station_weighted_m (m) is spatially weighted average of station elevations, 
	//Note: WeightedInputMap[MapPixel_ID]["Elevation_WeatherStations_weighted_m"] is computed for each map pixel, only at timeStep = 0
	double Elevation_Station_weighted_m = WeightedInputMap[MapPixel_ID]["Elevation_WeatherStations_weighted_m"][0];
	//Elevation_delta_km (km) is height of pixel above station, equals (DEM(MapPixel_ID) - Elevation_Station_weighted_m), convert from m to km
	double Elevation_delta_km = (input->Elevation_DEM_m[MapPixel_ID] - Elevation_Station_weighted_m) * Ratio_km_to_m;
	
	//Elevation_lcl_km (km) equals 1 / (Gamma_EnvironmentalLapseRate_C_p_km - Gamma_dew_C_p_km) * (Tair_weighted_K - Tdew_weighted_K)
	//Note: Elevation_lcl_km (km) is lifting condensation level, modified from from Eq 5.8 Stull (2000) use the environmental lapse rate
	double Elevation_lcl_km = 1 / (Gamma_EnvironmentalLapseRate_C_p_km - Gamma_dew_C_p_km) * (Tair_weighted_K - Tdew_weighted_K);

	//If Elevation_delta_km is below Elevation_lcl_km, then the lifting condensation level was not reached and air remains unsaturated
	if (Elevation_delta_km <= Elevation_lcl_km) {
		//Tair_weighted_K (K) equals Tair_weighted_K - Gamma_EnvironmentalLapseRate_C_p_km * Elevation_delta_km
		//Note: Tair_weighted_K (K) modified from Eq 5.12 Stull (2000) as explained above to use the environmental lapse rate
		//Note: Even if Flag_MultipleStations=1 and WeatherProcessor lapse rate adjusted Tair_weighted_K, HeatFluxCal lapse rate adjustment is needed 
		Tair_weighted_K = Tair_weighted_K - Gamma_EnvironmentalLapseRate_C_p_km * Elevation_delta_km;
	}
	//Else, Tair_K cools at adiabatic lapse rate dry Gamma_dry_C_p_km until Elevation_lcl_km, then becomes saturated and cools at Gamma_dew_C_p_km
	else {
		//Tair_weighted_K (K) equals Tair_weighted_K - Gamma_EnvironmentalLapseRate_C_p_km * Elevation_lcl_km;
		//Note: Tair_weighted_K (K) is initially updated using environmental lapse rate until Elevation_lcl_km, from Eq 5.12 Stull (2000)
		Tair_weighted_K = Tair_weighted_K - Gamma_EnvironmentalLapseRate_C_p_km * Elevation_lcl_km;

		//Tair_localCell_K (K) is subsequently updated using adiabatic lapse rate for dew point air, modified Eq 5.13 Stull (2000)
		//Note: Tair_weighted_K (K) modified from Eq 5.12 Stull (2000) as explained above to use the environmental lapse rate
		//Note: Even if Flag_MultipleStations=1 and WeatherProcessor lapse rate adjusted Tair_weighted_K, HeatFluxCal lapse rate adjustment is needed 
		Tair_weighted_K = Tair_weighted_K - Gamma_dew_C_p_km * (Elevation_delta_km - Elevation_lcl_km);
	}
	//Tdew_weighted_K (K) is updated based on the adiabatic lapse rate for saturated dew point air, modified Eq 5.13 Stull (2000)
	//Note: Even if Flag_MultipleStations=1 and WeatherProcessor lapse rate adjusted Tdew_weighted_K, HeatFluxCal lapse rate adjustment is needed 
	Tdew_weighted_K = Tdew_weighted_K - Gamma_dew_C_p_km * Elevation_delta_km;
}

//Distance_Calc_LatLong function computes distance between pixel center and weather station using latitude and longitude, not currently used by HydroPlus
//Note: Distance_Calc_LatLong function is not currently called by HydroPlus, and remains in code in case there is a need for this calculation
//Note: Theory uses Haversine Formula https://en.wikipedia.org/wiki/Haversine_formula, computing distance between points on a sphere using lat and long 
long double WeatherProcessor::Distance_Calc_LatLong(long double PixelCenter_Latitude_DD, long double PixelCenter_Longitude_DD, long double Station_Latitude_DD, long double Station_Longitude_DD)
{
	//PixelCenter_Latitude_rad (rad) is converted from decimal degree to radians
	long double PixelCenter_Latitude_rad = toRadians(PixelCenter_Latitude_DD);
	//PixelCenter_Longitude_rad (rad) is converted from decimal degree to radians
	long double PixelCenter_Longitude_rad = toRadians(PixelCenter_Longitude_DD);
	//station_lat_rad (rad) is converted from decimal degree to radians
	long double station_lat_rad = toRadians(Station_Latitude_DD);
	//station_long_rad (rad) is converted from decimal degree to radians
	long double station_long_rad = toRadians(Station_Longitude_DD);
	//Distance_Longitude_rad (radians) is distance along longitude between station and pixel center
	long double Distance_Longitude_rad = station_long_rad - PixelCenter_Longitude_rad;
	//Distance_Latitude_rad (radians) is distance along latitude between station and pixel center
	long double Distance_Latitude_rad = station_lat_rad - PixelCenter_Latitude_rad;

	//Distance_Pixel_to_Station_rad (rad) is partial result of latitude longitude points distance formula
	//Note: Theory based on Haversine Formula
	long double Distance_Pixel_to_Station_rad = pow(sin(Distance_Latitude_rad / 2), 2) + cos(PixelCenter_Latitude_rad) * cos(station_lat_rad) * pow(sin(Distance_Longitude_rad / 2), 2);
	//Distance_Pixel_to_Station_rad is completed result of distance formula
	Distance_Pixel_to_Station_rad = 2 * asin(sqrt(Distance_Pixel_to_Station_rad));

	//RadiusEarth_km (km) is average radius of Earth elipsoid, given as 6371 km
	long double RadiusEarth_km = 6371;
	//Distance_Pixel_to_Station_rad (km) is scaled from radians to km using radius of Earth
	double Distance_Pixel_to_Station_km = Distance_Pixel_to_Station_rad * RadiusEarth_km;
	return Distance_Pixel_to_Station_km;
}

//WeatherProcessor::toRadians function converts degrees to radians
long double WeatherProcessor::toRadians(const long double degree) {
	// M_PI as the value of pi accurate to 1e-30
	long double one_deg = (M_PI) / 180;
	return (one_deg * degree);
}

//function to calculate time difference in seconds
int WeatherProcessor::calculateTimeStep(int current_YYYYMMDD, int current_HMS, int previous_YYYYMMDD, int previous_HMS, int defaultStep_sec) {
	//SimulationTimeStep_Duration_sec initialized as defaultStep_sec, SimulationNumericalParams["TimeStep_sec"]
	int SimulationTimeStep_Duration_sec = defaultStep_sec;
	//If timeStepMetData not 0 (e.g., 1st timeStepMetData w 0 indexing) then reading 2nd row of data and enter with two pairs of time
	if (timeStepMetData != 0) {
		//struct tm has instance previous_tm and current_tm initialized to 0
		struct tm previous_tm = { 0 }, current_tm = { 0 };

		//previous_tm struct items are populated by parsing previous timestep variables
		previous_tm.tm_year = (previous_YYYYMMDD / 10000) - 1900;
		previous_tm.tm_mon = ((previous_YYYYMMDD / 100) % 100) - 1;
		previous_tm.tm_mday = previous_YYYYMMDD % 100;
		previous_tm.tm_hour = previous_HMS / 10000;
		previous_tm.tm_min = (previous_HMS / 100) % 100;
		previous_tm.tm_sec = previous_HMS % 100;

		//current_tm struct items are populated by parsing current timestep variables
		current_tm.tm_year = (current_YYYYMMDD / 10000) - 1900;
		current_tm.tm_mon = ((current_YYYYMMDD / 100) % 100) - 1;
		current_tm.tm_mday = current_YYYYMMDD % 100;
		current_tm.tm_hour = current_HMS / 10000;
		current_tm.tm_min = (current_HMS / 100) % 100;
		current_tm.tm_sec = current_HMS % 100;

		//previous_time is a time_t object derived from mktime function being sent &previous_tm
		time_t previous_time = mktime(&previous_tm);
		//current_time is a time_t object derived from mktime function being sent &current_tm
		time_t current_time = mktime(&current_tm);
		//SimulationTimeStep_Duration_sec (seconds) is derived from difftime function sent current_time and previous_time
		SimulationTimeStep_Duration_sec = static_cast<int>(difftime(current_time, previous_time));
		return SimulationTimeStep_Duration_sec;
	}
	//Else timeStepMetData is 0, 1st time step
	else {
		//return SimulationTimeStep_Duration_sec as defaultStep_sec, SimulationNumericalParams["TimeStep_sec"]
		return SimulationTimeStep_Duration_sec;
	}
}

//DegreeHour_CreateRatios_From_Sums_by_YearMonthHour function computes the degree hours using DegreeHour_statisticMethod, for StationID not equal to -9999
//Note: If StationID is passed as -9999, it indicates Flag_MultipleStations != 1, not multiple station
void WeatherProcessor::DegreeHour_CreateRatios_From_Sums_by_YearMonthHour(Inputs* input, string DegreeHour_statisticMethod, int StationID) {

	//Resize and reset vectors
	DegreeHour_Year_Cooling_Ratio_vec.clear();
	DegreeHour_Year_Heating_Ratio_vec.clear();
	DegreeHour_Month_Cooling_Ratio_vec.assign(12, 1.0);
	DegreeHour_Hour_Cooling_K_vec.assign(24, 1.0);
	DegreeHour_Month_Heating_Ratio_vec.assign(12, 1.0);
	DegreeHour_Hour_Heating_K_vec.assign(24, 1.0);

	//Ensure vectors for multi-station case exist
	if (StationID != -99999) {
		DegreeHour_Year_Cooling_Ratio_StationID_vec_map[StationID].clear();
		DegreeHour_Year_Heating_Ratio_StationID_vec_map[StationID].clear();
		DegreeHour_Month_Cooling_Ratio_StationID_vec_map[StationID].assign(12, 1.0);
		DegreeHour_Hour_Cooling_K_StationID_vec_map[StationID].assign(24, 1.0);
		DegreeHour_Month_Heating_Ratio_StationID_vec_map[StationID].assign(12, 1.0);
		DegreeHour_Hour_Heating_K_StationID_vec_map[StationID].assign(24, 1.0);
	}
	//If Inputs::Convert_String_to_Lower_Case(DegreeHour_statisticMethod) == "median"
	//Note: Yearly and Monthly Ratios should default to 1, Hourly degree adjustment should default to 0
	if (Inputs::Convert_String_to_Lower_Case(DegreeHour_statisticMethod) == "median") {

		//Ensure all months (0–11) have an index within the vector, e.g., DegreeHour_Cooling_Month_abs_K_vec
		for (int Month_MM = 0; Month_MM < 12; ++Month_MM) {
			DegreeHour_Cooling_Month_abs_K_vec[Month_MM];
			DegreeHour_Heating_Month_abs_K_vec[Month_MM];
		}

		//Ensure all months hours (0–23) have an index within the vector, e.g., DegreeHour_Cooling_Hour_abs_K_vec
		for (int Hour_HH = 0; Hour_HH < 24; ++Hour_HH) {
			DegreeHour_Cooling_Hour_abs_K_vec[Hour_HH];
			DegreeHour_Heating_Hour_abs_K_vec[Hour_HH];
		}

		//===== Part I: Computing median cooling/heating Degree Hour for all years =====
		//DegreeHour_Cooling_All_abs_median_K initialized to 0
		double DegreeHour_Cooling_All_abs_median_K = 0.0;
		//If !DegreeHour_Cooling_All_abs_K_vec.empty() then elements present to enter
		if (!DegreeHour_Cooling_All_abs_K_vec.empty()) {
			//nth_element(DegreeHour_Cooling_All_abs_K_vec.begin(), DegreeHour_Cooling_All_abs_K_vec.begin() + DegreeHour_Cooling_All_abs_K_vec.size() / 2, DegreeHour_Cooling_All_abs_K_vec.end())
			//Note: nth_element is rapid function to sort data about median value
			nth_element(DegreeHour_Cooling_All_abs_K_vec.begin(), DegreeHour_Cooling_All_abs_K_vec.begin() + DegreeHour_Cooling_All_abs_K_vec.size() / 2, DegreeHour_Cooling_All_abs_K_vec.end());
			//DegreeHour_Cooling_All_abs_median_K = DegreeHour_Cooling_All_abs_K_vec[DegreeHour_Cooling_All_abs_K_vec.size() / 2]
			//Note: nth_element sort length divided by 2 yields median value
			DegreeHour_Cooling_All_abs_median_K = DegreeHour_Cooling_All_abs_K_vec[DegreeHour_Cooling_All_abs_K_vec.size() / 2];
		}
		//DegreeHour_Heating_All_abs_median_K initialized to 0
		double DegreeHour_Heating_All_abs_median_K = 0.0;
		//If !DegreeHour_Heating_All_abs_K_vec.empty() then elements present to enter
		if (!DegreeHour_Heating_All_abs_K_vec.empty()) {
			//nth_element(DegreeHour_Heating_All_abs_K_vec.begin(), DegreeHour_Heating_All_abs_K_vec.begin() + DegreeHour_Heating_All_abs_K_vec.size() / 2, DegreeHour_Heating_All_abs_K_vec.end())
			//Note: nth_element is rapid function to sort data about median value
			nth_element(DegreeHour_Heating_All_abs_K_vec.begin(), DegreeHour_Heating_All_abs_K_vec.begin() + DegreeHour_Heating_All_abs_K_vec.size() / 2, DegreeHour_Heating_All_abs_K_vec.end());
			//Note: nth_element sort length divided by 2 yields median value
			DegreeHour_Heating_All_abs_median_K = DegreeHour_Heating_All_abs_K_vec[DegreeHour_Heating_All_abs_K_vec.size() / 2];
		}

		//===== Part II: Computing median cooling Degree Hour for each year, month and hour =====
		//If for (const auto& pair : DegreeHour_Cooling_Year_abs_K_vec) 
		//Note: DegreeHour_Year_Cooling_Ratio_vec for Yearly ratios of cooling ratio
		for (const auto& year_and_DH_vec_pair : DegreeHour_Cooling_Year_abs_K_vec) {
			// Extract the vector of cooling degree hour values for a specific year
			vector<double> DH_vec = year_and_DH_vec_pair.second;
			//initialize year_median_ratio to 1
			double year_median_ratio = 1;
			//If (!DH_vec.empty() && DegreeHour_Cooling_All_abs_median_K > 0) then cooling occurs
			if (!DH_vec.empty() && DegreeHour_Cooling_All_abs_median_K > 0) {
				// Make a copy of the yearly DH values to compute the median
				vector<double> sorted_DH_vec = DH_vec;
				//nth_element(sorted_DH_vec.begin(), sorted_DH_vec.begin() + sorted_DH_vec.size() / 2, sorted_DH_vec.end())
				//Note: nth_element is rapid function to sort data about median value
				nth_element(sorted_DH_vec.begin(), sorted_DH_vec.begin() + sorted_DH_vec.size() / 2, sorted_DH_vec.end());
				//Note: nth_element sort length divided by 2 yields median value
				year_median_ratio = sorted_DH_vec[sorted_DH_vec.size() / 2] / DegreeHour_Cooling_All_abs_median_K;
			}
			//If StationID is passed as -99999, it indicates Flag_MultipleStations != 1, not multiple station
			if (StationID == -99999) {
				DegreeHour_Year_Cooling_Ratio_vec.push_back(year_median_ratio);
			}
			else {
				DegreeHour_Year_Cooling_Ratio_StationID_vec_map[StationID].push_back(year_median_ratio);
			}
		}

		//median for monthly
		//Note: DegreeHour_Cooling_Month_abs_K_vec for Monthly ratios of cooling ratio
		for (const auto& month_and_DH_vec_pair : DegreeHour_Cooling_Month_abs_K_vec) {
			//Extract the month
			int Month_MM = month_and_DH_vec_pair.first;
			//Extract the vector of cooling degree hour values for a specific month
			vector<double> DH_vec = month_and_DH_vec_pair.second;
			//initialize month_median_ratio to 1
			double month_median_ratio = 1;
			//If DH_vec is not empty and overall cooling median is valid, compute median ratio
			if (!DH_vec.empty() && DegreeHour_Cooling_All_abs_median_K > 0) {
				// Create a copy to compute the median using nth_element
				vector<double> sorted_DH_vec = DH_vec;
				//nth_element(sorted_DH_vec.begin(), sorted_DH_vec.begin() + sorted_DH_vec.size() / 2, sorted_DH_vec.end())
				//Note: nth_element is rapid function to sort data about median value
				nth_element(sorted_DH_vec.begin(), sorted_DH_vec.begin() + sorted_DH_vec.size() / 2, sorted_DH_vec.end());
				//Note: nth_element sort length divided by 2 yields median value
				month_median_ratio = sorted_DH_vec[sorted_DH_vec.size() / 2] / DegreeHour_Cooling_All_abs_median_K;
			}
			//If StationID is passed as -99999, it indicates Flag_MultipleStations != 1, not multiple station
			if (StationID == -99999) {
				DegreeHour_Month_Cooling_Ratio_vec[Month_MM] = month_median_ratio;
			}
			else {
				DegreeHour_Month_Cooling_Ratio_StationID_vec_map[StationID][Month_MM] = month_median_ratio;
			}
		}

		//median for hourly
		//Note: DegreeHour_Cooling_Hour_abs_K_vec for hourly ratios of cooling ratio
		for (const auto& hour_and_DH_vec_pair : DegreeHour_Cooling_Hour_abs_K_vec) {
			//Extract the hour
			int Hour_HH = hour_and_DH_vec_pair.first;
			//Extract the vector of cooling degree hour values for a specific hour
			vector<double> DH_vec = hour_and_DH_vec_pair.second;
			//initialize median_val to 0
			double median_val = 0;
			//If DH_vec is not empty, compute median ratio
			if (!DH_vec.empty()) {
				// Create a copy to compute the median using nth_element
				vector<double> sorted_DH_vec = DH_vec;
				//nth_element(sorted_DH_vec.begin(), sorted_DH_vec.begin() + sorted_DH_vec.size() / 2, sorted_DH_vec.end())
				//Note: nth_element is rapid function to sort data about median value
				nth_element(sorted_DH_vec.begin(), sorted_DH_vec.begin() + sorted_DH_vec.size() / 2, sorted_DH_vec.end());
				//Note: nth_element sort length divided by 2 yields median value
				median_val = sorted_DH_vec[sorted_DH_vec.size() / 2];
			}
			//If StationID is passed as -99999, it indicates Flag_MultipleStations != 1, not multiple station
			if (StationID == -99999) {
				DegreeHour_Hour_Cooling_K_vec[Hour_HH] = median_val;
			}
			else {
				DegreeHour_Hour_Cooling_K_StationID_vec_map[StationID][Hour_HH] = median_val;
			}
		}

		//===== Part III: Computing median heating Degree Hour for each year, month and hour =====
		//median for yearly
		//Note: DegreeHour_Heating_Year_abs_K_vec for hourly ratios of heating ratio
		for (const auto& year_and_DH_vec_pair : DegreeHour_Heating_Year_abs_K_vec) {
			//Extract the vector of heating degree hour values for a specific year
			vector<double> DH_vec = year_and_DH_vec_pair.second;
			//initialize year_median_ratio to 1
			double year_median_ratio = 1;
			//If DH_vec is not empty and overall heating median is valid, compute median ratio
			if (!DH_vec.empty() && DegreeHour_Heating_All_abs_median_K > 0) {
				// Create a copy to compute the median using nth_element
				vector<double> sorted_DH_vec = DH_vec;
				//nth_element(sorted_DH_vec.begin(), sorted_DH_vec.begin() + sorted_DH_vec.size() / 2, sorted_DH_vec.end())
				//Note: nth_element is rapid function to sort data about median value
				nth_element(sorted_DH_vec.begin(), sorted_DH_vec.begin() + sorted_DH_vec.size() / 2, sorted_DH_vec.end());
				//Note: nth_element sort length divided by 2 yields median value
				year_median_ratio = sorted_DH_vec[sorted_DH_vec.size() / 2] / DegreeHour_Heating_All_abs_median_K;
			}
			//If StationID is passed as -99999, it indicates Flag_MultipleStations != 1, not multiple station
			if (StationID == -99999) {
				DegreeHour_Year_Heating_Ratio_vec.push_back(year_median_ratio);
			}
			else {
				DegreeHour_Year_Heating_Ratio_StationID_vec_map[StationID].push_back(year_median_ratio);
			}
		}

		//median for monthly
		//Note: DegreeHour_Heating_Month_abs_K_vec for monthly ratios of heating ratio
		for (const auto& month_and_DH_vec_pair : DegreeHour_Heating_Month_abs_K_vec) {
			//Extract the month
			int Month_MM = month_and_DH_vec_pair.first;
			//Extract the vector of heating degree hour values for a specific month
			vector<double> DH_vec = month_and_DH_vec_pair.second;
			//initialize month_median_ratio to 1
			double month_median_ratio = 1;
			//If DH_vec is not empty and overall heating median is valid, compute median ratio
			if (!DH_vec.empty() && DegreeHour_Heating_All_abs_median_K > 0) {
				// Create a copy to compute the median using nth_element
				vector<double> sorted_DH_vec = DH_vec;
				//nth_element(sorted_DH_vec.begin(), sorted_DH_vec.begin() + sorted_DH_vec.size() / 2, sorted_DH_vec.end())
				//Note: nth_element is rapid function to sort data about median value
				nth_element(sorted_DH_vec.begin(), sorted_DH_vec.begin() + sorted_DH_vec.size() / 2, sorted_DH_vec.end());
				//Note: nth_element sort length divided by 2 yields median value
				month_median_ratio = sorted_DH_vec[sorted_DH_vec.size() / 2] / DegreeHour_Heating_All_abs_median_K;
			}
			//If StationID is passed as -99999, it indicates Flag_MultipleStations != 1, not multiple station
			if (StationID == -99999) {
				DegreeHour_Month_Heating_Ratio_vec[Month_MM] = month_median_ratio;
			}
			else {
				DegreeHour_Month_Heating_Ratio_StationID_vec_map[StationID][Month_MM] = month_median_ratio;
			}
		}

		//median for hourly
		//Note: DegreeHour_Heating_Hour_abs_K_vec for hourly ratios of heating ratio
		for (const auto& hour_and_DH_vec_pair : DegreeHour_Heating_Hour_abs_K_vec) {
			//Extract the hour
			int Hour_HH = hour_and_DH_vec_pair.first;
			//Extract the vector of heating degree hour values for a specific hour
			vector<double> DH_vec = hour_and_DH_vec_pair.second;
			//initialize median_val to 0
			double median_val = 0;
			//If DH_vec is not empty, compute median ratio
			if (!DH_vec.empty()) {
				// Create a copy to compute the median using nth_element
				vector<double> sorted_DH_vec = DH_vec;
				//nth_element(sorted_DH_vec.begin(), sorted_DH_vec.begin() + sorted_DH_vec.size() / 2, sorted_DH_vec.end())
				//Note: nth_element is rapid function to sort data about median value
				nth_element(sorted_DH_vec.begin(), sorted_DH_vec.begin() + sorted_DH_vec.size() / 2, sorted_DH_vec.end());
				//Note: nth_element sort length divided by 2 yields median value
				median_val = sorted_DH_vec[sorted_DH_vec.size() / 2];
			}
			//If StationID is passed as -99999, it indicates Flag_MultipleStations != 1, not multiple station
			if (StationID == -99999) {
				DegreeHour_Hour_Heating_K_vec[Hour_HH] = median_val;
			}
			else {
				DegreeHour_Hour_Heating_K_StationID_vec_map[StationID][Hour_HH] = median_val;
			}
		}
	}

	//if Inputs::Convert_String_to_Lower_Case(DegreeHour_statisticMethod) == "mean"
	//Note: Yearly and Monthly Ratios should default to 1, Hourly degree adjustment should default to 0
	if (Inputs::Convert_String_to_Lower_Case(DegreeHour_statisticMethod) == "mean") {

		//Ensure all months (0–11) have an index within the vectors, e.g., DH_Cooling_Month_sum
		//Note: DH_Cooling_Month_sum and other vectors had real values added in DegreeHour_SumHeatingCooling_by_YearMonthHour
		for (int Month_MM = 0; Month_MM < 12; ++Month_MM) {
			DH_Cooling_Month_sum[Month_MM] += 0.0;
			DH_Heating_Month_sum[Month_MM] += 0.0;
			DH_Month_Cooling_count[Month_MM] += 0;
			DH_Month_Heating_count[Month_MM] += 0;
		}

		//Ensure all hour (0–23) have an index within the vectors, e.g., DH_Cooling_Hour_sum
		//Note: DH_Cooling_Hour_sum and other vectors had real values added in DegreeHour_SumHeatingCooling_by_YearMonthHour
		for (int Hour_HH = 0; Hour_HH < 24; ++Hour_HH) {
			DH_Cooling_Hour_sum[Hour_HH] += 0.0;
			DH_Heating_Hour_sum[Hour_HH] += 0.0;
			DH_Hour_Cooling_count[Hour_HH] += 0;
			DH_Hour_Heating_count[Hour_HH] += 0;
		}

		//cooling average for total
		double DH_Cooling_total_avg = Inputs::SafeDivide(DH_Cooling_total_sum, DH_data_Cooling_total_count, 0.0);
		//heating average for total
		double DH_Heating_total_avg = Inputs::SafeDivide(DH_Heating_total_sum, DH_data_Heating_total_count, 0.0);

		//average for yearly
		for (const auto& pair : DH_Cooling_Year_sum) {
			int Year_YYYY = pair.first;
			double DegreeHour_sum = pair.second;
			double DegreeHour_avg = Inputs::SafeDivide(Inputs::SafeDivide(DegreeHour_sum, DH_Year_Cooling_count[Year_YYYY], 1.0), DH_Cooling_total_avg, 1.0);
			//for single station, push back DegreeHour_avg to DegreeHour_Year_Cooling_Ratio_vec
			if (StationID == -99999) { DegreeHour_Year_Cooling_Ratio_vec.push_back(DegreeHour_avg); }
			//for multiple station, push back 
			else { DegreeHour_Year_Cooling_Ratio_StationID_vec_map[StationID].push_back(DegreeHour_avg); }
		}

		//average for monthly
		for (const auto& pair : DH_Cooling_Month_sum) {
			int Month_MM = pair.first;
			double DegreeHour_sum = pair.second;
			double DegreeHour_avg = Inputs::SafeDivide(Inputs::SafeDivide(DegreeHour_sum, DH_Month_Cooling_count[Month_MM], 1.0), DH_Cooling_total_avg, 1.0);
			if (StationID == -99999) { DegreeHour_Month_Cooling_Ratio_vec[Month_MM] = DegreeHour_avg; }
			else { DegreeHour_Month_Cooling_Ratio_StationID_vec_map[StationID][Month_MM] = DegreeHour_avg; }
		}
		//average for hourly
		for (const auto& pair : DH_Cooling_Hour_sum) {
			int Hour_HH = pair.first;
			double DegreeHour_sum = pair.second;
			double DegreeHour_avg = Inputs::SafeDivide(DegreeHour_sum, DH_Hour_Cooling_count[Hour_HH], 0.0);
			if (StationID == -99999) { DegreeHour_Hour_Cooling_K_vec[Hour_HH] = DegreeHour_avg; }
			else { DegreeHour_Hour_Cooling_K_StationID_vec_map[StationID][Hour_HH] = DegreeHour_avg; }
		}

		//average for yearly
		for (const auto& pair : DH_Heating_Year_sum) {
			int Year_YYYY = pair.first;
			double DegreeHour_sum = pair.second;
			double DegreeHour_avg = Inputs::SafeDivide(Inputs::SafeDivide(DegreeHour_sum, DH_Year_Heating_count[Year_YYYY], 1.0), DH_Heating_total_avg, 1.0);
			//If StationID is passed as -99999, it indicates Flag_MultipleStations != 1, not multiple station
			if (StationID == -99999) { DegreeHour_Year_Heating_Ratio_vec.push_back(DegreeHour_avg); }
			else { DegreeHour_Year_Heating_Ratio_StationID_vec_map[StationID].push_back(DegreeHour_avg); }
		}

		//average for monthly
		for (const auto& pair : DH_Heating_Month_sum) {
			int Month_MM = pair.first;
			double DegreeHour_sum = pair.second;
			double DegreeHour_avg = Inputs::SafeDivide(Inputs::SafeDivide(DegreeHour_sum, DH_Month_Heating_count[Month_MM], 1.0), DH_Heating_total_avg, 1.0);
			//If StationID is passed as -99999, it indicates Flag_MultipleStations != 1, not multiple station
			if (StationID == -99999) { DegreeHour_Month_Heating_Ratio_vec[Month_MM] = DegreeHour_avg; }
			else { DegreeHour_Month_Heating_Ratio_StationID_vec_map[StationID][Month_MM] = DegreeHour_avg; }
		}

		//average for hourly
		for (const auto& pair : DH_Heating_Hour_sum) {
			int Hour_HH = pair.first;
			double DegreeHour_sum = pair.second;
			double DegreeHour_avg = Inputs::SafeDivide(DegreeHour_sum, DH_Hour_Heating_count[Hour_HH], 0.0);
			//If StationID is passed as -99999, it indicates Flag_MultipleStations != 1, not multiple station
			if (StationID == -99999) { DegreeHour_Hour_Heating_K_vec[Hour_HH] = DegreeHour_avg; }
			else {
				DegreeHour_Hour_Heating_K_StationID_vec_map[StationID][Hour_HH] = DegreeHour_avg;
			}
		}
	}

	//cleaned up for next station
	DegreeHour_Cooling_All_abs_K_vec.clear();
	DegreeHour_Heating_All_abs_K_vec.clear();
	DegreeHour_Cooling_Year_abs_K_vec.clear();
	DegreeHour_Cooling_Month_abs_K_vec.clear();
	DegreeHour_Cooling_Hour_abs_K_vec.clear(); 
	DegreeHour_Heating_Year_abs_K_vec.clear();
	DegreeHour_Heating_Month_abs_K_vec.clear();
	DegreeHour_Heating_Hour_abs_K_vec.clear();

	DH_Year_Cooling_count.clear();
	DH_Month_Cooling_count.clear();
	DH_Hour_Cooling_count.clear();
	DH_data_Cooling_total_count = 0;

	DH_Cooling_Year_sum.clear();
	DH_Cooling_Month_sum.clear();
	DH_Cooling_Hour_sum.clear();
	DH_Cooling_total_sum = 0;

	DH_Year_Heating_count.clear();
	DH_Month_Heating_count.clear();
	DH_Hour_Heating_count.clear();
	DH_data_Heating_total_count = 0;

	DH_Heating_Year_sum.clear();
	DH_Heating_Month_sum.clear();
	DH_Heating_Hour_sum.clear();
	DH_Heating_total_sum = 0;
}

///ComputeSpatiallyWeighted_DegreeHour function to work when HydroPlusConfig.xml element Flag_DegreeHour_Dependent_AH_Flux_Qcr = 1
void WeatherProcessor::ComputeSpatiallyWeighted_DegreeHour(Inputs* input, int MapPixel_ID) {

	//DegreeHour_Year_Cooling_Ratio_vec and other vectors cleared and sized
	DegreeHour_Year_Cooling_Ratio_vec.clear();
	DegreeHour_Year_Heating_Ratio_vec.clear();
	DegreeHour_Month_Cooling_Ratio_vec.assign(12, 0.0);
	DegreeHour_Hour_Cooling_K_vec.assign(24, 0.0);
	DegreeHour_Month_Heating_Ratio_vec.assign(12, 0.0);
	DegreeHour_Hour_Heating_K_vec.assign(24, 0.0);

	//For (int i = 0; i < input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]; i++) to loop weather stations
	for (int i = 0; i < input->SimulationNumericalParams["MultipleStations_NumberOfStationsAveraged"]; i++) {
		
		int StationID = stoi(PixelNearestStation[MapPixel_ID][i][0]);
		double StationWeight_frac = stod(PixelNearestStation[MapPixel_ID][i][2]);

		// Yearly cooling
		for (size_t vec_length = 0; vec_length < DegreeHour_Year_Cooling_Ratio_StationID_vec_map[StationID].size(); ++vec_length) {
			//since yearly vecs don't have fixed length, pushed back 0 ifvec.size() <= vec_length before assignment
			if (DegreeHour_Year_Cooling_Ratio_vec.size() <= vec_length) DegreeHour_Year_Cooling_Ratio_vec.push_back(0.0);
			DegreeHour_Year_Cooling_Ratio_vec[vec_length] += DegreeHour_Year_Cooling_Ratio_StationID_vec_map[StationID][vec_length] * StationWeight_frac;
		}

		// Yearly heating
		for (size_t vec_length = 0; vec_length < DegreeHour_Year_Heating_Ratio_StationID_vec_map[StationID].size(); ++vec_length) {
			//since yearly vecs don't have fixed length, pushed back 0 ifvec.size() <= vec_length before assignment
			if (DegreeHour_Year_Heating_Ratio_vec.size() <= vec_length) DegreeHour_Year_Heating_Ratio_vec.push_back(0.0);
			DegreeHour_Year_Heating_Ratio_vec[vec_length] += DegreeHour_Year_Heating_Ratio_StationID_vec_map[StationID][vec_length] * StationWeight_frac;
		}

		// Monthly cooling
		for (int m = 0; m < 12; ++m) {
			DegreeHour_Month_Cooling_Ratio_vec[m] += DegreeHour_Month_Cooling_Ratio_StationID_vec_map[StationID][m] * StationWeight_frac;
		}

		// Monthly heating
		for (int m = 0; m < 12; ++m) {
			DegreeHour_Month_Heating_Ratio_vec[m] += DegreeHour_Month_Heating_Ratio_StationID_vec_map[StationID][m] * StationWeight_frac;
		}

		// Hourly cooling
		for (int h = 0; h < 24; ++h) {
			DegreeHour_Hour_Cooling_K_vec[h] += DegreeHour_Hour_Cooling_K_StationID_vec_map[StationID][h] * StationWeight_frac;
		}

		// Hourly heating
		for (int h = 0; h < 24; ++h) {
			DegreeHour_Hour_Heating_K_vec[h] += DegreeHour_Hour_Heating_K_StationID_vec_map[StationID][h] * StationWeight_frac;
		}
	}
}

//| Year Range   | DST Start (2:00 a.m. local)   | DST End (2:00 a.m. local)     | Notes                           |
//|--------------|-------------------------------|-------------------------------|---------------------------------|
//| 1967–1973    | Last Sunday in April          | Last Sunday in October        | Uniform Time Act                |
//| 1974         | January 6                     | October 27                    | Energy crisis (special case)    |
//| 1975         | February 23                   | October 26                    | Energy crisis (special case)    |
//| 1976–1986    | Last Sunday in April          | Last Sunday in October        | Post-crisis standard            |
//| 1987–2006    | First Sunday in April         | Last Sunday in October        | Amended Uniform Time Act (1986) |
//| 2007–present | Second Sunday in March        | First Sunday in November      | Energy Policy Act of 2005       |
 
//Get_DaylightSavingTime_Periods function creates vectors pairs of DST transition dates, by reference DaylightSaving_StartStopPairs_YYYYMMDDHH
void WeatherProcessor::Get_DaylightSavingTime_Periods(int SimulationDateStart_YYYYMMDDHH, int SimulationDateStop_YYYYMMDDHH, vector<pair<int, int>>& DaylightSaving_StartStopPairs_YYYYMMDDHH) {
	
	// Extract start and end years from the simulation dates
	int StartYear = SimulationDateStart_YYYYMMDDHH / 1000000;
	int EndYear = SimulationDateStop_YYYYMMDDHH / 1000000;
	int DST_Transition_Hour_of_Day_HH = 2;

	//if (StartYear < 1976) StartYear = 1976
	//Note: Consider refactor to handle 1976 - 1986
	if (StartYear < 1976) {
		//Note: Consider refactor to send more detailed warning if year is earlier it is treated as 1976
		cout << "Notice: Daylight Saving Time will use 1976 start and stop dates" << endl;
		StartYear = 1976;
	}
	//For (int Year_YYYY = StartYear; Year_YYYY <= EndYear; Year_YYYY++) 
	//Note: Loop Year_YYYY from StartYear to EndYear
	for (int Year_YYYY = StartYear; Year_YYYY <= EndYear; Year_YYYY++) {
		//int Month_DST_Start_MM, Week_of_Month_DST_Start_WW, Month_DST_Stop_MM, Week_of_Month_DST_Stop_WW
		int Month_DST_Start_MM, Week_of_Month_DST_Start_WW, Month_DST_Stop_MM, Week_of_Month_DST_Stop_WW;
		//If (Year_YYYY >= 2007) then
		if (Year_YYYY >= 2007) {
			//Month_DST_Start_MM = 3, March
			Month_DST_Start_MM = 3;
			//Week_of_Month_DST_Start_WW = 2, 2nd Sunday
			Week_of_Month_DST_Start_WW = 2;
			//Month_DST_Stop_MM = 11, November
			Month_DST_Stop_MM = 11;
			//Week_of_Month_DST_Stop_WW = 1, 1st Sunday
			Week_of_Month_DST_Stop_WW = 1;
		}
		//else if (Year_YYYY >= 1987)
		else if (Year_YYYY >= 1987) {
			//Month_DST_Start_MM = 4, April
			Month_DST_Start_MM = 4;
			//Week_of_Month_DST_Start_WW = 1, 1st Sunday
			Week_of_Month_DST_Start_WW = 1;
			//Month_DST_Stop_MM = 10, October
			Month_DST_Stop_MM = 10;
			//Week_of_Month_DST_Stop_WW = -1, Last Sunday
			Week_of_Month_DST_Stop_WW = -1;
		}
		//else if (Year_YYYY >= 1976)
		else if (Year_YYYY >= 1976) {
			//Month_DST_Start_MM = 4, April
			Month_DST_Start_MM = 4;
			//Week_of_Month_DST_Start_WW = -1, Last Sunday
			Week_of_Month_DST_Start_WW = -1;
			//Month_DST_Stop_MM = 10, October
			Month_DST_Stop_MM = 10;
			//Week_of_Month_DST_Stop_WW = -1, Last Sunday
			Week_of_Month_DST_Stop_WW = -1;
		}
		//Else skip years before 1976
		else {
			//continue returns to for loop
			continue; 
		}

		//Find_DST_Transition_DayOfMonth = [&](int Year_YYYY, int Month_MM, int Week_of_Month_int) -> int 
		//Note: [&] Capture clause to captures variables by reference
		//Note: -> int Explicitly declares that the lambda returns an int
		//Note: Helper lambda function assumes Sunday is DST Transition by using 0 and 7 values below
		auto Find_DST_Transition_DayOfMonth = [&](int Year_YYYY, int Month_MM, int Week_of_Month_int) -> int {

			int DST_Transition_Day_of_Month_DD;
			//If (Week_of_Month_int > 0) then DST is not last Sunday of month
			if (Week_of_Month_int > 0) {
				//DateParts_Gregorian[3] = { Year_YYYY, Month_MM, 1 } gets 1st day of month
				int DateParts_Gregorian[3] = { Year_YYYY, Month_MM, 1 };
				//Date_JulianDate_base_4713BCE_int = Inputs::GregorianToJulianDate(DateParts_Gregorian) is Julian Date of 1st day
				int Date_JulianDate_base_4713BCE_int = Inputs::GregorianToJulianDate(DateParts_Gregorian);
				//Date_WeekDayCount_Code_0Sunday = Inputs::JulianDateToWeekday(Date_JulianDate_base_4713BCE_int) is weekday ID of 1st day
				int Date_WeekDayCount_Code_0Sunday = Inputs::JulianDateToWeekday(Date_JulianDate_base_4713BCE_int);
				//DST_Transition_Day_of_Month_DD = (7 - Date_WeekDayCount_Code_0Sunday) % 7 + (Week_of_Month_int - 1) * 7
				//Note: Function determines how many days to add to reach the next Sunday, combining part of week and whole week(s)
				//Note: (7 - Date_WeekDayCount_Code_0Sunday) % 7 yields 0 if Sunday, or 6 if Monday, ..., 1 if Saturday 
				//Note: (Week_of_Month_int - 1) * 7 yields 7 days per week of month, not counting the 1st week of month
				DST_Transition_Day_of_Month_DD = (7 - Date_WeekDayCount_Code_0Sunday) % 7 + (Week_of_Month_int - 1) * 7;
				//return DST_Transition_Day_of_Month_DD + 1
				//Note: 1 added to correct for indexing given, where day of month cannot start at 0
				return DST_Transition_Day_of_Month_DD + 1;
			}
			//Else DST is the last Sunday of Month
			else { 
				//daysInMonth_DD = Inputs::Month_Year_to_Days_in_Month(Year_YYYY, Month_MM) determines days in month
				int daysInMonth_DD = Inputs::Month_Year_to_Days_in_Month(Year_YYYY, Month_MM);
				//for (int dayToSearch_DD = daysInMonth_DD; dayToSearch_DD >= daysInMonth_DD - 6; --dayToSearch_DD)
				//Note: Reverse for loop from daysInMonth_DD to dayToSearch_DD >= daysInMonth_DD - 6, which contains the final 7 days of the month
				//Note: dayToSearch_DD range: if daysInMonth_DD=31: 25-31, if daysInMonth_DD=30: 24-30, if daysInMonth_DD=28: 22-28
				for (int dayToSearch_DD = daysInMonth_DD; dayToSearch_DD >= daysInMonth_DD - 6; --dayToSearch_DD) {
					//DateParts_Gregorian[3] = { Year_YYYY, Month_MM, dayToSearch_DD } gets dayToSearch_DD of month
					int DateParts_Gregorian[3] = { Year_YYYY, Month_MM, dayToSearch_DD };
					//Date_JulianDate_base_4713BCE_int = Inputs::GregorianToJulianDate(DateParts_Gregorian) is Julian Date of dayToSearch_DD
					int Date_JulianDate_base_4713BCE_int = Inputs::GregorianToJulianDate(DateParts_Gregorian);
					//Date_WeekDayCount_Code_0Sunday = Inputs::JulianDateToWeekday(Date_JulianDate_base_4713BCE_int) is weekday ID of dayToSearch_DD
					int Date_WeekDayCount_Code_0Sunday = Inputs::JulianDateToWeekday(Date_JulianDate_base_4713BCE_int);
					//If (Date_WeekDayCount_Code_0Sunday == 0) is 0 then 
					if (Date_WeekDayCount_Code_0Sunday == 0) { 
						//DST_Transition_Day_of_Month_DD = dayToSearch_DD
						//Note: dayToSearch_DD range: if daysInMonth_DD=31: 25-31, if daysInMonth_DD=30: 24-30, if daysInMonth_DD=28: 22-28
						DST_Transition_Day_of_Month_DD = dayToSearch_DD;
						//return DST_Transition_Day_of_Month_DD
						//Note: Nothing added to correct indexing due to it was derived from days in month, which indexes from 1
						return DST_Transition_Day_of_Month_DD;
					}
				}
			}
		};

		//Day_DST_Start_DD = Find_DST_Transition_DayOfMonth(Year_YYYY, Month_DST_Start_MM, Week_of_Month_DST_Start_WW)
		//Note: Call Find_DST_Transition_DayOfMonth to find Day_DST_Start_DD, day of month with DST start
		int Day_DST_Start_DD = Find_DST_Transition_DayOfMonth(Year_YYYY, Month_DST_Start_MM, Week_of_Month_DST_Start_WW);
		//Day_DST_Stop_DD = Find_DST_Transition_DayOfMonth(Year_YYYY, Month_DST_Stop_MM, Week_of_Month_DST_Stop_WW)
		//Note: Call Find_DST_Transition_DayOfMonth to find Day_DST_Stop_DD, day of month with DST stop
		int Day_DST_Stop_DD = Find_DST_Transition_DayOfMonth(Year_YYYY, Month_DST_Stop_MM, Week_of_Month_DST_Stop_WW);

		//DST_Start_DateTime = Year_YYYY * 1000000 + Month_DST_Start_MM * 10000 + Day_DST_Start_DD * 100 + DST_Transition_Hour_of_Day_HH 
		//Note: DST_Start_DateTime (YYYYMMDDHH) includes DST_Transition_Hour_of_Day_HH for hour 2 am at DST transition
		int DST_Start_DateTime = Year_YYYY * 1000000 + Month_DST_Start_MM * 10000 + Day_DST_Start_DD * 100 + DST_Transition_Hour_of_Day_HH;
		//DST_End_DateTime = Year_YYYY * 1000000 + Month_DST_Stop_MM * 10000 + Day_DST_Stop_DD * 100 + DST_Transition_Hour_of_Day_HH 
		//Note: DST_End_DateTime (YYYYMMDDHH) includes DST_Transition_Hour_of_Day_HH for hour 2 am at DST transition
		int DST_End_DateTime = Year_YYYY * 1000000 + Month_DST_Stop_MM * 10000 + Day_DST_Stop_DD * 100 + DST_Transition_Hour_of_Day_HH;

		//DaylightSaving_StartStopPairs_YYYYMMDDHH.emplace_back(DST_Start_DateTime, DST_End_DateTime)
		DaylightSaving_StartStopPairs_YYYYMMDDHH.emplace_back(DST_Start_DateTime, DST_End_DateTime);
	}
}
//bool IsDuring_DaylightSavingTime function will check if dateTime_YYYYMMDDHH is within known Daylight Saving Time based on DaylightSavingTime_Pairs 
bool WeatherProcessor::IsDuring_DaylightSavingTime(int dateTime_YYYYMMDDHH, const vector<pair<int, int>>& DaylightSavingTime_Pairs) {
	//For loop of DST_Pair in DaylightSavingTime_Pairs
	for (const auto& DST_Pair : DaylightSavingTime_Pairs) {
		//If dateTime_YYYYMMDDHH is witin DST_Pair.first && DST_Pair.second then return true
		if (dateTime_YYYYMMDDHH >= DST_Pair.first && dateTime_YYYYMMDDHH < DST_Pair.second) {
			//Return true
			return true;
		}
	}
	//Else return false
	return false;
}

//Apply_DaylightSavingTime_Adjustment function will determine if DST adjustment needs to roll time to new day, month, year
//Note: C++ std tm libary fully handles converting standard time to local time, observing Daylight Saving Time, but ...
//Note: ... this will not work if local machine is in wrong setting for DST, so function comes with too much vulnerability
void WeatherProcessor::Apply_DaylightSavingTime_Adjustment(int Date_StdT_YYYYMMDD, int Time_StdT_HH, int& Date_DST_YYYYMMDD, int& Date_DST_YYYYMMDDHH, string& Time_DST_HH) {

	//Date_StdT_YYYY = Date_StdT_YYYYMMDD / 10000; creates year by division (10,000)
	int Date_StdT_YYYY = Date_StdT_YYYYMMDD / 10000;
	//Date_StdT_MM = (Date_StdT_YYYYMMDD % 10000) / 100; creates month by modulus (10,000) then division by (100)
	int Date_StdT_MM = (Date_StdT_YYYYMMDD % 10000) / 100;
	//Date_StdT_DD = Date_StdT_YYYYMMDD % 100; creates day by modulus (100)
	int Date_StdT_DD = Date_StdT_YYYYMMDD % 100;
	//Date_DST_HH = Time_StdT_HH + 1; creates time advancement for Daylight Saving Time
	int Date_DST_HH = Time_StdT_HH + 1;

	//If (Date_DST_HH >= 24) then beyond a day; reset to 0 and increment the day
	if (Date_DST_HH >= 24) {
		//Date_DST_HH = 0
		Date_DST_HH = 0;
		//Date_StdT_DD = Date_StdT_DD + 1
		Date_StdT_DD = Date_StdT_DD + 1;

		//DaysInMonth = Inputs::Month_Year_to_Days_in_Month(Date_StdT_YYYY, Date_StdT_MM); get days in month to check if new month is needed
		int DaysInMonth = Inputs::Month_Year_to_Days_in_Month(Date_StdT_YYYY, Date_StdT_MM);
		//If (Date_StdT_DD > DaysInMonth) then advance month
		if (Date_StdT_DD > DaysInMonth) {
			//Date_StdT_DD = 1
			Date_StdT_DD = 1;
			//Date_StdT_MM = Date_StdT_MM + 1
			Date_StdT_MM = Date_StdT_MM + 1;

			//If (Date_StdT_MM > 12) then advance year
			if (Date_StdT_MM > 12) {
				///Date_StdT_MM = 1
				Date_StdT_MM = 1;
				//Date_StdT_YYYY = Date_StdT_YYYY + 1
				Date_StdT_YYYY = Date_StdT_YYYY + 1;
			}
		}
	}

	//Date_DST_YYYYMMDD = Date_StdT_YYYY * 10000 + Date_StdT_MM * 100 + Date_StdT_DD; create Daylight Saving Time
	Date_DST_YYYYMMDD = Date_StdT_YYYY * 10000 + Date_StdT_MM * 100 + Date_StdT_DD;
	//Date_DST_YYYYMMDDHH = Date_DST_YYYYMMDD * 100 + Date_DST_HH; create Daylight Saving Time
	Date_DST_YYYYMMDDHH = Date_DST_YYYYMMDD * 100 + Date_DST_HH;

	// Format the time as HH:MM:SS
	char HHMMSS_char[9];
	///snprintf(HHMMSS_char, sizeof(HHMMSS_char), "%02d:00:00", Date_DST_HH)
	//Note: Function populates HHMMSS_char based on formatting guidance from variables DateTime_DST_structure tm_hour, tm_min, tm_sec
	snprintf(HHMMSS_char, sizeof(HHMMSS_char), "%02d:00:00", Date_DST_HH);
	//Time_DST_HH = string(HHMMSS_char); create DST Hour
	Time_DST_HH = string(HHMMSS_char);
}

///Update_SimulationTime_Vectors will update the simulation time vectors for all but one case, which is with SimulationDate_Hottest_GD
//Note: When DST starts time jumps from 2 am to 3 am, skipping the hour of 2 am, which creates a gap in the output file time series
//Note: When DST ends time jumps from 2 am to 1 am, repeating the hour of 2 am, which creates a duplicate in the output file time series 
void WeatherProcessor::Update_SimulationTime_Vectors(Inputs* input, int Date_YYYYMMDD, string Time_HHMMSS_str, int date_InputDataRow_HH, int date_InputDataRow_GDH, bool Flag_OutputTime_DST, const vector<pair<int, int>>& DaylightSaving_StartStopPairs_YYYYMMDDHH, int Flag_Update_SimulationDate_GD_GDH) {

	//Time_HHMMSS_adjusted_HMS initialized for output with DST
	string Time_HHMMSS_adjusted_HMS;
	//Date_YYYYMMDD_adjusted_GD, Date_YYYYMMDDHH_adjusted_GDH initialized for output with DST
	int Date_YYYYMMDD_adjusted_GD, Date_YYYYMMDDHH_adjusted_GDH;

	//If (Flag_OutputTime_DST == 1 && IsDuring_DaylightSavingTime(date_InputDataRow_GDH, DaylightSaving_StartStopPairs_YYYYMMDDHH))
	//Note: Define Date_YYYYMMDD_adjusted_GD, Date_YYYYMMDDHH_adjusted_GDH, and Time_HHMMSS_adjusted_HMS based on Daylight Saving Time
	if (Flag_OutputTime_DST == 1 && IsDuring_DaylightSavingTime(date_InputDataRow_GDH, DaylightSaving_StartStopPairs_YYYYMMDDHH)) {

		//Apply_DaylightSavingTime_Adjustment(Date_YYYYMMDD, date_InputDataRow_HH, Date_YYYYMMDD_adjusted_GD, Date_YYYYMMDDHH_adjusted_GDH, Time_HHMMSS_adjusted_HMS)
		Apply_DaylightSavingTime_Adjustment(Date_YYYYMMDD, date_InputDataRow_HH, Date_YYYYMMDD_adjusted_GD, Date_YYYYMMDDHH_adjusted_GDH, Time_HHMMSS_adjusted_HMS);
	}
	//Else Define terms based on Standard Time, not Daylight Saving Time
	else {
		//Time_HHMMSS_adjusted_HMS is Time_HHMMSS_str
		Time_HHMMSS_adjusted_HMS = Time_HHMMSS_str;
		//Date_YYYYMMDD_adjusted_GD is Date_YYYYMMDD 
		Date_YYYYMMDD_adjusted_GD = Date_YYYYMMDD;
		//Date_YYYYMMDDHH_adjusted_GDH is date_InputDataRow_GDH
		Date_YYYYMMDDHH_adjusted_GDH = date_InputDataRow_GDH;
	}

	//SimulationTime_Output_HMS vector appended with Time_HHMMSS_adjusted_HMS for Standard Time or Daylight Saving Time
	input->SimulationTime_Output_HMS.push_back(Time_HHMMSS_adjusted_HMS);
	//SimulationDate_Output_GD created here for efficiency, and not of While loop with SimulationDate_Hottest_GD
	input->SimulationDate_Output_GD.push_back(Date_YYYYMMDD_adjusted_GD);
	//SimulationDate_Output_GDH created here for efficiency, and not of While loop with SimulationDate_Hottest_GD * 100 + SimulationDate_HH[timeStep]
	input->SimulationDate_Output_GDH.push_back(Date_YYYYMMDDHH_adjusted_GDH);

	//SimulationTime_HMS vector appended with Time_HHMMSS_str for Standard Time
	//Note: SimulationDate_YYYYMMDD and SimulationDate_GDH are not created or updated in this conditional
	input->SimulationTime_HMS.push_back(Time_HHMMSS_str);

	//If Flag_Update_SimulationDate_GD_GDH is true then append to standard-time vectors
	if (Flag_Update_SimulationDate_GD_GDH == true) {
		//SimulationDate_YYYYMMDD push_back Date_YYYYMMDD with Standard Time
		input->SimulationDate_YYYYMMDD.push_back(Date_YYYYMMDD);
		//SimulationDate_GDH push_back date_InputDataRow_GDH with Standard Time
		input->SimulationDate_GDH.push_back(date_InputDataRow_GDH);
	}
}
//isMissing will search if variable is present within input file during read
bool WeatherProcessor::isMissing(const string& variable_input) {
	return variable_input.empty() || variable_input == "NA" || variable_input == "null" || variable_input == "NULL";
}

//safe_stod function will check if it is safe to convert from string to double
double WeatherProcessor::safe_stod(const string& variable_input, const string& var_name, int line_num) {
	if (WeatherProcessor::isMissing(variable_input)) {
		cerr << "Error: Missing value for " << var_name << " on line " << line_num << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cerr << "Explanation: Expected columns in WeatherStationAttributeTable.csv:" << endl;
		cerr << WeatherStationAttributeTable_header;
		abort();
	}
	try {
		return stod(variable_input);
	}
	catch (const exception& e) {
		cerr << "Error: Converting " << var_name << " = '" << variable_input << "' to double on line " << line_num << ": " << e.what() << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cerr << "Explanation: Expected columns in WeatherStationAttributeTable.csv:" << endl;
		cerr << WeatherStationAttributeTable_header;
		abort();
	}
}

//safe_stod function will check if it is safe to convert from string to int
int WeatherProcessor::safe_stoi(const string& variable_input, const string& var_name, int line_num) {
	if (WeatherProcessor::isMissing(variable_input)) {
		cerr << "Error: Missing value for " << var_name << " on line " << line_num << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cerr << "Explanation: Expected columns in WeatherStationAttributeTable.csv:" << endl;
		cerr << WeatherStationAttributeTable_header;
		abort();
	}
	try {
		return stoi(variable_input);
	}
	catch (const exception& e) {
		cerr << "Error: Converting " << var_name << " = '" << variable_input << "' to int on line " << line_num << ": " << e.what() << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cerr << "Explanation: Expected columns in WeatherStationAttributeTable.csv:" << endl;
		cerr << WeatherStationAttributeTable_header;
		abort();
	}
}

//check_column_count function will check length of lines in WeatherStationAttributeTable.csv
bool WeatherProcessor::check_column_count(const string& line_of_input, int expected_columns, int line_num) {
	istringstream ss(line_of_input);
	string cell;
	int count = 0;

	while (getline(ss, cell, ',')) {
		++count;
	}

	if (count != expected_columns) {
		cerr << "Warning: WeatherStationAttributeTable.csv line " << line_num << " has " << count << " columns, but expected "
			<< expected_columns << " columns.\n";
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cerr << "Explanation: Line contents:\n" << line_of_input << "\n\n";
		cerr << "Explanation: Expected columns in WeatherStationAttributeTable.csv:" << endl;
		cerr << WeatherStationAttributeTable_header;
		//abort simulation
		abort();
	}
	return true;
}

//check_range function will check if WeatherStationAttributeTable.csv are reasonable values
void WeatherProcessor::check_range(double value_of_input, const string& var_name, int line_num, double min_val, double max_val) {
	if (value_of_input < min_val || value_of_input > max_val) {
		cerr << "Error: " << var_name << " = " << value_of_input << " is out of valid range (" << min_val << " to " << max_val << ") on line " << line_num << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cerr << "Explanation: Expected columns in WeatherStationAttributeTable.csv:" << endl;
		cerr << WeatherStationAttributeTable_header;
		//abort simulation
		abort();
	}
}

void WeatherProcessor::TrimWeatherMapToAdjustedTimesteps(int adjustedTimesteps, map<string, vector<double>>& StationDataMap) {
	for (auto& pair : StationDataMap) {
		vector<double>& vec = pair.second;
		if (vec.size() > adjustedTimesteps) {
			// Erase from the beginning, keeping only the last 'adjustedTimesteps' values
			vec.erase(vec.begin(), vec.end() - adjustedTimesteps);
		}
	}
}
