Updated a number of functions including improved hilltop flow routing…

…, a swath profile tool, and fixed a bug in geojson printing

src/LSDFlowInfo.cpp

@@ -1306,60 +1306,99 @@ LSDRaster LSDFlowInfo::find_nodes_not_influenced_by_edge_draining_to_nodelist(ve
   // The node list is loaded from another function. This can be a channel network or a
   // main stream channel. That channel can be generated using the LSDFlowInfo method LSDFlowInfo::get_flow_path
   int n_nodes = int(node_list.size());
+  bool problem_node = false;        // for debugging
   cout << "I am trying to isolate nodes draining to a fixed channel that are not influenced by nodata." << endl;
   cout << "I am programmed in a rather rudimentary, brute force way. It may take some time. " << endl;
   cout << "I need to process " << n_nodes << " nodes." << endl;
+
+  if(problem_node)
+  {
+    cout << "I am in debugging mode! To switch me back you need to go into FlowInfo line 1319" << endl;
+  }
+
   for(int node = 0; node < n_nodes; node++)
   {
-    if (node % 500 == 0)
+    if ( node_list[node] != NoDataValue)
     {
-      cout << "Node " << node << " of " << n_nodes << endl;
-    }
+      if (node % 500 == 0)
+      {
+        cout << "Node " << node << " of " << n_nodes << endl;
+      }
 
-    // first get the current node
-    retrieve_current_row_and_col(node_list[node],row,col);
-    new_topography_data[row][col] = topography.get_data_element(row,col);
+      // first get the current node
+      retrieve_current_row_and_col(node_list[node],row,col);
+      new_topography_data[row][col] = topography.get_data_element(row,col);
 
-    // now get the donors
-    vector<int> donors = retrieve_donors_to_node(  node_list[node]);
 
-    // Loop through all the donors to this particular node. We will look through the donors to
-    // find upslope influence of nodata.
-    int n_donors = int(donors.size());
-    for (int donor = 0; donor<n_donors; donor++)
-    {
-      int this_d_node = donors[donor];
 
-      // check to see if the donor is in the fixed channel. If it is, get rid of it.
-      bool donor_in_node_list = false;
-      for (int nl = 0; nl<n_nodes; nl++)
+
+      // now get the donors
+      vector<int> donors = retrieve_donors_to_node(  node_list[node]);
+      int n_donors = int(donors.size());
+
+      if( problem_node)
       {
-        if (node_list[nl] == this_d_node)
-        {
-          donor_in_node_list = true;
-        }
+        cout << "=====================================================================" << endl;
+        cout << "Inspecting problem nodes: " << row << ", " << col << "," << new_topography_data[row][col] << endl;
       }
 
-      if (donor_in_node_list == false)
+      // Loop through all the donors to this particular node. We will look through the donors to
+      // find upslope influence of nodata.
+      for (int donor = 0; donor<n_donors; donor++)
       {
-        bool is_influenced = is_upstream_influenced_by_nodata(this_d_node, topography);
+        int this_d_node = donors[donor];
+
+        if (problem_node)
+        {
+          retrieve_current_row_and_col(this_d_node,row,col);
+          cout << "donor " << this_d_node << " r: " << row << ", " << col << " ";
+        }
 
-        // If this donor is not influenced by nodata anywhere upstream, add all the pixels
-        // draining to this node to the data.
-        if (is_influenced == false)
+        // check to see if the donor is in the fixed channel. If it is, get rid of it.
+        bool donor_in_node_list = false;
+        for (int nl = 0; nl<n_nodes; nl++)
         {
-          vector<int> influenced_nodes = get_upslope_nodes_include_outlet(this_d_node);
-          int n_in_basin = int(influenced_nodes.size());
-          for (int in = 0; in < n_in_basin; in++)
+          if (node_list[nl] == this_d_node)
           {
-            current_node = influenced_nodes[in];
-            retrieve_current_row_and_col(current_node,row,col);
-            new_topography_data[row][col] = topography.get_data_element(row,col);
+            donor_in_node_list = true;
+
+            if (problem_node)
+            {
+              cout << "This donor is in the channel" << endl;
+            }
           }
         }
-      }
 
+        if (donor_in_node_list == false)
+        {
+          bool is_influenced = is_upstream_influenced_by_nodata(this_d_node, topography);
+
+          if (problem_node)
+          {
+            cout << "This donor is not in the channel, Is it influenced by nodata? " << is_influenced << endl;
+          }
+
+          // If this donor is not influenced by nodata anywhere upstream, add all the pixels
+          // draining to this node to the data.
+          if (is_influenced == false)
+          {
+            vector<int> influenced_nodes = get_upslope_nodes_include_outlet(this_d_node);
+            int n_in_basin = int(influenced_nodes.size());
+            for (int in = 0; in < n_in_basin; in++)
+            {
+              current_node = influenced_nodes[in];
+              retrieve_current_row_and_col(current_node,row,col);
+              new_topography_data[row][col] = topography.get_data_element(row,col);
+            }
+          }
+        }
+      }
     }
+    else
+    {
+      cout << "This nodeindex has a nodata value" << endl;
+    }
+
   }
 
   LSDRaster temp_topo(NRows,NCols,XMinimum,YMinimum,DataResolution,NoDataValue,new_topography_data,GeoReferencingStrings);
@@ -2257,6 +2296,7 @@ vector<int> LSDFlowInfo::Ingest_Channel_Heads(string filename, int input_switch)
   return Sources;
 }
 
+
 //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
 // Method to ingest sources from OS MasterMap Water Network Layer (csv)
 // into a vector of source nodes so that an LSDJunctionNetwork can be created easily
@@ -4111,12 +4151,13 @@ vector<int> LSDFlowInfo::sort_node_list_based_on_raster(vector<int> node_vec, LS
 int LSDFlowInfo::get_node_index_of_coordinate_point(float X_coordinate, float Y_coordinate)
 {
   // Shift origin to that of dataset
-  float X_coordinate_shifted_origin = X_coordinate - XMinimum - DataResolution*0.5;
-  float Y_coordinate_shifted_origin = Y_coordinate - YMinimum - DataResolution*0.5;
+  float X_coordinate_shifted_origin = X_coordinate - XMinimum;
+  float Y_coordinate_shifted_origin = Y_coordinate - YMinimum;
 
   // Get row and column of point
   int col_point = int(X_coordinate_shifted_origin/DataResolution);
-  int row_point = (NRows-1) - int(ceil(Y_coordinate_shifted_origin/DataResolution));
+  int row_point = (NRows - 1) - int(ceil(Y_coordinate_shifted_origin/DataResolution)-0.5);
+
 
   // Get node of point
   int CurrentNode;

src/LSDFlowInfo.hpp

@@ -699,7 +699,7 @@ class LSDFlowInfo
 
   ///@brief This function tests whether one node is upstream of another node
   ///@param current_node
-  ///@param test_node
+  ///@param test_node The node which we want to test whether it is upstream or not
   ///@return Boolean indicating whether node is upstream or not
   ///@author FC
   ///@date 08/10/13

src/LSDGeometry.cpp

@@ -87,6 +87,7 @@
 using namespace std;
 using namespace TNT;
 
+
 //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
 // Creates an LSDGeometry from an LSDRaster
 //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

src/LSDGeometry.hpp

@@ -88,7 +88,7 @@ using namespace std;
 using namespace TNT;
 
 
-/// @brief This object packages a number of tools for chi analysis
+/// @brief This object has some routines for dealing with points and lines
 class LSDGeometry
 {
   public:
@@ -398,7 +398,7 @@ class LSDPolyline: public LSDGeometry
     vector<int> get_flowinfo_nodes_of_line(LSDRasterInfo& RI, LSDFlowInfo& FlowInfo);
 
 /*
-    /// @detail THis traces to the next pixel ensureing no nodes are missed
+    /// @detail THis traces to the next pixel ensuring no nodes are missed
     void trace_to_next_pixel(LSDRasterInfo& RI, double StartEasting,double StartNorthing, int start_row, int start_col,
                          int end_row, int end_col, int current_row, int current_col,
                          double slope,

src/LSDRaster.cpp

@@ -10420,8 +10420,180 @@ vector<LSDRaster> LSDRaster::get_nearest_distance_and_value_masks()
 }
 
 
+
+//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+/// This takes a list of points. Then, for every pixel in the 
+///  raster it finds the point amongst that list that is closest to the given pixel. 
+//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
+vector< LSDRaster > LSDRaster::find_nearest_point_from_list_of_points(vector<float> Eastings, vector<float> Northings,
+                                              vector<float> values, float swath_width)
+{
+  Array2D<float> Distances(NRows, NCols, NoDataValue);
+  Array2D<float> Values(NRows, NCols, NoDataValue);
+  Array2D<float> Nodes(NRows, NCols, NoDataValue);
+
+  // Get the easting and northing vectors from the raster
+  vector<float> raster_eastings;
+  vector<float> raster_northings;
+  get_easting_and_northing_vectors(raster_eastings, raster_northings);
+
+  int n_nodes = int(Eastings.size());
+
+  float this_min_distance, this_value, this_node, this_distance; 
+  float x_0,y_0,x_1,y_1; 
+
+  // Now we loop through all pixels in the raster
+  cout << "I am going to find the closest point in a list of points to every pixel in your raster." << endl;
+  cout << "This has not been programmed in an efficient way," << endl;
+  cout << "so if this is a big DEM it will take a while" << endl;
+  for(int row = 0; row<NRows; row++)
+  {
+    for (int col = 0; col<NCols; col++)
+    {
+      if(RasterData[row][col] != NoDataValue)
+      {
+
+        // reset the minimum distance
+        this_min_distance = 10000000000;
+        this_value = NoDataValue;
+        this_node = NoDataValue;
+
+        // now loop through all the points getting the minimum point
+        x_0 = raster_eastings[col];
+        y_0 = raster_northings[row];
+
+        for(int i = 0; i< n_nodes; i++)
+        {
+          x_1 = Eastings[i];
+          y_1 = Northings[i];
+
+          this_distance = distance_between_two_points(x_0,y_0,x_1,y_1);
+
+          if (this_distance < this_min_distance)
+          {
+            this_min_distance = this_distance;
+            this_value = values[i];
+            this_node = float(i);
+          }
+        }
+
+        //cout << "r: " << row << ", c: "<< col << " min dist: " << this_min_distance << endl;
+
+        // Finished looping through node, update the data 
+        // only record the nodes that are withing the swath width
+        if (this_min_distance<0.5*swath_width)
+        {
+          Distances[row][col] = this_min_distance;
+          Values[row][col] = this_value;
+          Nodes[row][col] = this_node;
+        }
+      }
+    }
+  }
+
+  LSDRaster DistRaster(NRows,NCols, XMinimum, YMinimum, DataResolution,
+                            NoDataValue, Distances, GeoReferencingStrings);
+  LSDRaster ValuesRaster(NRows,NCols, XMinimum, YMinimum, DataResolution,
+                            NoDataValue, Values, GeoReferencingStrings);                            
+  LSDRaster NodesRaster(NRows,NCols, XMinimum, YMinimum, DataResolution,
+                            NoDataValue, Nodes, GeoReferencingStrings);
+
+  vector<LSDRaster> return_rasters;
+  return_rasters.push_back(DistRaster);
+  return_rasters.push_back(ValuesRaster);
+  return_rasters.push_back(NodesRaster);
+
+  return return_rasters;
+
+}
+
+
+void LSDRaster::make_swath(vector<float> Eastings, vector<float> Northings,
+                                              vector<float> values, float swath_width, float bin_width,
+                                              string swath_data_prefix, bool print_swath_rasters)
+{
+  // first we get the swath rasters
+  vector< LSDRaster > swath_rasther_vec = find_nearest_point_from_list_of_points(Eastings, Northings,values, swath_width);
+
+  string swath_data_name = swath_data_prefix+"_swath.csv";
+  if (print_swath_rasters)
+  {
+    cout << "I am now going to print your swath rasters." << endl;
+    string dist_fname = swath_data_prefix+"_swathdist";
+    string val_fname = swath_data_prefix+"_swathval";
+    string node_fname = swath_data_prefix+"_swathnode";
+
+    swath_rasther_vec[0].write_raster(dist_fname,"bil");
+    swath_rasther_vec[1].write_raster(val_fname,"bil");
+    swath_rasther_vec[2].write_raster(node_fname,"bil");
+  }
+
+  // now we need to vectorize the data
+  vector<float> distance_vector;
+  vector<float> elevation_vector;
+  float this_distance;
+  for(int row = 0; row<NRows; row++)
+  {
+    for(int col = 0; col< NCols; col++)
+    {
+      this_distance = swath_rasther_vec[1].get_data_element(row,col);
+      if(this_distance != NoDataValue)
+      {
+        distance_vector.push_back(this_distance);
+        elevation_vector.push_back(RasterData[row][col]);
+      }
+    }
+  }
+
+
+
+
+
+  // the million data vectors that go into the binning
+  vector<float> midpoints_output;
+  vector<float>  MeanX_output;
+  vector<float>  MedianX_output;
+  vector<float> StandardDeviationX_output;
+  vector<float> StandardErrorX_output;
+  vector<float> MADX_output;
+  vector<float> MeanY_output;
+  vector<float> MinimumY_output;
+  vector<float> FirstQuartileY_output;
+  vector<float> MedianY_output;
+  vector<float> ThirdQuartileY_output;
+  vector<float> MaximumY_output;
+  vector<float> StandardDeviationY_output;
+  vector<float> StandardErrorY_output;
+  vector<float> MADY_output;
+  vector<int> number_observations_output;
+
+
+  bin_data(distance_vector, elevation_vector, bin_width, midpoints_output, MeanX_output,
+           MedianX_output, StandardDeviationX_output, StandardErrorX_output, MADX_output,
+           MeanY_output, MinimumY_output, FirstQuartileY_output, MedianY_output,
+           ThirdQuartileY_output, MaximumY_output, StandardDeviationY_output, StandardErrorY_output,
+           MADY_output, number_observations_output, NoDataValue);
+
+
+  ofstream bin_data_out(swath_data_name);
+  bin_data_out << "distance,mean_elevation,minimum_z,first_quartile_z,median_z,third_quartile_z,max_z,n_samples"<< endl;
+  int n_bins = (midpoints_output.size());
+  // For some reason the last bin is always empty so I just don't print it
+  for(int i = 0; i<n_bins-1; i++)
+  {
+    bin_data_out << midpoints_output[i] << "," << MeanY_output[i] << "," << MinimumY_output[i] <<","
+                  << FirstQuartileY_output[i] << "," <<  MedianY_output[i] << "," << ThirdQuartileY_output[i] 
+                  << "," << MaximumY_output[i] << "," << number_observations_output[i] << endl;
+  }
+  bin_data_out.close();
+
+}
+
+
+
+
 //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
-// THis pads a smaller index raster to the same extent as a bigger raster by adding
+// This pads a smaller index raster to the same extent as a bigger raster by adding
 // no data values
 //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
 LSDIndexRaster LSDRaster::PadSmallerRaster(LSDIndexRaster& smaller_raster)