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

//.........这里部分代码省略.........
    LSDBasin Basin(basin_junctions[w], FlowInfo, ChanNetwork);
    Basin.set_FlowLength(StreamNetwork, FlowInfo);
    Basin.set_DrainageDensity();
    Basin.set_all_HillslopeLengths(FlowInfo, HFR_LH, slope, DinfArea, log_bin_width, SplineResolution, bin_threshold);
    Basin.set_SlopeMean(FlowInfo, slope);
    Basin.set_AspectMean(FlowInfo, aspect);
    Basin.set_ElevationMean(FlowInfo, FilledDEM);
    Basin.set_ReliefMean(FlowInfo, relief);
    Basin.set_CHTMean(FlowInfo, CHT);
    Basin.set_EStar_RStar(CriticalSlope);
    
    Basins.push_back(Basin);
                             
  }
 
  //create a filestream to write the output data
  // use the input arguments to generate a path and filename for the output file
  ofstream WriteData;                 
  stringstream ss;
  ss << path << filename << "_dreich_PaperData.txt";                
  WriteData.open(ss.str().c_str());

  //write headers
  WriteData << "BasinID HFR_mean HFR_median HFR_stddev HFR_stderr HFR_Nvalues HFR_range HFR_min HFR_max SA_binned_LH SA_Spline_LH LH_Density Area Basin_Slope_mean Basin_Slope_median Basin_Slope_stddev Basin_Slope_stderr Basin_Slope_Nvalues Basin_Slope_range Basin_Slope_min Basin_Slope_max Basin_elev_mean Basin_elev_median Basin_elev_stddev Basin_elev_stderr Basin_elev_Nvalues Basin_elev_Range Basin_elev_min Basin_elev_max Aspect_mean CHT_mean CHT_median CHT_stddev CHT_stderr CHT_Nvalues CHT_range CHT_min CHT_max EStar RStar HT_Slope_mean HT_Slope_median HT_Slope_stddev HT_Slope_stderr HT_Slope_Nvalues HT_Slope_range HT_Slope_min HT_Slope_max HT_relief_mean HT_relief_median HT_relief_stddev HT_relief_stderr HT_relief_Nvalues HT_relief_range HT_relief_min HT_relief_max" << endl;

  cout << "\nWriting data to file\n" << endl;

  //write all data to the opened file, ensuring that there are data points to be written in each basin                                         
  for (int q = 0; q < int(Basins.size()); ++q){
    // only work where we have data points
    if (Basins[q].CalculateNumDataPoints(FlowInfo, HFR_LH) != 0 && Basins[q].CalculateNumDataPoints(FlowInfo, slope) != 0 && Basins[q].CalculateNumDataPoints(FlowInfo, FilledDEM) != 0 && Basins[q].CalculateNumDataPoints(FlowInfo, CHT) != 0 && Basins[q].CalculateNumDataPoints(FlowInfo, HFR_Slope) != 0 && Basins[q].CalculateNumDataPoints(FlowInfo, relief) != 0){
      
      // BasinID
      WriteData << Basins[q].get_Junction()<< " ";
      
      //HFR
      WriteData << Basins[q].get_HillslopeLength_HFR() << " " << Basins[q].CalculateBasinMedian(FlowInfo, HFR_LH) << " " << Basins[q].CalculateBasinStdDev(FlowInfo, HFR_LH) << " " << Basins[q].CalculateBasinStdError(FlowInfo, HFR_LH) << " " << Basins[q].CalculateNumDataPoints(FlowInfo, HFR_LH) << " " << Basins[q].CalculateBasinRange(FlowInfo, HFR_LH) << " " << Basins[q].CalculateBasinMin(FlowInfo, HFR_LH) << " " << Basins[q].CalculateBasinMax(FlowInfo, HFR_LH)<< " ";         
      
      //SA_Bins
      WriteData << Basins[q].get_HillslopeLength_Binned()<< " ";
      
      //SA_Spline
      WriteData << Basins[q].get_HillslopeLength_Spline()<< " ";
      
      //Density
      WriteData << Basins[q].get_HillslopeLength_Density()<< " ";
      
      //Area
      WriteData << Basins[q].get_Area()<< " ";
      
      //Slope_Basin
      WriteData << Basins[q].get_SlopeMean() << " " << Basins[q].CalculateBasinMedian(FlowInfo, slope) << " " << Basins[q].CalculateBasinStdDev(FlowInfo, slope) << " " << Basins[q].CalculateBasinStdError(FlowInfo, slope) << " " << Basins[q].CalculateNumDataPoints(FlowInfo, slope) << " " << Basins[q].CalculateBasinRange(FlowInfo, slope) << " " << Basins[q].CalculateBasinMin(FlowInfo, slope) << " " << Basins[q].CalculateBasinMax(FlowInfo, slope)<< " ";
      
      //Elev_Basin
      WriteData << Basins[q].get_ElevationMean() << " " << Basins[q].CalculateBasinMedian(FlowInfo, FilledDEM) << " " << Basins[q].CalculateBasinStdDev(FlowInfo, FilledDEM) << " " << Basins[q].CalculateBasinStdError(FlowInfo, FilledDEM) << " " << Basins[q].CalculateNumDataPoints(FlowInfo, FilledDEM) << " " << Basins[q].CalculateBasinRange(FlowInfo, FilledDEM) << " " << Basins[q].CalculateBasinMin(FlowInfo, FilledDEM) << " " << Basins[q].CalculateBasinMax(FlowInfo, FilledDEM)<< " ";
      
      //Aspect_Basin
      WriteData << Basins[q].get_AspectMean()<< " ";
      
      //CHT
      WriteData << Basins[q].get_CHTMean() << " " << Basins[q].CalculateBasinMedian(FlowInfo, CHT) << " " << Basins[q].CalculateBasinStdDev(FlowInfo, CHT) << " " << Basins[q].CalculateBasinStdError(FlowInfo, CHT) << " " << Basins[q].CalculateNumDataPoints(FlowInfo, CHT) << " " << Basins[q].CalculateBasinRange(FlowInfo, CHT) << " " << Basins[q].CalculateBasinMin(FlowInfo, CHT) << " " << Basins[q].CalculateBasinMax(FlowInfo, CHT) << " ";
      
      //EStar
      WriteData << Basins[q].get_EStar()<< " ";
      
      //RStar
      WriteData << Basins[q].get_RStar()<< " ";
      
      //Slope_mean
      WriteData << Basins[q].CalculateBasinMean(FlowInfo, HFR_Slope) << " " << Basins[q].CalculateBasinMedian(FlowInfo, HFR_Slope) << " " << Basins[q].CalculateBasinStdDev(FlowInfo, HFR_Slope) << " " << Basins[q].CalculateBasinStdError(FlowInfo, HFR_Slope) << " " << Basins[q].CalculateNumDataPoints(FlowInfo, HFR_Slope) << " " << Basins[q].CalculateBasinRange(FlowInfo, HFR_Slope) << " " << Basins[q].CalculateBasinMin(FlowInfo, HFR_Slope) << " " << Basins[q].CalculateBasinMax(FlowInfo, HFR_Slope)<< " ";
      
      //Relief_mean
      WriteData << Basins[q].get_ReliefMean() << " " << Basins[q].CalculateBasinMedian(FlowInfo, relief) << " " << Basins[q].CalculateBasinStdDev(FlowInfo, relief) << " " << Basins[q].CalculateBasinStdError(FlowInfo, relief) << " " << Basins[q].CalculateNumDataPoints(FlowInfo, relief) << " " << Basins[q].CalculateBasinRange(FlowInfo, relief) << " " << Basins[q].CalculateBasinMin(FlowInfo, relief) << " " << Basins[q].CalculateBasinMax(FlowInfo, relief)<< "\n";
  
    }    
  } 
   
  // close the output file
  WriteData.close();

  //if the user requests the raster to be written, write the rasters
  if (WriteRasters == 1){
    cout << "Writing Rasters\n" << endl;                                   
   // FilledDEM.write_raster((path+filename+"_Fill"), "flt");
    Surfaces[1].write_raster((path+filename+"_dreich_Slope"),"flt");
    //Surfaces[2].write_raster((path+filename+"_Aspect"),"flt");
    //Surfaces[3].write_raster((path+filename+"_Curvature"),"flt");
    //StreamNetwork.write_raster((path+filename+"_STNET"), "flt"); 
    //Basin_Raster.write_raster((path+filename+"_Basins"), "flt"); 
    CHT.write_raster((path+filename+"_dreich_CHT"),"flt");
    HFR_LH.write_raster((path+filename+"_dreich_HFR_LH"),"flt"); 
    HFR_Slope.write_raster((path+filename+"_dreich_HFR_SLP"),"flt");
    relief.write_raster((path+filename+"_dreich_Relief"),"flt");
    
    //perform a hillshade
    //LSDRaster Hillshade = FilledDEM.hillshade(45.0,315.0,1.0);
    //Hillshade.write_raster((path+filename+"_HS"),"flt");
  
  }
}

//.........这里部分代码省略.........
       << "minimum_segment_length: " << minimum_segment_length << endl
       << "sigma: " << sigma << endl
       << "start_movern " << start_movern << endl
       << "d_movern: " << d_movern << endl
       << "n_movern: " << n_movern << endl
       << "target_nodes: " << target_nodes << endl
       << "n_iterarions: " << n_iterations << endl
       << "fraction_dchi_for_variation: " << fraction_dchi_for_variation << endl
       << "vertical interval: " << vertical_interval << endl
       << "horizontal interval: " << horizontal_interval << endl
       << "area thinning fraction for SA analysis: " << area_thin_frac << endl
       << "target_skip is: " << target_skip << endl;
  
  cout << "=================================================================" << endl
       << "You are also loading a precipitation raster to get the discharge" << endl
       << "The precipitation raster should: " << endl
       << "1) Have units of m/yr" << endl
       << "Have the same prefix as the DEM, but will have the extension _PRECIP"<< endl
       << "For example, if the DEM is my_DEM.bil, the precip file will be" << endl
       << "my_DEM_precip.bil" << endl
       << "=================================================================" << endl;


  string jn_name = itoa(junction_number);
  string uscore = "_";
  jn_name = uscore+jn_name;
  file_info_in.close();

  string DEM_f_name = path_name+DEM_name+fill_ext;
  string DEM_bil_extension = "bil";
  
  string precip_ext = "_PRECIP";
  string Precip_f_name = path_name+DEM_name+precip_ext;

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

  // load the filled DEM
  LSDRaster filled_topo_test((path_name+DEM_name+fill_ext), DEM_bil_extension);

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

  // calcualte the distance from outlet
  LSDRaster DistanceFromOutlet = FlowInfo.distance_from_outlet();
  LSDIndexRaster ContributingPixels = FlowInfo.write_NContributingNodes_to_LSDIndexRaster();

  // calculate the discharge
  // note: not discharge yet, need to multiply by cell area
  LSDRaster VolumePrecipitation(Precip_f_name, DEM_bil_extension);
  float dx = VolumePrecipitation.get_DataResolution();
  
  // volume precipitation per time precipitation times the cell areas
  VolumePrecipitation.raster_multiplier(dx*dx);
  
  // discharge accumulates this precipitation
  LSDRaster Discharge = FlowInfo.upslope_variable_accumulator(VolumePrecipitation);

  string Q_ext = "_Q";
  string Q_f_name = path_name+DEM_name+Q_ext;
  Discharge.write_raster(Q_f_name,DEM_bil_extension);

  // get the sources
  vector<int> sources;
  sources = FlowInfo.get_sources_index_threshold(ContributingPixels, threshold);

  // now get the junction network
  LSDJunctionNetwork ChanNetwork(sources, FlowInfo);

  // now get a junction and look for the longest channel upstream
  cout << "creating main stem" << endl;
  LSDIndexChannel main_stem = ChanNetwork.generate_longest_index_channel_in_basin(junction_number, FlowInfo, DistanceFromOutlet);
        cout << "got main stem channel, with n_nodes " << main_stem.get_n_nodes_in_channel() <<  endl;

  string Basin_name = "_basin";
  LSDIndexRaster BasinArray = ChanNetwork.extract_basin_from_junction(junction_number,junction_number,FlowInfo);
  BasinArray.write_raster((path_name+DEM_name+Basin_name+jn_name),DEM_bil_extension);

  // now get the best fit m over n for all the tributaries
  int organization_switch = 1;
  int pruning_switch = 1;
  LSDIndexChannelTree ChannelTree(FlowInfo, ChanNetwork, junction_number, organization_switch,
                                        DistanceFromOutlet, pruning_switch, pruning_threshold);

  LSDRaster DrainageArea = FlowInfo.write_DrainageArea_to_LSDRaster();

  // print a file that can be ingested by the chi fitting algorithm
  string Chan_fname = "_ChanNet";
  string Chan_ext = ".chan";
  string Chan_for_chi_ingestion_fname = path_name+DEM_name+Chan_fname+jn_name+Chan_ext;
  ChannelTree.print_LSDChannels_for_chi_network_ingestion(FlowInfo,
                             filled_topo_test, DistanceFromOutlet, Chan_for_chi_ingestion_fname,
                             Discharge);
  ChannelTree.convert_chan_file_for_ArcMap_ingestion(Chan_for_chi_ingestion_fname,DrainageArea,Discharge);

}

//.........这里部分代码省略.........
    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[])
{
    // the driver version
    string driver_version = "Driver_version: 0.01";

    // some paramters
    //Test for correct input arguments
    if (nNumberofArgs!=5)
    {
        cout << endl;
        cout << "=====================================================================" << endl;
        cout << "|| Welcome to the Topographic Shielding tool!                      ||" << endl;
        cout << "|| This program is used to calculate topographic shielding_driver  ||" << endl;
        cout << "|| rasters following the method of Codilean (2006).                ||" << endl;
        cout << "=====================================================================" << endl;
        cout << "This program requires four inputs: " << endl;
        cout << "* First the path to the DEM files." << endl;
        cout << "   The path must have a slash at the end." << endl;
        cout << "  (Either \\ or / depending on your operating system.)" << endl;
        cout << "* Second the prefix of the DEM (that is, without the .bil)." << endl;
        cout << "   For example if you DEM is Ladakh.bil, you should enter Ladakh" << endl;
        cout << "   (Note that the DEM should be in *.bil format)" << endl;
        cout << "* Third the increment in azimuth (in degrees) over which you" << endl;
        cout << "   want shielding calculated. Recommended values is 5" << endl;
        cout << "* Fourth the increment in inclination (in degrees) over which you" << endl;
        cout << "   want shielding calculated. Recommended values is 5" << endl;
        cout << "=====================================================================" << endl;
        cout << "For more documentation, see readme and online documentation" << endl;
        cout << "=====================================================================" << endl;
        cout << endl;
        exit(EXIT_SUCCESS);
    }

    cout << endl;
    cout << "===============================================================" << endl;
    cout << "Welcome to the Topographic Shielding tool" << endl;
    cout << "This software was developed at the University of Edinburgh," << endl;
    cout << "by the Land Surface Dynamics group. For questions email" << endl;
    cout << "simon.m.mudd _at_ ed.ac.uk" << endl;
    cout << "This software is released under a GNU public license." << endl;
    cout << "You are using " << driver_version << endl;
    cout << "================================================================" << endl;
    cout << "++IMPORTANT++ The DEM must be an ENVI bil format file" << endl;
    cout << "ENVI bil files are required because, unlike asc or flt files, " << endl;
    cout << "they use georeferencing information." << endl;
    cout << "For more information about changing DEM formatting, see: " << endl;
    cout << "http://lsdtopotools.github.io/LSDTT_book/#_gdal_2" << endl;
    cout << "================================================================" << endl;


    //Assign values from input arguments
    string path_name = argv[1];
    string file_name = argv[2];
    string azimuth_str = argv[3];
    string inclination_str = argv[4];

    //Set the DEM filename
    file_name = path_name+file_name;

    //Convert the angle increments to integers
    int azimuth_step =  atoi(azimuth_str.c_str());
    int inclination_step =  atoi(inclination_str.c_str());

    // check the parameter values
    check_azimuth_and_inclination(azimuth_step,inclination_step);

    //load dem
    LSDRaster DEM(file_name, "bil");

    //launch toposhielding
    LSDRaster TopoShielding = DEM.TopographicShielding(azimuth_step,inclination_step);
    TopoShielding.write_raster(file_name+"_TopoShield","bil");

    cout << "Done!" << endl << endl;
}