diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt
index 32378a93..a1fa9349 100644
--- a/apps/CMakeLists.txt
+++ b/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)
\ No newline at end of file
diff --git a/apps/rd_depression_hierarchy.cpp b/apps/rd_depression_hierarchy.cpp
new file mode 100644
index 00000000..4e1d17c3
--- /dev/null
+++ b/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;
+}
diff --git a/include/richdem/common/disjoint_dense_int_set.hpp b/include/richdem/common/disjoint_dense_int_set.hpp
new file mode 100644
index 00000000..632c091c
--- /dev/null
+++ b/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
diff --git a/include/richdem/common/radix_heap.hpp b/include/richdem/common/radix_heap.hpp
new file mode 100644
index 00000000..774f41f4
--- /dev/null
+++ b/include/richdem/common/radix_heap.hpp
@@ -0,0 +1,394 @@
+#pragma once
+#ifndef HEADER_RADIX_HEAP
+#define HEADER_RADIX_HEAP
+
+#include <algorithm>
+#include <array>
+#include <cassert>
+#include <climits>
+#include <cstdint>
+#include <iostream>
+#include <limits>
+#include <type_traits>
+#include <utility>
+#include <vector>
+#include <tuple>
+
+namespace radix_heap {
+namespace internal {
+template<bool Is64bit> class find_bucket_impl;
+
+template<>
+class find_bucket_impl<false> {
+ public:
+  static inline constexpr size_t find_bucket(uint32_t x, uint32_t last) {
+    return x == last ? 0 : 32 - __builtin_clz(x ^ last);
+  }
+};
+
+template<>
+class find_bucket_impl<true> {
+ public:
+  static inline constexpr size_t find_bucket(uint64_t x, uint64_t last) {
+    return x == last ? 0 : 64 - __builtin_clzll(x ^ last);
+  }
+};
+
+template<typename T>
+inline constexpr size_t find_bucket(T x, T last) {
+  return find_bucket_impl<sizeof(T) == 8>::find_bucket(x, last);
+}
+
+template<typename KeyType, bool IsSigned> class encoder_impl_integer;
+
+template<typename KeyType>
+class encoder_impl_integer<KeyType, false> {
+ public:
+  typedef KeyType key_type;
+  typedef KeyType unsigned_key_type;
+
+  inline static constexpr unsigned_key_type encode(key_type x) {
+    return x;
+  }
+
+  inline static constexpr key_type decode(unsigned_key_type x) {
+    return x;
+  }
+};
+
+template<typename KeyType>
+class encoder_impl_integer<KeyType, true> {
+ public:
+  typedef KeyType key_type;
+  typedef typename std::make_unsigned<KeyType>::type unsigned_key_type;
+
+  inline static constexpr unsigned_key_type encode(key_type x) {
+    return static_cast<unsigned_key_type>(x) ^
+        (unsigned_key_type(1) << unsigned_key_type(std::numeric_limits<unsigned_key_type>::digits - 1));
+  }
+
+  inline static constexpr key_type decode(unsigned_key_type x) {
+    return static_cast<key_type>
+        (x ^ (unsigned_key_type(1) << (std::numeric_limits<unsigned_key_type>::digits - 1)));
+  }
+};
+
+template<typename KeyType, typename UnsignedKeyType>
+class encoder_impl_decimal {
+public:
+  typedef KeyType key_type;
+  typedef UnsignedKeyType unsigned_key_type;
+
+  inline static constexpr unsigned_key_type encode(key_type x) {
+    return raw_cast<key_type, unsigned_key_type>(x) ^
+        ((-(raw_cast<key_type, unsigned_key_type>(x) >> (std::numeric_limits<unsigned_key_type>::digits - 1))) |
+         (unsigned_key_type(1) << (std::numeric_limits<unsigned_key_type>::digits - 1)));
+  }
+
+  inline static constexpr key_type decode(unsigned_key_type x) {
+    return raw_cast<unsigned_key_type, key_type>
+        (x ^ (((x >> (std::numeric_limits<unsigned_key_type>::digits - 1)) - 1) |
+              (unsigned_key_type(1) << (std::numeric_limits<unsigned_key_type>::digits - 1))));
+  }
+
+ private:
+  template<typename T, typename U>
+  union raw_cast {
+   public:
+    constexpr raw_cast(T t) : t_(t) {}
+    operator U() const { return u_; }
+
+   private:
+    T t_;
+    U u_;
+  };
+};
+
+template<typename KeyType>
+class encoder : public encoder_impl_integer<KeyType, std::is_signed<KeyType>::value> {};
+template<>
+class encoder<float> : public encoder_impl_decimal<float, uint32_t> {};
+template<>
+class encoder<double> : public encoder_impl_decimal<double, uint64_t> {};
+
+template <typename UIntDecType>
+class bucket_flags {
+ public:
+  typedef UIntDecType unsigned_key_type;
+
+  bucket_flags() {clear();}
+
+  void set_empty(size_t bucket) {
+    assert(bucket <= num_buckets);
+
+    flags_ &= ~((unsigned_key_type(1) << (bucket - 1)) * !!bucket);
+  }
+
+  void set_non_empty(unsigned int bucket) {
+    assert(bucket <= num_buckets);
+
+    flags_ |= (unsigned_key_type(1) << (bucket - 1)) * !!bucket;
+  }
+
+  void clear() {
+    flags_ = 0;
+  }
+
+  size_t find_first_non_empty() const {
+    if (sizeof(flags_) == 8)
+      return __builtin_ffsll(flags_);
+
+    return __builtin_ffs(flags_);
+  }
+
+  void swap(bucket_flags& a) {
+    std::swap(flags_, a.flags_);
+  }
+
+private:
+  unsigned_key_type flags_;
+  constexpr static size_t num_buckets = (sizeof(flags_) == 8 ? 64 : 32);
+};
+}  // namespace internal
+
+
+
+template<typename KeyType, typename EncoderType = internal::encoder<KeyType>>
+class radix_heap {
+ public:
+  typedef KeyType key_type;
+  typedef EncoderType encoder_type;
+  typedef typename encoder_type::unsigned_key_type unsigned_key_type;
+
+  radix_heap() : size_(0), last_(), buckets_() {
+    buckets_min_.fill(std::numeric_limits<unsigned_key_type>::max());
+  }
+
+  void push(key_type key) {
+    const unsigned_key_type x = encoder_type::encode(key);
+    assert(last_ <= x);
+    ++size_;
+    const size_t k = internal::find_bucket(x, last_);
+    buckets_[k].emplace_back(x);
+    bucket_flags_.set_non_empty(k);
+    buckets_min_[k] = std::min(buckets_min_[k], x);
+  }
+
+  key_type top() {
+    pull();
+    return encoder_type::decode(last_);
+  }
+
+  key_type min() const {
+    assert(size_ > 0);
+
+    const size_t i = bucket_flags_.find_first_non_empty() * buckets_[0].empty();
+    return encoder_type::decode(buckets_min_[i]);
+  }
+
+  void pop() {
+    pull();
+    buckets_[0].pop_back();
+    --size_;
+  }
+
+  size_t size() const {
+    return size_;
+  }
+
+  bool empty() const {
+    return size_ == 0;
+  }
+
+  void clear() {
+    size_ = 0;
+    last_ = key_type();
+    for (auto &b : buckets_) b.clear();
+    bucket_flags_.clear();
+    buckets_min_.fill(std::numeric_limits<unsigned_key_type>::max());
+  }
+
+  void swap(radix_heap<KeyType, EncoderType> &a) {
+    std::swap(size_, a.size_);
+    std::swap(last_, a.last_);
+    buckets_.swap(a.buckets_);
+    buckets_min_.swap(a.buckets_min_);
+    bucket_flags_.swap(a.bucket_flags_);
+  }
+
+ private:
+  size_t size_;
+  unsigned_key_type last_;
+  std::array<std::vector<unsigned_key_type>,
+             std::numeric_limits<unsigned_key_type>::digits + 1> buckets_;
+  std::array<unsigned_key_type,
+             std::numeric_limits<unsigned_key_type>::digits + 1> buckets_min_;
+
+  internal::bucket_flags<unsigned_key_type> bucket_flags_;
+
+  void pull() {
+    assert(size_ > 0);
+    if (!buckets_[0].empty()) return;
+
+    const size_t i = bucket_flags_.find_first_non_empty();
+    last_ = buckets_min_[i];
+
+    for (unsigned_key_type x : buckets_[i]) {
+      const size_t k = internal::find_bucket(x, last_);
+      buckets_[k].emplace_back(x);
+      bucket_flags_.set_non_empty(k);
+      buckets_min_[k] = std::min(buckets_min_[k], x);
+    }
+
+    buckets_[i].clear();
+    bucket_flags_.set_empty(i);
+    buckets_min_[i] = std::numeric_limits<unsigned_key_type>::max();
+  }
+};
+
+template<typename KeyType, typename ValueType, typename EncoderType = internal::encoder<KeyType>>
+class pair_radix_heap {
+ public:
+  typedef KeyType key_type;
+  typedef ValueType value_type;
+  typedef EncoderType encoder_type;
+  typedef typename encoder_type::unsigned_key_type unsigned_key_type;
+  typedef typename std::vector<std::pair<key_type, value_type> >::const_iterator value_iterator;
+
+  pair_radix_heap() : size_(0), last_(), buckets_() {
+    buckets_min_.fill(std::numeric_limits<unsigned_key_type>::max());
+  }
+
+  void push(key_type key, const value_type &value) {
+    const unsigned_key_type x = encoder_type::encode(key);
+    assert(last_ <= x);
+    ++size_;
+    const size_t k = internal::find_bucket(x, last_);
+    buckets_[k].emplace_back(x, value);
+    bucket_flags_.set_non_empty(k);
+    buckets_min_[k] = std::min(buckets_min_[k], x);
+  }
+
+  void push(key_type key, value_type &&value) {
+    const unsigned_key_type x = encoder_type::encode(key);
+    assert(last_ <= x);
+    ++size_;
+    const size_t k = internal::find_bucket(x, last_);
+    buckets_[k].emplace_back(x, std::move(value));
+    bucket_flags_.set_non_empty(k);
+    buckets_min_[k] = std::min(buckets_min_[k], x);
+  }
+
+  template <class... Args>
+  void emplace(key_type key, Args&&... args) {
+    const unsigned_key_type x = encoder_type::encode(key);
+    assert(last_ <= x);
+    ++size_;
+    const size_t k = internal::find_bucket(x, last_);
+    buckets_[k].emplace_back(std::piecewise_construct,
+                             std::forward_as_tuple(x), std::forward_as_tuple(args...));
+    bucket_flags_.set_non_empty(k);
+    buckets_min_[k] = std::min(buckets_min_[k], x);
+  }
+
+  key_type min_key() const {
+    assert(size_ > 0);
+    const size_t i = buckets_[0].empty() ? bucket_flags_.find_first_non_empty() : 0;
+    assert(last_ <= buckets_min_[i]);
+    return encoder_type::decode(buckets_min_[i]);
+  }
+
+  key_type top_key() {
+    pull();
+    return encoder_type::decode(last_);
+  }
+
+  value_type &top_value() {
+    pull();
+    return buckets_[0].back().second;
+  }
+
+  std::tuple<value_iterator, value_iterator> top_values() {
+    pull();
+    return std::forward_as_tuple(buckets_[0].cbegin(), buckets_[0].cend());
+  }
+
+  void pop() {
+    pull();
+    buckets_[0].pop_back();
+    --size_;
+  }
+
+  void pop_top_values() {
+      pull();
+      size_ -= buckets_[0].size();
+      buckets_[0].clear();
+  }
+
+
+  std::vector<std::pair<unsigned_key_type, value_type> > extract_top_values() {
+    pull();
+    std::vector<std::pair<unsigned_key_type, value_type> > vec;
+    size_ -= buckets_[0].size();
+    buckets_[0].swap(vec);
+    return vec;
+  }
+
+  size_t size() const {
+    return size_;
+  }
+
+  bool empty() const {
+    return size_ == 0;
+  }
+
+  void clear() {
+    size_ = 0;
+    last_ = key_type();
+    for (auto &b : buckets_) b.clear();
+    buckets_min_.fill(std::numeric_limits<unsigned_key_type>::max());
+    bucket_flags_.clear();
+  }
+
+  void swap(pair_radix_heap<KeyType, ValueType, EncoderType> &a) {
+    std::swap(size_, a.size_);
+    std::swap(last_, a.last_);
+    bucket_flags_.swap(a.bucket_flags_);
+    buckets_.swap(a.buckets_);
+    buckets_min_.swap(a.buckets_min_);
+  }
+
+ private:
+  size_t size_;
+  unsigned_key_type last_;
+  std::array<std::vector<std::pair<unsigned_key_type, value_type>>,
+             std::numeric_limits<unsigned_key_type>::digits + 1> buckets_;
+  std::array<unsigned_key_type,
+             std::numeric_limits<unsigned_key_type>::digits + 1> buckets_min_;
+
+  internal::bucket_flags<unsigned_key_type> bucket_flags_;
+
+  void pull() {
+    assert(size_ > 0);
+    if (!buckets_[0].empty()) return;
+    buckets_min_[0] = std::numeric_limits<unsigned_key_type>::max();
+
+    const size_t i = bucket_flags_.find_first_non_empty();
+    last_ = buckets_min_[i];
+
+    for (size_t j = 0; j < buckets_[i].size(); ++j) {
+      const unsigned_key_type x = buckets_[i][j].first;
+      const size_t k = internal::find_bucket(x, last_);
+      buckets_[k].emplace_back(std::move(buckets_[i][j]));
+      bucket_flags_.set_non_empty(k);
+      buckets_min_[k] = std::min(buckets_min_[k], x);
+    }
+
+    buckets_[i].clear();
+    bucket_flags_.set_empty(i);
+    buckets_min_[i] = std::numeric_limits<unsigned_key_type>::max();
+  }
+};
+}  // namespace radix_heap
+
+#endif // HEADER_RADIX_HEAP
\ No newline at end of file
diff --git a/include/richdem/depressions/depression_hierarchy.hpp b/include/richdem/depressions/depression_hierarchy.hpp
new file mode 100644
index 00000000..e202e55f
--- /dev/null
+++ b/include/richdem/depressions/depression_hierarchy.hpp
@@ -0,0 +1,859 @@
+#pragma once
+
+#include <richdem/common/Array2D.hpp>
+#include <richdem/common/constants.hpp>
+#include <richdem/common/disjoint_dense_int_set.hpp>
+#include <richdem/common/grid_cell.hpp>
+#include <richdem/common/logger.hpp>
+#include <richdem/common/math.hpp>
+#include <richdem/common/ProgressBar.hpp>
+#include <richdem/common/radix_heap.hpp>
+#include <richdem/common/timer.hpp>
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstdint>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <queue>
+#include <stdexcept>
+#include <string>
+#include <unordered_map>
+#include <utility>
+
+namespace richdem::dephier {
+
+//We use a 32-bit integer for labeling depressions. This allows for a maximum of
+//2,147,483,647 depressions. This should be enough for most practical purposes.
+typedef uint32_t dh_label_t;
+typedef uint32_t flat_c_idx;
+
+//Some special valuess
+const dh_label_t NO_PARENT = std::numeric_limits<dh_label_t>::max();
+const dh_label_t NO_VALUE  = std::numeric_limits<dh_label_t>::max();
+
+//This class holds information about a depression. Its pit cell and outlet cell
+//(in flat-index form) as well as the elevations of these cells. It also notes                                                   //what is flat-index form?
+//the depression's parent. The parent of the depression is the outlet through
+//which it must flow in order to reach the ocean. If a depression has more than
+//one outlet at the same level one of them is arbitrarily chosen; hopefully this                                                  //so, everything should have a parent except for the ocean, right?
+//happens only rarely in natural environments.
+template<class elev_t>
+struct Depression {
+  //Flat index of the pit cell, the lowest cell in the depression. If more than
+  //one cell shares this lowest elevation, then one is arbitrarily chosen.
+  flat_c_idx pit_cell = NO_VALUE;
+  //Flat index of the outlet cell. If there is more than one outlet cell at this
+  //cell's elevation, then one is arbitrarily chosen.
+  flat_c_idx out_cell = NO_VALUE;
+  //Parent depression. If both this depression and its neighbour fill up, this
+  //parent depression is the one which will contain the overflow.
+  dh_label_t parent   = NO_PARENT;
+  //Outlet depression. The metadepression into which this one overflows. Usually
+  //its neighbour depression, but sometimes the ocean.
+  dh_label_t odep     = NO_VALUE;
+  //When a metadepression overflows it does so into the metadepression indicated
+  //by `odep`. However, odep must flood from the bottom up. Therefore, we keep
+  //track of the `geolink`, which indicates what leaf depression the overflow is
+  //initially routed into.
+  dh_label_t geolink  = NO_VALUE;
+  //Elevation of the pit cell. Since the pit cell has the lowest elevation of
+  //any cell in the depression, we initialize this to infinity.
+  elev_t  pit_elev = std::numeric_limits<elev_t>::infinity();
+  //Elevation of the outlet cell. Since the outlet cell has the lowest elevation
+  //of any path leading from a depression, we initialize this to infinity.
+  elev_t  out_elev = std::numeric_limits<elev_t>::infinity();
+  //The depressions form a binary tree. Each depression has two child
+  //depressions: one left and one right.
+  dh_label_t lchild = NO_VALUE;
+  dh_label_t rchild = NO_VALUE;
+  //Indicates whether the parent link is to either the ocean or a depression
+  //that links to the ocean
+  bool ocean_parent = false;
+  //Indicates depressions which link to the ocean through this depression, but
+  //are not subdepressions. That is, these ocean-linked depressions may be at
+  //the top of high cliffs and spilling into this depression.
+  std::vector<dh_label_t> ocean_linked;
+  //the label of the depression, for calling it up again
+  dh_label_t dep_label = 0;
+  //Number of cells contained within the depression and its children
+  uint32_t cell_count = 0;
+  //Total of elevations within the depression, used in the WLE. Because I think I need to start adding up total elevations before I know the outlet of the depression.
+  //double dep_sum_elevations = 0;
+  //Volume of the depression and its children. Used in the Water Level Equation (see below).
+  double   dep_vol    = 0;
+  //Water currently contained within the depression. Used in the Water Level
+  //Equation (see below).
+  double   water_vol  = 0;
+
+  //Total elevation of cells contained with the depression and its children
+  double total_elevation = 0;
+};
+
+
+
+//A key part of the algorithm is keeping track of the outlets which connect
+//depressions. While each depression has only one outlet, a depression may have
+//many inlets. All of the outlets and inlets taken together form a graph which
+//we can traverse to determine which way water flows. This class keeps track of
+//which cell links two depressions, as well as the elevation of that cell.
+
+//The OutletLink is used as a key for a hashtable which stores information about
+//the outlets.
+struct OutletLink {
+  dh_label_t depa;
+  dh_label_t depb;
+  OutletLink() = default;
+  OutletLink(dh_label_t depa0, dh_label_t depb0) : depa(depa0), depb (depb0) {}
+  //This is used to compare two outlets. The outlets are the same regardless of
+  //the order in which they store depressions
+  bool operator==(const OutletLink &o) const {
+    return depa==o.depa && depb==o.depb;
+  }
+};
+
+//The outlet class stores, again, the depressions being linked as well as
+//information about the link
+template<class elev_t>
+struct Outlet {
+  dh_label_t depa;                //Depression A
+  dh_label_t depb;                //Depression B
+  flat_c_idx out_cell = NO_VALUE; //Flat-index of cell at which A and B meet.
+  //Elevation of the cell linking A and B
+  elev_t  out_elev = std::numeric_limits<elev_t>::infinity();
+
+  Outlet() = default;
+
+  //Standard issue constructor
+  Outlet(dh_label_t depa0, dh_label_t depb0, flat_c_idx out_cell0, elev_t out_elev0){
+    depa       = depa0;
+    depb       = depb0;
+    if(depa>depb)           //Create a preferred ordering so that comparisons and hashing are faster
+      std::swap(depa,depb);
+    out_cell   = out_cell0;
+    out_elev   = out_elev0;
+  }
+
+  //Determines whether one outlet is the same as another. Note that we do not
+  //check elevation for this! This is because we'll be using this operator to
+  //determine if we've already found an outlet for a depression. We'll look at
+  //outlets from lowest to highest, so if an outlet already exists for a
+  //depression, it is that depression's lowest outlet.
+  bool operator==(const Outlet &o) const {                                                                                              //so beyond just checking, is this somehow preventing it from being recorded if one already exists? How does this work?
+    //Outlets are the same if they link two depressions, regardless of the
+    //depressions' labels storage order within this class.
+    return depa==o.depa && depb==o.depb;
+  }
+};
+
+//We'll initially keep track of outlets using a hash table. Hash tables require
+//that every item they contain be reducible to a numerical "key". We provide
+//such a key for outlets here.
+template<class elev_t>
+struct OutletHash {
+  std::size_t operator()(const OutletLink &out) const {
+    //Since depa and depb are sorted on construction, we don't have to worry
+    //about which order the invoking code called them in and our hash function
+    //doesn't need to be symmetric with respect to depa and depb.
+
+    //Boost hash combine function: https://stackoverflow.com/q/35985960/752843
+    return out.depa^(out.depb + 0x9e3779b9 + (out.depa << 6) + (out.depa >> 2));
+  }
+};
+
+
+
+//The regular mod function allows negative numbers to stay negative. This mod
+//function wraps negative numbers around. For instance, if a=-1 and n=100, then
+//the result is 99.
+int ModFloor(int a, int n) {
+  return ((a % n) + n) % n;
+}
+
+
+template<class elev_t>
+using PriorityQueue = radix_heap::pair_radix_heap<elev_t,uint64_t>;
+// using PriorityQueue = GridCellZk_high_pq<elev_t> pq;
+
+
+//Cell is not part of a depression
+const dh_label_t NO_DEP = std::numeric_limits<dh_label_t>::max();
+//Cell is part of the ocean and a place from which we begin searching for
+//depressions.
+const dh_label_t OCEAN  = 0;
+
+template<typename elev_t>
+using DepressionHierarchy = std::vector<Depression<elev_t>>;
+
+template<class elev_t>
+void CalculateMarginalVolumes(DepressionHierarchy<elev_t> &deps, const Array2D<elev_t> &dem, const Array2D<dh_label_t> &label);
+
+template<class elev_t>
+void CalculateTotalVolumes(DepressionHierarchy<elev_t> &deps);
+
+
+
+template<class elev_t>
+std::ostream& operator<<(std::ostream &out, const DepressionHierarchy<elev_t> &deps){
+  std::vector<int> child_count;
+  std::vector<size_t> stack;
+
+  const std::function<void(const size_t root, const size_t depth)> print_helper = [&](const size_t root, const size_t depth) -> void {
+    const auto &dep = deps.at(root);
+    stack.push_back(root);
+    child_count.push_back( (dep.lchild!=NO_VALUE) + (dep.rchild!=NO_VALUE) + dep.ocean_linked.size() );
+
+    for(int i=0;i<depth;i++){
+      if(child_count.at(i)>1 && i==depth-1)
+        out<<(dep.ocean_parent?"╠═":"├─");
+      else if(child_count.at(i)==1 && i==depth-1)
+        out<<(dep.ocean_parent?"╚═":"└─");
+      else if(child_count.at(i)>2)
+        out<<"║ ";
+      else if(child_count.at(i)>1)
+        out<<"│ ";
+      else
+        out<<"  ";
+    }
+
+    out<<"Id="<<root
+       <<", dep_vol="<<dep.dep_vol
+       <<", water_vol="<<dep.water_vol
+       <<", pit_cell="<<dep.pit_cell
+       <<", out_cell="<<dep.out_cell
+       <<", out_elev="<<dep.out_elev
+       <<", parent="<<dep.parent
+       <<", odep="<<dep.odep
+       <<"\n";
+
+    if(dep.lchild!=NO_VALUE) { print_helper(dep.lchild,depth+1); child_count.back()--; }
+    if(dep.rchild!=NO_VALUE) { print_helper(dep.rchild,depth+1); child_count.back()--; }
+    for(const auto &x: dep.ocean_linked){
+      print_helper(x,depth+1);
+      child_count.back()--;
+    }
+    stack.pop_back();
+    child_count.pop_back();
+  };
+
+  print_helper(0, 0);
+
+  return out;
+}
+
+
+
+//Calculate the hierarchy of depressions. Takes as input a digital elevation
+//model and a set of labels. The labels should have `OCEAN` for cells
+//representing the "ocean" (the place to which depressions drain) and `NO_DEP`
+//for all other cells.
+//
+//@param  dem   - 2D array of elevations. May be in any data format.
+//
+//@return label - A label indiciate which depression the cell belongs to.
+//                The indicated label is always the leaf of the depression
+//                hierarchy, or the OCEAN.
+//
+//        flowdirs - A value [0,7] indicated which direction water from the cell
+//                   flows in order to go "downhill". All cells have a flow
+//                   direction (even flats) except for pit cells.
+template<class elev_t, Topology topo>
+DepressionHierarchy<elev_t> GetDepressionHierarchy(
+  const Array2D<elev_t> &dem,
+  Array2D<dh_label_t>   &label,
+  Array2D<int8_t>       &flowdirs
+){
+  ProgressBar progress;
+  Timer timer_overall;
+  Timer timer_dephier;
+  timer_overall.start();
+
+  timer_dephier.start();
+
+  RDLOG_ALG_NAME<<"DepressionHierarchy";
+
+  //A D4 or D8 topology can be used.
+  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 dinverse = get_dinverse_for_topology<topo>();
+  constexpr auto neighbours = get_nmax_for_topology<topo>();
+
+  //Depressions are identified by a number [0,*). The ocean is always
+  //"depression" 0. This vector holds the depressions.
+  DepressionHierarchy<elev_t> depressions;
+
+  //This keeps track of the outlets we find. Each pair of depressions can only
+  //be linked once and the lowest link found between them is the one which is
+  //retained.
+  typedef std::unordered_map<OutletLink, Outlet<elev_t>, OutletHash<elev_t>> outletdb_t;
+  outletdb_t outlet_database;
+
+  //Places to seed depression growth from. These vectors are used to make the
+  //search for seeds parallel, yet deterministic.
+  std::vector<flat_c_idx> ocean_seeds;
+  std::vector<flat_c_idx> land_seeds;
+  //Reduce reallocations by assuming 2.5% of the map is seeds
+  ocean_seeds.reserve(dem.width()*dem.height()/40);
+  land_seeds.reserve(dem.width()*dem.height()/40);
+
+  #pragma omp declare reduction(merge : std::vector<flat_c_idx> : omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end()))
+
+  RDLOG_PROGRESS<<"p Adding ocean cells to priority-queue...";
+  //We assume the user has already specified a few ocean cells from which to
+  //begin looking for depressions. We add all of these ocean cells to the
+  //priority queue now.
+  uint64_t ocean_cells = 0;
+  #pragma omp parallel for collapse(2) reduction(+:ocean_cells) reduction(merge:ocean_seeds)
+  for(int y=0;y<dem.height();y++)
+  for(int x=0;x<dem.width();x++){
+    //Ensure the input only has OCEAN and NO_DEP labels.
+    if(label(x,y)!=OCEAN){
+      if(label(x,y)!=NO_DEP){
+        throw std::runtime_error("Label array given to GetDepressionHierarchy must contain only NO_DEP and OCEAN labels!");
+      }
+      continue;
+    }
+
+    //We'll only add ocean cells to the PQ if they border a non-ocean cell
+    bool has_non_ocean = false;
+    for(int n=1;n<=neighbours;n++){
+      if(label.inGrid(x+dx[n],y+dy[n]) && label(x+dx[n],y+dy[n])!=OCEAN){
+        has_non_ocean = true;
+        break;
+      }
+    }
+    if(has_non_ocean){       //If they are ocean cells, put them in the priority queue
+      ocean_seeds.emplace_back(dem.xyToI(x,y));
+      ocean_cells++;
+    }
+  }
+
+  if(ocean_cells==0){
+    throw std::runtime_error("No OCEAN cells found, could not make a DepressionHierarchy!");
+  }
+
+  //The 0th depression is the ocean. We add it to the list of depressions now
+  //that we're sure there is an ocean!
+  { //Use a little scope to avoid having `oceandep` linger around
+    auto &oceandep    = depressions.emplace_back();
+    //The ocean is deep
+    oceandep.pit_elev = -std::numeric_limits<elev_t>::infinity();
+    //It's so deep we can't find its bottom
+    oceandep.pit_cell = NO_VALUE;
+    oceandep.dep_label = 0;
+  }
+
+
+  RDLOG_PROGRESS<<"p Finding pit cells...";
+
+  //Here we find the pit cells of internally-draining regions. We define these
+  //to be cells without any downstream neighbours. Note that this means we will
+  //identify all flat cells as being pit cells. For DEMs with many flat cells,
+  //this will bloat the priortiy queue slightly. If your DEM includes extensive,
+  //predictably located flat regions, you may wish to add these in some special
+  //way. Alternatively, you could use Barnes (2014, "An Efficient Assignment of
+  //Drainage Direction Over Flat Surfaces") as a way of reducing the number of
+  //flat cells. Regardless, the algorithm will deal gracefully with the flats it
+  //finds and this shouldn't slow things down too much!
+  int pit_cell_count = 0;
+  progress.start(dem.size());
+  #pragma omp parallel for collapse(2) reduction(+:pit_cell_count) reduction(merge:land_seeds)
+  for(int y=0;y<dem.height();y++)  //Look at all the cells
+  for(int x=0;x<dem.width() ;x++){ //Yes, all of them
+    ++progress;
+    if(label(x,y)==OCEAN)          //Already in priority queue
+      continue;
+    const auto my_elev = dem(x,y); //Focal cell's elevation
+    bool has_lower     = false;    //Pretend we have no lower neighbours
+    for(int n=1;n<=neighbours;n++){ //Check out our neighbours
+      //Use offset to get neighbour x coordinate, wrapping as needed
+      // const int nx = ModFloor(x+dx[n],dem.width());
+      const int nx = x+dx[n];
+      //Use offset to get neighbour y coordinate
+      const int ny = y+dy[n];
+      if(!dem.inGrid(nx,ny))  //Is cell outside grid (too far North/South)?
+        continue;             //Yup: skip it.
+      if(dem(nx,ny)<my_elev){ //Is this neighbour lower than focal cell?
+        has_lower = true;     //Make a note of it
+        break;                //Don't need to look at additional neighbours
+      }
+    }
+    if(!has_lower){           //The cell can't drain, so it is a pit cell
+      land_seeds.emplace_back(dem.xyToI(x,y));
+      pit_cell_count++;       //Add to pit cell count. Parallel safe because of reduction.
+    }
+  }
+  progress.stop();
+  RDLOG_TIME_USE<<"t Pit cells found in = "<<progress.time_it_took()<<" s";
+
+  //Since the above runs in parallel, the ordering of the seed cells is
+  //nondeterministic. Let's fix that.
+  std::sort(ocean_seeds.begin(), ocean_seeds.end());
+  std::sort(land_seeds.begin(), land_seeds.end());
+
+  //The priority queue ensures that cells are visited in order from lowest to
+  //highest. If two or more cells are of equal elevation then the one added last
+  //(most recently) is returned from the queue first. This ensures that a single
+  //depression gets all the cells within a flat area.
+  PriorityQueue<elev_t> pq;
+
+  //Add all the seed cells to the PQ
+  for(const auto &x: ocean_seeds)
+    pq.emplace(dem(x), x);
+  ocean_seeds.clear();
+  ocean_seeds.shrink_to_fit();
+
+  for(const auto &x: land_seeds)
+    pq.emplace(dem(x), x);
+  land_seeds.clear();
+  land_seeds.shrink_to_fit();
+
+  //The priority queue now contains all of the ocean cells as well as all of the
+  //pit cells. We will now proceed through the cells by always pulling the cell
+  //of lowest elevation from the priority queue. This ensures that, when we find
+  //an outlet between two depressions, it is always the lowest outlet linking
+  //them.
+
+  //Once two depressions meet, they do not battle for dominance. Rather, we
+  //build an invisible, and solely conceptual, wall between them. Each
+  //depression continues to grow, by encompassing cells which have not yet been
+  //assigned to a depression, until it is surrounded by other depressions on all
+  //sides. The depression then contains its pit cell, its outlet, every cell
+  //below its outlet, and possibly many cells above its outlet which ultimately
+  //drain into its pit cell. (Though note that if the depression has a flat
+  //bottom the pit cell is chosen arbitrarily as one of the flat cells.) Later,
+  //when we construct the depression hierarchy, we will separate the cells above
+  //a depression's outlet into new meta-depressions.
+
+  //In the following we'll temporarily relax our definition of an outlet to mean
+  //"the lowest connection between two depressions" rather than "the lowest
+  //connection out of a depression". This means depressions may have outlets at
+  //many elevations. Later on we'll fix this and some of those outlets will
+  //become inlets or the outlets of meta-depressions.
+
+  //The hash table of outlets will dynamically resize as we add elements to it.
+  //However, this slows things down a bit. Therefore, we presize the hash set to
+  //be equal to be 3x the number of pit cells plus the ocean cell. 3 is just a
+  //guess as to how many neighbouring depressions each depression will have. If
+  //we get this value too small we lose a little speed due to rehashing. If we
+  //get this value too large then we waste space. If we get this value far too
+  //large then we run out of RAM.
+  outlet_database.reserve(3*(pit_cell_count+1));
+
+  //Visit cells in order of elevation from lowest to highest. If two or more
+  //cells are of the same elevation then we visit the one added last (most
+  //recently) first.
+
+  RDLOG_PROGRESS<<"p Searching for outlets...";
+
+  progress.start(dem.size());
+  while(!pq.empty()){
+    ++progress;
+
+    const auto ci    = pq.top_value();     //Copy cell with lowest elevation from priority queue
+    const auto celev = pq.top_key();       //Elevation of focal cell
+    pq.pop();                              //Remove the copied cell from the priority queue
+    auto clabel = label(ci);               //Nominal label of cell
+    int cx,cy;
+    dem.iToxy(ci,cx,cy);
+
+    if(clabel==OCEAN){
+      //This cell is an ocean cell or a cell that flows into the ocean without
+      //encountering any depressions on the way. Upon encountering it we do not
+      //need to do anything special.
+    } else if(clabel==NO_DEP){
+      //Since cells label their neighbours and ocean cells are labeled in the
+      //initialization, the only way to get to a cell that is still labeled as
+      //not being part of a depression is if that cell were added as a pit cell.
+      //For each pit cell we find, we make a new depression and label it
+      //accordingly. Not all the pit cells originally added will form new
+      //depressions as flat cells will relabel their neighbours and the first
+      //cell found in a flat determines the label for the entirety of that flat.
+      clabel            = depressions.size();         //In a 0-based indexing system, size is equal to the id of the next flat
+      auto &newdep      = depressions.emplace_back(); //Add the next flat (increases size by 1)
+      newdep.pit_cell   = dem.xyToI(cx,cy);           //Make a note of the pit cell's location
+      newdep.pit_elev   = celev;                      //Make a note of the pit cell's elevation
+      newdep.dep_label  = clabel;                     //I am storing the label in the object so that I can find it later and call up the number of cells and volume (better way of doing this?) -- I have since realised I can use the index in the depressions array. So perhaps the label is no longer needed?
+      label(ci)         = clabel;                     //Update cell with new label
+    } else {
+
+      //Cell has already been assigned to a depression. In this case, one of two
+      //things is true. (1) This cell is on the frontier of our search, in which
+      //case the cell has neighbours which have not yet been seen. (2) This cell
+      //was part of a flat which has previously been processed by a wavefront
+      //beginning at some other cell. In this case, all of this cell's
+      //neighbours will have already been seen and added to the priority queue.
+      //However, it is harmless to check on them again.
+    }
+
+    //TODO: Update the appropriate depression's cell_count and dep_vol variables                        I did this in the else if above, and then in the if and the else below. I add a cell to the count whenever it is added to the depression and add its elevation to the total elevations.
+    //here.                                                                                             Then I calculate the total volume only when we find an outlet (Good way to test this? Print values of volumes only of those that make the outlet queue? I get some negative values sometimes so I may have done something wrong, but what if it's an 'outlet' at the highest point of the depression?)
+
+    //Consider the cell's neighbours
+    for(int n=1;n<=neighbours;n++){
+      // const int nx = ModFloor(cx+dx[n],dem.width()); //Get neighbour's x-coordinate using an offset and wrapping
+      const int nx = cx + dx[n];                      //Get neighbour's y-coordinate using an offset
+      const int ny = cy + dy[n];                      //Get neighbour's y-coordinate using an offset
+      if(!dem.inGrid(nx,ny))                          //Is this cell in the grid?
+        continue;                                     //Nope: out of bounds.
+      const auto ni     = dem.xyToI(nx,ny);           //Flat index of neighbour
+      const auto nlabel = label(ni);                  //Label of neighbour
+
+      if(nlabel==NO_DEP){                             //Neighbour has not been visited yet
+        label(ni) = clabel;                           //Give the neighbour my label
+        pq.emplace(dem(ni), dem.xyToI(nx,ny));        //Add the neighbour to the priority queue
+        flowdirs(nx,ny) = dinverse[n];                //Neighbour flows in the direction of this cell
+      } else if (nlabel==clabel) {
+        //Skip because we are not interested in ourself. That would be vain.
+        //Note that this case will come up frequently as we traverse flats since
+        //the first cell to be visited in a flat labels all the cells in the
+        //flat like itself. All of the other cells in the flat will come off of
+        //the priority queue later and, in looking at their neighbours, reach
+        //this point.
+      } else {
+        //We've found a neighbouring depression!
+
+        //Determine whether the focal cell or this neighbour is the outlet of
+        //the depression. The outlet is the higher of the two.
+        auto out_cell = ci;    //Pretend focal cell is the outlet
+        auto out_elev = celev; //Note its height
+
+        if(dem(ni)>out_elev){  //Check to see if we were wrong and the neighbour cell is higher.
+          out_cell = ni;       //Neighbour cell was higher. Note it.
+          out_elev = dem(ni);  //Note neighbour's elevation
+        }
+
+        //We've found an outlet between two depressions. Now we need to
+        //determine if it is the lowest outlet between the two.
+
+        //Even though we pull cells off of the priority queue in order of
+        //increasing elevation, we can still add a link between depressions that
+        //is not as low as it could be. This can happen at saddle points, for
+        //instance, consider the cells A-H and their corresponding elevations.
+        //Cells in parantheses are in a neighbouring depression
+        //     (B) (C) (D)   (256) (197) (329)
+        //      A   X   E     228    X    319
+        //      H   G   F     255   184   254
+
+        //In this case we are at Cell X. Cells B, C, and D have previously been
+        //added. Cell G has added X and X has just been popped. X considers its
+        //neighbours in order from A to H. It finds B, at elevation 256, and
+        //makes a note that its depression links with B's depression. It then
+        //sees Cell C, which is in the same depression as B, and has to update
+        //the outlet information between the two depressions.
+
+        const OutletLink olink(clabel,nlabel);      //Create outlet link (order of clabel and nlabel doesn't matter)
+        if(outlet_database.count(olink)!=0){        //Determine if the outlet is already present
+          auto &outlet = outlet_database.at(olink); //It was. Use the outlet link to get the outlet information
+          if(outlet.out_elev>out_elev){             //Is the previously stored link higher than the new one?
+            outlet.out_cell = out_cell;             //Yes. So update the link with new outlet cell
+            outlet.out_elev = out_elev;             //Also, update the outlet's elevation
+          }
+        } else {                                    //No preexisting link found; create a new one
+          outlet_database[olink] = Outlet<elev_t>(clabel,nlabel,out_cell,out_elev);
+        }
+      }
+
+    }
+  }
+  progress.stop();
+  RDLOG_TIME_USE<<"t Outlets found in = "<<progress.time_it_took()<<" s";
+
+  //At this point every cell is associated with the label of a depression. Each
+  //depression contains the cells lower than its outlet elevation as well as all
+  //cells whose flow ultimately terminates somewhere within the depression. The
+  //next order of business is to determine which depressions flow into which
+  //other depressions. That is, we need to build a hierarchy of depressions.
+
+  //Since outlets link two depressions, any time we have an outlet we can form a
+  //meta-depression whose two children most both fill before the meta-depression
+  //itself can spill. This meta-depression has an outlet which differs from
+  //either of its children.
+
+  //We can identify the cells belonging to a meta-depression as those which are
+  //less than or equal to its outlet elevation, but greater than the elevation
+  //of the outlet linking its two children.
+
+  //Since each depression has one and only one parent and at most two children,
+  //the depressions will form a binary tree.
+
+  //In order to build the depression hierarchy, it is convenient to visit
+  //outlets from lowest to highest.
+
+  //Since the `unordered_set` is, well, unordered. We create a vector into which
+  //we will move all of the outlets so we can sort them by elevation. Note that
+  //this temporarily doubles the memory required by the program. TODO: Is there
+  //a way to avoid this doubling?
+  std::vector<Outlet<elev_t>> outlets;
+
+  //Pre-size the vector to avoid expensive copy operations as we expand it
+  outlets.reserve(outlet_database.size());
+
+  //Copy the database into the vector
+  for(const auto &o: outlet_database)
+    outlets.push_back(o.second);
+
+  //It's a little difficult to free memory, but this should do it by replacing
+  //the outlet database with an empty database.
+  outlet_database = outletdb_t();
+
+  //Sort outlets in order from lowest to highest. Takes O(N log N) time.
+  std::sort(outlets.begin(), outlets.end(), [](const Outlet<elev_t> &a, const Outlet<elev_t> &b){
+    return a.out_elev<b.out_elev;
+  });
+
+  //TODO: For debugging
+  if(outlets.size()>0){
+    for(unsigned int i=0;i<outlets.size()-1;i++)
+      assert(outlets.at(i).out_elev<=outlets.at(i+1).out_elev);
+  }
+
+  //Now that we have the outlets in order, we'll visit them from lowest to
+  //highest. If two outlets are at the same elevation we visit them in an
+  //arbitrary order. Each outlet we find is the unique lowest connection between
+  //two depressions. We join these depressions to make a meta-depression. The
+  //problem is, once we've formed a meta-depression, there may still be many
+  //outlets which believe they link to one of the child depressions.
+
+  //To deal with this, we use a Disjoint-Set/Union-Find data structure. This
+  //data structure, when passed a depression label as a query, returns the label
+  //of the upper-most meta-depression in the chain of parent depressions
+  //starting at the query label. The Disjoint-Set data structure has some nice
+  //caching properties which, *roughly speaking*, ensure that all queries
+  //execute in O(1) time.
+
+  //Presize the DisjointDenseIntSet to twice the number of depressions. Since we
+  //are building a binary tree the number of leaf nodes is about equal to the
+  //number of non-leaf nodes. The data structure will expand dynamically as
+  //needed.
+  DisjointDenseIntSet djset(depressions.size());
+
+  RDLOG_PROGRESS<<"p Constructing hierarchy from outlets...";
+
+  //Visit outlets in order of elevation from lowest to highest. If two outlets
+  //are at the same elevation, choose one arbitrarily.
+  progress.start(outlets.size());
+  for(auto &outlet: outlets){
+    ++progress;
+    auto depa_set = djset.findSet(outlet.depa); //Find the ultimate parent of Depression A
+    auto depb_set = djset.findSet(outlet.depb); //Find the ultimate parent of Depression B
+
+    //If the depressions are already part of the same meta-depression, then
+    //nothing needs to be done.
+    if(depa_set==depb_set)
+      continue; //Do nothing, move on to the next highest outlet
+
+
+    if(depa_set==OCEAN || depb_set==OCEAN){
+      //If we're here then both depressions cannot link to the ocean, since we
+      //would have used `continue` above. Therefore, one and only one of them
+      //links to the ocean. We swap them to ensure that `depb` is the one which
+      //links to the ocean.
+      if(depa_set==OCEAN){
+        std::swap(outlet.depa, outlet.depb);
+        std::swap(depa_set, depb_set);
+      }
+
+      //We now have four values, the Depression A Label, the Depression B Label,
+      //the Depression A MetaLabel, and the Depression B MetaLabel. We know that
+      //the Depression B MetaLabel is OCEAN. Depression B Label is the label of
+      //the actual depression this outlet links to, not the meta-depressions of
+      //which it is a part. Depression A MetaLabel is the meta-depression that
+      //has just found a path to the ocean via Depression B. Depression A Label
+      //is some value we don't care about.
+
+      //What we will do is link Depression A MetaLabel to Depression B.
+      //Depression B ultimately terminates in the ocean, but the only way to get
+      //there in real-life is to crawl into Depression B, not into its meta-
+      //depression. At this point its meta-depression is the ocean, so crawling
+      //into the meta-depression would form a direct link to the ocean, which is
+      //not realistic.
+
+      //Get a reference to Depression A MetaLabel.
+      auto &dep = depressions.at(depa_set);
+
+      //If this depression has already found the ocean then don't merge it
+      //again. (TODO: Richard)
+      // if(dep.out_cell==OCEAN)
+        // continue;
+
+      //Ensure we don't modify depressions that have already found their paths
+      assert(dep.out_cell==NO_VALUE);
+      assert(dep.odep==NO_VALUE);
+
+      //Point this depression to the ocean through Depression B Label
+      dep.parent       = outlet.depb;        //Set Depression Meta(A) parent
+      dep.out_elev     = outlet.out_elev;    //Set Depression Meta(A) outlet elevation
+      dep.out_cell     = outlet.out_cell;    //Set Depression Meta(A) outlet cell index
+      dep.odep         = NO_VALUE;           //Since this is an ocean link, A has now overflow depression
+      dep.ocean_parent = true;
+      dep.geolink      = outlet.depb;        //Metadepression(A) overflows, geographically, into Depression B
+      depressions.at(outlet.depb).ocean_linked.emplace_back(depa_set);
+      djset.mergeAintoB(depa_set,OCEAN); //Make a note that Depression A MetaLabel has a path to the ocean
+    } else {
+      //Neither depression has found the ocean, so we merge the two depressions
+      //into a new depression.
+      auto &depa          = depressions.at(depa_set); //Reference to Depression A MetaLabel
+      auto &depb          = depressions.at(depb_set); //Reference to Depression B MetaLabel
+
+      //Ensure we haven't already given these depressions outlet information
+      assert(depa.odep==NO_VALUE);
+      assert(depb.odep==NO_VALUE);
+
+      const auto newlabel = depressions.size();       //Label of A and B's new parent depression
+      depa.parent   = newlabel;        //Set Meta(A)'s parent to be the new meta-depression
+      depb.parent   = newlabel;        //Set Meta(B)'s parent to be the new meta-depression
+      depa.out_cell = outlet.out_cell; //Note that this is Meta(A)'s outlet
+      depb.out_cell = outlet.out_cell; //Note that this is Meta(B)'s outlet
+      depa.out_elev = outlet.out_elev; //Note that this is Meta(A)'s outlet's elevation
+      depb.out_elev = outlet.out_elev; //Note that this is Meta(B)'s outlet's elevation
+      depa.odep     = depb_set;        //Note that Meta(A) overflows, logically, into Meta(B)
+      depb.odep     = depa_set;        //Note that Meta(B) overflows, logically, into Meta(A)
+      depa.geolink  = outlet.depb;     //Meta(A) overflows, geographically, into B
+      depb.geolink  = outlet.depa;     //Meta(B) overflows, geographically, into A
+
+      //Be sure that this happens AFTER we are done using the `depa` and `depb`
+      //references since they will be invalidated if `depressions` has to
+      //resize!
+      const auto depa_pitcell_temp = depa.pit_cell;
+
+      auto &newdep     = depressions.emplace_back();                                                                       //is it right to create a new depression for the metadepression like this?
+      newdep.lchild    = depa_set;
+      newdep.rchild    = depb_set;
+      newdep.dep_label = newlabel;
+      newdep.pit_cell  = depa_pitcell_temp;
+
+      djset.mergeAintoB(depa_set, newlabel); //A has a parent now
+      djset.mergeAintoB(depb_set, newlabel); //B has a parent now
+    }
+  }
+  progress.stop();
+
+  RDLOG_TIME_USE<<"t Time to construct Depression Hierarchy = "<<timer_dephier.stop()<<" s";
+
+
+  //At this point we have a 2D array in which each cell is labeled. This label
+  //corresponds to either the root node (the ocean) or a leaf node of a binary
+  //tree representing the hierarchy of depressions.
+
+  //The labels array has been modified in place. The depression hierarchy is
+  //returned.
+
+  Timer timer_volumes;
+  timer_volumes.start();
+
+  CalculateMarginalVolumes(depressions, dem, label);
+
+  CalculateTotalVolumes(depressions);
+
+  RDLOG_TIME_USE<<"t Time to calculate volumes = "<<timer_volumes.stop()<<" s";
+  RDLOG_TIME_USE<<"t Total time in depression hierarchy calculations = "<<timer_overall.stop()<<" s";
+
+  return depressions;
+}
+
+
+
+template<class elev_t>
+void CalculateMarginalVolumes(
+  DepressionHierarchy<elev_t> &deps,
+  const Array2D<elev_t>       &dem,
+  const Array2D<dh_label_t>   &label
+){
+  ProgressBar progress;
+
+  RDLOG_PROGRESS<<"p Calculating depression marginal volumes...";
+
+  //Get the marginal depression cell counts and total elevations
+  progress.start(dem.size());
+  #pragma omp parallel default(none) shared(progress,deps,dem,label)
+  {
+    std::vector<uint32_t> cell_counts     (deps.size(), 0);
+    std::vector<double>   total_elevations(deps.size(), 0);
+
+    #pragma omp for
+    for(unsigned int i=0;i<dem.size();i++){
+      ++progress;
+      const auto my_elev = dem(i);
+      auto clabel        = label(i);
+
+      while(clabel!=OCEAN && my_elev>deps.at(clabel).out_elev)
+        clabel = deps[clabel].parent;
+
+      if(clabel==OCEAN)
+        continue;
+
+      cell_counts[clabel]++;
+      total_elevations[clabel] += dem(i);
+    }
+
+    #pragma omp critical
+    for(unsigned int i=0;i<deps.size();i++){
+      deps[i].cell_count      += cell_counts[i];
+      deps[i].total_elevation += total_elevations[i];
+    }
+  }
+  progress.stop();
+}
+
+
+
+template<class elev_t>
+void CalculateTotalVolumes(
+  DepressionHierarchy<elev_t> &deps
+){
+  ProgressBar progress;
+
+  RDLOG_PROGRESS<<"p Calculating depression total volumes...";
+  //Calculate total depression volumes and cell counts
+  progress.start(deps.size());
+  for(unsigned int d=0;d<deps.size();d++){
+    ++progress;
+
+    auto &dep = deps.at(d);
+    if(dep.lchild!=NO_VALUE){
+      assert(dep.rchild!=NO_VALUE); //Either no children or two children
+      assert(dep.lchild<d);         //ID of child must be smaller than parent's
+      assert(dep.rchild<d);         //ID of child must be smaller than parent's
+      dep.cell_count      += deps.at(dep.lchild).cell_count;
+      dep.total_elevation += deps.at(dep.lchild).total_elevation;
+      dep.cell_count      += deps.at(dep.rchild).cell_count;
+      dep.total_elevation += deps.at(dep.rchild).total_elevation;
+    }
+    //This has to be after the foregoing because the cells added by the if-
+    //clauses have additional volume above their spill elevations that cannot be
+    //counted simply by adding their volumes to their parent depression.
+    dep.dep_vol = dep.cell_count*static_cast<double>(dep.out_elev)-dep.total_elevation;
+
+    assert(dep.lchild==NO_VALUE || fp_le(deps.at(dep.lchild).dep_vol+deps.at(dep.rchild).dep_vol,dep.dep_vol));
+  }
+  progress.stop();
+}
+
+
+
+//Utility function for doing various relabelings based on the depression
+//hierarchy.
+template<class elev_t>
+void LastLayer(Array2D<dh_label_t> &label, const Array2D<elev_t> &dem, const DepressionHierarchy<elev_t> &depressions){
+  #pragma omp parallel for collapse(2)
+  for(int y=0;y<label.height();y++)
+  for(int x=0;x<label.width();x++){
+    auto mylabel = label(x,y);
+    while(true){
+      if(dem(x,y)>=depressions.at(mylabel).out_elev)
+        mylabel = depressions.at(mylabel).parent;
+      else {
+        if(mylabel!=0)
+          mylabel = -3;
+        break;
+      }
+    }
+    label(x,y) = mylabel;
+  }
+}
+
+}