Fix highest peaks being filtered from bathymetry (#44)

* Fix highest peaks being filtered from bathymetry.

We were assuming that prominence=elevation for the highest peak
in a tree.  When the peak's elevation is negative, as in bathymetry,
this causes the prominence of the highest peak to be negative,
which means the highest peak will always be filtered out.

* Update README.md

Author: akirmse

README.md

@@ -355,6 +355,16 @@ The input divide tree must be free of runoffs (see the options to ```merge_divid
 parent, and its line parent on each line.  Landmass high points (where the prominence is equal to the elevation) are not included.
 Their key saddles are the ocean, and there isn't a well-defined way to connect such peaks to other land masses through the divide tree.
 
+## Bathymetry
+
+For "regular" terrain on the Earth, we set the prominence of the highest peak equal to its elevation by convention.  This obviously 
+doesn't work for bathymetry (terrain below sea level).  It's not clear what the prominence of such underwater peaks should be,
+but clearly we don't want to filter them away by assigning them negative prominence values.
+
+If the ```-b``` flag is set on the ```prominence``` and ```merge_divide_trees```, we compute the prominence of the highest peak to
+be its elevation, minus the elevation of the lowest saddle in the divide tree.  This guarantees that the highest peak will have
+the highest prominence.  Setting the ```--bathymetry``` flag on ```run_prominence.py``` will also turn on this behavior.
+
 ## Anti-prominence
 
 The "anti-prominence" of low points can be computed by the same algorithm, simply by changing

code/compute_parents.cpp

@@ -112,7 +112,7 @@ int main(int argc, char **argv) {
 
   // Build island tree to get prominence values
   IslandTree *islandTree = new IslandTree(*divideTree);
-  islandTree->build();
+  islandTree->build(false);  // TODO: Flag for bathymetry?
 
   // Build line parent tree
   LineTree *lineTree = new LineTree(*divideTree);

code/island_tree.cpp

@@ -36,7 +36,7 @@ IslandTree::IslandTree(const DivideTree &divideTree) :
     mDivideTree(divideTree) {
 }
 
-void IslandTree::build() { 
+void IslandTree::build(bool isBathymetry) { 
   // Initialize topology
   mNodes.resize(mDivideTree.nodes().size());
   int index = 1;  // Peaks are 1-indexed; put in a dummy node 0
@@ -54,7 +54,7 @@ void IslandTree::build() {
   uninvertPeaks();
   uninvertSaddles();
 
-  computeProminences();
+  computeProminences(isBathymetry);
 }
 
 // Sort peaks so that parent is always higher elevation.
@@ -152,7 +152,7 @@ void IslandTree::uninvertSaddle(int nodeId) {
   }
 }
 
-void IslandTree::computeProminences() {
+void IslandTree::computeProminences(bool isBathymetry) {
   for (int i = 1; i < (int) mNodes.size(); ++i) {
     Elevation elev = getPeak(i).elevation;
     int childNodeId = i;
@@ -177,8 +177,17 @@ void IslandTree::computeProminences() {
     }
 
     if (parentNodeId == Node::Null) {
-      // This is the highest point in the tree
-      mNodes[i].prominence = elev;
+      // This is the highest point in the tree.  In "normal" circumstances where
+      // all the elevations are nonnegative, we can safely define the prominence as equal
+      // to the elevation.  But for bathymetry the prominence is not so clear.  We
+      // use the elevation of the minimum saddle as "sea level".  This should agree
+      // with the traditional definition (prominence = elevation) for regular terrain,
+      // and for bathymetry, it will guarantee that the highest point has the
+      // highest prominence, which matches intuition.
+      //
+      // Other definitions of the prominence of the highest point are possible,
+      // for example, height above the minimum sample value.
+      mNodes[i].prominence = elev - getSeaLevelValue(isBathymetry);
     } else {
       // Prominence = peak elevation - key col elevation
       int saddlePeakId = mNodes[childNodeId].saddlePeakId;
@@ -189,6 +198,19 @@ void IslandTree::computeProminences() {
   }
 }
 
+Elevation IslandTree::getSeaLevelValue(bool isBathymetry) {
+  // If this is bathymetry, find the minimum saddle elevation in the tree.
+  if (isBathymetry && !mNodes.empty()) {
+    Elevation elevation = std::numeric_limits<Elevation>::max();
+    for (auto const &saddle : mDivideTree.saddles()) {
+      elevation = std::min(elevation, saddle.elevation);
+    }
+    return elevation;
+  }
+  
+  return 0;
+}
+
 const Peak &IslandTree::getPeak(int peakId) const {
   return mDivideTree.peaks()[peakId - 1];  // 1-indexed
 }

code/island_tree.h

@@ -49,7 +49,8 @@ class IslandTree {
 
   explicit IslandTree(const DivideTree &divideTree);
 
-  void build();
+  // If isBathymetry is true, don't assume sea level = 0.
+  void build(bool isBathymetry);
 
   bool writeToFile(const std::string &filename) const;
 
@@ -66,11 +67,15 @@ class IslandTree {
   // Sort peaks by increasing divide tree saddle elevation
   void uninvertSaddles();
 
-  void computeProminences();
+  void computeProminences(bool isBathymetry);
 
   void uninvertPeak(int nodeId);
   void uninvertSaddle(int nodeId);
 
+  // Return the elevation of "sea level", i.e. the elevation used as the
+  // base of the highest peak in the tree.
+  Elevation getSeaLevelValue(bool isBathymetry);
+
   // Indices start at 1; use these helper functions to deal with offset.
   const Peak &getPeak(int peakId) const;
   const Saddle &getSaddle(int saddleId) const;

code/merge_divide_trees.cpp

@@ -48,6 +48,7 @@ static void usage() {
   printf("                    in same units as divide tree, default = 100\n");
   printf("  -s scale_factor   Scale output elevations by this (e.g. meters to feet), default = 1\n");
   printf("  -t num_threads    Number of threads to use, default = 1\n");
+  printf("  -b                Input data is bathymetric (do not use sea level)\n");
   exit(1);
 }
 
@@ -83,6 +84,7 @@ int main(int argc, char **argv) {
   Elevation minProminence = 100;
   bool finalize = false;
   bool flipElevations = false;
+  bool bathymetry = false;
   int numThreads = 1;
   float elevationScale = 1.0f;
 
@@ -96,12 +98,16 @@ int main(int argc, char **argv) {
     {"v", required_argument, nullptr, 0},
     {nullptr, 0, 0, 0},
   };
-  while ((ch = getopt_long(argc, argv, "afm:s:t:", long_options, nullptr)) != -1) {
+  while ((ch = getopt_long(argc, argv, "abfm:s:t:", long_options, nullptr)) != -1) {
     switch (ch) {
     case 'a':
       flipElevations = true;
       break;
-      
+
+    case 'b':
+      bathymetry = true;
+      break;
+
     case 'f':
       finalize = true;
       break;
@@ -183,7 +189,7 @@ int main(int argc, char **argv) {
   VLOG(1) << "Building prominence island tree";
 
   IslandTree *unprunedIslandTree = new IslandTree(*divideTree);
-  unprunedIslandTree->build();
+  unprunedIslandTree->build(bathymetry);
 
   if (finalize) {
     divideTree->deleteRunoffs();  // Does not affect island tree
@@ -211,7 +217,7 @@ int main(int argc, char **argv) {
 
   // Build new island tree on pruned divide tree to get final prominence values
   IslandTree *prunedIslandTree = new IslandTree(*divideTree);
-  prunedIslandTree->build();
+  prunedIslandTree->build(bathymetry);
 
   // Write final prominence value table
   string filename = outputFilename + ".txt";

code/prominence.cpp

@@ -62,6 +62,7 @@ static void usage() {
   printf("  -t num_threads    Number of threads, default = 1\n");
   printf("  -z                UTM zone (if input data is in UTM)\n");
   printf("  -a                Compute anti-prominence instead of prominence\n");
+  printf("  -b                Input DEM is bathymetric (do not use sea level)\n");
   printf("  --kml             Generate KML output of divide tree\n"); 
   exit(1);
 }
@@ -80,6 +81,7 @@ int main(int argc, char **argv) {
   int ch;
   string str;
   bool antiprominence = false;
+  bool bathymetry = false;
   int writeKml = false;
   int utmZone = NO_UTM_ZONE;
 
@@ -91,12 +93,16 @@ int main(int argc, char **argv) {
     {"kml", no_argument, &writeKml, 1},
     {nullptr, 0, 0, 0},
   };
-  while ((ch = getopt_long(argc, argv, "af:i:k:m:o:t:z:", long_options, nullptr)) != -1) {
+  while ((ch = getopt_long(argc, argv, "abf:i:k:m:o:t:z:", long_options, nullptr)) != -1) {
     switch (ch) {
     case 'a':
       antiprominence = true;
       break;
-      
+
+    case 'b':
+      bathymetry = true;
+      break;
+
     case 'f': {
       auto format = FileFormat::fromName(optarg);
       if (format == nullptr) {
@@ -183,6 +189,7 @@ int main(int argc, char **argv) {
 
   // Don't write out unpruned divide tree--it's too large and slow
   ProminenceOptions options = {output_directory, minProminence, false, antiprominence,
+                               bathymetry,
                                static_cast<bool>(writeKml)};
 
   // Caching doesn't do anything for our calculation and the tiles are huge

code/prominence_task.cpp

@@ -86,7 +86,7 @@ bool ProminenceTask::run(double lat, double lng, const CoordinateSystem &coordin
   VLOG(1) << "Pruning divide tree to " << mOptions.minProminence << " prominence";
 
   IslandTree islandTree(*divideTree);
-  islandTree.build();
+  islandTree.build(mOptions.bathymetry);
 
   divideTree->prune(mOptions.minProminence, islandTree);
 

code/prominence_task.h

@@ -43,6 +43,9 @@ struct ProminenceOptions {
   // or anti-prominence, which is the "prominence" of low points.
   bool antiprominence;
 
+  // Data is bathymetry; do not assume sea level = 0 for prom calculations.
+  bool bathymetry;
+
   // Write KML of pruned divide tree?
   bool writeKml;
 };

scripts/run_prominence.py

@@ -219,6 +219,8 @@ def main():
                         help='Size of reprojected tiles to generate')
     parser.add_argument('--samples_per_tile', type=int, default=10000,
                         help='Number of samples per edge of a reprojected tile')
+    parser.add_argument('--bathymetry', action=argparse.BooleanOptionalAction,
+                        help='Is the input DEM bathymetry?')
     args = parser.parse_args()
 
     gdal.UseExceptions()
@@ -304,9 +306,13 @@ def main():
     # Run prominence and merge_divide_trees
     print("Running prominence")
     prominence_binary = os.path.join(args.binary_dir, "prominence")
+    extra_args = ""
+    if args.bathymetry:
+        extra_args += " -b"
     prom_command = f"{prominence_binary} --v=1 -f CUSTOM-{args.degrees_per_tile}-{args.samples_per_tile}" + \
         f" -i {args.tile_dir} -o {args.output_dir}" + \
         f" -t {args.threads} -m {args.min_prominence}" + \
+        extra_args + \
         f" -- {ymin} {ymax} {xmin} {xmax}"
     run_command(prom_command)
 
@@ -317,8 +323,12 @@ def main():
     units_multiplier = 1
     if args.output_units == "feet":
         units_multiplier = 3.28084
+    extra_args = ""
+    if args.bathymetry:
+        extra_args += " -b"
     merge_command = f"{merge_binary} --v=1 -f" + \
         f" -t {args.threads} -m {args.min_prominence} -s {units_multiplier}" + \
+        extra_args + \
         f" {merge_output} {merge_inputs}"
     run_command(merge_command)