#include "MapResampler.h"
#include "Inputs.h"

//ResampleAnthropogenicHeatFlux_Map function to resample Qtot, Qcr, Qncr maps from their AH4GUC grid resolution of ~900m or 0.00833 deg
//Note: Downscaling algorithms developed by T. Endreny for smoothing the blocky nature of AH4GUC pixels originally at 0.0083 degree resolution
void MapResampler::ResampleAnthropogenicHeatFlux_Map(Inputs* input, vector<double>& AnthropogenicHeat_Flux_Qtot_Avg_Wpm2, vector<double>& AnthropogenicHeat_Flux_Qcr_Avg_Wpm2, vector<double>& AnthropogenicHeat_Flux_Qncr_Avg_Wpm2, const vector<double>& ImperviousCover_frac, int nCols, int nRows, double lowerLeftX_m, double lowerLeftY_m, double Map_Resample_Grid_m, double NODATA_code, int Flag_AH_Flux_Qcr_Qncr_not_Qtot, int Flag_AH_Flux_Maps_Resample, int Map_OriginalScaleGrid_m, double AnthropogenicHeat_Flux_Qtot_Map_Sum_Wpm2, double AnthropogenicHeat_Flux_Qcr_Map_Sum_Wpm2, double AnthropogenicHeat_Flux_Qncr_Map_Sum_Wpm2, double ImperviousCover_Sum_frac, string OutputFolder_Path, double RefWeatherLocationAH_Flux_Qtot_Avg_Wpm2)
{
    //Coefficient_Map_Sampler is 1.5 to ensure pixel search radius is more than ratio of coarse to base map to properly smooth AH Flux
    double Coefficient_Map_Sampler = 1.5;
    //Resample_influence_radius_pixels = round(Coefficient_Map_Sampler * (Map_OriginalScaleGrid_m / Map_Resample_Grid_m))
    //Note: The coarse map averages at 900 m for the Qtot, Qcr, Qncr products derived from AH4GUC grid resolution of 0.00833 deg
    //Note: Coefficient_Map_Sampler adjusts Map_OriginalScaleGrid_m / Map_Resample_Grid_m ratio
    int Resample_influence_radius_pixels = round(Coefficient_Map_Sampler * (Map_OriginalScaleGrid_m / Map_Resample_Grid_m));

    //AnthropogenicHeat_Flux_Avg_temp_Wpm2 is assigned AnthropogenicHeat_Flux_Qtot_Avg_Wpm2
    vector<double> AnthropogenicHeat_Flux_Avg_temp_Wpm2 = AnthropogenicHeat_Flux_Qtot_Avg_Wpm2;
    vector<double> AnthropogenicHeat_Flux_2nd_Avg_temp_Wpm2;
    //targetVector1 assigned reference AnthropogenicHeat_Flux_Qtot_Avg_Wpm2, so AnthropogenicHeat_Flux_Qtot_Avg_Wpm2 can be wiped clean
    vector<double>* targetVector1 = &AnthropogenicHeat_Flux_Qtot_Avg_Wpm2;
    //targetVector2 not needed in function, so set to nullptr
    vector<double>* targetVector2 = nullptr;

    //If Flag_AH_Flux_Qcr_Qncr_not_Qtot is 1 then reassign values for two AHE vectors
    //Note: Flag_AH_Flux_Qcr_Qncr_not_Qtot with Qcr is commercial residential AH flux, and Qncr is non commercial residential AH flux
    if (Flag_AH_Flux_Qcr_Qncr_not_Qtot == 1) {
        //AnthropogenicHeat_Flux_Avg_temp_Wpm2 is assigned AnthropogenicHeat_Flux_Qcr_Avg_Wpm2
        AnthropogenicHeat_Flux_Avg_temp_Wpm2 = AnthropogenicHeat_Flux_Qcr_Avg_Wpm2;
        //AnthropogenicHeat_Flux_2nd_Avg_temp_Wpm2 is assigned AnthropogenicHeat_Flux_Qcr_Avg_Wpm2
        AnthropogenicHeat_Flux_2nd_Avg_temp_Wpm2 = AnthropogenicHeat_Flux_Qncr_Avg_Wpm2;
        //targetVector1 assigned reference AnthropogenicHeat_Flux_Qcr_Avg_Wpm2, so AnthropogenicHeat_Flux_Qcr_Avg_Wpm2 can be wiped clean
        targetVector1 = &AnthropogenicHeat_Flux_Qcr_Avg_Wpm2;
        //targetVector1 assigned reference AnthropogenicHeat_Flux_Qncr_Avg_Wpm2, so AnthropogenicHeat_Flux_Qncr_Avg_Wpm2 can be wiped clean
        targetVector2 = &AnthropogenicHeat_Flux_Qncr_Avg_Wpm2;
    }

    //For loop of MapPixel_ID through AnthropogenicHeat_Flux_Avg_temp_Wpm2.size to process each pixel
    for (int MapPixel_ID = 0; MapPixel_ID < AnthropogenicHeat_Flux_Avg_temp_Wpm2.size(); MapPixel_ID++) {
        //If LandCover_NLCD_Class is NODATA_code then define all maps as NODATA_code
        if (input->LandCover_NLCD_Class[MapPixel_ID] == NODATA_code) {
            AnthropogenicHeat_Flux_Avg_temp_Wpm2[MapPixel_ID] = NODATA_code;
            AnthropogenicHeat_Flux_2nd_Avg_temp_Wpm2[MapPixel_ID] = NODATA_code;
        }
        //If AnthropogenicHeat_Flux_Avg_temp_Wpm2[MapPixel_ID] is NODATA_code then continue to next item in for loop
        if (AnthropogenicHeat_Flux_Avg_temp_Wpm2[MapPixel_ID] == NODATA_code) {
            continue;
        }
        //If RefWeatherLocationAH_Flux_Qtot_Avg_Wpm2 not -1 and isReferenceStationPixel(MapPixel_ID) is true then exit
        //Note: For reference station pixels, when RefWeatherLocationAH_Flux_Qtot_Avg_Wpm2 is not -1 then use HydroPlusConfig.xml value, not map value
        if (RefWeatherLocationAH_Flux_Qtot_Avg_Wpm2 != -1 && input->isReferenceStationPixel(MapPixel_ID)) {
            continue;
        }
        //row = MapPixel_ID / nCols; converting MapPixel_ID to row 
        //Note: Conversion with MapPixel_ID starting at 0, row & col starting at 0 for internal processing
        //Note: Eqs: row=MapPixel_ID / Inputs::nCols; col=MapPixel_ID % Inputs::nCols; MapPixel_ID = (row * nCols) + col
        int row = MapPixel_ID / nCols ;
        //col = MapPixel_ID % nCols; converting MapPixel_ID to col
        int col = MapPixel_ID % nCols;

        //If Flag_AH_Flux_Maps_Resample is 1 then resample AH Flux coarse map values so AH Flux is scaled across nearby impervious cover based on its fraction
        //Note: Cost: algorithm does not necessarily distribute AH Flux from original map to all impervious cover across domain due to search radius
        //Note: Benefit: algorithm maintains original map energy and distributes energy in pattern matching original map but focused on impervious cover
        //Note: Finding: Setting Flag_AH_Flux_Maps_Resample to 1 vs 2 performed worse for Washington, DC 20180828 against Shandas observations
        if (Flag_AH_Flux_Maps_Resample == 1) {
            //AnthropogenicHeat_Flux_Avg_mean_Wpm2 defined by calling ComputeMeanMapPixelValue function within the circular radius
            //Note: AnthropogenicHeat_Flux_Avg_temp_Wpm2 is Qtot or Qcr
            double AnthropogenicHeat_Flux_Avg_mean_Wpm2 = ComputeMeanMapPixelValue(AnthropogenicHeat_Flux_Avg_temp_Wpm2, row, col, Resample_influence_radius_pixels, nCols, nRows, NODATA_code);

            //AnthropogenicHeat_Flux_2nd_Avg_mean_Wpm2 will be defined by calling ComputeMedianMapPixelValue for Qcnr
            double AnthropogenicHeat_Flux_2nd_Avg_mean_Wpm2 = 0;
            //If Flag_AH_Flux_Qcr_Qncr_not_Qtot is 1 then call ComputeMedianMapPixelValue second to create AnthropogenicHeat_Flux_2nd_Avg_mean_Wpm2
            if (Flag_AH_Flux_Qcr_Qncr_not_Qtot == 1) {
                //Note: AnthropogenicHeat_Flux_2nd_Avg_temp_Wpm2 is Qcnr
                AnthropogenicHeat_Flux_2nd_Avg_mean_Wpm2 = ComputeMeanMapPixelValue(AnthropogenicHeat_Flux_2nd_Avg_temp_Wpm2, row, col, Resample_influence_radius_pixels, nCols, nRows, NODATA_code);
            }
            //ImperCoverData_median_frac defined by calling ComputeMedianMapPixelValue function within the circular radius, where ImperviousCover_frac is fraction values
            double ImperCoverData_median_frac = ComputeMedianMapPixelValue(ImperviousCover_frac, row, col, Resample_influence_radius_pixels, nCols, nRows, NODATA_code);

            //if ImperCoverData_median_frac > 0 then compute resampled value
            if (Inputs::isGreaterThan(ImperCoverData_median_frac, 0, Epsilon_Tolerance_1E_negative15)) {
                //*targetVector1[MapPixel_ID] equals (ImperviousCover_frac[MapPixel_ID] / ImperCoverData_median_frac) * AnthropogenicHeat_Flux_Avg_mean_Wpm2
                //Note: downscaling weight is ImperviousCover_frac/ImperCoverData_median_frac, where ImperCoverData_median_frac is for non-zero values
                (*targetVector1)[MapPixel_ID] = (ImperviousCover_frac[MapPixel_ID] / ImperCoverData_median_frac) * AnthropogenicHeat_Flux_Avg_mean_Wpm2;

                //If targetVector2 is true then
                if (targetVector2) {
                    //*targetVector1[MapPixel_ID] equals (ImperviousCover_frac[MapPixel_ID] / ImperCoverData_median_frac) * AnthropogenicHeat_Flux_2nd_Avg_mean_Wpm2
                    (*targetVector2)[MapPixel_ID] = (ImperviousCover_frac[MapPixel_ID] / ImperCoverData_median_frac) * AnthropogenicHeat_Flux_2nd_Avg_mean_Wpm2;
                }
            }
        }
        //If Flag_AH_Flux_Maps_Resample is 2 then resample AH Flux coarse map sum so AH Flux is scaled across all impervious cover based on its fraction
        //Note: Cost: algorithm does not necessarily maintain AH Flux pattern from original map when distributed to all impervious cover across domain 
        //Note: Benefit: algorithm maintains original map energy and distributes energy across all IC pixels in pattern regulated by impervious cover
        //Note: Finding: Setting Flag_AH_Flux_Maps_Resample to 2 vs 1 performed best for Washington, DC 20180828 against Shandas observations
        if (Flag_AH_Flux_Maps_Resample == 2) {
            //If ImperviousCover_Sum_frac <=0 then set to 1 for division
            if (ImperviousCover_Sum_frac <= 0) { ImperviousCover_Sum_frac = 1;}
            if (Flag_AH_Flux_Qcr_Qncr_not_Qtot == 0) {
                //*targetVector1[MapPixel_ID] is ImperviousCover_frac[MapPixel_ID] (frac) * (AnthropogenicHeat_Flux_Qtot_Map_Sum_Wpm2 / ImperviousCover_Sum_frac)
                //Note: Multiplication by ImperviousCover_frac (frac) distributes the area average AnthropogenicHeat_Flux_Wpm2 across domain 
                (*targetVector1)[MapPixel_ID] = ImperviousCover_frac[MapPixel_ID] * (AnthropogenicHeat_Flux_Qtot_Map_Sum_Wpm2 / ImperviousCover_Sum_frac);
            }
            else {
                //*targetVector1[MapPixel_ID] is ImperviousCover_frac[MapPixel_ID] (frac) * (AnthropogenicHeat_Flux_Qcr_Map_Sum_Wpm2 / ImperviousCover_Sum_frac)
                //Note: Multiplication by ImperviousCover_frac (frac) distributes the area average AnthropogenicHeat_Flux_Wpm2 across domain 
                (*targetVector1)[MapPixel_ID] = ImperviousCover_frac[MapPixel_ID] * (AnthropogenicHeat_Flux_Qcr_Map_Sum_Wpm2 / ImperviousCover_Sum_frac);
            }
            //If targetVector2 is true then
            if (targetVector2) {
                //*targetVector1[MapPixel_ID] is ImperviousCover_frac[MapPixel_ID] (frac) * (AnthropogenicHeat_Flux_Qncr_Map_Sum_Wpm2 / ImperviousCover_Sum_frac)
                //Note: Multiplication by ImperviousCover_frac (frac) distributes the area average AnthropogenicHeat_Flux_Wpm2 across domain 
                (*targetVector2)[MapPixel_ID] = ImperviousCover_frac[MapPixel_ID] * (AnthropogenicHeat_Flux_Qncr_Map_Sum_Wpm2 / ImperviousCover_Sum_frac);
            }
        }
    }

    //call WriteResampledMaps function to organize writing of resampled Qtot, Qcr, Qncr maps
    MapResampler::WriteResampledMaps(input, AnthropogenicHeat_Flux_Qtot_Avg_Wpm2, AnthropogenicHeat_Flux_Qcr_Avg_Wpm2, AnthropogenicHeat_Flux_Qncr_Avg_Wpm2, OutputFolder_Path, nCols, nRows, lowerLeftX_m, lowerLeftY_m, Map_Resample_Grid_m, NODATA_code, Flag_AH_Flux_Qcr_Qncr_not_Qtot);
}

//ComputeMeanMapPixelValue function to compute Mean within the search_radius
double MapResampler::ComputeMeanMapPixelValue(const vector<double>& map_data_to_search, int centerRow, int centerCol, int search_radius, int nCols, int nRows, double NODATA_code) {

    //map_data_sum and map_data_count initialized to 0
    double map_data_sum = 0.0;
    int map_data_count = 0;

    //For loop of dx from -search_radius to search_radius distance for computing median
    for (int dx = -search_radius; dx <= search_radius; dx++) {
        //For loop of dy from -search_radius to search_radius distance for computing median
        for (int dy = -search_radius;  dy <= search_radius; dy++) {
            //newRow is centerRow + dx, where centerRow is x location of pixel about which the median is computed
            int newRow = centerRow + dx;
            //newCol is centerRow + dy, where centerCol is y location of pixel about which the median is computed
            int newCol = centerCol + dy;

            //If newRow >= 0 and new Row < nRows and newCol >= 0 and newCol < nCols then enter
            if (newRow >= 0 && newRow < nRows && newCol >= 0 && newCol < nCols) {
                //MapPixel_ID is newRow * nCols + newCol, formula to convert row and column to MapPixel_ID
                int MapPixel_ID = newRow * nCols + newCol;
                //If map_data_to_search[MapPixel_ID] not NODATA_code and sqrt(dx * dx + dy * dy) <= search_radius then data is valid and within search search_radius
                if (map_data_to_search[MapPixel_ID] != NODATA_code 
                   && (dx * dx + dy * dy) <= search_radius * search_radius){
                    //map_data_sum increased by map_data_to_search[MapPixel_ID]
                    map_data_sum += map_data_to_search[MapPixel_ID];
                    //map_data_count incremented
                    map_data_count++;
                }
            }
        }
    }

    //mean_value is map_data_sum / map_data_count
    double mean_value = map_data_sum / map_data_count;
    //return based on ternary, mean_value or NODATA_code
    return (map_data_count > 0) ? map_data_sum / map_data_count : NODATA_code;
}

//ComputeMedianMapPixelValue function to compute Median within the search_radius of a search area and for non-zero data values
//Note: Median of non-zero values is critical for getting reasonable ratios of ic/ic_median
double MapResampler::ComputeMedianMapPixelValue(const vector<double>& map_data_to_search, int centerRow, int centerCol, int search_radius, int nCols, int nRows, double NODATA_code) {

    //map_data_to_analyze takes input magnitude and sign for computing median
    vector<double> map_data_to_analyze;

    //For loop of dx from -search_radius to search_radius distance for computing median
    for (int dx = -search_radius;  dx <= search_radius; dx++) {
        //For loop of dy from -search_radius to search_radius distance for computing median
        for (int dy = -search_radius; dy <= search_radius; dy++) {
            //newRow is centerRow + dx, where centerRow is x location of pixel about which the median is computed
            int newRow = centerRow + dx;
            //newCol is centerRow + dy, where centerCol is y location of pixel about which the median is computed
            int newCol = centerCol + dy;

            //If newRow >= 0 and new Row < nRows and newCol >= 0 and newCol < nCols then enter
            if (newRow >= 0 && newRow < nRows && newCol >= 0 && newCol < nCols) {
                //MapPixel_ID is newRow * nCols + newCol, formula to convert row and column to MapPixel_ID
                //Note: Conversion with MapPixel_ID starting at 0, row & col starting at 0 for internal processing
                //Note: Eqs: row=MapPixel_ID / Inputs::nCols; col=MapPixel_ID % Inputs::nCols; MapPixel_ID = (row * nCols) + col
                int MapPixel_ID = newRow * nCols + newCol;
                //If map_data_to_search > 0, not NODATA_code, and sqrt(dx * dx + dy * dy) <= search_radius then data is valid and within search search_radius
                //Note: sqrt(dx * dx + dy * dy) <= search_radius is conditional to ensure search is circular in shape, not rectangular
                if (map_data_to_search[MapPixel_ID] > 0 &&
                    map_data_to_search[MapPixel_ID] != NODATA_code &&
                    (dx * dx + dy * dy) <= search_radius * search_radius) {
                    //map_data_to_analyze vector has map_data_to_search added with push_back function
                    map_data_to_analyze.push_back(map_data_to_search[MapPixel_ID]);
                }
            }
        }
    }

    //If map_data_to_analyzeis empty return NODATA_code
    if (map_data_to_analyze.empty()) {
        return NODATA_code;
    }

    //else
    //nth_element modifies map_data_to_analyze, only ensures that the middle element is in the correct position and not sorting other data
    nth_element(map_data_to_analyze.begin(), map_data_to_analyze.begin() + map_data_to_analyze.size() / 2, map_data_to_analyze.end());
    //median_value = map_data_to_analyze[map_data_to_analyze.size() / 2], which is the 50% index in partially sorted map_data_to_analyze
    double median_value = map_data_to_analyze[map_data_to_analyze.size() / 2];
    //return median_value
    return median_value;
}

//WriteResampledMaps function to organize writing of resampled Qtot, Qcr, Qncr maps
void MapResampler::WriteResampledMaps(Inputs* input, const vector<double>& AnthropogenicHeat_Flux_Qtot_Avg_Wpm2, const vector<double>& AnthropogenicHeat_Flux_Qcr_Avg_Wpm2, const vector<double>& AnthropogenicHeat_Flux_Qncr_Avg_Wpm2, const string& outputFolderPath, int nCols, int nRows, double lowerLeftX_m, double lowerLeftY_m, double cellSize_m, double NODATA_code, int Flag_AH_Flux_Qcr_Qncr_not_Qtot)
{
    //define resampled AnthropogenicHeat_Flux_Qtot_Avg_Wpm2_str, Qcr, Qcnr ASCII map filenames with _new extension
    string filename_Qtot = outputFolderPath + "/" + input->AnthropogenicHeat_Flux_Qtot_Avg_Wpm2_str.substr(0, input->AnthropogenicHeat_Flux_Qtot_Avg_Wpm2_str.find(".asc")) + "_new.asc";
    string filename_Qcr = outputFolderPath + "/" + input->AnthropogenicHeat_Flux_Qcr_Avg_Wpm2_str.substr(0, input->AnthropogenicHeat_Flux_Qcr_Avg_Wpm2_str.find(".asc")) + "_new.asc";
    string filename_Qncr = outputFolderPath + "/" + input->AnthropogenicHeat_Flux_Qncr_Avg_Wpm2_str.substr(0, input->AnthropogenicHeat_Flux_Qncr_Avg_Wpm2_str.find(".asc")) + "_new.asc";


    //If Flag_AH_Flux_Qcr_Qncr_not_Qtot == 0 then write Qtot
    if (Flag_AH_Flux_Qcr_Qncr_not_Qtot == 0) {
        // Write each file if the vector is not empty
        if (!AnthropogenicHeat_Flux_Qtot_Avg_Wpm2.empty()) {
            //call WriteRasterToFile function to write resampled Qtot map for future use
            WriteRasterToFile(filename_Qtot, AnthropogenicHeat_Flux_Qtot_Avg_Wpm2, nCols, nRows, lowerLeftX_m, lowerLeftY_m, cellSize_m, NODATA_code);
        }
    }
    //Else If Flag_AH_Flux_Qcr_Qncr_not_Qtot not 0 then write Qcr and Qncr
    else {
        if (!AnthropogenicHeat_Flux_Qcr_Avg_Wpm2.empty()) {
            //call WriteRasterToFile function to write resampled Qcr map for future use
            WriteRasterToFile(filename_Qcr, AnthropogenicHeat_Flux_Qcr_Avg_Wpm2, nCols, nRows, lowerLeftX_m, lowerLeftY_m, cellSize_m, NODATA_code);
        }
        if (!AnthropogenicHeat_Flux_Qncr_Avg_Wpm2.empty()) {
            //call WriteRasterToFile function to write resampled Qncr map for future use
            WriteRasterToFile(filename_Qncr, AnthropogenicHeat_Flux_Qncr_Avg_Wpm2, nCols, nRows, lowerLeftX_m, lowerLeftY_m, cellSize_m, NODATA_code);
        }
    }
}

//WriteRasterToFile helper function to write an individual raster to ASCII map
void MapResampler::WriteRasterToFile(const string& filename, const vector<double>& aha_data_value_Wpm2, int nCols, int nRows, double lowerLeftX_m, double lowerLeftY_m, double cellSize_m, double NODATA_code)
{
    //ofstream file_output_map from filename
    ofstream file_output_map(filename);
    //If not file_output_map then notify user 
    if (!file_output_map) {
        cerr << "Error: Cannot open file_output_map " << filename << " for writing.\n";
        return;
    }

    //file_output_map write header for ASCII map
    file_output_map << "nCols         " << nCols << "\n";
    file_output_map << "nRows         " << nRows << "\n";
    file_output_map << "xllcorner     " << lowerLeftX_m << "\n";
    file_output_map << "yllcorner     " << lowerLeftY_m << "\n";
    file_output_map << "cellsize      " << cellSize_m << "\n";
    file_output_map << "NODATA_value  " << NODATA_code << "\n";

    //For loop of row through map until nRows
    for (int row = 0; row < nRows; row++) {
        //For loop of col through map until nCols
        for (int col = 0; col < nCols; col++) {
            //MapPixel_ID equals row * nCols + col;
            int MapPixel_ID = row * nCols + col;

            //AnthropogenicHeat_Flux_Wpm2 = (aha_data_value_Wpm2[MapPixel_ID] != NODATA_code) ternary of actual value or NODATA_code
            double AnthropogenicHeat_Flux_Wpm2 = (aha_data_value_Wpm2[MapPixel_ID] != NODATA_code)
                ? aha_data_value_Wpm2[MapPixel_ID]
                : NODATA_code;

            //file_output_map is fixed setprecision(2)
            file_output_map << fixed << setprecision(2) << AnthropogenicHeat_Flux_Wpm2 << " ";
        }

        //file_output_map has line return at end of row
        file_output_map << "\n"; 
    }

    //file_output_map closed
    file_output_map.close();
}
