Add Depression Hierarchy to RichDEM proper

apps/CMakeLists.txt

@@ -1,43 +1,45 @@
 cmake_minimum_required (VERSION 3.9)
 
-add_executable(rd_compare.exe                 rd_compare.cpp) # -Wno-sign-compare
+add_executable(rd_compare.exe                 rd_compare.cpp)
 add_executable(rd_d8_flowdirs.exe             rd_d8_flowdirs.cpp)
-add_executable(rd_raster_display.exe          rd_raster_display.cpp)
+add_executable(rd_depression_hierarchy.exe    rd_depression_hierarchy.cpp)
 add_executable(rd_depressions_breach.exe      rd_depressions_breach.cpp)
+add_executable(rd_depressions_flood.exe       rd_depressions_flood.cpp)
+add_executable(rd_depressions_has.exe         rd_depressions_has.cpp)
+add_executable(rd_depressions_mask.exe        rd_depressions_mask.cpp)
 add_executable(rd_expand_dimensions.exe       rd_expand_dimensions.cpp)
 add_executable(rd_flood_for_flowdirs.exe      rd_flood_for_flowdirs.cpp)
 add_executable(rd_flow_accumulation.exe       rd_flow_accumulation.cpp)
 add_executable(rd_geotransform.exe            rd_geotransform.cpp)
-add_executable(rd_depressions_has.exe         rd_depressions_has.cpp)
 add_executable(rd_hist.exe                    rd_hist.cpp)
-add_executable(rd_depressions_mask.exe        rd_depressions_mask.cpp)
-add_executable(rd_raster_inspect.exe          rd_raster_inspect.cpp)
 add_executable(rd_loop_check.exe              rd_loop_check.cpp)
 add_executable(rd_merge_rasters_by_layout.exe rd_merge_rasters_by_layout.cpp)
-add_executable(rd_depressions_flood.exe       rd_depressions_flood.cpp)
 add_executable(rd_no_data.exe                 rd_no_data.cpp)
-add_executable(rd_taudem_d8_to_richdem_d8.exe rd_taudem_d8_to_richdem_d8.cpp)
-add_executable(rd_projection.exe              rd_projection.cpp)
 add_executable(rd_processing_history.exe      rd_processing_history.cpp)
+add_executable(rd_projection.exe              rd_projection.cpp)
+add_executable(rd_raster_display.exe          rd_raster_display.cpp)
+add_executable(rd_raster_inspect.exe          rd_raster_inspect.cpp)
+add_executable(rd_taudem_d8_to_richdem_d8.exe rd_taudem_d8_to_richdem_d8.cpp)
 add_executable(rd_terrain_property.exe        rd_terrain_property.cpp)
 
 target_link_libraries(rd_compare.exe                  richdem)
 target_link_libraries(rd_d8_flowdirs.exe              richdem)
-target_link_libraries(rd_raster_display.exe           richdem)
+target_link_libraries(rd_depression_hierarchy.exe     richdem)
 target_link_libraries(rd_depressions_breach.exe       richdem)
+target_link_libraries(rd_depressions_flood.exe        richdem)
+target_link_libraries(rd_depressions_has.exe          richdem)
+target_link_libraries(rd_depressions_mask.exe         richdem)
 target_link_libraries(rd_expand_dimensions.exe        richdem)
 target_link_libraries(rd_flood_for_flowdirs.exe       richdem)
 target_link_libraries(rd_flow_accumulation.exe        richdem)
 target_link_libraries(rd_geotransform.exe             richdem)
-target_link_libraries(rd_depressions_has.exe          richdem)
 target_link_libraries(rd_hist.exe                     richdem)
-target_link_libraries(rd_depressions_mask.exe         richdem)
-target_link_libraries(rd_raster_inspect.exe           richdem)
 target_link_libraries(rd_loop_check.exe               richdem)
 target_link_libraries(rd_merge_rasters_by_layout.exe  richdem)
-target_link_libraries(rd_depressions_flood.exe        richdem)
 target_link_libraries(rd_no_data.exe                  richdem)
-target_link_libraries(rd_taudem_d8_to_richdem_d8.exe  richdem)
-target_link_libraries(rd_projection.exe               richdem)
 target_link_libraries(rd_processing_history.exe       richdem)
+target_link_libraries(rd_projection.exe               richdem)
+target_link_libraries(rd_raster_display.exe           richdem)
+target_link_libraries(rd_raster_inspect.exe           richdem)
+target_link_libraries(rd_taudem_d8_to_richdem_d8.exe  richdem)
 target_link_libraries(rd_terrain_property.exe         richdem)

apps/rd_depression_hierarchy.cpp

@@ -0,0 +1,73 @@
+#include <richdem/common/Array2D.hpp>
+#include <richdem/common/gdal.hpp>
+#include <richdem/depressions/depression_hierarchy.hpp>
+
+#include <iostream>
+#include <string>
+#include <stdexcept>
+
+namespace rd = richdem;
+namespace dh = richdem::dephier;
+
+int main(int argc, char **argv){
+  if(argc!=4){
+    std::cout<<"Syntax: "<<argv[0]<<" <Input> <Output Prefix> <Ocean Level>"<<std::endl;
+    return -1;
+  }
+
+  const std::string in_name     = argv[1];
+  const std::string out_name    = argv[2];
+  const float       ocean_level = std::stod(argv[3]);
+
+  std::cout<<"m Processing    = "<<argv[1]<<std::endl;
+  std::cout<<"m Output prefix = "<<argv[2]<<std::endl;
+  std::cout<<"m Ocean level   = "<<argv[3]<<std::endl;
+
+  rd::Timer timer_io;
+  timer_io.start();
+  rd::Array2D<float> topo(in_name);   //Recharge (Percipitation minus Evapotranspiration)
+  timer_io.stop();
+
+  std::cout<<"m Data width  = "<<topo.width ()<<std::endl;
+  std::cout<<"m Data height = "<<topo.height()<<std::endl;
+  std::cout<<"m Data cells  = "<<topo.numDataCells()<<std::endl;
+
+  rd::Array2D<dh::dh_label_t> label   (topo.width(), topo.height(), dh::NO_DEP ); //No cells are part of a depression
+  rd::Array2D<rd::flowdir_t>  flowdirs(topo.width(), topo.height(), rd::NO_FLOW); //No cells flow anywhere
+
+  //Label the ocean cells. This is a precondition for using
+  //`GetDepressionHierarchy()`.
+  #pragma omp parallel for
+  for(unsigned int i=0;i<label.size();i++){
+    if(topo.isNoData(i) || topo(i)==ocean_level){ //Ocean Level is assumed to be lower than any other cells (even Death Valley)
+      label(i) = dh::OCEAN;
+    }
+  }
+
+  //Generate flow directions, label all the depressions, and get the hierarchy
+  //connecting them
+  auto deps = dh::GetDepressionHierarchy<float,rd::Topology::D8>(topo, label, flowdirs);
+
+  timer_io.start();
+
+  //GraphViz dot-style output for drawing depression hierarchy graphs.
+  {
+    std::ofstream fgraph(out_name+"-graph.csv");
+    fgraph<<"parent,child,oceanlink\n";
+    for(unsigned int i=1;i<deps.size();i++){
+      fgraph<<deps[i].parent<<","<<i<<",";
+      //Is it an ocean link?
+      fgraph<<(deps[i].parent!=dh::NO_VALUE && (deps[i].parent==dh::OCEAN || !(deps[deps[i].parent].lchild==i || deps[deps[i].parent].rchild==i)))<<"\n";
+    }
+  }
+
+  label.projection = topo.projection;
+  label.saveGDAL(out_name+"-label.tif");
+
+  timer_io.stop();
+
+  std::cout<<"Finished"<<std::endl;
+  std::cout<<"IO time   = "<<timer_io.accumulated()<<" s"<<std::endl;
+
+  return 0;
+}

include/richdem/common/disjoint_dense_int_set.hpp

@@ -0,0 +1,150 @@
+//DisjointDenseIntSet
+//Author: Richard Barnes (rijard.barnes@gmail.com)
+#ifndef _disjoint_dense_int_sets_
+#define _disjoint_dense_int_sets_
+
+#include <vector>
+#include <cassert>
+
+
+//A disjoint-set/union-find class. Starting from a collection of sets, this data
+//structure efficiently keeps track of which sets have been merged. It is
+//assumed that every integer between 0 and some maximum value N is a set, so the
+//data structure takes O(N) space. If only `findSet()` and `unionSet()` are
+//used, then all accesses are in O(a(N)) time, where `a()` is the inverse
+//Ackermann function. For all practical purposes, this is `(1). If
+//`mergeAintoB()` is used then `findSet()` can have a worst-case of `O(N)`.
+class DisjointDenseIntSet {
+ private:
+  //How many sets are children of this set. Initially 0.
+  std::vector<unsigned int> rank;
+  //Which set is this set's parent. May be the set itself.
+  std::vector<unsigned int> parent;
+
+  //When a set value X is passed to the data structure, this checks to see if
+  //the set exists. If not, the set is created, along with all the sets between
+  //the old maximum set id and X.
+  void checkSize(unsigned int newN){
+    if(newN<rank.size())               //Does the set exist?
+      return;                          //Yup.
+                                       //Nope.
+    const auto old_size = rank.size(); //Get old maximum set value
+
+    //Resize so that `newN` is a valid value. None of the new sets have
+    //children.
+    rank.resize(newN+1,0);
+
+    //Resize so that `newN` is a valid value
+    parent.resize(newN+1);
+
+    //Ensure that each new set is its own parent since they have not yet been
+    //merged.
+    for(auto i=old_size;i<newN+1;i++)
+      parent[i] = i;
+  }
+
+ public:
+  //Construct a DisjointDenseIntSet without any sets. Sets will be dynamically
+  //created as the data structure is used.
+  DisjointDenseIntSet(){}
+
+  //Create a DisjointDenseIntSet with `N` initial sets preallocated. More sets
+  //can be dynamically allocated as the data structure is used.
+  DisjointDenseIntSet(unsigned int N){
+    rank.resize(N,0);
+    parent.resize(N);
+    for(unsigned int i=0;i<N;i++)
+      parent[i] = i;
+  }
+
+  //Explicitly creates a set. May incidentally create several intermediate sets
+  //of `n` is more than one larger than the maximum set id previously seen.
+  void makeSet(unsigned int n){
+    checkSize(n);
+  }
+
+  //Returns the highest set id.
+  unsigned int maxElement() const {
+    return rank.size()-1;
+  }
+
+  //Follows a set's chain of parents until a set which is its own parent is
+  //reached. This ultimate parent's id is returned as the representative id of
+  //the set in question. Note that this collapses the chain of parents so that
+  //after this method has run every set between the one in question and the
+  //ultimate parent points to the ultimate parent. This means that while the
+  //first call to this function may take `O(N)` lookups in the worst-case (less
+  //due to the use of ranks, as explained below), subsequent calls to any set in
+  //the chain will take `O(1)` time. This technique is known as "path
+  //compression".
+  unsigned int findSet(unsigned int n){
+    if(n>=parent.size()){
+      throw std::runtime_error("DisjointDenseIntSet::findSet(" + std::to_string(n) + ") is looking for a set outside the valid range, which is [0," + std::to_string(parent.size()) + ")!");
+    }
+    if(parent[n]==n){                  //Am I my own parent?
+      return n;                        //Yes: I represent the set in question.
+    } else {                           //No.
+      //Who is my parent's ultimate parent? Make them my parent.
+      return parent[n] = findSet(parent[n]);
+    }
+  }
+
+  //Join two sets into a single set. Note that we "cannot" predict the `id` of
+  //the resulting set ahead of time.
+  void unionSet(unsigned int a, unsigned int b){
+    const auto roota = findSet(a); //Find the ultimate parent of A
+    const auto rootb = findSet(b); //Find the ultimate parent of B
+    //Note that the foregoing collapses any chain of parents so that each set in
+    //the chain points to the ultimate parent. Therefore, any subsequent call to
+    //`findSet` involving any set in the chain will take `O(1)` time.
+
+    //If A and B already share a parent, then they do not need merging.
+    if(roota==rootb)
+      return;
+
+    //If we always naively tacked A onto B then we could develop a worst-case
+    //scenario in which each set pointed to just one other set in a long, linear
+    //chain. If this happened then calls to `findSet()` would take `O(N)` time.
+    //Instead, we keep track of how many child sets each set has and ensure that
+    //the shorter tree of sets becomes part of the taller tree of sets. This
+    //ensures that the tree does not grow taller unless the two trees were of
+    //equal height in which case the resultant tree is taller by 1. In essence,
+    //this bounds the depth of any query to being `log_2(N)`. However, due to
+    //the use of path compression above, the query path is actually less than
+    //this.
+
+    if(rank[roota]<rank[rootb]){          //Is A shorter?
+      parent[roota] = rootb;              //Make A a child of B.
+    } else if(rank[roota]>rank[rootb]) {  //Is B shorter?
+      parent[rootb] = roota;              //Make B a child of A.
+    } else {                              //A and B are the same height
+      parent[rootb] = roota;              //Arbitrarily make B a child of A
+      rank[roota]++;                      //Increase A's height.
+    }
+  }
+
+  //Using `unionSet` merges two sets in a way which does not allow us to decide
+  //which set is the parent; however, `unionSet` helps guarantee fast queries.
+  //`mergeAintoB` sacrifices speed but preserves parenthood by always making A a
+  //child of B, regardless of the height of `B`.
+  void mergeAintoB(unsigned int a, unsigned int b){
+    checkSize(a);
+    checkSize(b);
+    parent[a] = b;
+    if(rank[a]==rank[b]){
+      rank[b]++;
+    } else if(rank[a]>rank[b]){
+      rank[b] = rank[a]+1;
+    } else {
+      //If `rank[b]>rank[a]` then making A a child of B does not increase B's
+      //height.
+    }
+  }
+
+  //Returns true if A and B belong to the same set.
+  bool sameSet(unsigned int a, unsigned int b){
+    return findSet(a)==findSet(b);
+  }
+};
+
+#endif