C++ LSDRaster类代码示例

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

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

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

}

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

//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
// 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;

}

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

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

}