separate C++ classes for each variant of neighbor pixel iteration

src/natcap/invest/managed_raster/ManagedRaster.h

@@ -433,8 +433,18 @@ class DownslopeNeighborIterator {
             return next();
         }
     }
+};
+
+class DownslopeNeighborIteratorNoSkip: public DownslopeNeighborIterator {
+public:
+
+    DownslopeNeighborIteratorNoSkip() {}
 
-    NeighborTuple next_no_skip() {
+    DownslopeNeighborIteratorNoSkip(ManagedFlowDirRaster managed_raster, int x, int y)
+        : DownslopeNeighborIterator(managed_raster, x, y)   // Call the superclass constructor in the subclass' initialization list.
+        { }
+
+    NeighborTuple next() {
         NeighborTuple n;
         long xj, yj;
         int flow_ij;
@@ -455,7 +465,7 @@ class DownslopeNeighborIterator {
             return n;
         } else {
             n_dir += 1;
-            return next_no_skip();
+            return next();
         }
     }
 };
@@ -518,8 +528,18 @@ class UpslopeNeighborIterator {
             return next();
         }
     }
+};
+
+class UpslopeNeighborIteratorNoDivide: public UpslopeNeighborIterator {
+public:
+
+    UpslopeNeighborIteratorNoDivide() { }
 
-    NeighborTuple next_no_divide() {
+    UpslopeNeighborIteratorNoDivide(ManagedFlowDirRaster managed_raster, int x, int y)
+        : UpslopeNeighborIterator(managed_raster, x, y)
+    {}
+
+    NeighborTuple next() {
 
         NeighborTuple n;
         long xj, yj;
@@ -537,7 +557,7 @@ class UpslopeNeighborIterator {
         if (xj < 0 or xj >= raster.raster_x_size or
                 yj < 0 or yj >= raster.raster_y_size) {
             n_dir += 1;
-            return next_no_divide();
+            return next();
         }
 
         flow_dir_j = raster.get(xj, yj);
@@ -549,11 +569,25 @@ class UpslopeNeighborIterator {
             return n;
         } else {
             n_dir += 1;
-            return next_no_divide();
+            return next();
         }
     }
+};
 
-    NeighborTuple next_skip(int skip) {
+
+class UpslopeNeighborIteratorSkip: public UpslopeNeighborIterator {
+public:
+    int skip;
+
+    UpslopeNeighborIteratorSkip() { }
+
+    UpslopeNeighborIteratorSkip(ManagedFlowDirRaster managed_raster, int x, int y, int skip)
+        : UpslopeNeighborIterator(managed_raster, x, y)
+    {
+        this->skip = skip;
+    }
+
+    NeighborTuple next() {
 
         NeighborTuple n;
         long xj, yj;
@@ -573,7 +607,7 @@ class UpslopeNeighborIterator {
                 yj < 0 or yj >= raster.raster_y_size or
                 INFLOW_OFFSETS[n_dir] == skip) {
             n_dir += 1;
-            return next_skip(skip);
+            return next();
         }
 
         flow_dir_j = raster.get(xj, yj);
@@ -590,11 +624,12 @@ class UpslopeNeighborIterator {
             return n;
         } else {
             n_dir += 1;
-            return next_skip(skip);
+            return next();
         }
     }
 };
 
+
 inline bool is_close(double x, double y) {
     if (isnan(x) and isnan(y)) {
         return true;

src/natcap/invest/managed_raster/managed_raster.pxd

@@ -87,7 +87,18 @@ cdef extern from "ManagedRaster.h":
         DownslopeNeighborIterator()
         DownslopeNeighborIterator(ManagedFlowDirRaster, int, int)
         NeighborTuple next()
-        NeighborTuple next_no_skip()
+
+    cdef cppclass DownslopeNeighborIteratorNoSkip:
+        ManagedFlowDirRaster raster
+        int col
+        int row
+        int n_dir
+        int flow_dir
+        int flow_dir_sum
+
+        DownslopeNeighborIteratorNoSkip()
+        DownslopeNeighborIteratorNoSkip(ManagedFlowDirRaster, int, int)
+        NeighborTuple next()
 
     cdef cppclass UpslopeNeighborIterator:
         ManagedFlowDirRaster raster
@@ -99,8 +110,28 @@ cdef extern from "ManagedRaster.h":
         UpslopeNeighborIterator()
         UpslopeNeighborIterator(ManagedFlowDirRaster, int, int)
         NeighborTuple next()
-        NeighborTuple next_no_divide()
-        NeighborTuple next_skip(int skip)
+
+    cdef cppclass UpslopeNeighborIteratorNoDivide:
+        ManagedFlowDirRaster raster
+        int col
+        int row
+        int n_dir
+        int flow_dir
+
+        UpslopeNeighborIteratorNoDivide()
+        UpslopeNeighborIteratorNoDivide(ManagedFlowDirRaster, int, int)
+        NeighborTuple next()
+
+    cdef cppclass UpslopeNeighborIteratorSkip:
+        ManagedFlowDirRaster raster
+        int col
+        int row
+        int n_dir
+        int flow_dir
+
+        UpslopeNeighborIteratorSkip()
+        UpslopeNeighborIteratorSkip(ManagedFlowDirRaster, int, int, int)
+        NeighborTuple next()
 
     bint is_close(double, double)
 

src/natcap/invest/ndr/ndr_core.pyx

@@ -13,7 +13,7 @@ from libcpp.stack cimport stack
 from libc.math cimport exp
 from ..managed_raster.managed_raster cimport ManagedRaster
 from ..managed_raster.managed_raster cimport ManagedFlowDirRaster
-from ..managed_raster.managed_raster cimport DownslopeNeighborIterator
+from ..managed_raster.managed_raster cimport DownslopeNeighborIteratorNoSkip
 from ..managed_raster.managed_raster cimport UpslopeNeighborIterator
 from ..managed_raster.managed_raster cimport NeighborTuple
 from ..managed_raster.managed_raster cimport is_close
@@ -185,15 +185,15 @@ def ndr_eff_calculation(
             else:
                 working_retention_eff = 0.0
 
-                dn_iterator = DownslopeNeighborIterator(
+                dn_iterator = DownslopeNeighborIteratorNoSkip(
                     mfd_flow_direction_raster, global_col, global_row)
                 has_outflow = False
-                neighbor = dn_iterator.next_no_skip()
+                neighbor = dn_iterator.next()
                 while neighbor.direction < 8:
                     has_outflow = True
                     if (neighbor.x < 0 or neighbor.x >= n_cols or
                         neighbor.y < 0 or neighbor.y >= n_rows):
-                        neighbor = dn_iterator.next_no_skip()
+                        neighbor = dn_iterator.next()
                         continue
                     if neighbor.direction % 2 == 1:
                         step_size = <float>(cell_size * 1.41421356237)
@@ -227,7 +227,7 @@ def ndr_eff_calculation(
 
                     working_retention_eff += (
                         intermediate_retention * neighbor.flow_proportion)
-                    neighbor = dn_iterator.next_no_skip()
+                    neighbor = dn_iterator.next()
 
                 if has_outflow:
                     working_retention_eff /= dn_iterator.flow_dir_sum

src/natcap/invest/sdr/sdr_core.pyx

@@ -13,6 +13,7 @@ from ..managed_raster.managed_raster cimport ManagedFlowDirRaster
 from ..managed_raster.managed_raster cimport is_close
 from ..managed_raster.managed_raster cimport DownslopeNeighborIterator
 from ..managed_raster.managed_raster cimport UpslopeNeighborIterator
+from ..managed_raster.managed_raster cimport UpslopeNeighborIteratorSkip
 
 
 cdef extern from "time.h" nogil:
@@ -118,6 +119,8 @@ def calculate_sediment_deposition(
     cdef time_t last_log_time = ctime(NULL)
     cdef float f_j_weighted_sum
     cdef NeighborTuple neighbor
+    cdef UpslopeNeighborIterator up_iterator
+    cdef UpslopeNeighborIteratorSkip up_iterator_skip
 
     for offset_dict in pygeoprocessing.iterblocks(
             (mfd_flow_direction_path, 1), offset_only=True, largest_block=0):
@@ -205,16 +208,16 @@ def calculate_sediment_deposition(
                         # completed
                         upslope_neighbors_processed = 1
                         # iterate over each neighbor-of-neighbor
-                        up_iterator = UpslopeNeighborIterator(
-                            mfd_flow_direction_raster, neighbor.x, neighbor.y)
-                        neighbor_of_neighbor = up_iterator.next_skip(neighbor.direction)
+                        up_iterator_skip = UpslopeNeighborIteratorSkip(
+                            mfd_flow_direction_raster, neighbor.x, neighbor.y, neighbor.direction)
+                        neighbor_of_neighbor = up_iterator_skip.next()
                         while neighbor_of_neighbor.direction < 8:
                             if is_close(sediment_deposition_raster.get(
                                 neighbor_of_neighbor.x, neighbor_of_neighbor.y
                             ), target_nodata):
                                 upslope_neighbors_processed = 0
                                 break
-                            neighbor_of_neighbor = up_iterator.next_skip(neighbor.direction)
+                            neighbor_of_neighbor = up_iterator_skip.next()
                         # if all upslope neighbors of neighbor j are
                         # processed, we can push j onto the stack.
                         if upslope_neighbors_processed:

src/natcap/invest/seasonal_water_yield/seasonal_water_yield_core.pyx

@@ -21,6 +21,7 @@ from ..managed_raster.managed_raster cimport ManagedRaster
 from ..managed_raster.managed_raster cimport ManagedFlowDirRaster
 from ..managed_raster.managed_raster cimport DownslopeNeighborIterator
 from ..managed_raster.managed_raster cimport UpslopeNeighborIterator
+from ..managed_raster.managed_raster cimport UpslopeNeighborIteratorNoDivide
 from ..managed_raster.managed_raster cimport NeighborTuple
 from ..managed_raster.managed_raster cimport is_close
 
@@ -217,8 +218,8 @@ cpdef calculate_local_recharge(
                     # initialize to 0 so we indicate we haven't tracked any
                     # mfd values yet
                     l_sum_avail_i = 0.0
-                    up_iterator = UpslopeNeighborIterator(flow_raster, xi, yi)
-                    neighbor = up_iterator.next_no_divide()
+                    up_iterator = UpslopeNeighborIteratorNoDivide(flow_raster, xi, yi)
+                    neighbor = up_iterator.next()
                     mfd_dir_sum = 0
                     while neighbor.direction < 8:
                         # pixel flows inward, check upslope
@@ -233,7 +234,7 @@ cpdef calculate_local_recharge(
                         l_sum_avail_i += (
                             l_sum_avail_j + l_avail_j) * neighbor.flow_proportion
                         mfd_dir_sum += <int>neighbor.flow_proportion
-                        neighbor = up_iterator.next_no_divide()
+                        neighbor = up_iterator.next()
                     # calculate l_sum_avail_i by summing all the valid
                     # directions then normalizing by the sum of the mfd
                     # direction weights (Equation 8)