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

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=---=-=-=-=-=-=-=-=-=-=-=-=-
// This function prints a chan file for assimilation into the chi analysis, 
// but in this cases uses a discharge rather than a drainage area
//
// the file format is
// channel_number node_index node_on_reciever row column flow_dist elevation drainage_area
//
// SMM 07/05/2015
//
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
void LSDIndexChannelTree::print_LSDChannels_for_chi_network_ingestion(LSDFlowInfo& FlowInfo,
                             LSDRaster& Elevation_Raster, LSDRaster& FlowDistance, string fname,
                             LSDRaster& Discharge)
{
  if (organization_switch != 1)
  {
    cout << "LSDIndexChannelTree you can't run LSDIndexChannelTree::retrieve_LSDChannels_from_tree" << endl;
    cout << "with this channel organization, organization switch: " << organization_switch << endl;
    exit(EXIT_FAILURE);
  }

  // open the outfile
  ofstream channelfile_out;
  channelfile_out.open(fname.c_str());

  channelfile_out.precision(10);

  float m_over_n = 0.5;
  float A_0 = 1;

  // get the vector of channels
  vector<LSDChannel> vector_of_channels = retrieve_LSDChannels_from_tree(m_over_n, A_0, FlowInfo,Elevation_Raster);

  int n_channels = vector_of_channels.size();
  int n_nodes_in_channel;
  int node,row,col;
  float elev,chi,drain_area,flow_dist,this_discharge;

  // first print out some data about the dem
  channelfile_out << get_NRows() << endl;
  channelfile_out << get_NCols() << endl;
  channelfile_out << get_XMinimum() << endl;
  channelfile_out << get_YMinimum() << endl;
  channelfile_out << get_DataResolution() << endl;
  channelfile_out << get_NoDataValue() << endl;

  //loop through the channels
  for (int i = 0; i< n_channels; i++)
  {
    // get the number of nodes in the channel
    n_nodes_in_channel =IndexChannelVector[i].get_n_nodes_in_channel();

    // now loop through the channel, printing out the data.
    for(int ch_node= 0; ch_node<n_nodes_in_channel; ch_node++)
    {
      IndexChannelVector[i].get_node_row_col_in_channel(ch_node, node, row, col);
      vector_of_channels[i].retrieve_node_information(ch_node, elev, chi, drain_area);
      flow_dist = FlowDistance.get_data_element(row,col);
      this_discharge = Discharge.get_data_element(row,col);

      // print data to file
      channelfile_out << i << " " << receiver_channel[i] << " " << node_on_receiver_channel[i] << " "
                      << node << " " << row << " " << col << " " << flow_dist << " "
                      << " " << elev << " " << this_discharge << endl;
    }
  }

  channelfile_out.close();

}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
// Calculate w, a hyrdological index, used in the factor of safety equation.
// Call with the ratio of recharge to transmissivity.
// SWDG 13/6/16
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
LSDSoilHydroRaster LSDSoilHydroRaster::Calculate_w(LSDRaster& Slope, LSDRaster& DrainageArea){

  Array2D<float> w(NRows, NCols, NoDataValue);

  for (int i = 1; i < NRows - 1; ++i){
    for (int j = 1; j < NCols - 1; ++j){

      if (RasterData[i][j] != NoDataValue){

        float value = RasterData[i][j] * (DrainageArea.get_data_element(i,j)/sin(Slope.get_data_element(i,j)));
        if (value < 1.0){
          w[i][j] = value;
        }
        else{
          w[i][j] = 1.0;
        }

      }
    }
  }

  LSDSoilHydroRaster output(NRows,NCols,XMinimum,YMinimum,DataResolution,NoDataValue,w,GeoReferencingStrings);
  return output;

}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
// This is an incredibly rudimentary function used to modify landslide raster
// It takes a few rasters from the
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
void LSDSoilHydroRaster::NaiveLandslide(LSDRaster& FilledElevation, int initiationPixels,
                                      int MinPixels, float landslide_thickness)
{
  // Get a flow info object
  // set no flux boundary conditions
  vector<string> boundary_conditions(4);
  boundary_conditions[0] = "No";
  boundary_conditions[1] = "no flux";
  boundary_conditions[2] = "no flux";
  boundary_conditions[3] = "No flux";


  // some values from the rasters
  float local_elev;
  float local_mask;

  // get a flow info object
  LSDFlowInfo FlowInfo(boundary_conditions,FilledElevation);

  // get the contributing pixels
  LSDIndexRaster ContributingPixels = FlowInfo.write_NContributingNodes_to_LSDIndexRaster();
  vector<int> sources = FlowInfo.get_sources_index_threshold(ContributingPixels, initiationPixels);

  // get a value vector for the landslides
  vector<float> landslide_thicknesses;
  for (int i = 0; i< int(sources.size()); i++)
  {
    landslide_thicknesses.push_back(landslide_thickness);
  }


  // get the mask
  LSDRaster Mask = FlowInfo.get_upslope_node_mask(sources,landslide_thicknesses);

  // now set all points that have elevation data but not landslide data to
  // the value of the landslide thickness, removing data that is below the minium
  // pixel area
  for (int row = 0; row<NRows; row++)
  {
    for (int col = 0; col<NCols; col++)
    {
      local_elev =  FilledElevation.get_data_element(row,col);
      local_mask =  Mask.get_data_element(row,col);

      RasterData[row][col] = local_mask;

      // Turn nodata points into 0s
      if( local_mask == NoDataValue)
      {
        RasterData[row][col] = 0.0;
      }

      // remove data where there is no topographic information
      if( local_elev == NoDataValue)
      {
        RasterData[row][col] = NoDataValue;
      }
    }
  }
}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=---=-=-=-=-=-=-=-=-=-=-=-=-
// Same as above but also reports the discharge
//
// SMM 06/05/2015
//
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
void LSDIndexChannelTree::convert_chan_file_for_ArcMap_ingestion(string fname, LSDRaster& DrainageArea, LSDRaster& Discharge)
{

  // open the outfile
  ifstream channelfile_in;
  channelfile_in.open(fname.c_str());

  unsigned dot = fname.find_last_of(".");

  string prefix = fname.substr(0,dot);
  //string suffix = str.substr(dot);
    string insert = "_for_Arc.csv";
    string outfname = prefix+insert;

    cout << "the Arc channel filename is: " << outfname << endl;

    ofstream ArcChan_out;
    ArcChan_out.open(outfname.c_str());
    ArcChan_out.precision(10);

  // print the first line of the arcchan. This is going to be comma seperated!
  ArcChan_out << "id,x,y,channel,reciever_channel,node_on_reciever_channel,node,row,col,flow_distance_m,elevation_m,drainage_area_m2,discharge_m2_times_precipunits" << endl;

  // now go throught the file, collecting the data
  int id,ch,rc,norc,n,r,c;
  float fd,elev,da;
  float x,y;

  float xll;
  float yll;
  float datares;
  float ndv;
  int nrows;
  int ncols;
  float this_discharge;
  float this_da;

  // read in the first lines with DEM information
  channelfile_in >> nrows >> ncols >> xll >> yll >> datares >> ndv;
  id = 0;

  // now loop through the file, calculating x and y locations as you go
  while(channelfile_in >> ch >> rc >> norc >> n >> r >> c >> fd >> elev >> da)
  {
    id++;
    x = xll + float(c)*datares + 0.5*datares;
    y = yll + float(nrows-r)*datares - 0.5*datares;		// this is because the DEM starts from the top corner

    this_discharge = Discharge.get_data_element(r,c);
    this_da = DrainageArea.get_data_element(r,c);

    ArcChan_out << id << "," << x << "," << y << "," << ch << "," << rc 
                << "," << norc << "," << n << "," << r << "," << c << ","
                << fd << "," << elev << "," << this_da << "," << this_discharge << endl;
  }

  channelfile_in.close();
  ArcChan_out.close();
}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
// 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();
}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
// Calculate the the sinmap Stability Index (SI).
// This is a wrapper around a lightly modified port of the original sinmap 2.0 implementation.
// call with any SoilHydroRaster, it's values are used for identification of NoDataValues.
//
// SWDG 15/6/16
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
LSDSoilHydroRaster LSDSoilHydroRaster::Calculate_sinmap_SI(LSDRaster Slope, LSDRaster DrainageArea, LSDSoilHydroRaster lo_C, LSDSoilHydroRaster hi_C, LSDSoilHydroRaster lo_phi, LSDSoilHydroRaster hi_phi, LSDSoilHydroRaster lo_RoverT, LSDSoilHydroRaster hi_RoverT, LSDSoilHydroRaster lo_r, LSDSoilHydroRaster hi_r, LSDSoilHydroRaster lo_FS, LSDSoilHydroRaster hi_FS){

  Array2D<float> SI(NRows, NCols, NoDataValue);

  for (int i = 1; i < NRows - 1; ++i){
    for (int j = 1; j < NCols - 1; ++j){

      if (RasterData[i][j] != NoDataValue){
        SI[i][j] = StabilityIndex(Slope.get_data_element(i, j), DrainageArea.get_data_element(i, j), lo_C.get_data_element(i, j), hi_C.get_data_element(i, j), lo_phi.get_data_element(i, j), hi_phi.get_data_element(i, j), lo_RoverT.get_data_element(i, j), hi_RoverT.get_data_element(i, j), lo_r.get_data_element(i, j), hi_r.get_data_element(i, j), lo_FS.get_data_element(i, j), hi_FS.get_data_element(i, j));
      }
    }
  }

  LSDSoilHydroRaster output(NRows,NCols,XMinimum,YMinimum,DataResolution,NoDataValue,SI,GeoReferencingStrings);
  return output;

}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
// This function calculates a snow thickenss (effective, in g cm^-2 for cosmogenic
// applications)  based on a bilinear model such as that of P Kirchner: http://escholarship.org/uc/item/9zn1c1mk#page-8
// The paper is here: http://www.hydrol-earth-syst-sci.net/18/4261/2014/hess-18-4261-2014.html
// This paper also agrees withy this general trend:
// http://www.the-cryosphere.net/8/2381/2014/tc-8-2381-2014.pdf
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
void LSDSoilHydroRaster::SetSnowEffDepthBilinear(float SlopeAscend, float SlopeDescend,
                          float PeakElevation, float PeakSnowpack, LSDRaster& Elevation)
{
  float LocalElevation;
  float ascendEffDepth;
  float descendEffDepth;
  float thisEffDepth = 0;

  for (int row = 0; row <NRows; row++)
  {
    for (int col = 0; col<NCols; col++)
    {
      if (RasterData[row][col] != NoDataValue)
      {
        LocalElevation = Elevation.get_data_element(row,col);

        if (LocalElevation != NoDataValue)
        {
          // get the effective depth on both the ascending and descending limb
          ascendEffDepth = SlopeAscend*(LocalElevation-PeakElevation)+PeakSnowpack;
          descendEffDepth = SlopeDescend*(LocalElevation-PeakElevation)+PeakSnowpack;

          // the correct depth is the lesser of the two
          if (ascendEffDepth < descendEffDepth)
          {
            thisEffDepth =  ascendEffDepth;
          }
          else
          {
            thisEffDepth = descendEffDepth;
          }

          // if the depth is less than zero, then set to zero
          if(thisEffDepth <0)
          {
            thisEffDepth = 0;
          }

          RasterData[row][col] =  thisEffDepth;

        }
        else        // if there ins't any elevation data, set the snow data to NoData
        {
          RasterData[row][col] = NoDataValue;
        }
      }
    }
  }
}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
// This function calculates a snow thickenss (effective, in g cm^-2 for cosmogenic
// applications)  based on a richard's equsion sigmoidal growth model
// It was propoesd to represent peak SWE so we cruedly apply it to average annual SWE
// see
// http://onlinelibrary.wiley.com/doi/10.1002/2015GL063413/epdf
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
void LSDSoilHydroRaster::SetSnowEffDepthRichards(float MaximumEffDepth, float MaximumSlope, float v,
                          float lambda, LSDRaster& Elevation)
{
  // Don't let V be less than or equal to zero
  if (v <= 0)
  {
    v = 0.001;
  }

  // some variables to speed up compuation
  float exp_term;
  float thisEffDepth = 0;
  float elev_mulitplier = (MaximumSlope/MaximumEffDepth)*pow((1+v),1+(1/v));
  float LocalElevation;

  for (int row = 0; row <NRows; row++)
  {
    for (int col = 0; col<NCols; col++)
    {
      if (RasterData[row][col] != NoDataValue)
      {
        LocalElevation = Elevation.get_data_element(row,col);

        if (LocalElevation != NoDataValue)
        {
          // get the effective depth using the richards sigmoidal gorth function
          exp_term = 1+v*exp(elev_mulitplier*(lambda-LocalElevation));
          thisEffDepth = MaximumEffDepth*pow(exp_term,-(1/v));

          // if the depth is less than zero, then set to zero
          if(thisEffDepth <0)
          {
            thisEffDepth = 0;
          }

          // update the data
          RasterData[row][col] =  thisEffDepth;

        }
        else        // if there ins't any elevation data, set the snow data to NoData
        {
          RasterData[row][col] = NoDataValue;
        }
      }
    }
  }
}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
// Calculate the factor of safety using the sinmap definition.
// Call with the dimensionless cohesion (C) raster.
// SWDG 13/6/16
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
LSDSoilHydroRaster LSDSoilHydroRaster::Calculate_sinmap_Fs(LSDRaster& Slope, LSDSoilHydroRaster& w, LSDSoilHydroRaster& r, LSDSoilHydroRaster& phi){

  Array2D<float> Fs(NRows, NCols, NoDataValue);

  for (int i = 1; i < NRows - 1; ++i){
    for (int j = 1; j < NCols - 1; ++j){

      if (RasterData[i][j] != NoDataValue){
        Fs[i][j] = ( RasterData[i][j] + cos(Slope.get_data_element(i,j)) * (1.0-w.get_data_element(i,j)*r.get_data_element(i,j)) * tan(phi.get_data_element(i,j)) ) / sin(Slope.get_data_element(i,j));
      }
    }
  }

  LSDSoilHydroRaster output(NRows,NCols,XMinimum,YMinimum,DataResolution,NoDataValue,Fs,GeoReferencingStrings);
  return output;

}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
// Calculate h, the soil depth normal to the slope, used in the factor of safety equation.
// Call with the soil thickness raster.
// SWDG 13/6/16
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
LSDSoilHydroRaster LSDSoilHydroRaster::Calculate_h(LSDRaster& Slope){

  Array2D<float> h(NRows, NCols, NoDataValue);

  for (int i = 1; i < NRows - 1; ++i){
    for (int j = 1; j < NCols - 1; ++j){

      if (RasterData[i][j] != NoDataValue){
        h[i][j] = RasterData[i][j] * cos(Slope.get_data_element(i,j));

      }
    }
  }

  LSDSoilHydroRaster output(NRows,NCols,XMinimum,YMinimum,DataResolution,NoDataValue,h,GeoReferencingStrings);
  return output;

}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
// this function prints chi values. It is used on the channel tree when channels are organized by links
// (that is organization switch == 0
//
// SMM 01/09/2012
//
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
void LSDIndexChannelTree::print_chi_vs_elevation_from_channel_tree(LSDRaster& Elevation, LSDFlowInfo& FlowInfo, LSDJunctionNetwork& ChannelNetwork,
                                                     float m_over_n, float A_0, string chi_vs_elev_fname)
{

	if(organization_switch != 0)
	{
		cout << "LSDIndexChannelTree you can't run LSDIndexChannelTree::print_chi_vs_elevation_from_channel_tree with this channel organization" << endl;
		exit(EXIT_FAILURE);
	}

	int n_channels = IndexChannelVector.size();
	int n_nodes_in_link,current_node_index;
	vector< vector<float> > chi_vectors = calculate_chi_from_channel_tree(FlowInfo, ChannelNetwork,
                                                     m_over_n, A_0);

	// the chi iterator
  	vector< vector<float> >::iterator chi_iter;
    chi_iter =  chi_vectors.begin();
	float current_chi;
	float current_elev;
	int curr_row,curr_col;

	ofstream chi_elev_out;
	chi_elev_out.open(chi_vs_elev_fname.c_str());

	for (int i = 0; i<n_channels; i++)
	{
		n_nodes_in_link = IndexChannelVector[i].get_n_nodes_in_channel();

	  	for (int this_node = 0; this_node<n_nodes_in_link; this_node++)
		{
			current_node_index = IndexChannelVector[i].get_node_in_channel(this_node);
			current_chi  =  (*chi_iter)[this_node];
			FlowInfo.retrieve_current_row_and_col(current_node_index,curr_row,curr_col);
			current_elev = Elevation.get_data_element(curr_row,curr_col);

			//cout << "current_node: " << current_node_index << " and chi: " << current_chi << endl;
			//chi_elev_out << current_chi << " " << current_elev << endl;
		}
		chi_iter++;
	}
	chi_elev_out.close();

}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=---=-=-=-=-=-=-=-=-=-=-=-=-
// this function prints chi and elevation, along with flow distance and the 
// number of the tributary it all goes to one file
//
// the file format is
// channel_number node_index row column flow_dist chi elevation drainage_area
//
// SMM 01/09/2012
//
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
void LSDIndexChannelTree::print_LSDChannels_from_tree(float m_over_n, float A_0, LSDFlowInfo& FlowInfo,
                             LSDRaster& Elevation_Raster, LSDRaster& FlowDistance, string fname)
{
  if (organization_switch != 1)
  {
    cout << "LSDIndexChannelTree you can't run LSDIndexChannelTree::retrieve_LSDChannels_from_tree" << endl;
    cout << "with this channel organization, organization switch: " << organization_switch << endl;
    exit(EXIT_FAILURE);
  }

  // open the outfile
  ofstream channelfile_out;
  channelfile_out.open(fname.c_str());

  // get the vector of channels
  vector<LSDChannel> vector_of_channels = retrieve_LSDChannels_from_tree(m_over_n, A_0, FlowInfo,Elevation_Raster);

  int n_channels = vector_of_channels.size();
  int n_nodes_in_channel;
  int node,row,col;
  float elev,chi,drain_area,flow_dist;
  //loop through the channels
  for (int i = 0; i< n_channels; i++)
  {
    // get the number of nodes in the channel
    n_nodes_in_channel =IndexChannelVector[i].get_n_nodes_in_channel();

    // now loop through the channel, printing out the data.
    for(int ch_node= 0; ch_node<n_nodes_in_channel; ch_node++)
    {
      IndexChannelVector[i].get_node_row_col_in_channel(ch_node, node, row, col);
      vector_of_channels[i].retrieve_node_information(ch_node, elev, chi, drain_area);
      flow_dist = FlowDistance.get_data_element(row,col);

      // print data to file
      channelfile_out << i << " " << node << " " << row << " " << col << " " << flow_dist << " "
                      << chi << " " << elev << " " << drain_area << endl;
    }
  }

  channelfile_out.close();

}

//.........这里部分代码省略.........
  vector< vector<int> > chi_node_indices;
  
  // these are for the individual channels
  vector<float> these_chi_m_means;
  vector<float> these_chi_b_means;
  vector<float> these_chi_coordinates;
  vector<int> these_chi_node_indices;
  
  // these are maps that will store the data
  map<int,float> m_means_map;
  map<int,float> b_means_map;
  map<int,float> chi_coord_map;
  map<int,float> elev_map;
  map<int,float> area_map;
  map<int,float> flow_distance_map;
  vector<int> node_sequence_vec;

  // these are for working with the FlowInfo object
  int this_node,row,col;

  // get the number of channels
  int n_channels = int(source_nodes.size());
  
  for(int chan = 0; chan<n_channels; chan++)
  {
    cout << "Sampling channel " << chan+1 << " of " << n_channels << endl;
    
    // get this particualr channel (it is a chi network with only one channel)
    LSDChiNetwork ThisChiChannel(FlowInfo, source_nodes[chan], outlet_nodes[chan], 
                                Elevation, FlowDistance, DrainageArea,chi_coordinate);
    
    // split the channel
    //cout << "Splitting channels" << endl;
    ThisChiChannel.split_all_channels(A_0, m_over_n, n_iterations, skip, target_nodes, minimum_segment_length, sigma);
    
    // monte carlo sample all channels
    //cout << "Entering the monte carlo sampling" << endl;
    ThisChiChannel.monte_carlo_sample_river_network_for_best_fit_after_breaks(A_0, m_over_n, n_iterations, skip, minimum_segment_length, sigma);
  
    // okay the ChiNetwork has all the data about the m vales at this stage. 
    // Get these vales and print them to a raster
    chi_m_means = ThisChiChannel.get_m_means();
    chi_b_means = ThisChiChannel.get_b_means();
    chi_coordinates = ThisChiChannel.get_chis();
    chi_node_indices = ThisChiChannel.get_node_indices();
    
    // now get the number of channels. This should be 1!
    int n_channels = int(chi_m_means.size());
    if (n_channels != 1)
    {
      cout << "Whoa there, I am trying to make a chi map but something seems to have gone wrong with the channel extraction."  << endl;
      cout << "I should only have one channel per look but I have " << n_channels << " channels." << endl;
    }
    
    // now get the m_means out
    these_chi_m_means = chi_m_means[0];
    these_chi_b_means = chi_b_means[0];
    these_chi_coordinates = chi_coordinates[0];
    these_chi_node_indices = chi_node_indices[0];
    
    //cout << "I have " << these_chi_m_means.size() << " nodes." << endl;
    
    int n_nodes_in_channel = int(these_chi_m_means.size());
    for (int node = 0; node< n_nodes_in_channel; node++)
    {
      
      this_node =  these_chi_node_indices[node];
      //cout << "This node is " << this_node << endl;
      
      // only take the nodes that have not been found
      if (m_means_map.find(this_node) == m_means_map.end() )
      {
        FlowInfo.retrieve_current_row_and_col(this_node,row,col);
      
        //cout << "This is a new node; " << this_node << endl;
        m_means_map[this_node] = these_chi_m_means[node];
        b_means_map[this_node] = these_chi_b_means[node];
        chi_coord_map[this_node] = these_chi_coordinates[node];
        elev_map[this_node] = Elevation.get_data_element(row,col);
        area_map[this_node] = DrainageArea.get_data_element(row,col);
        flow_distance_map[this_node] = FlowDistance.get_data_element(row,col);
        node_sequence_vec.push_back(this_node);
      }
      else
      {
        //cout << "I already have node: " << this_node << endl;
      }
    }
  }
  
  // set the opject data members
  M_chi_data_map =m_means_map; 
  b_chi_data_map = b_means_map;
  elev_data_map = elev_map;
  chi_data_map = chi_coord_map;
  flow_distance_data_map = flow_distance_map;
  drainage_area_data_map = area_map;
  node_sequence = node_sequence_vec;
  
}

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
// This function is a much more rudimentary version that mimics the
// channel steepness caluclations.
// chi needs tobe calculated outside of the function 
//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
void LSDChiTools::chi_map_automator_rudimentary(LSDFlowInfo& FlowInfo, 
                                    vector<int> source_nodes,
                                    vector<int> outlet_nodes,
                                    LSDRaster& Elevation, LSDRaster& FlowDistance, 
                                    LSDRaster& DrainageArea, LSDRaster& chi_coordinate, 
                                    int regression_nodes)
{

  // the data is stored in maps, for easier testing if a node has been
  // visited. 
  // You might consider having these as data elements in the object so you don't 
  // have to pass them
  map<int,float> gradient_data_map;
  map<int,float> intercept_data_map;
  map<int,float> R2_data_map;
  map<int,float> chi_coordinate_data_map;
  map<int,float> elevation_data_map;
  map<int,float> flow_distance_map;
  map<int,float> area_map;
  vector<int> node_order;
  
  // check if the number of nodes are odd .If not add 1
  if (regression_nodes % 2 == 0)
  {
    cout << "Hello user. You need an odd number of regression nodes." << endl;
    regression_nodes = regression_nodes+1;
    cout << " Changing your regression nodes to " << regression_nodes << endl;
  }
  
  // now get the midpoint
  int mp_nodes = (regression_nodes-1)/2;
  
  //cout << "The number of mp nodes is: " << mp_nodes << endl;
  
  // these keep track of the beginning and ending nodes of a given channel
  int channel_start_node;
  int channel_end_node;
  float channel_end_elevation;
  
  // vectors for holding the chi elevation data
  vector<float> chi_vec;
  vector<float> elev_vec;
  vector<float> empty_vec;
  
  // these are extracted from the channel segments using linear regression
  float intercept,gradient,R_squared;
  
  //float this_chi;
  //float this_elev;
  int this_mp_node;
  int this_end_node;
  int this_start_node;
  
  // these are for getting information out of the FlowInfo object
  int row,col, this_node;
  int r_node, r_row,r_col;          // reciever row and column. 

  // The way this works is that it starts at the top of a channel. It then works 
  // its way down and find the node that is the midpoint and the node that is the
  // end point. The midpoint node is where the data will be recorded.
  // It then puts the data from the start node to the end node into a vector 
  // and performs a linear regression of this vector. The regression data from these
  // vectors are recorded at the nodes.
  // We then want to cover all the nodes with data so what happens if some nodes
  // do not become midpoints? 
  // We could start at the top and get the first midpoint. 
  // From there we can work our way down checking if the top of the regression segment
  // is more than one node down from the end point...

  // get the number of channels
  int n_channels = int(source_nodes.size());
  // now loop through the channels
  for(int chan = 0; chan<n_channels; chan++)
  {
    channel_start_node = source_nodes[chan];
    channel_end_node = outlet_nodes[chan];
    
    // Get the elevation of the end node as a secondary check of the ending of the channel
    // segment
    FlowInfo.retrieve_current_row_and_col(channel_end_node,row,col);
    channel_end_elevation = Elevation.get_data_element(row,col);
    
    // reset the flag for ending the channel
    bool is_end_of_channel = false;

    // set the segment start node to the channel start node
    this_start_node = channel_start_node;

    // now retrieve the midpoint node
    this_node = channel_start_node;
    for(int n = 0; n<mp_nodes; n++)
    {
      FlowInfo.retrieve_receiver_information(this_node,r_node,r_row,r_col);
      this_node = r_node;
    }
//.........这里部分代码省略.........