Added upslope cell finding

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/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/tests.cpp

@@ -5,6 +5,8 @@
 #include <richdem/common/loaders.hpp>
 #include <richdem/misc/misc_methods.hpp>
 #include <richdem/richdem.hpp>
+#include <richdem/methods/upslope_cells_functions.hpp>
+#include <richdem/methods/upslope_cells.hpp>
 
 #include <filesystem>
 #include <queue>
@@ -426,3 +428,183 @@ TEST_CASE("BucketFill"){
     CHECK(set_raster==good_raster);
   }
 }
+
+TEST_CASE("Checking Catchment Delineation Generic"){
+  SUBCASE("mask"){
+    SUBCASE("single"){
+    Array2D<uint8_t> mask = {
+      {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},
+    };
+    mask.setNoData(0);
+    Array2D<uint8_t> upslope(mask,0);
+
+    std::queue<GridCell> q_expected;
+    q_expected.emplace(2,2);
+    const Array2D<uint8_t> upslope_e = {
+      {0, 0, 0, 0, 0, 0},
+      {0, 0, 0, 0, 0, 0},
+      {0, 0, 2, 0, 0, 0},
+      {0, 0, 0, 0, 0, 0},
+      {0, 0, 0, 0, 0, 0},
+    };
+
+    std::queue<GridCell> q = queue_from_mask(mask, upslope);
+    std::cout<< mask.noData() << std::endl;
+    CHECK(q.front().x == q_expected.front().x);
+    CHECK(q.front().y == q_expected.front().y);
+    CHECK(q.size() == 1);
+    CHECK(upslope(2,2) == 2);
+    }
+    SUBCASE("multi"){
+    Array2D<uint8_t> mask = {
+      {0, 0, 0, 0, 0, 0},
+      {0, 1, 0, 0, 0, 0},
+      {0, 0, 1, 0, 0, 0},
+      {0, 1, 0, 0, 0, 0},
+      {0, 0, 0, 0, 0, 0},
+    };
+    mask.setNoData(0);
+    Array2D<uint8_t> upslope(mask,0);
+
+    std::queue<GridCell> q_expected;
+    q_expected.emplace(1,1);
+    q_expected.emplace(1,3);
+    q_expected.emplace(2,2);
+
+    const Array2D<uint8_t> upslope_e = {
+      {0, 0, 0, 0, 0, 0},
+      {0, 2, 0, 0, 0, 0},
+      {0, 0, 2, 0, 0, 0},
+      {0, 2, 0, 0, 0, 0},
+      {0, 0, 0, 0, 0, 0},
+    };
+
+    std::queue q = queue_from_mask(mask, upslope);
+    CHECK(q.size() == 3);
+    CHECK(q == q_expected);
+    CHECK(upslope == upslope_e);
+    }
+    // TODO: make more tests?
+  }
+
+  SUBCASE("line"){
+    SUBCASE("horizontal"){
+      int x0=0;
+      int y0=0;
+      int x1=2;
+      int y1=0;
+
+      std::queue<GridCell> q_expected;
+      q_expected.emplace(0,0);
+      q_expected.emplace(1,0);
+      q_expected.emplace(2,0);
+
+      Array2D<uint8_t> upslope_e = {
+        {2, 2, 2},
+        {0, 0, 0},
+        {0, 0, 0},
+      };
+      upslope_e.setNoData(0);
+      Array2D<uint8_t> upslope(upslope_e,0);
+      upslope.setNoData(0);
+
+      std::queue<GridCell> q = queue_from_linepoints(x0,y0,x1,y1,upslope);
+
+      CHECK(upslope == upslope_e);
+      CHECK(q == q_expected);
+    }
+    //TODO: add more tests. definately for slope>1, so that swapping is
+    //routinely tested. (I now did it by inspection, and it works)
+  }
+
+  Array2D<int> dem = {
+    {1,2,3,4,5},
+    {1,2,3,4,5},
+    {1,2,3,4,5},
+    {1,2,3,4,5},
+    {1,2,3,4,5},
+  };
+  dem.setNoData(0);
+
+  SUBCASE("mflow"){
+
+    Array2D<uint8_t> result(dem,0);
+    std::queue<GridCell> expansion;
+    expansion.emplace(0,2);
+
+    Array2D<uint8_t> expected = {
+      {0, 0, 1, 1, 1},
+      {0, 1, 1, 1, 1},
+      {0, 1, 1, 1, 1},
+      {0, 1, 1, 1, 1},
+      {0, 0, 1, 1, 1}
+    };
+    result.setNoData(0);
+    expected.setNoData(0); //the function sets nodata to 0
+
+    upslope_cells_mflow(expansion, dem, result);
+    CHECK(result == expected);
+  }
+
+  SUBCASE("props (take quinn)"){
+    //more methods is not necessary
+    Array2D<uint8_t> result(dem,0);
+    result.setNoData(0);
+    std::queue<GridCell> expansion;
+    expansion.emplace(0,2);
+
+    //construct the proportions array
+    Array3D<float> props(dem);
+    FM_Quinn(dem,props);
+
+    Array2D<uint8_t> expected = {
+      {0, 0, 0, 0, 0},
+      {0, 1, 1, 1, 0},
+      {0, 1, 1, 1, 0},
+      {0, 1, 1, 1, 0},
+      {0, 0, 0, 0, 0}      // TODO: maybe add edge cells if they are higher? is slightly complex though, cuz Array3D does not support that...
+    };                     // so flow always goes off the edge of the grid, hence the 0s.
+    expected.setNoData(0); //the function sets nodata to 0
+
+    upslope_cells_props(expansion, props, result);
+    CHECK(result == expected);
+  }
+}
+
+TEST_CASE("Checking upslope cells") {
+  for(auto p: fs::directory_iterator("upslope_cells")){
+    fs::path this_path = p.path();
+    if(this_path.extension()==".dem"){
+        Array2D<int>     dem      (this_path,                            false);
+        Array2D<uint8_t> mask     (this_path.replace_extension("mask"),  false);
+        Array2D<int32_t> ans_mflow(this_path.replace_extension("mflow"), false);
+        Array2D<int32_t> ans_quinn(this_path.replace_extension("quinn"), false);
+        Array2D<int32_t> ans_d8   (this_path.replace_extension("d8"),    false);
+        Array2D<int32_t> result;
+        Array3D<float>   props(dem);
+
+        UC_mask_mflow(mask,dem,result);
+        CHECK(result == ans_mflow);
+        result.printAll("mflow result");
+        ans_mflow.printAll("expected result");
+
+        FM_Quinn(dem, props);
+        UC_mask_props(mask,props,result);
+        CHECK(result == ans_quinn);
+        result.printAll("Quinn upslope result");
+        ans_quinn.printAll("expected result");
+
+        FM_D8(dem, props);
+        UC_mask_props(mask,props,result);
+        CHECK(result == ans_d8);
+        result.printAll("D8 upslope result");
+        ans_d8.printAll("expected result");
+
+
+      }
+    }
+  }

tests/upslope_cells/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/upslope_cells/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/upslope_cells/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/upslope_cells/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
+2 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 0 1 1 1 1 1 1
+0 0 0 0 0 1 1 1 1 1