C++ LSDRaster::get_NoDataValue方法代码示例

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
// Create function that takes the dimensions and georeferencing of a raster
// but then sets all data to value, setting the NoDataValues to
// the NoData of the raster
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
void LSDSoilHydroRaster::create(LSDRaster& OtherRaster, float value)
{
  NRows = OtherRaster.get_NRows();
  NCols = OtherRaster.get_NCols();
  XMinimum = OtherRaster.get_XMinimum();
  YMinimum = OtherRaster.get_YMinimum();
  DataResolution = OtherRaster.get_DataResolution();
  NoDataValue = OtherRaster.get_NoDataValue();
  GeoReferencingStrings = OtherRaster.get_GeoReferencingStrings();

  // set the raster data to be a certain value
  Array2D<float> data(NRows,NCols,NoDataValue);

  for (int row = 0; row <NRows; row++)
  {
    for (int col = 0; col<NCols; col++)
    {

      if (OtherRaster.get_data_element(row,col) != NoDataValue)
      {
        data[row][col] = value;
        //cout << value << endl;
      }
    }
  }

  RasterData = data.copy();
}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
// This prints a chi map to csv with an area threshold in m^2
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
void LSDChiTools::chi_map_to_csv(LSDFlowInfo& FlowInfo, string chi_map_fname, 
                                 float A_0, float m_over_n, float area_threshold)
{
  
  ofstream chi_map_csv_out;
  chi_map_csv_out.open(chi_map_fname.c_str());
  
  chi_map_csv_out.precision(9);
  
  float chi_coord;
  double latitude,longitude;
  LSDCoordinateConverterLLandUTM Converter;
  
  chi_map_csv_out << "latitude,longitude,chi" << endl;
  
  LSDRaster Chi = FlowInfo.get_upslope_chi_from_all_baselevel_nodes(m_over_n, A_0, area_threshold);
  
  float NDV = Chi.get_NoDataValue();

  for(int row = 0; row<NRows; row++)
  {
    for(int col = 0; col<NCols; col++)
    {
      chi_coord =  Chi.get_data_element(row,col);
      
      if (chi_coord != NDV)
      {
        get_lat_and_long_locations(row, col, latitude, longitude, Converter);
        chi_map_csv_out << latitude << "," << longitude  << "," << chi_coord << endl;
      }
    }
  }
  
  chi_map_csv_out.close();
}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
// Creates an LSDChiTools from an LSDRaster 
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
void LSDChiTools::create(LSDRaster& ThisRaster)
{
  NRows = ThisRaster.get_NRows();
  NCols = ThisRaster.get_NCols();
  XMinimum = ThisRaster.get_XMinimum();
  YMinimum = ThisRaster.get_YMinimum();
  DataResolution = ThisRaster.get_DataResolution();
  NoDataValue = ThisRaster.get_NoDataValue();
  GeoReferencingStrings = ThisRaster.get_GeoReferencingStrings();
}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
//
// Create function that reads from a raster
//
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
void LSDRasterInfo::create(LSDRaster& Raster)
{
      ///Number of rows.
  NRows = Raster.get_NRows();
  NCols = Raster.get_NCols();
  XMinimum = Raster.get_XMinimum();
  YMinimum = Raster.get_YMinimum();
  DataResolution = Raster.get_DataResolution();
  NoDataValue = Raster.get_NoDataValue();
  GeoReferencingStrings = Raster.get_GeoReferencingStrings();
}

//.........这里部分代码省略.........
    raster_selection.push_back(0);
    raster_selection.push_back(0);

    int CurrentWindowSize = window_radius;

    vector<LSDRaster> Surfaces = FilledDEM.calculate_polyfit_surface_metrics(CurrentWindowSize, raster_selection);
    LSDRaster slope = Surfaces[1];
    LSDRaster aspect = Surfaces[2];

    cout << "\nGetting drainage network and basins\n" << endl;

    // get a flow info object
  	LSDFlowInfo FlowInfo(BoundaryConditions,FilledDEM);

    //get stream net from channel heads
    stringstream ss_ch;
    ss_ch << Path << "Pelletier/" << Prefix << "_" << Resolutions[a] << "_CH";

    vector<int> sources = FlowInfo.Ingest_Channel_Heads(ss_ch.str(), "bil");
    LSDJunctionNetwork ChanNetwork(sources, FlowInfo);
    LSDIndexRaster StreamNetwork = ChanNetwork.StreamOrderArray_to_LSDIndexRaster();

    //Extract basins based on input stream order
    vector< int > basin_junctions = ChanNetwork.ExtractBasinJunctionOrder(BasinOrder, FlowInfo);
    LSDIndexRaster Basin_Raster = ChanNetwork.extract_basins_from_junction_vector(basin_junctions, FlowInfo);

    cout << "\nExtracting hilltops and hilltop curvature" << endl;

    // extract hilltops - no critical slope filtering is performed here
    LSDRaster hilltops = ChanNetwork.ExtractRidges(FlowInfo);

    //get hilltop curvature using filter to remove positive curvatures
    LSDRaster cht_raster = FilledDEM.get_hilltop_curvature(Surfaces[3], hilltops);
    LSDRaster CHT = FilledDEM.remove_positive_hilltop_curvature(cht_raster);

    cout << "Starting hilltop flow routing\n" << endl;

    // these params do not need changed during normal use of the HFR algorithm
    bool print_paths_switch = false;
    int thinning = 1;
    string trace_path = "";
    bool basin_filter_switch = false;
    vector<int> Target_Basin_Vector;

    //run HFR
    vector< Array2D<float> > HFR_Arrays = FlowInfo.HilltopFlowRouting(FilledDEM, hilltops, slope, StreamNetwork, aspect, Filename_variable, Basin_Raster, Surfaces[4], print_paths_switch, thinning, trace_path, basin_filter_switch, Target_Basin_Vector);

    LSDRaster HFR_LH = hilltops.LSDRasterTemplate(HFR_Arrays[1]);
    LSDRaster relief = hilltops.LSDRasterTemplate(HFR_Arrays[3]);

    //Filter Relief and LH data to remove any values < 2 pixels, as in Grieve et al (2015)
    LSDRaster LH1 = HFR_LH.RemoveBelow(2.0*Resolutions[a]);
    LSDRaster Relief1 = relief.RemoveBelow(2.0*Resolutions[a]);

    //Filter Relief and LH data to remove any values > 10000
    //These are created due to using 1m channel heads on low res data, and inadvertantly
    //sampling nodata values, which gives us massive relief values
    LSDRaster LH = LH1.RemoveAbove(10000);
    LSDRaster Relief = Relief1.RemoveAbove(10000);

    //go through the lh raster and get every value into a 1D vector
    vector<float> LH_vec = Flatten_Without_Nodata(LH.get_RasterData(), LH.get_NoDataValue());
    vector<float> Boxplot = BoxPlot(LH_vec);

    //write the values to the output file
    WriteLHData << Resolutions[a] << " " << Boxplot[0] << " " << Boxplot[1] << " " << Boxplot[2] << " " << Boxplot[3] << " " << Boxplot[4] << " " << Boxplot[5] << " " << Boxplot[6] << " " << Boxplot[7];

    WriteLHData << endl;

    stringstream ss3;
    ss3 << Path << Prefix << "_" << Resolutions[a] << "_Hist_LH_variable.txt";
    print_histogram(LH_vec, 1, ss3.str());

    //go through the relief raster and get every value into a 1D vector
    vector<float> R_vec = Flatten_Without_Nodata(Relief.get_RasterData(), Relief.get_NoDataValue());
    vector<float> Boxplot_R = BoxPlot(R_vec);

    //write the values to the output file
    WriteRData << Resolutions[a] << " " << Boxplot_R[0] << " " << Boxplot_R[1] << " " << Boxplot_R[2] << " " << Boxplot_R[3] << " " << Boxplot_R[4] << " " << Boxplot_R[5] << " " << Boxplot_R[6] << " " << Boxplot_R[7];

    WriteRData << endl;

    stringstream ss4;
    ss4 << Path << Prefix << "_" << Resolutions[a] << "_Hist_R_variable.txt";
    print_histogram(R_vec, 1, ss4.str());

    //if the user requests the rasters to be written, write the rasters
    if (WriteRasters == 1){
      cout << "Writing Rasters\n" << endl;

      Surfaces[1].write_raster((Filename + "_Slope_variable"), DEM_Format);
      CHT.write_raster((Filename + "_CHT_variable"), DEM_Format);
      LH.write_raster((Filename + "_HFR_LH_variable"), DEM_Format);
      Relief.write_raster((Filename + "_Relief_variable"), DEM_Format);

    }
  }
    WriteLHData.close();
    WriteRData.close();
}

int main(int nNumberofArgs, char *argv[])
{
  //Test for correct input arguments
	if (nNumberofArgs!=5)
	{
		cout << "FATAL ERROR: wrong number of inputs. The program needs the path (with trailing slash), the filename prefix, the DEM file format, the window size in spatial units.";
		exit(EXIT_FAILURE);
	}

  //get input args
  string path = argv[1];
  string Prefix = argv[2];
  string DEM_Format = argv[3];
  int WindowSize = atoi(argv[4]);

  //surface fitting
  vector<int> raster_selection;
  raster_selection.push_back(0);
  raster_selection.push_back(1); //slope
  raster_selection.push_back(0);
	raster_selection.push_back(1); //curvature

  //set up a writer to write the output data
  ofstream WriteData;

  //create an output filename based on the dem name
  stringstream ss;

  ss << path << Prefix << "_ChtResData.txt";
  WriteData.open(ss.str().c_str());

  //write headers
  WriteData << "resoulution 2pc 25pc median mean 75pc 98pc minimum maximum" << endl;

  //array of resolutions to load
  int Resolutions[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};

  // vectors to hold the stats about the fitted surface
  vector<float> Curv_vec;

	//set boundary conditions
  vector<string> BoundaryConditions(4, "No Flux");

  for (int a = 0; a < 10; ++a){

	  cout << "Processing DEM " << a+1 << " of " << "10" << endl;

		//load the DEM
    //build the string of the filename to load
		stringstream ss2;
		ss2 << path << Prefix << "_" << Resolutions[a] << "_DEM";
		LSDRaster DEM(ss2.str(), DEM_Format);

		//Fill
		float MinSlope = 0.0001;
		LSDRaster FilledDEM = DEM.fill(MinSlope);

	  int CurrentWindowSize = WindowSize;
	  vector<LSDRaster> Surfaces = FilledDEM.calculate_polyfit_surface_metrics(CurrentWindowSize, raster_selection);

		// get a flow info object
	  LSDFlowInfo FlowInfo(BoundaryConditions,FilledDEM);

	  //get stream net from channel heads
	  vector<int> sources = FlowInfo.Ingest_Channel_Heads((path+Prefix+"_CH"), "csv", 2);

	  LSDJunctionNetwork ChanNetwork(sources, FlowInfo);

		// extract hilltops
    LSDRaster hilltops = ChanNetwork.ExtractRidges(FlowInfo);
    LSDRaster Hilltops = ChanNetwork.ExtractHilltops(hilltops, Surfaces[1], 0.4);


    //get hilltop curvature using filter to remove positive curvatures
    LSDRaster cht_raster = FilledDEM.get_hilltop_curvature(Surfaces[3], Hilltops);
    LSDRaster CHT = FilledDEM.remove_positive_hilltop_curvature(cht_raster);

	  //go through the landscape and get every curvature value into a 1D vector
    Curv_vec = Flatten_Without_Nodata(CHT.get_RasterData(), CHT.get_NoDataValue());

    stringstream ss3;
    ss3 << path << Prefix << "_" << Resolutions[a] << "_Hist_CHT.txt";

    print_histogram(Curv_vec, 0.01, ss3.str());

    vector<float> Boxplot = BoxPlot(Curv_vec);

	  //write the values to the output file
	  WriteData << Resolutions[a] << " " << Boxplot[0] << " " << Boxplot[1] << " " << Boxplot[2] << " " << Boxplot[3] << " " << Boxplot[4] << " " << Boxplot[5] << " " << Boxplot[6] << " " << Boxplot[7];

		WriteData << endl;

  }
  WriteData.close();

}