BUGFIX: Stop MV search at ocean link

When calculating marginal volumes for depressions we sum the total elevation of the cells contained within the depression.

If a cell's elevation is higher than the outlet for the depression it's in, we go to the parent depression and make the same check there. This continues until we reach the ocean. But that means that a cell's elevation can contribute to a depression we reach through an ocean link! That is, if a cell's elevation exceeds that of any outlet of a depression it's part of, then we continue looking for a place to add its elevation in downhill depressions. But this cell cannot contribute marginal volume to any such depression because no such depression contains it.

This bug likely doesn't manifest much (if ever?) because the outlets of downhill depressions are almost certainly lower than the elevation of the cell - after all, we're looking in downhill depressions because the cell's elevation was higher than that of its outlet depression and the downhill depression is downhill because its outlet is lower.

Maybe the only way this would ever be an issue would be if outlets were somehow equal? But in that case I think (it's been a while since I've run this through in my head) that they would be part of the same meta-depression.

So maybe this never actually manifests as an issue, but it does mean that we can stop looking for places to add marginal volume early! That'll be a perf improvement.

An example of this issue is here:

21(36.000000) -> 30(36.000000)  -> 31(37.000000) -> 34(41.000000)  ->
                 38(43.000000) o-> 16(41.000000) -> 36(42.000000) o->
                 10(35.000000) o-> 0(inf)

A cell with an elevation of 78 doesn't contribute marginal volume to depression 21 because 21's outlet elevation is 36.

So we look at 36's parent, Dep 30, but the outlet elevation there is also 36. We continue until we reach Dep 38, with an outlet elevation of 43. Now we follow an oceanlink to Dep 16, whose outlet elevation is still too low. We expect this because Dep 16 is downhill of Dep 38. At this point we can continue following the hierarchy to the ocean, but that's unnecessary both because it would be incorrect to add the cell's marginal volume to any depression past 38 but also because all of the downhill outlets will be too low.

include/richdem/depressions/depression_hierarchy.hpp

@@ -746,11 +746,34 @@ void CalculateMarginalVolumes(
       const auto my_elev = dem(i);
       auto clabel        = label(i);
 
-      while (clabel != OCEAN && my_elev > deps.at(clabel).out_elev)
+      // This cell contributes volume to some depression or meta-depression.
+      // Let's find out which. To do so, we start with the leaf depression the
+      // cell is part of. If it's not contained by that depression, we look at
+      // the meta-depression and work our way up until we find a depression in
+      // the hierarchy whose outlet elevation is higher than the cell's
+      // elevation. If at any point we pass through an oceanlink, our search
+      // terminates: any water falling on the cell will ultimately end up in the
+      // ocean and the cell is not contained by any depression.
+      while (clabel != OCEAN && my_elev > deps.at(clabel).out_elev) {
+        if (deps[clabel].ocean_parent){
+          // There is no meta-depression: the next depression in the hierarchy
+          // is downhill of the one we started in. Therefore, our marginal
+          // volume is contributed to the ocean.
+          clabel = OCEAN;
+          break;
+        }
+
+        // We weren't contained in this depression, try its parent (the
+        // meta-depression)
         clabel = deps[clabel].parent;
+      }
 
-      if (clabel == OCEAN)
+      // This cell was the above the outlet of all the metadepressions that
+      // could contain it, so its flow goes to the ocean. Therefore, we ignore
+      // it because it doesn't contribute marginal volume to any depression.
+      if(clabel == OCEAN){
         continue;
+      }
 
       cell_counts[clabel]++;
       total_elevations[clabel] += dem(i);