Added catchment delineation

include/richdem/common/grid_cell.hpp

@@ -23,6 +23,7 @@ class GridCell {
   GridCell(){}
   /// Initiate the grid cell to the coordinates (x0,y0)
   GridCell(int x, int y):x(x),y(y){}
+  bool operator == (const GridCell &a) const {return x == a.x && y == a.y;}
 };
 
 

include/richdem/common/random.cpp

@@ -0,0 +1,69 @@
+#include <richdem/common/random.hpp>
+
+#include <algorithm>
+#include <cassert>
+#include <cstdint>
+#include <functional>
+#include <iostream>
+#include <limits>
+#include <random>
+#include <sstream>
+
+namespace richdem {
+
+our_random_engine& rand_engine(){
+  static our_random_engine e[PRNG_THREAD_MAX];
+  return e[omp_get_thread_num()];
+}
+
+
+//Be sure to read: http://www.pcg-random.org/posts/cpp-seeding-surprises.html
+//and http://www.pcg-random.org/posts/cpps-random_device.html
+void seed_rand(unsigned long seed){
+  #pragma omp parallel      //All threads must come here
+  {
+    #pragma omp critical    //But only one at a time
+    if(seed==0){
+      std::uint_least32_t seed_data[std::mt19937::state_size];
+      std::random_device r;
+      std::generate_n(seed_data, std::mt19937::state_size, std::ref(r));
+      std::seed_seq q(std::begin(seed_data), std::end(seed_data));
+      rand_engine().seed(q);
+    } else
+      rand_engine().seed( seed*omp_get_thread_num() );
+  }
+}
+
+
+int uniform_rand_int(int from, int thru){
+  static std::uniform_int_distribution<> d[PRNG_THREAD_MAX];
+  using parm_t = std::uniform_int_distribution<>::param_type;
+  return d[omp_get_thread_num()]( rand_engine(), parm_t{from, thru} );
+}
+
+
+double uniform_rand_real(double from, double thru){
+  static std::uniform_real_distribution<> d[PRNG_THREAD_MAX];
+  using parm_t = std::uniform_real_distribution<>::param_type;
+  return d[omp_get_thread_num()]( rand_engine(), parm_t{from, thru} );
+}
+
+
+double normal_rand(double mean, double stddev){
+  static std::normal_distribution<double> d[PRNG_THREAD_MAX];
+  using parm_t = std::normal_distribution<double>::param_type;
+  return d[omp_get_thread_num()]( rand_engine(), parm_t{mean, stddev} );
+}
+
+RandomEngineState SaveRandomState(){
+  std::ostringstream oss;
+  oss<<rand_engine();
+  return oss.str();
+}
+
+void SetRandomState(const RandomEngineState &res){
+  std::istringstream iss(res);
+  iss>>rand_engine();
+}
+
+}

include/richdem/methods/catchment_delineation.hpp

@@ -0,0 +1,17 @@
+#ifndef _richdem_catchment_delineation_hpp_
+#define _richdem_catchment_delineation_hpp_
+
+#include <richdem/common/Array2D.hpp>
+#include <richdem/common/Array3D.hpp>
+#include <richdem/methods/catchment_delineation_generic.hpp>
+
+namespace richdem {
+
+// TODO: do we have to resize  and initialize the results array?
+template<class mask_t,               class result_t> void DC_mask_props(const Array2D<mask_t> &mask,    const Array3D<float> &props, Array2D<result_t> &upslope){std::queue<GridCell> expansion = queue_from_mask(mask,upslope)             ; upslope_cells_props(expansion, props, upslope);}
+template<class mask_t, class elev_t, class result_t> void DC_mask_mflow(const Array2D<mask_t> &mask,    const Array2D<elev_t> &elev, Array2D<result_t> &upslope){std::queue<GridCell> expansion = queue_from_mask(mask,upslope)             ; upslope_cells_mflow(expansion, elev,  upslope);}
+template<                            class result_t> void DC_line_props(int x0, int y0, int x1, int y1, const Array3D<float> &props, Array2D<result_t> &upslope){std::queue<GridCell> expansion = queue_from_linepoints(x0,y0,x1,y1,upslope); upslope_cells_props(expansion, props, upslope);}
+template<              class elev_t, class result_t> void DC_line_mflow(int x0, int y0, int x1, int y1, const Array2D<elev_t> &elev, Array2D<result_t> &upslope){std::queue<GridCell> expansion = queue_from_linepoints(x0,y0,x1,y1,upslope); upslope_cells_mflow(expansion, elev,  upslope);}
+
+}
+#endif

include/richdem/methods/catchment_delineation_generic.hpp

@@ -0,0 +1,165 @@
+#ifndef _richdem_catchment_delineation_generic_hpp_
+#define _richdem_catchment_delineation_generic_hpp_
+
+#include <richdem/common/Array2D.hpp>
+#include <richdem/common/Array3D.hpp>
+#include <richdem/common/logger.hpp>
+#include <richdem/common/ProgressBar.hpp>
+#include <richdem/methods/d8_methods.hpp> //only for sgn() function: TODO
+#include <queue>
+
+namespace richdem {
+
+template<class T, class U>
+std::queue<GridCell> queue_from_mask(
+  const Array2D<U>  &mask,
+  Array2D<T>              &upslope_cells
+){
+  std::queue<GridCell> expansion;
+  // this cannot be parallelized, because it appends everythign to the same queue
+  for (auto x=0;x<mask.width();x++){
+    for (auto y=0;y<mask.height();y++){
+      if(mask(x,y)!=mask.noData()){
+        upslope_cells(x,y)=2;
+        expansion.emplace(x,y);
+      }
+    }
+  }
+
+  return expansion;
+}
+
+template<class T>
+std::queue<GridCell> queue_from_linepoints(
+  int x0,
+  int y0,
+  int x1,
+  int y1,
+  Array2D<T>       &upslope_cells
+){
+  std::queue<GridCell> expansion;
+
+  if(x0>x1){
+    std::swap(x0,x1);
+    std::swap(y0,y1);
+  }
+
+  //Modified Bresenham Line-Drawing Algorithm
+  int deltax     = x1-x0;
+  int deltay     = y1-y0;
+  float error    = 0;
+  float deltaerr = (float)deltay/(float)deltax;
+
+  if (deltaerr<0)
+    deltaerr = -deltaerr;
+
+  RDLOG_MISC<<"Line slope is "<<deltaerr;
+  int y=y0;
+  for(int x=x0;x<=x1;x++){
+    expansion.push(GridCell(x,y));
+    upslope_cells(x,y)=2;
+    error+=deltaerr;
+    if (error>=0.5) {
+      expansion.push(GridCell(x+1,y));
+      upslope_cells(x+1,y) = 2;
+      y                   += sgn(deltay);
+      error               -= 1;
+    }
+  }
+
+  return expansion;
+}
+
+//upslope_from_props
+/**
+ * @brief Calculates uplsope cells for a given set of input cells, given a
+ * proportions matrix of flow. All cells that have flow into input cells will be added
+ * 
+ * @param[in]  &expansion      queue with starting cells
+ * @param[in]  &elevation      DEM
+ * @param[out] &upslope_cells  A grid of 1/2/NoData
+ */
+template<class T, class U>
+void upslope_cells_props(
+  std::queue<GridCell> &expansion,
+  const Array3D<T>     &props,
+  Array2D<U>           &upslope_cells
+){
+  //TODO: ALG NAME?
+  RDLOG_PROGRESS<<"Setting up the upslope_cells matrix..."<<std::flush;
+  // upslope_cells.resize(props.width(),props.height());
+  // upslope_cells.setAll(0);
+  // upslope_cells.setNoData(0);
+
+  ProgressBar progress;
+
+  progress.start(props.size());
+  long int ccount=0;
+  while(expansion.size()>0){
+    GridCell c=expansion.front();
+    expansion.pop();
+
+    progress.update(ccount++);
+
+    for(int n=1;n<=8;n++)
+      if(!props.inGrid(c.x+dx[n],c.y+dy[n]))
+        continue;
+      else if(upslope_cells.isNoData(c.x+dx[n],c.y+dy[n]) && props(c.x+dx[n],c.y+dy[n],d8_inverse[n])>0) {
+        expansion.emplace(c.x+dx[n],c.y+dy[n]);
+        upslope_cells(c.x+dx[n],c.y+dy[n])=1;
+      }
+  }
+  RDLOG_TIME_USE<<"Succeeded in "<<progress.stop();
+  RDLOG_MISC<<"Found "<<ccount<<" up-slope cells."; //TODO
+}
+
+//upslope_cells multiple flow implementation
+/**
+ * @brief Calculates uplsope cells for a given set of input cells, assuming
+ * multiple flow. That is, all cells that are higher than the cell being
+ * processed will be added, regardless of whether the current cell is the lowest
+ * neighbour.
+ * 
+ * 
+ * @param[in]  &expansion      queue with starting cells
+ * @param[in]  &elevation      DEM
+ * @param[out] &upslope_cells  A grid of 1/2/NoData
+ */
+template<class T, class U>
+void upslope_cells_mflow(
+  std::queue<GridCell> &expansion,
+  const Array2D<T>     &elevation,
+  Array2D<U>           &upslope_cells
+){
+  //TODO: ALG NAME?
+  RDLOG_PROGRESS<<"Setting up the upslope_cells matrix..."<<std::flush;
+  // upslope_cells.resize(elevation);
+  // upslope_cells.setAll(0);
+  // upslope_cells.setNoData(0);
+
+  ProgressBar progress;
+
+  progress.start(elevation.size());
+  long int ccount=0;
+  while(expansion.size()>0){
+    GridCell c=expansion.front();
+    expansion.pop();
+
+    progress.update(ccount++);
+
+    for(int n=1;n<=8;n++)
+      if(!elevation.inGrid(c.x+dx[n],c.y+dy[n]))
+        continue;
+      else if(elevation.isNoData(c.x+dx[n],c.y+dy[n]))
+        continue;
+      else if(upslope_cells.isNoData(c.x+dx[n],c.y+dy[n]) && elevation(c.x,c.y)<elevation(c.x+dx[n],c.y+dy[n])) {
+        expansion.emplace(c.x+dx[n],c.y+dy[n]);
+        upslope_cells(c.x+dx[n],c.y+dy[n])=1;
+      }
+  }
+  RDLOG_TIME_USE<<"Succeeded in "<<progress.stop();
+  RDLOG_MISC<<"Found "<<ccount<<" up-slope cells."; //TODO
+}
+
+}
+#endif

include/richdem/richdem.cpp

@@ -0,0 +1,46 @@
+#include <richdem/common/logger.hpp>
+
+#include <map>
+#include <string>
+
+namespace richdem {
+
+//TODO: Could use vector for this since the enum just compiles to an int. But
+//need to make sure things are in the right order.
+std::map<LogFlag, std::string> log_flag_chars_begin = {
+  {ALG_NAME, "\nA"},
+  {CITATION, "C"},
+  {CONFIG,   "c"},
+  {DEBUG,    "\033[95md"},
+  {ERROR,    "E"},
+  {MEM_USE,  " "},
+  {MISC,     "m"},
+  {PROGRESS, "p"},
+  {TIME_USE, "t"},
+  {WARN,     "\033[91mW"}
+};
+
+std::map<LogFlag, std::string> log_flag_chars_end = {
+  {ALG_NAME, ""},
+  {CITATION, "\n"},
+  {CONFIG,   ""},
+  {DEBUG,    ""},
+  {ERROR,    ""},
+  {MEM_USE,  ""},
+  {MISC,     ""},
+  {PROGRESS, ""},
+  {TIME_USE, ""},
+  {WARN,     ""}
+};
+
+void RDLOGfunc(LogFlag flag, const char* file, const char* func, unsigned line, std::string msg) {
+  std::cerr<<log_flag_chars_begin.at(flag)<<" "<<msg
+  #ifdef RICHDEM_DEBUG
+  <<" ("<<file<<" : "<<func<<" : "<<line<<")"
+  #endif
+  <<"\033[39m"
+  <<log_flag_chars_end.at(flag)
+  <<std::endl;
+}
+
+}

tests/catch_delineation/test_east.d8

@@ -0,0 +1,16 @@
+ncols         10
+nrows         10
+xllcorner     421568
+yllcorner     4872699
+cellsize      3
+NODATA_value  0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0
+2 1 1 1 1 1 1 1 1 0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0

tests/catch_delineation/test_east.dem

@@ -0,0 +1,16 @@
+ncols         10
+nrows         10
+xllcorner     421568
+yllcorner     4872699
+cellsize      3
+NODATA_value  -1
+1 2 3 4 5 6 7 8 9 10
+1 2 3 4 5 6 7 8 9 10
+1 2 3 4 5 6 7 8 9 10
+1 2 3 4 5 6 7 8 9 10
+1 2 3 4 5 6 7 8 9 10
+1 2 3 4 5 6 7 8 9 10
+1 2 3 4 5 6 7 8 9 10
+1 2 3 4 5 6 7 8 9 10
+1 2 3 4 5 6 7 8 9 10
+1 2 3 4 5 6 7 8 9 10

tests/catch_delineation/test_east.mask

@@ -0,0 +1,16 @@
+ncols         10
+nrows         10
+xllcorner     421568
+yllcorner     4872699
+cellsize      3
+NODATA_value  0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0
+1 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0
+0 0 0 0 0 0 0 0 0 0

tests/catch_delineation/test_east.mflow

@@ -0,0 +1,16 @@
+ncols         10
+nrows         10
+xllcorner     421568
+yllcorner     4872699
+cellsize      3
+NODATA_value  0
+0 0 0 0 1 1 1 1 1 1
+0 0 0 1 1 1 1 1 1 1
+0 0 1 1 1 1 1 1 1 1
+0 1 1 1 1 1 1 1 1 1
+0 1 1 1 1 1 1 1 1 1
+0 1 1 1 1 1 1 1 1 1
+0 0 1 1 1 1 1 1 1 1
+0 0 0 1 1 1 1 1 1 1
+0 0 0 1 1 1 1 1 1 1
+0 0 0 0 1 1 1 1 1 1