Add FSM to RichDEM

apps/CMakeLists.txt

@@ -8,6 +8,7 @@ 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_fill_spill_merge.exe        rd_fill_spill_merge.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)
@@ -30,6 +31,7 @@ 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_fill_spill_merge.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)

apps/rd_fill_spill_merge.cpp

@@ -0,0 +1,104 @@
+#include <richdem/common/Array2D.hpp>
+#include <richdem/common/gdal.hpp>
+#include <richdem/depressions/fill_spill_merge.hpp>
+#include <richdem/misc/misc_methods.hpp>
+#include <richdem/ui/cli_options.hpp>
+
+#include <iostream>
+#include <string>
+#include <stdexcept>
+
+namespace rd = richdem;
+namespace dh = richdem::dephier;
+
+int main(int argc, char **argv){
+  CLI::App app("Fill-Spill-Merge Example Program");
+
+  std::string topography_filename;
+  std::string output_prefix;
+  double      surface_water_level = std::numeric_limits<double>::quiet_NaN();
+  std::string surface_water_filename;
+  double      ocean_level;
+
+  size_t num;
+  size_t len;
+  app.add_option("topography", topography_filename, "Topography to run FSM on")->required();
+  app.add_option("output",     output_prefix,       "Path of GeoTiff output file")->required();
+  const auto swl_ptr = app.add_option("--swl", surface_water_level, "Surface water level as a numeric constant");
+  app.add_option("--swf", surface_water_filename, "File containing surface water levels")->excludes(swl_ptr);
+  app.add_option("ocean_level",         ocean_level, "Elevation of the ocean")->required();
+
+  CLI11_PARSE(app, argc, argv);
+
+  std::cout<<"m Input DEM           = "<<topography_filename   <<std::endl;
+  std::cout<<"m Output prefix       = "<<output_prefix         <<std::endl;
+  std::cout<<"m Surface water level = "<<surface_water_level   <<std::endl;
+  std::cout<<"m Surface water file  = "<<surface_water_filename<<std::endl;
+  std::cout<<"m Ocean level         = "<<ocean_level           <<std::endl;
+
+  rd::Timer timer_io;
+  timer_io.start();
+  rd::Array2D<double> topo(topography_filename);
+  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<double> wtd;
+  if(surface_water_filename.empty()){
+    //All cells have the same amount of water
+    wtd.resize(topo.width(), topo.height(), surface_water_level);
+  } else {
+    wtd = rd::Array2D<double>(surface_water_filename);
+  }
+
+  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
+
+  wtd.setNoData(topo.noData());
+
+  //Label the ocean cells. This is a precondition for using
+  //`GetDepressionHierarchy()`.
+  rd::BucketFillFromEdges<rd::Topology::D8>(topo, label, ocean_level, dh::OCEAN);
+
+  //Make NoData cells also ocean cells. Ocean has no water on it to begin with.
+  #pragma omp parallel for
+  for(unsigned int i=0;i<label.size();i++){
+    if(topo.isNoData(i) || label(i)==dh::OCEAN){
+      label(i) = dh::OCEAN;
+      wtd  (i) = 0;
+    }
+  }
+
+  rd::Timer timer_calc;
+  timer_calc.start();
+
+  //Generate flow directions, label all the depressions, and get the hierarchy
+  //connecting them
+  auto deps = dh::GetDepressionHierarchy<double,rd::Topology::D8>(topo, label, flowdirs);
+
+  dh::FillSpillMerge(topo, label, flowdirs, deps, wtd);
+
+  timer_calc.stop();
+
+  timer_io.start();
+
+  //Output the water table depth
+  wtd.saveGDAL(output_prefix+"-wtd.tif");
+
+  for(unsigned int i=0;i<topo.size();i++)
+    if(!topo.isNoData(i))
+      wtd(i) += topo(i);
+
+  //Output the new height of the hydraulic surface
+  wtd.saveGDAL(output_prefix+"-hydrologic-surface-height.tif");
+
+  timer_io.stop();
+
+  std::cout<<"Finished."<<std::endl;
+  std::cout<<"IO time   = "<<timer_io.accumulated()  <<" s"<<std::endl;
+  std::cout<<"Calc time = "<<timer_calc.accumulated()<<" s"<<std::endl;
+
+  return 0;
+}

include/doctest.h

@@ -652,7 +652,7 @@ namespace assertType {
         DT_WARN_THROWS_WITH    = is_throws_with | is_warn,
         DT_CHECK_THROWS_WITH   = is_throws_with | is_check,
         DT_REQUIRE_THROWS_WITH = is_throws_with | is_require,
-        
+
         DT_WARN_THROWS_WITH_AS    = is_throws_with | is_throws_as | is_warn,
         DT_CHECK_THROWS_WITH_AS   = is_throws_with | is_throws_as | is_check,
         DT_REQUIRE_THROWS_WITH_AS = is_throws_with | is_throws_as | is_require,
@@ -828,9 +828,9 @@ namespace detail {
     template<class T> struct remove_reference<T&>  { typedef T type; };
     template<class T> struct remove_reference<T&&> { typedef T type; };
 
-    template<typename T, typename U = T&&> U declval(int); 
+    template<typename T, typename U = T&&> U declval(int);
 
-    template<typename T> T declval(long); 
+    template<typename T> T declval(long);
 
     template<typename T> auto declval() DOCTEST_NOEXCEPT -> decltype(declval<T>(0)) ;
 
@@ -1683,7 +1683,7 @@ DOCTEST_CLANG_SUPPRESS_WARNING_POP
     DOCTEST_INTERFACE void toStream(std::ostream* s, int long long in);
     DOCTEST_INTERFACE void toStream(std::ostream* s, int long long unsigned in);
 
-    // ContextScope base class used to allow implementing methods of ContextScope 
+    // ContextScope base class used to allow implementing methods of ContextScope
     // that don't depend on the template parameter in doctest.cpp.
     class DOCTEST_INTERFACE ContextScopeBase : public IContextScope {
     protected:
@@ -1742,7 +1742,7 @@ DOCTEST_CLANG_SUPPRESS_WARNING_POP
         bool log();
         void react();
     };
-    
+
     template <typename L>
     ContextScope<L> MakeContextScope(const L &lambda) {
         return ContextScope<L>(lambda);
@@ -3192,7 +3192,7 @@ namespace detail {
 
 namespace timer_large_integer
 {
-    
+
 #if defined(DOCTEST_PLATFORM_WINDOWS)
     typedef ULONGLONG type;
 #else // DOCTEST_PLATFORM_WINDOWS
@@ -3949,7 +3949,7 @@ namespace detail {
             if(matchesAny(m_signature.m_name.c_str(), s->filters[7], false, s->case_sensitive))
                 return;
         }
-        
+
         // if a Subcase on the same level has already been entered
         if(s->subcasesStack.size() < size_t(s->subcasesCurrentMaxLevel)) {
             s->should_reenter = true;
@@ -4754,7 +4754,7 @@ namespace detail {
             m_string = tlssPop();
             logged = true;
         }
-        
+
         DOCTEST_ITERATE_THROUGH_REPORTERS(log_message, *this);
 
         const bool isWarn = m_severity & assertType::is_warn;
@@ -5303,7 +5303,7 @@ namespace {
             test_case_start_impl(in);
             xml.ensureTagClosed();
         }
-        
+
         void test_case_reenter(const TestCaseData&) override {}
 
         void test_case_end(const CurrentTestCaseStats& st) override {
@@ -6019,7 +6019,7 @@ namespace {
             subcasesStack.clear();
             currentSubcaseLevel = 0;
         }
-        
+
         void test_case_reenter(const TestCaseData&) override {
             subcasesStack.clear();
         }
@@ -6703,7 +6703,7 @@ int Context::run() {
             DOCTEST_ITERATE_THROUGH_REPORTERS(test_case_start, tc);
 
             p->timer.start();
-            
+
             bool run_test = true;
 
             do {
@@ -6743,7 +6743,7 @@ DOCTEST_MSVC_SUPPRESS_WARNING_POP
                     run_test = false;
                     p->failure_flags |= TestCaseFailureReason::TooManyFailedAsserts;
                 }
-                
+
                 if(p->should_reenter && run_test)
                     DOCTEST_ITERATE_THROUGH_REPORTERS(test_case_reenter, tc);
                 if(!p->should_reenter)

tests/CMakeLists.txt

@@ -1,12 +1,12 @@
 cmake_minimum_required (VERSION 3.9)
 
-project (dephier)
+project (richdem_tests)
 
 add_executable(richdem_tests.exe
+  fsm_tests.cpp
+  test_main.cpp
   tests.cpp
 )
 
-target_link_libraries(richdem_tests.exe
-  PRIVATE 
-    richdem
-)
+target_link_libraries(richdem_tests.exe PRIVATE richdem)
+target_compile_definitions(richdem_tests.exe PRIVATE RICHDEM_NO_PROGRESS)