Update to version 0.8

src/CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 2.6 FATAL_ERROR)
+cmake_minimum_required(VERSION 2.8.12 FATAL_ERROR)
 project(LSDTopoTools)
 
 # Dependencies
@@ -38,6 +38,7 @@ if(OpenCV_FOUND)
     LSDMostLikelyPartitionsFinder.hpp
     LSDParameterParser.hpp
     LSDParticle.hpp
+    LSDParticleColumn.hpp
     LSDRaster.hpp
     LSDRasterAggregator.hpp
     LSDRasterInfo.hpp
@@ -106,6 +107,7 @@ else()
     LSDMostLikelyPartitionsFinder.hpp
     LSDParameterParser.hpp
     LSDParticle.hpp
+    LSDParticleColumn.hpp
     LSDRaster.hpp
     LSDRasterAggregator.hpp
     LSDRasterInfo.hpp
@@ -152,4 +154,4 @@ else()
     ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
   )
   endif()
-
+

src/LSDBasin.cpp

@@ -1849,7 +1849,7 @@ void LSDBasin::add_basin_to_LSDIndexRaster(LSDIndexRaster& basin_raster,
   if (RNR != NRows || RNC != NCols)
   {
     cout << "LSDBasin::Add_basin_to_LSDIndexRaster, failed to add basins" << endl
-         << "since LSDIndexRaster does not have same georeferncing as basin." << endl;
+         << "since LSDIndexRaster does not have same georeferencing as basin." << endl;
   }
   else
   {
@@ -1864,16 +1864,16 @@ void LSDBasin::add_basin_to_LSDIndexRaster(LSDIndexRaster& basin_raster,
       }
       else
       {
-        // if there is already data there, the smaller basin overrwites the larger
+        // if there is already data there, the smaller basin overwrites the larger
         // basin. You will be able to see the larger basins since the smaller
         // basins are nested within the larger basins
         existing_basin_index = basin_raster.get_data_element(row,col);
 
         if(drainage_of_other_basins.find(existing_basin_index) == drainage_of_other_basins.end())
         {
-          cout << "LSDBasin::Add_basin_to_LSDIndexRaster, Something has gone wrong." << endl
-               << "The existing basin is not in the basin map. " << endl
-               << "Overwriting data." << endl;
+          //cout << "LSDBasin::Add_basin_to_LSDIndexRaster, Something has gone wrong." << endl
+          //     << "The existing basin is not in the basin map. " << endl
+          //     << "Overwriting data." << endl;
           basin_raster.set_data_element(row,col,this_basin_index);
         }
         else if (this_basin_pixel_area < drainage_of_other_basins[existing_basin_index])

src/LSDCRNParameters.cpp

@@ -271,8 +271,8 @@ void LSDCRNParameters::load_parameters_for_atmospheric_scaling(string path_to_da
 //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
 
 //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
-// This function sets a number of paramters that are used for acalucaltion
-// of scaling and error propigation
+// This function sets a number of paramters that are used for calculation
+// of scaling and error propagation
 //
 // They are constants used in the CRONUS caluclator, and have been ported
 // from make_al_be_consts_v22.m written by Greg Balco
@@ -1152,6 +1152,10 @@ void LSDCRNParameters::set_BraucherBorchers_parameters()
 {
   //S_t = 1;
 
+  // Notes:
+  // The F values for 10Be, 26Al and 36Cl come from Braucher et al EPSL 309 1-9 (2011)
+  //
+
   // from Vermeesh 2007
   // 10Be from Chmeleff/Korschinek 10Be decay constant;
   lambda_10Be = 500e-9;    // in yr-1
@@ -1161,13 +1165,19 @@ void LSDCRNParameters::set_BraucherBorchers_parameters()
 
   // from the Braucher and Borchers papers
   // data compiled by Mirjam Schaller in personal communication
-  //
-  // All but 10Be are calibrated to the Stone scaling
+  // for 10Be and 26Al
+  // These come from dividing the Borchers et al 2016 stone spallation
+  // number by F0, the spallation production from the Braucher 2011 paper
+  // I did this myself for 36C
+  // All are calibrated to the Stone scaling
+  // The 14C uses the lupker spallation scaling and the Borchers stone spallation
+  // 12.24/0.788
   // Also linked to the nishizumii standards
+  // 
   P0_10Be = 4.061;          // in a/g/yr
   P0_26Al = 28.851;         // in a/g/yr
-  P0_14C = 15.21;          // in a/g/yr
-  P0_36Cl = 58.95;         // in a/g/yr
+  P0_14C = 15.533;          // in a/g/yr
+  P0_36Cl = 55.705;         // in a/g/yr
   P0_21Ne = 18.23;         // in a/g/yr
   P0_3He = 121.59;         // in a/g/yr
 
@@ -1190,16 +1200,19 @@ void LSDCRNParameters::set_BraucherBorchers_parameters()
   F_26Al[3] = 0.0028;
 
   // dimensionless
-  F_14C[0] = 0.83;
-  F_14C[1] = 0.15;
+  // This data comes from the Lupker et al QG 28 80-87 2015 paper
+  // There are no fast muons in that paper (F[3])
+  F_14C[0] = 0.788;
+  F_14C[1] = 0.212;
   F_14C[2] = 0.0;
-  F_14C[3] = 0.02;
+  F_14C[3] = 0.0;
 
   // dimensionless
-  F_36Cl[0] = 0.9456;
-  F_36Cl[1] = 0.0324;
+  // calculated from table 5 in Borchers et al 2011
+  F_36Cl[0] = 0.9396;
+  F_36Cl[1] = 0.0369;
   F_36Cl[2] = 0.00;
-  F_36Cl[3] = 0.022;
+  F_36Cl[3] = 0.0235;
 }
 //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
 

src/LSDChiTools.cpp

@@ -2163,6 +2163,29 @@ LSDIndexRaster LSDChiTools::get_basin_raster(LSDFlowInfo& FlowInfo, LSDJunctionN
 
 }
 
+
+//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+// Get the basin key given a outlet node 
+// SMM 12/10/2022
+//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+int LSDChiTools::get_basin_key_from_outlet_node(int outlet_node)
+{
+  // This uses the baselevel_keys_map, where the key is the node index of a baselevel node and the
+  // value is the baselevel index
+  int basin_key;
+  if ( baselevel_keys_map.find(outlet_node) == baselevel_keys_map.end() )
+  {
+    basin_key = NoDataValue;
+  }
+  else
+  {
+    basin_key = baselevel_keys_map[outlet_node];
+  }
+  return basin_key;
+
+}
+
+
 //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
 // extract lithology data per basins from a litho map, a FlowInfo and a junction networks
 // ongoing work
@@ -2882,7 +2905,7 @@ LSDRaster LSDChiTools::prefilter_river_topography(LSDRaster& filled_topography,
 
     // Going through the vector of nodes
     size_t rit_a = vecnode.size() - 1;
-    int last_node; 
+    //int last_node; 
     float last_elev;
     for(vector<int>::reverse_iterator ti2 = vecnode.rbegin(); ti2 != vecnode.rend(); ++ti2)
     {
@@ -2903,7 +2926,7 @@ LSDRaster LSDChiTools::prefilter_river_topography(LSDRaster& filled_topography,
         // cout << "D:" << detrend[rit_a] << endl;
         // cout << "LE:" << last_elev << endl;
       }
-      last_node = this_node;
+      //last_node = this_node;
       last_elev = filled_topography.get_data_element(row,col);
       rit_a--;
     }
@@ -2925,7 +2948,7 @@ LSDRaster LSDChiTools::prefilter_river_topography(LSDRaster& filled_topography,
     for(int rit_b = int(vecnode.size()) - 1 ; rit_b >= 0; --rit_b)
     {
 
-      if(rit_b == vecnode.size() - 1)
+      if(rit_b == int(vecnode.size()) - 1)
       {
         // cout << "Premier node: " << this_base_elevation << endl;
         int tnode = vecnode[rit_b], trow=0,tcol=0;
@@ -4032,7 +4055,7 @@ void  LSDChiTools::TVD_on_segelev(LSDFlowInfo& Flowinfo)
 
 
         // Now I just have to save the results into global maps
-        int last_node = -9999; // forcing segfault in case it happens 
+        //int last_node = -9999; // forcing segfault in case it happens 
         for(pferd=0; pferd<this_vec.size(); pferd++)
         {
           int this_node = this_vec[pferd];
@@ -4599,7 +4622,7 @@ void LSDChiTools::print_bandwidth_ksn_knickpoint(string filename)
   file_out << "source_key,basin_key,length,chi" << endl;
 
   int this_source_key, this_basin_key, this_node = 0;
-  float this_bandwidth = 0;
+  float this_bandwidth;
   map<int,vector<int> >::iterator OL;
 
   for(OL = map_node_source_key.begin(); OL !=  map_node_source_key.end() ; OL++)
@@ -4966,9 +4989,9 @@ void LSDChiTools::generate_knickpoint_overview(LSDFlowInfo& FlowInfo, LSDRaster&
   Array2D<float> tCx_TVD_map(NRows,NCols,NoDataValue);
   Array2D<float> tMx_TVD_map(NRows,NCols,NoDataValue);
 
-  for(size_t h=0;h<NRows;h++)
+  for(int h=0;h<NRows;h++)
   {
-    for(size_t g=0;g<NCols;g++)
+    for(int g=0;g<NCols;g++)
     {
       if(Elevation.get_data_element(h,g) != NoDataValue)
       {

src/LSDChiTools.hpp

@@ -1084,6 +1084,14 @@ class LSDChiTools
     LSDIndexRaster get_basin_raster(LSDFlowInfo& FlowInfo, LSDJunctionNetwork& JunctionNetwork,
                                vector<int> Juntions);
 
+
+    /// @brief This gets the basin key given an outlet node
+    /// @param outlet_node 
+    /// @return The basin key
+    /// @author SMM
+    /// @date 12/10/2022
+    int get_basin_key_from_outlet_node(int outlet_node);
+
     /// @brief This return a map with the basin ID as key, and map of the count of the different lithology as key/value.
     /// @param FlowInfo
     /// @param JunctionNetwork