From c51ebafddfd7780b8d9a0a968bcbb90c302349e7 Mon Sep 17 00:00:00 2001
From: Joeperdefloep <ikbendeleukstethuis@gmail.com>
Date: Wed, 25 May 2022 14:48:37 +0200
Subject: [PATCH] Added upslope cell finding
going to do some things
fixed?
---
include/richdem/common/grid_cell.hpp | 1 +
include/richdem/methods/upslope_cells.hpp | 29 +++
.../methods/upslope_cells_functions.hpp | 215 ++++++++++++++++++
tests/tests.cpp | 183 ++++++++++++++-
tests/upslope_cells/test_east.d8 | 16 ++
tests/upslope_cells/test_east.dem | 16 ++
tests/upslope_cells/test_east.mask | 16 ++
tests/upslope_cells/test_east.mflow | 16 ++
tests/upslope_cells/test_east.quinn | 16 ++
wrappers/pyrichdem/richdem/__init__.py | 63 +++++
wrappers/pyrichdem/src/pywrapper.hpp | 15 ++
11 files changed, 585 insertions(+), 1 deletion(-)
create mode 100644 include/richdem/methods/upslope_cells.hpp
create mode 100644 include/richdem/methods/upslope_cells_functions.hpp
create mode 100644 tests/upslope_cells/test_east.d8
create mode 100644 tests/upslope_cells/test_east.dem
create mode 100644 tests/upslope_cells/test_east.mask
create mode 100644 tests/upslope_cells/test_east.mflow
create mode 100644 tests/upslope_cells/test_east.quinn
diff --git a/include/richdem/common/grid_cell.hpp b/include/richdem/common/grid_cell.hpp
index 6a8dfafe..3aff1f44 100644
--- a/include/richdem/common/grid_cell.hpp
+++ b/include/richdem/common/grid_cell.hpp
@@ -22,6 +22,7 @@ struct GridCell {
GridCell() = default;
/// Initiate the grid cell to the coordinates (x0,y0)
GridCell(int x0, int y0):x(x0),y(y0){}
+ bool operator == (const GridCell &a) const {return x == a.x && y == a.y;}
};
diff --git a/include/richdem/methods/upslope_cells.hpp b/include/richdem/methods/upslope_cells.hpp
new file mode 100644
index 00000000..0fc58d18
--- /dev/null
+++ b/include/richdem/methods/upslope_cells.hpp
@@ -0,0 +1,29 @@
+#ifndef _richdem_upslope_cells_hpp_
+#define _richdem_upslope_cells_hpp_
+
+#include <richdem/common/Array2D.hpp>
+#include <richdem/common/Array3D.hpp>
+#include <richdem/methods/upslope_cells_functions.hpp>
+
+namespace richdem {
+
+template<class T, class U>
+void init(
+ Array2D<T> &to_init,
+ const Array2D<U> &example,
+ const int &val
+){
+ to_init.resize(example);
+ to_init.setAll(val);
+ to_init.setNoData(val);
+}
+
+
+// TODO: do we have to resize and initialize the results array?
+template<class mask_t, class result_t> void UC_mask_props(const Array2D<mask_t> &mask, const Array3D<float> &props, Array2D<result_t> &upslope){init(upslope,mask,0); std::queue<GridCell> expansion = queue_from_mask(mask,upslope) ; upslope_cells_props<Topology::D8>(expansion, props, upslope);}
+template<class mask_t, class elev_t, class result_t> void UC_mask_mflow(const Array2D<mask_t> &mask, const Array2D<elev_t> &elev, Array2D<result_t> &upslope){init(upslope,elev,0); std::queue<GridCell> expansion = queue_from_mask(mask,upslope) ; upslope_cells_mflow<Topology::D8>(expansion, elev, upslope);}
+template< class result_t> void UC_line_props(int x0, int y0, int x1, int y1, const Array3D<float> &props, Array2D<result_t> &upslope){init(upslope,props,0); std::queue<GridCell> expansion = queue_from_linepoints(x0,y0,x1,y1,upslope); upslope_cells_props<Topology::D8>(expansion, props, upslope);}
+template< class elev_t, class result_t> void UC_line_mflow(int x0, int y0, int x1, int y1, const Array2D<elev_t> &elev, Array2D<result_t> &upslope){init(upslope,elev,0); std::queue<GridCell> expansion = queue_from_linepoints(x0,y0,x1,y1,upslope); upslope_cells_mflow<Topology::D8>(expansion, elev, upslope);}
+
+}
+#endif //_richdem_upslope_cells_
\ No newline at end of file
diff --git a/include/richdem/methods/upslope_cells_functions.hpp b/include/richdem/methods/upslope_cells_functions.hpp
new file mode 100644
index 00000000..cb8a1ac8
--- /dev/null
+++ b/include/richdem/methods/upslope_cells_functions.hpp
@@ -0,0 +1,215 @@
+#ifndef _richdem_upslope_cells_functions_hpp_
+#define _richdem_upslope_cells_functions_hpp_
+
+#include <richdem/common/Array2D.hpp>
+#include <richdem/common/Array3D.hpp>
+#include <richdem/common/logger.hpp>
+#include <richdem/common/ProgressBar.hpp>
+#include <queue>
+
+namespace richdem {
+
+/**
+ * @brief Add all non-nodata cells to the return queue and sets the
+ * ``upslope_cells`` value to 2 there.
+ *
+ * @tparam T upslope cells dtype
+ * @tparam U mask dtype
+ * @param[in] mask The mask
+ * @param[out] upslope_cells grid for indicating upslope cells
+ * @return std::queue<GridCell>
+ */
+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 everything 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;
+}
+
+/**
+ * @brief Draws a line between the points (x0,y0) and (x1,y1). Adds visited
+ * cells to the return queue and sets them to 2 in ``upslope_cells``.
+ *
+ * @tparam T upslope cells dtype
+ * @param[in] x0 x coordinate of first point
+ * @param[in] y0 y coordinate of first point
+ * @param[in] x1 x coordinate of second point
+ * @param[in] y1 y coordinate of second point
+ * @param[out] upslope_cells grid for indicating upslope cells
+ * @return std::queue<GridCell> A queue with gridcells that were visited
+ */
+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 the slope>1, swap x and y coordinates.
+ // It works exatly the same, but points
+ // are drawn with swapped coordinates
+ bool swapxy = false;
+ if(std::abs(float(y0-y1)/float(x0-x1))>1){
+ std::cout << "swapped x and y" << std::endl;
+ std::swap(x0,y0);
+ std::swap(x1,y1);
+ swapxy= true;
+ }
+
+ if(x0>x1){
+ std::swap(x0,x1);
+ std::swap(y0,y1);
+ }
+
+ // modified Bresenham Line-Drawing Algorithm from
+ // https://www.geeksforgeeks.org/bresenhams-line-generation-algorithm/
+ int deltaerr = 2*(y1-y0);
+ int error = deltaerr - (x1-x0);
+
+ RDLOG_MISC<<"Line slope is "<<deltaerr;
+ int y=y0;
+ for(int x=x0;x<=x1;x++){
+ if (!swapxy){
+ expansion.emplace(x,y);
+ upslope_cells(x,y)=2;
+ } else {
+ expansion.emplace(y,x);
+ upslope_cells(y,x)=2;
+ }
+ error += deltaerr;
+ if (error>=0 && x!=x1) { //safeguard for segfault if we draw to the edge of the grid
+ // Draw the point above if we move to te next pixel
+ // the difference between a "solid" and non-solid line:
+ // o o X x x <--the X
+ // x x x 0 0
+ if(!swapxy){
+ expansion.emplace(x+1,y);
+ upslope_cells(x+1,y) = 2;
+ } else {
+ expansion.emplace(y,x+1);
+ upslope_cells(y,x+1) = 2;
+ }
+ y++;
+ error -= 2*(x1-x0);
+ }
+ }
+ 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<Topology topo, 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);
+
+ static_assert(topo==Topology::D8 || topo==Topology::D4);
+ constexpr auto dx = get_dx_for_topology<topo>();
+ constexpr auto dy = get_dy_for_topology<topo>();
+ constexpr auto nmax = get_nmax_for_topology<topo>();
+
+ 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<Topology topo, 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);
+
+ static_assert(topo==Topology::D8 || topo==Topology::D4);
+ constexpr auto dx = get_dx_for_topology<topo>();
+ constexpr auto dy = get_dy_for_topology<topo>();
+ constexpr auto nmax = get_nmax_for_topology<topo>();
+
+ 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
\ No newline at end of file
diff --git a/tests/tests.cpp b/tests/tests.cpp
index e3bb4263..5eb77ef7 100644
--- a/tests/tests.cpp
+++ b/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>
@@ -455,4 +457,183 @@ TEST_CASE("DirectionsMatchExpectations"){
CHECK(get_nmax_for_topology<Topology::D4>() == 4);
}
-}
\ No newline at end of file
+}
+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<Topology::D8>(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<Topology::D8>(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");
+
+
+ }
+ }
+ }
diff --git a/tests/upslope_cells/test_east.d8 b/tests/upslope_cells/test_east.d8
new file mode 100644
index 00000000..a62ba2a5
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/tests/upslope_cells/test_east.dem b/tests/upslope_cells/test_east.dem
new file mode 100644
index 00000000..7e763b9b
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/tests/upslope_cells/test_east.mask b/tests/upslope_cells/test_east.mask
new file mode 100644
index 00000000..5d5543fb
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/tests/upslope_cells/test_east.mflow b/tests/upslope_cells/test_east.mflow
new file mode 100644
index 00000000..8f12c7aa
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/tests/upslope_cells/test_east.quinn b/tests/upslope_cells/test_east.quinn
new file mode 100644
index 00000000..1839a2d2
--- /dev/null
+++ b/tests/upslope_cells/test_east.quinn
@@ -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 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
+2 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 0
+0 0 0 0 0 0 0 0 0 0
\ No newline at end of file
diff --git a/wrappers/pyrichdem/richdem/__init__.py b/wrappers/pyrichdem/richdem/__init__.py
index fadc80a5..7e5b5399 100644
--- a/wrappers/pyrichdem/richdem/__init__.py
+++ b/wrappers/pyrichdem/richdem/__init__.py
@@ -897,3 +897,66 @@ def fill_spill_merge(dem: rdarray, labels: rdarray, flowdirs: rdarray, deps: Lis
wtdw = wtd.wrap()
dhret = depression_hierarchy.fill_spill_merge(demw, labelsw, flowdirsw, deps, wtdw)
+
+def UpslopeCells(
+ dem,
+ x0: int = None,
+ y0: int = None,
+ x1: int = None,
+ y1: int = None,
+ mask: rdarray = None,
+ method: str = None
+):
+ """
+ Delineates a catchment, based on input points, or a true/false raster.
+
+ Args:
+ dem (rdarray): A digital elevation model
+ x0 (int): x coordinate of first point
+ y0 (int): y coordinate of first point
+ x1 (int): x coordinate of second point
+ y1 (int): y coordinate of second point
+ init_mask (rdarray): Mask array with nonzero values for input cells
+ method (string): The method to use. Can be "mflow", "D8" or any flow metric, such as "Dinf"
+ """
+ if type(dem) is not rdarray:
+ raise Exception("A richdem.rdarray or numpy.ndarray is required!")
+
+ result = rdarray(np.zeros(shape=dem.shape,dtype="uint8"),meta_obj=dem,no_data=0)
+ resultw = result.wrap()
+
+ if method == "mflow":
+ demw = dem.wrap()
+ if x0 is not None and y0 is not None and x1 is not None and y1 is not None:
+ _richdem.UC_line_mflow(x0,y0,x1,y1,demw,resultw)
+ elif x0 is not None and y0 is not None:
+ _richdem.UC_line_mflow(x0,y0,x0,y0,demw,resultw)
+ elif mask is not None:
+ maskw = mask.wrap()
+ _richdem.UC_mask_mflow(maskw,demw,resultw)
+ else:
+ propsw = FlowProportions(dem,method).wrap()
+ if x0 is not None and y0 is not None and x1 is not None and y1 is not None:
+ _richdem.UC_line_props(x0,y0,x1,y1,propsw,resultw)
+ elif x0 is not None and y0 is not None:
+ _richdem.UC_line_props(x0,y0,x0,y0,propsw,resultw)
+ elif mask is not None:
+ maskw = mask.wrap()
+ _richdem.UC_mask_props(maskw,propsw,resultw)
+
+ # init_maskw = init_mask.wrap()
+ # initqueue = _richdem.InitQueueFromMask(demw,init_maskw)
+ # elif x0 is not None and y0 is not None and x1 is not None and y1 is not None:
+ # initqueue = _richdem.InitQueueFromLinePoints(x0,y0,x1,y1,demw)
+ # elif x0 is not None and y0 is not None:
+ # initqueue = _richdem.InitQueueFromLinePoints(x0,y0,x0,y0,demw)
+
+ # if a method is specified, calculate proportions
+ if method == "mflow":
+ method = "Quinn" #TODO: make this specialized
+ # else:
+
+ result.copyFromWrapped(resultw)
+
+ return result
+
diff --git a/wrappers/pyrichdem/src/pywrapper.hpp b/wrappers/pyrichdem/src/pywrapper.hpp
index c6419ef5..3cb8d9b1 100644
--- a/wrappers/pyrichdem/src/pywrapper.hpp
+++ b/wrappers/pyrichdem/src/pywrapper.hpp
@@ -79,6 +79,21 @@ void TemplatedFunctionsWrapper(pybind11::module &m, std::string tname){
m.def("FM_OCallaghanD4", &FM_OCallaghan <Topology::D4,T>, "TODO");
m.def("FM_D8", &FM_D8 <T>, "TODO");
m.def("FM_D4", &FM_D4 <T>, "TODO");
+
+ m.def("UC_mask_props", &UC_mask_props <T,uint8_t>, "TODO");
+ // TODO: make this work for different datatypes. Do we just make a big list?
+ m.def("UC_mask_mflow", &UC_mask_mflow <float ,T,uint8_t>, "TODO");
+ m.def("UC_mask_mflow", &UC_mask_mflow <double ,T,uint8_t>, "TODO");
+ m.def("UC_mask_mflow", &UC_mask_mflow <int8_t ,T,uint8_t>, "TODO");
+ m.def("UC_mask_mflow", &UC_mask_mflow <int16_t ,T,uint8_t>, "TODO");
+ m.def("UC_mask_mflow", &UC_mask_mflow <int32_t ,T,uint8_t>, "TODO");
+ m.def("UC_mask_mflow", &UC_mask_mflow <int64_t ,T,uint8_t>, "TODO");
+ m.def("UC_mask_mflow", &UC_mask_mflow <uint8_t ,T,uint8_t>, "TODO");
+ m.def("UC_mask_mflow", &UC_mask_mflow <uint16_t,T,uint8_t>, "TODO");
+ m.def("UC_mask_mflow", &UC_mask_mflow <uint32_t,T,uint8_t>, "TODO");
+ m.def("UC_mask_mflow", &UC_mask_mflow <uint64_t,T,uint8_t>, "TODO");
+ m.def("UC_line_props", &UC_line_props <uint8_t>, "TODO");
+ m.def("UC_line_mflow", &UC_line_mflow <T,uint8_t>, "TODO");
}