Updating to version 9. See history for changes

LSDtopotools · Jun 23, 2023 · f3ffbe8 · f3ffbe8
1 parent 2540659
commit f3ffbe8
Show file tree

Hide file tree

Showing 13 changed files with 1,315 additions and 261 deletions.
diff --git a/history.asc b/history.asc
@@ -35,19 +35,27 @@
 
 == 0.7 (2022-09-13)
 
-* This includes a number of improvements to the valley width extraction algorithms including those used in the paper https://doi.org/10.5194/esurf-10-437-2022
-* This release has a new dependency on the OpenCV library. The programs will compile without this library but the valley metrics component will only compile with OpenCV.
-* This version removed dependency on the Point Cloud Library.
-* It also includes updated functionality for the cosmogenic tools, so the analyses from the CAIRN (Mudd et al 2016) method are streamlined.
+* This includes a number of improvements to the valley width extraction algorithms including those used in the paper https://doi.org/10.5194/esurf-10-437-2022 (SMM and FJC)
+* This release has a new dependency on the OpenCV library. The programs will compile without this library but the valley metrics component will only compile with OpenCV. (FJC)
+* This version removed dependency on the Point Cloud Library. (FJC)
+* It also includes updated functionality for the cosmogenic tools, so the analyses from the CAIRN (Mudd et al 2016) method are streamlined. (SMM)
 
 == 0.8 (2023-05-03)
 
-* New contributions from SMM and FJC
-* Additions to channel extraction include getting longest channel, including more information in the output files, and an option to print channels as shapefiles. 
-* A number of changes to the cosmo tool including better accounting for a column, an interface to bring in transient erosion rates, better support for CRN, smoother operation of the nesting functions, and fixes to the shielding calculations that allow you to compute erosion rates in one step. 
-* Added a function for creating steady state fluvial landscapes
-* Included some new metrics, most notably the normalised concavity index or NCI that was used in Chen et al Nature paper. 
-* Added a routine that facilitates tagging of major drainage divides.
-* Added perona-malik filter to hillslope analysis
-* added a small function to analyse channel tips for a channel extraction routine
-* Other minor bug and typo fixes. 
+* Additions to channel extraction include getting longest channel, including more information in the output files, and an option to print channels as shapefiles. (SMM)
+* A number of changes to the cosmo tool including better accounting for a column, an interface to bring in transient erosion rates, better support for CRN, smoother operation of the nesting functions, and fixes to the shielding calculations that allow you to compute erosion rates in one step. (SMM)
+* Added a function for creating steady state fluvial landscapes (SMM)
+* Included some new metrics, most notably the normalised concavity index or NCI that was used in Chen et al Nature paper. (SMM)
+* Added a routine that facilitates tagging of major drainage divides. (SMM)
+* Added perona-malik filter to hillslope analysis (FJC)
+* added a small function to analyse channel tips for a channel extraction routine (SMM)
+* Other minor bug and typo fixes. (SMM and FJC)
+
+== 0.9 (2023-06-23)
+
+* Added function in lsdtt-basic-metrics to extract basin outlines (SMM)
+* Added function in lsdtt-basic-metrics to get the bearing of channel segments (SMM)
+* Added function in lsdtt-basic-metrics to print flow direction codes in ArcMap format (SMM)
+* More testing of lsdtt-valley-metrics to reduce bugs and give more informative error messages. (SMM)
+* Added swath output to the terrace routine (SMM)
+* Some testing of ridge extraction and minor bug fixes (FJC)
diff --git a/src/LSDBasin.cpp b/src/LSDBasin.cpp
@@ -873,7 +873,7 @@ void LSDBasin::set_Perimeter(LSDFlowInfo& FlowInfo)
   {
 
     FlowInfo.retrieve_current_row_and_col(BasinNodes[q], i, j);
-      BasinData[i][j] = BasinNodes[q];
+    BasinData[i][j] = BasinNodes[q];
 
   }
 
@@ -884,40 +884,38 @@ void LSDBasin::set_Perimeter(LSDFlowInfo& FlowInfo)
     NDVCount = 0;
     // cout << i << "||" << j << endl;
 
-      if (i != 0 && j != 0 && i<NRows-1 && j < NCols-1)
-      {
-        //count border cells that are NDV
-
-        if (BasinData[i-1][j-1] == NoDataValue){ ++NDVCount; }
-        if (BasinData[i][j-1] == NoDataValue){ ++NDVCount; }
-        if (BasinData[i+1][j-1] == NoDataValue){ ++NDVCount; }
-        if (BasinData[i-1][j] == NoDataValue){ ++NDVCount; }
-        if (BasinData[i+1][j] == NoDataValue){ ++NDVCount; }
-        if (BasinData[i-1][j+1] == NoDataValue){ ++NDVCount; }
-        if (BasinData[i][j+1] == NoDataValue){ ++NDVCount; }
-        if (BasinData[i+1][j+1] == NoDataValue){ ++NDVCount; }
-
-        if (NDVCount >= 1 && NDVCount < 8)
-        {  //increase the first value to get a simpler polygon  (changed to 1 by FJC 23/03/15 to get only internal hilltops.
-          //edge pixel                       // Otherwise not all external ridges were being excluded from the analysis).
-          I.push_back(i);
-          J.push_back(j);
-          B.push_back(BasinNodes[q]);
-        }
-      }
-      else
-      {
-        ++i;
+    if (i != 0 && j != 0 && i<NRows-1 && j < NCols-1)
+    {
+      //count border cells that are NDV
+
+      if (BasinData[i-1][j-1] == NoDataValue){ ++NDVCount; }
+      if (BasinData[i][j-1] == NoDataValue){ ++NDVCount; }
+      if (BasinData[i+1][j-1] == NoDataValue){ ++NDVCount; }
+      if (BasinData[i-1][j] == NoDataValue){ ++NDVCount; }
+      if (BasinData[i+1][j] == NoDataValue){ ++NDVCount; }
+      if (BasinData[i-1][j+1] == NoDataValue){ ++NDVCount; }
+      if (BasinData[i][j+1] == NoDataValue){ ++NDVCount; }
+      if (BasinData[i+1][j+1] == NoDataValue){ ++NDVCount; }
+
+      if (NDVCount >= 1 && NDVCount < 8)
+      {  //increase the first value to get a simpler polygon  (changed to 1 by FJC 23/03/15 to get only internal hilltops.
+        //edge pixel                       // Otherwise not all external ridges were being excluded from the analysis).
+        I.push_back(i);
+        J.push_back(j);
+        B.push_back(BasinNodes[q]);
       }
+    }
+    else
+    {
+      ++i;
+    }
 
   }
 
   //now have 2 vectors of i and j indexes of every point
   Perimeter_i = I;
   Perimeter_j = J;
   Perimeter_nodes = B;
-
-
 }
 
 //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
@@ -967,10 +965,6 @@ void LSDBasin::print_perimeter_hypsometry_to_csv(LSDFlowInfo& FlowInfo, string p
   set_Perimeter(FlowInfo);
   clean_perimeter(FlowInfo);
 
-  // TESTING
-  string perim_test = "/home/clubb/Data_for_papers/drainage_capture/Santa_Cruz/FUCK_THIS.csv";
-  print_perimeter_to_csv(FlowInfo, perim_test);
-
   //  open the file
   ofstream perim_out;
   perim_out.open(perimeter_fname.c_str());
@@ -1100,10 +1094,10 @@ vector<int> LSDBasin::order_perimeter_nodes(LSDFlowInfo& FlowInfo)
     next_j = Outlet_j+1;
     cout << "The closest node is to the SE" << endl;
   }
-  else if (PerimeterNodes[Outlet_i-1][Outlet_j+1] == 1) // southwest
+  else if (PerimeterNodes[Outlet_i+1][Outlet_j-1] == 1) // southwest
   {
-    next_i = Outlet_i-1;
-    next_j = Outlet_j+1;
+    next_i = Outlet_i+1;
+    next_j = Outlet_j-1;
     cout << "The closest node is to the SW" << endl;
   }
   else
@@ -1357,6 +1351,209 @@ vector<int> LSDBasin::order_perimeter_nodes(LSDFlowInfo& FlowInfo)
   return sorted_nodes;
 }
 
+
+
+
+
+
+
+
+//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+// Order perimeter nodes from the outlet
+// Must CLEAN the perimeter first using clean_perimeter function
+// FJC 16/01/18
+//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+vector<int> LSDBasin::order_perimeter_nodes_SMM(LSDFlowInfo& FlowInfo)
+{
+  // first clean the perimeter nodes
+  //clean_perimeter(FlowInfo);
+  cout << "Ordering the perimeter nodes..." << endl;
+  vector<int> sorted_nodes;
+  Array2D<int> PerimeterNodes(NRows,NCols,0);
+
+  // first get the outlet node
+  int outlet_node = get_Outlet_node();
+  int outlet_row, outlet_col;
+  FlowInfo.retrieve_current_row_and_col(outlet_node, outlet_row, outlet_col);
+  Outlet_i = outlet_row;
+  Outlet_j = outlet_col;
+  cout << "The outlet node is: " << outlet_node << endl;
+
+  // get an array of perimeter nodes
+  for (int i = 0; i < int(Perimeter_nodes.size()); i++)
+  {
+    int this_row, this_col;
+    FlowInfo.retrieve_current_row_and_col(Perimeter_nodes[i], this_row, this_col);
+    PerimeterNodes[this_row][this_col] = 1;
+  }
+  // The outlet always needs to be in the perimeter
+  // The outlet is already added so it gets an increment to 2 
+  PerimeterNodes[Outlet_i][Outlet_j] = 2;
+  cout << "Got the array of perimeter nodes" << endl;
+  cout << "The outlet row and column are: " << Outlet_i << "," << Outlet_j << endl;
+
+  // push back the outlet node, node 0
+  sorted_nodes.push_back(outlet_node);
+  int next_i, next_j;
+  bool first_node = true; // bool to check if you're at the first node. This makes a difference because if you're at the first node we don't want to go back to the outlet.
+
+
+  // N, S, E, and W will always be the shortest distances, so do these first
+  if (PerimeterNodes[Outlet_i][Outlet_j-1] == 1) // West
+  {
+    next_i = Outlet_i;
+    next_j = Outlet_j-1;
+    //cout << "The closest node is to the west" << endl;
+  }
+  else if (PerimeterNodes[Outlet_i-1][Outlet_j] == 1) // North
+  {
+    next_i = Outlet_i-1;
+    next_j = Outlet_j;
+    //cout << "The closest node is to the north" << endl;
+  }
+  else if (PerimeterNodes[Outlet_i][Outlet_j+1] == 1)  // east
+  {
+    next_i = Outlet_i;
+    next_j = Outlet_j+1;
+    //cout << "The closest node is to the east" << endl;
+  }
+  else if (PerimeterNodes[Outlet_i+1][Outlet_j] == 1)  // south
+  {
+    next_i = Outlet_i+1;
+    next_j = Outlet_j;
+    //cout << "The closest node is to the south" << endl;
+  }
+  else if (PerimeterNodes[Outlet_i-1][Outlet_j-1] == 1) // northwest
+  {
+    next_i = Outlet_i-1;
+    next_j = Outlet_j-1;
+    //cout << "The closest node is to the NW" << endl;
+  }
+  else if (PerimeterNodes[Outlet_i-1][Outlet_j+1] == 1) // northeast
+  {
+    next_i = Outlet_i-1;
+    next_j = Outlet_j+1;
+    //cout << "The closest node is to the NE" << endl;
+  }
+  else if (PerimeterNodes[Outlet_i+1][Outlet_j+1] == 1) // southeast
+  {
+    next_i = Outlet_i+1;
+    next_j = Outlet_j+1;
+    //cout << "The closest node is to the SE" << endl;
+  }
+  else if (PerimeterNodes[Outlet_i+1][Outlet_j-1] == 1) // southwest
+  {
+    next_i = Outlet_i+1;
+    next_j = Outlet_j-1;
+    //cout << "The closest node is to the SW" << endl;
+  }
+  else
+  {
+    //cout << "None of these were perimeter nodes, oops" << endl;
+  }
+
+  cout << "The first node is at: " << next_i << "," << next_j << endl;
+
+  // push back the next node to the sorted node vector
+  int next_node = FlowInfo.retrieve_node_from_row_and_column(next_i, next_j);
+  sorted_nodes.push_back(next_node);
+
+  bool reached_outlet = false;
+  int this_i, this_j;
+  // now start at the outlet node and find the nearest perimeter node.
+  while (reached_outlet == false)
+  {
+    // start at the next node and find the one with the closest distance that
+    // hasn't already been visited
+    this_i = next_i;
+    this_j = next_j;
+    PerimeterNodes[this_i][this_j] = 2;
+
+    // N, S, E, and W will always be the shortest distances, so do these first
+    if (PerimeterNodes[this_i][this_j-1] == 1) // West
+    {
+      next_i = this_i;
+      next_j = this_j-1;
+      //cout << "The closest node is to the west" << endl;
+    }
+    else if (PerimeterNodes[this_i-1][this_j] == 1) // North
+    {
+      next_i = this_i-1;
+      next_j = this_j;
+      //cout << "The closest node is to the north" << endl;
+    }
+    else if (PerimeterNodes[this_i][this_j+1] == 1)  // east
+    {
+      next_i = this_i;
+      next_j = this_j+1;
+      //cout << "The closest node is to the east" << endl;
+    }
+    else if (PerimeterNodes[this_i+1][this_j] == 1)  // south
+    {
+      next_i = this_i+1;
+      next_j = this_j;
+      //cout << "The closest node is to the south" << endl;
+    }
+    else if (PerimeterNodes[this_i-1][this_j-1] == 1) // northwest
+    {
+      next_i = this_i-1;
+      next_j = this_j-1;
+      //cout << "The closest node is to the NW" << endl;
+    }
+    else if (PerimeterNodes[this_i-1][this_j+1] == 1) // northeast
+    {
+      next_i = this_i-1;
+      next_j = this_j+1;
+      //cout << "The closest node is to the NE" << endl;
+    }
+    else if (PerimeterNodes[this_i+1][this_j+1] == 1) // southeast
+    {
+      next_i = this_i+1;
+      next_j = this_j+1;
+      //cout << "The closest node is to the SE" << endl;
+    }
+    else if (PerimeterNodes[this_i+1][this_j-1] == 1) // southwest
+    {
+      next_i = this_i+1;
+      next_j = this_j-1;
+      //cout << "The closest node is to the SW" << endl;
+    }
+    else
+    {
+      reached_outlet = true;
+    }
+
+    //cout << "The next node is at: " << next_i << "," << next_j << endl;
+
+    if (!reached_outlet)
+    {
+      next_node = FlowInfo.retrieve_node_from_row_and_column(next_i, next_j);
+      sorted_nodes.push_back(next_node);
+    }
+  }
+
+  return sorted_nodes;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
 //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
 // Set the four different hillslope length measurements for the basin.
@@ -2527,7 +2724,7 @@ void LSDBasin::organise_perimeter(LSDFlowInfo& flowpy)
 
   // TESTING FUNCTION, DELETE IT AFTERWARDS BORIS!!!!
   Perimeter_nodes = Perimeter_nodes_sorted;
-    print_perimeter_to_csv(flowpy, "/home/boris/Desktop/LSD/capture/sorbas/peritest_AFTER.csv");
+  print_perimeter_to_csv(flowpy, "/home/boris/Desktop/LSD/capture/sorbas/peritest_AFTER.csv");
 
 
 

diff --git a/src/LSDBasin.hpp b/src/LSDBasin.hpp
@@ -464,6 +464,12 @@ class LSDBasin
   /// @date 16/01/18
   vector<int> order_perimeter_nodes(LSDFlowInfo& FlowInfo);
 
+  /// @brief Orders perimeter nodes from the outlet. Uses a different approach
+  /// @param FlowInfo the LSDFlowInfo object
+  /// @author SMM
+  /// @date 29/05/23
+  vector<int> order_perimeter_nodes_SMM(LSDFlowInfo& FlowInfo);
+
   /// @brief Set the four different hillslope length measurements for the basin.
   /// @param FlowInfo Flowinfo object.
   /// @param HillslopeLengths LSDRaster of hillslope lengths from the hilltop flow routing method.