Improve readability of calculate_subwatershed_boundary

src/pygeoprocessing/extensions.pxd

@@ -39,6 +39,7 @@ cdef extern from "extensions/ManagedRaster.h":
         int block_ybits
         long raster_x_size
         long raster_y_size
+        long n_pixels
         int block_nx
         int block_ny
         int write_mode
@@ -53,6 +54,9 @@ cdef extern from "extensions/ManagedRaster.h":
         double get(long xi, long yi)
         void _load_block(int block_index) except *
         void close()
+        bint is_out_of_bounds(long x, long y)
+        bint is_nodata(double val)
+        bint is_nodata(long x, long y)
 
     cdef cppclass ManagedFlowDirRaster[T]:
         LRUCache[int, double*]* lru_cache
@@ -80,6 +84,9 @@ cdef extern from "extensions/ManagedRaster.h":
         void set(long xi, long yi, double value)
         double get(long xi, long yi)
         void close()
+        bint is_out_of_bounds(long x, long y)
+        bint is_nodata(double val)
+        bint is_nodata(long x, long y)
 
     cdef cppclass D8
     cdef cppclass MFD

src/pygeoprocessing/extensions/ManagedRaster.h

@@ -73,6 +73,17 @@ static void log_msg(LogLevel level, string msg)
   Py_DECREF(pyString);
 }
 
+// Note: I was concerned that checking each value for nan would be too slow, but
+// I compared its performance against another implementation where we first reclassify
+// nans to a regular float, and then skip the nan check, and that was much slower:
+// https://github.com/natcap/invest/issues/1714#issuecomment-2762134419
+inline bool is_close(double x, double y) {
+  if (isnan(x) and isnan(y)) {
+    return true;
+  }
+  return abs(x - y) <= (pow(10, -8) + pow(10, -05) * abs(y));
+}
+
 class NeighborTuple {
 public:
   int direction, x, y;
@@ -104,6 +115,7 @@ class ManagedRaster {
     long raster_y_size;
     int block_nx;
     int block_ny;
+    long n_pixels;
     char* raster_path;
     int band_id;
     GDALDataset* dataset;
@@ -138,6 +150,7 @@ class ManagedRaster {
 
       raster_x_size = dataset->GetRasterXSize();
       raster_y_size = dataset->GetRasterYSize();
+      n_pixels = raster_x_size * raster_y_size;
 
       if (band_id < 1 or band_id > dataset->GetRasterCount()) {
         throw std::invalid_argument(
@@ -380,6 +393,21 @@ class ManagedRaster {
       delete lru_cache;
       free(actualBlockWidths);
     }
+
+    bool is_out_of_bounds(long x, long y) {
+      if (x < 0 or x >= raster_x_size or y < 0 or y >= raster_y_size) {
+        return true;
+      }
+      return false;
+    }
+
+    bool is_nodata(double val) {
+      return is_close(val, nodata);
+    }
+
+    bool is_nodata(long x, long y) {
+      return is_close(get(x, y), nodata);
+    }
 };
 
 // Represents a flow direction raster, which may be of type MFD or D8
@@ -851,15 +879,4 @@ class UpslopeNeighborsNoDivide: public Neighbors<T> {
   UpslopeNeighborNoDivideIterator<T> end() { return UpslopeNeighborNoDivideIterator<T>(&endVal); }
 };
 
-// Note: I was concerned that checking each value for nan would be too slow, but
-// I compared its performance against another implementation where we first reclassify
-// nans to a regular float, and then skip the nan check, and that was much slower:
-// https://github.com/natcap/invest/issues/1714#issuecomment-2762134419
-inline bool is_close(double x, double y) {
-  if (isnan(x) and isnan(y)) {
-    return true;
-  }
-  return abs(x - y) <= (pow(10, -8) + pow(10, -05) * abs(y));
-}
-
 #endif  // NATCAP_INVEST_MANAGEDRASTER_H_