diff --git a/pysheds/__init__.py b/pysheds/__init__.py
index 6a35e85..260c070 100644
--- a/pysheds/__init__.py
+++ b/pysheds/__init__.py
@@ -1 +1 @@
-__version__ = "0.3"
+__version__ = "0.3.1"
diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index e096b5e..b60c5dc 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -269,7 +269,7 @@ def _d8_catchment_recursion(ix, catch, fdir, offsets, r_dirmap):
 
 @njit(boolean[:,:](int64[:,:], UniTuple(int64, 2), UniTuple(int64, 8)),
       cache=True)
-def _d8_catchment_numba(fdir, pour_point, dirmap):
+def _d8_catchment_recur_numba(fdir, pour_point, dirmap):
     catch = np.zeros(fdir.shape, dtype=np.bool8)
     offset = fdir.shape[1]
     i, j = pour_point
@@ -282,6 +282,33 @@ def _d8_catchment_numba(fdir, pour_point, dirmap):
     _d8_catchment_recursion(ix, catch, fdir, offsets, r_dirmap)
     return catch
 
+@njit(boolean[:,:](int64[:,:], UniTuple(int64, 2), UniTuple(int64, 8)))
+def _d8_catchment_iter_numba(fdir, pour_point, dirmap):
+    catch = np.zeros(fdir.shape, dtype=np.bool8)
+    offset = fdir.shape[1]
+    i, j = pour_point
+    ix = (i * offset) + j
+    offsets = np.array([-offset, 1 - offset, 1, 1 + offset,
+                        offset, - 1 + offset, - 1, - 1 - offset])
+    r_dirmap = np.array([dirmap[4], dirmap[5], dirmap[6],
+                         dirmap[7], dirmap[0], dirmap[1],
+                         dirmap[2], dirmap[3]])
+    queue = [ix]
+    while queue:
+        parent = queue.pop()
+        catch.flat[parent] = True
+        neighbors = offsets + parent
+        for k in range(8):
+            neighbor = neighbors[k]
+            visited = catch.flat[neighbor]
+            if visited:
+                continue
+            else:
+                points_to = (fdir.flat[neighbor] == r_dirmap[k])
+                if points_to:
+                    queue.append(neighbor)
+    return catch
+
 @njit(void(int64, boolean[:,:], int64[:,:], int64[:,:], int64[:], int64[:]),
       cache=True)
 def _dinf_catchment_recursion(ix, catch, fdir_0, fdir_1, offsets, r_dirmap):
@@ -299,7 +326,7 @@ def _dinf_catchment_recursion(ix, catch, fdir_0, fdir_1, offsets, r_dirmap):
 
 @njit(boolean[:,:](int64[:,:], int64[:,:], UniTuple(int64, 2), UniTuple(int64, 8)),
       cache=True)
-def _dinf_catchment_numba(fdir_0, fdir_1, pour_point, dirmap):
+def _dinf_catchment_recur_numba(fdir_0, fdir_1, pour_point, dirmap):
     catch = np.zeros(fdir_0.shape, dtype=np.bool8)
     dirmap = np.array(dirmap)
     offset = fdir_0.shape[1]
@@ -314,6 +341,38 @@ def _dinf_catchment_numba(fdir_0, fdir_1, pour_point, dirmap):
     _dinf_catchment_recursion(ix, catch, fdir_0, fdir_1, offsets, r_dirmap)
     return catch
 
+@njit(boolean[:,:](int64[:,:], int64[:,:], UniTuple(int64, 2), UniTuple(int64, 8)),
+      cache=True)
+def _dinf_catchment_iter_numba(fdir_0, fdir_1, pour_point, dirmap):
+    catch = np.zeros(fdir_0.shape, dtype=np.bool8)
+    dirmap = np.array(dirmap)
+    offset = fdir_0.shape[1]
+    i, j = pour_point
+    ix = (i * offset) + j
+    offsets = np.array([-offset, 1 - offset, 1,
+                        1 + offset, offset, - 1 + offset,
+                        - 1, - 1 - offset])
+    r_dirmap = np.array([dirmap[4], dirmap[5], dirmap[6],
+                         dirmap[7], dirmap[0], dirmap[1],
+                         dirmap[2], dirmap[3]])
+    queue = [ix]
+    while queue:
+        parent = queue.pop()
+        catch.flat[parent] = True
+        neighbors = offsets + parent
+        for k in range(8):
+            neighbor = neighbors[k]
+            visited = catch.flat[neighbor]
+            if visited:
+                continue
+            else:
+                points_to_0 = (fdir_0.flat[neighbor] == r_dirmap[k])
+                points_to_1 = (fdir_1.flat[neighbor] == r_dirmap[k])
+                points_to = points_to_0 or points_to_1
+                if points_to:
+                    queue.append(neighbor)
+    return catch
+
 # Functions for 'accumulation'
 
 @njit(void(int64, int64, float64[:,:], int64[:,:], uint8[:]),
@@ -328,7 +387,7 @@ def _d8_accumulation_recursion(startnode, endnode, acc, fdir, indegree):
 
 @njit(float64[:,:](float64[:,:], int64[:,:], uint8[:], int64[:]),
       cache=True)
-def _d8_accumulation_numba(acc, fdir, indegree, startnodes):
+def _d8_accumulation_recur_numba(acc, fdir, indegree, startnodes):
     n = startnodes.size
     for k in range(n):
         startnode = startnodes[k]
@@ -336,6 +395,20 @@ def _d8_accumulation_numba(acc, fdir, indegree, startnodes):
         _d8_accumulation_recursion(startnode, endnode, acc, fdir, indegree)
     return acc
 
+@njit(float64[:,:](float64[:,:], int64[:,:], uint8[:], int64[:]),
+      cache=True)
+def _d8_accumulation_iter_numba(acc, fdir, indegree, startnodes):
+    n = startnodes.size
+    for k in range(n):
+        startnode = startnodes[k]
+        endnode = fdir.flat[startnode]
+        while(indegree[startnode] == 0):
+            acc.flat[endnode] += acc.flat[startnode]
+            indegree[endnode] -= 1
+            startnode = endnode
+            endnode = fdir.flat[startnode]
+    return acc
+
 @njit(void(int64, int64, float64[:,:], int64[:,:], uint8[:], float64[:,:]),
       cache=True)
 def _d8_accumulation_eff_recursion(startnode, endnode, acc, fdir, indegree, eff):
@@ -348,7 +421,7 @@ def _d8_accumulation_eff_recursion(startnode, endnode, acc, fdir, indegree, eff)
 
 @njit(float64[:,:](float64[:,:], int64[:,:], uint8[:], int64[:], float64[:,:]),
       cache=True)
-def _d8_accumulation_eff_numba(acc, fdir, indegree, startnodes, eff):
+def _d8_accumulation_eff_recur_numba(acc, fdir, indegree, startnodes, eff):
     n = startnodes.size
     for k in range(n):
         startnode = startnodes[k]
@@ -356,6 +429,20 @@ def _d8_accumulation_eff_numba(acc, fdir, indegree, startnodes, eff):
         _d8_accumulation_eff_recursion(startnode, endnode, acc, fdir, indegree, eff)
     return acc
 
+@njit(float64[:,:](float64[:,:], int64[:,:], uint8[:], int64[:], float64[:,:]),
+      cache=True)
+def _d8_accumulation_eff_iter_numba(acc, fdir, indegree, startnodes, eff):
+    n = startnodes.size
+    for k in range(n):
+        startnode = startnodes[k]
+        endnode = fdir.flat[startnode]
+        while(indegree[startnode] == 0):
+            acc.flat[endnode] += (acc.flat[startnode] * eff.flat[startnode])
+            indegree[endnode] -= 1
+            startnode = endnode
+            endnode = fdir.flat[startnode]
+    return acc
+
 @njit(void(int64, int64, float64[:,:], int64[:,:], int64[:,:], uint8[:], float64,
            boolean[:,:], float64[:,:], float64[:,:]),
       cache=True)
@@ -378,8 +465,8 @@ def _dinf_accumulation_recursion(startnode, endnode, acc, fdir_0, fdir_1,
 @njit(float64[:,:](float64[:,:], int64[:,:], int64[:,:], uint8[:], int64[:],
                    float64[:,:], float64[:,:]),
       cache=True)
-def _dinf_accumulation_numba(acc, fdir_0, fdir_1, indegree, startnodes,
-                             props_0, props_1):
+def _dinf_accumulation_recur_numba(acc, fdir_0, fdir_1, indegree, startnodes,
+                                   props_0, props_1):
     n = startnodes.size
     visited = np.zeros(acc.shape, dtype=np.bool8)
     for k in range(n):
@@ -396,6 +483,35 @@ def _dinf_accumulation_numba(acc, fdir_0, fdir_1, indegree, startnodes,
         visited.flat[startnode] = True
     return acc
 
+@njit(float64[:,:](float64[:,:], int64[:,:], int64[:,:], uint8[:], int64[:],
+                   float64[:,:], float64[:,:]),
+      cache=True)
+def _dinf_accumulation_iter_numba(acc, fdir_0, fdir_1, indegree, startnodes,
+                                  props_0, props_1):
+    n = startnodes.size
+    queue = [0]
+    _ = queue.pop()
+    for k in range(n):
+        startnode = startnodes.flat[k]
+        queue.append(startnode)
+        while queue:
+            startnode = queue.pop()
+            endnode_0 = fdir_0.flat[startnode]
+            endnode_1 = fdir_1.flat[startnode]
+            prop_0 = props_0.flat[startnode]
+            prop_1 = props_1.flat[startnode]
+            acc.flat[endnode_0] += (prop_0 * acc.flat[startnode])
+            acc.flat[endnode_1] += (prop_1 * acc.flat[startnode])
+            indegree.flat[endnode_0] -= 1
+            indegree.flat[endnode_1] -= 1
+            if (indegree.flat[endnode_0] == 0):
+                queue.append(endnode_0)
+            if (indegree.flat[endnode_1] == 0):
+                # Account for cases where both fdirs point in same direction
+                if (endnode_0 != endnode_1):
+                    queue.append(endnode_1)
+    return acc
+
 @njit(void(int64, int64, float64[:,:], int64[:,:], int64[:,:], uint8[:], float64,
            boolean[:,:], float64[:,:], float64[:,:], float64[:,:]),
       cache=True)
@@ -436,6 +552,36 @@ def _dinf_accumulation_eff_numba(acc, fdir_0, fdir_1, indegree, startnodes,
         visited.flat[startnode] = True
     return acc
 
+@njit(float64[:,:](float64[:,:], int64[:,:], int64[:,:], uint8[:], int64[:],
+                   float64[:,:], float64[:,:], float64[:,:]),
+      cache=True)
+def _dinf_accumulation_eff_iter_numba(acc, fdir_0, fdir_1, indegree, startnodes,
+                                      props_0, props_1, eff):
+    n = startnodes.size
+    queue = [0]
+    _ = queue.pop()
+    for k in range(n):
+        startnode = startnodes.flat[k]
+        queue.append(startnode)
+        while queue:
+            startnode = queue.pop()
+            endnode_0 = fdir_0.flat[startnode]
+            endnode_1 = fdir_1.flat[startnode]
+            prop_0 = props_0.flat[startnode]
+            prop_1 = props_1.flat[startnode]
+            transfer = acc.flat[startnode] * eff.flat[startnode]
+            acc.flat[endnode_0] += (prop_0 * transfer)
+            acc.flat[endnode_1] += (prop_1 * transfer)
+            indegree.flat[endnode_0] -= 1
+            indegree.flat[endnode_1] -= 1
+            if (indegree.flat[endnode_0] == 0):
+                queue.append(endnode_0)
+            if (indegree.flat[endnode_1] == 0):
+                # Account for cases where both fdirs point in same direction
+                if (endnode_0 != endnode_1):
+                    queue.append(endnode_1)
+    return acc
+
 # Functions for 'flow_distance'
 
 @njit(void(int64, int64[:,:], boolean[:,:], float64[:,:], float64[:,:],
@@ -458,7 +604,7 @@ def _d8_flow_distance_recursion(ix, fdir, visits, dist, weights, r_dirmap,
 
 @njit(float64[:,:](int64[:,:], float64[:,:], UniTuple(int64, 2), UniTuple(int64, 8)),
       cache=True)
-def _d8_flow_distance_numba(fdir, weights, pour_point, dirmap):
+def _d8_flow_distance_recur_numba(fdir, weights, pour_point, dirmap):
     visits = np.zeros(fdir.shape, dtype=np.bool8)
     dist = np.full(fdir.shape, np.inf, dtype=np.float64)
     r_dirmap = np.array([dirmap[4], dirmap[5], dirmap[6],
@@ -474,6 +620,38 @@ def _d8_flow_distance_numba(fdir, weights, pour_point, dirmap):
                                 r_dirmap, 0., offsets)
     return dist
 
+@njit(float64[:,:](int64[:,:], float64[:,:], UniTuple(int64, 2), UniTuple(int64, 8)),
+      cache=True)
+def _d8_flow_distance_iter_numba(fdir, weights, pour_point, dirmap):
+    visits = np.zeros(fdir.shape, dtype=np.bool8)
+    dist = np.full(fdir.shape, np.inf, dtype=np.float64)
+    r_dirmap = np.array([dirmap[4], dirmap[5], dirmap[6],
+                         dirmap[7], dirmap[0], dirmap[1],
+                         dirmap[2], dirmap[3]])
+    m, n = fdir.shape
+    offsets = np.array([-n, 1 - n, 1,
+                        1 + n, n, - 1 + n,
+                        - 1, - 1 - n])
+    i, j = pour_point
+    ix = (i * n) + j
+    dist.flat[ix] = 0.
+    queue = [ix]
+    while queue:
+        parent = queue.pop()
+        visits.flat[parent] = True
+        neighbors = offsets + parent
+        for k in range(8):
+            neighbor = neighbors[k]
+            visited = visits.flat[neighbor]
+            if visited:
+                continue
+            else:
+                points_to = (fdir.flat[neighbor] == r_dirmap[k])
+                if points_to:
+                    dist.flat[neighbor] = dist.flat[parent] + weights.flat[neighbor]
+                    queue.append(neighbor)
+    return dist
+
 @njit(void(int64, int64[:,:], int64[:,:], boolean[:,:], float64[:,:],
            float64[:,:], float64[:,:], int64[:], float64, int64[:]),
       cache=True)
@@ -501,7 +679,7 @@ def _dinf_flow_distance_recursion(ix, fdir_0, fdir_1, visits, dist,
 @njit(float64[:,:](int64[:,:], int64[:,:], float64[:,:], float64[:,:],
                    UniTuple(int64, 2), UniTuple(int64, 8)),
       cache=True)
-def _dinf_flow_distance_numba(fdir_0, fdir_1, weights_0, weights_1,
+def _dinf_flow_distance_recur_numba(fdir_0, fdir_1, weights_0, weights_1,
                               pour_point, dirmap):
     visits = np.zeros(fdir_0.shape, dtype=np.bool8)
     dist = np.full(fdir_0.shape, np.inf, dtype=np.float64)
@@ -518,6 +696,44 @@ def _dinf_flow_distance_numba(fdir_0, fdir_1, weights_0, weights_1,
                                   weights_0, weights_1, r_dirmap, 0., offsets)
     return dist
 
+@njit(float64[:,:](int64[:,:], int64[:,:], float64[:,:], float64[:,:],
+                   UniTuple(int64, 2), UniTuple(int64, 8)),
+      cache=True)
+def _dinf_flow_distance_iter_numba(fdir_0, fdir_1, weights_0, weights_1,
+                                   pour_point, dirmap):
+    dist = np.full(fdir_0.shape, np.inf, dtype=np.float64)
+    r_dirmap = np.array([dirmap[4], dirmap[5], dirmap[6],
+                         dirmap[7], dirmap[0], dirmap[1],
+                         dirmap[2], dirmap[3]])
+    m, n = fdir_0.shape
+    offsets = np.array([-n, 1 - n, 1,
+                        1 + n, n, - 1 + n,
+                        - 1, - 1 - n])
+    i, j = pour_point
+    ix = (i * n) + j
+    dist.flat[ix] = 0.
+    queue = [ix]
+    while queue:
+        parent = queue.pop()
+        parent_dist = dist.flat[parent]
+        neighbors = offsets + parent
+        for k in range(8):
+            neighbor = neighbors[k]
+            current_neighbor_dist = dist.flat[neighbor]
+            points_to_0 = (fdir_0.flat[neighbor] == r_dirmap[k])
+            points_to_1 = (fdir_1.flat[neighbor] == r_dirmap[k])
+            if points_to_0:
+                neighbor_dist_0 = parent_dist + weights_0.flat[neighbor]
+                if (neighbor_dist_0 < current_neighbor_dist):
+                    dist.flat[neighbor] = neighbor_dist_0
+                    queue.append(neighbor)
+            elif points_to_1:
+                neighbor_dist_1 = parent_dist + weights_1.flat[neighbor]
+                if (neighbor_dist_1 < current_neighbor_dist):
+                    dist.flat[neighbor] = neighbor_dist_1
+                    queue.append(neighbor)
+    return dist
+
 @njit(void(int64, int64, int64[:,:], int64[:,:], float64[:,:],
            int64[:,:], uint8[:], float64[:,:]),
       cache=True)
@@ -536,8 +752,8 @@ def _d8_reverse_distance_recursion(startnode, endnode, min_order, max_order,
 @njit(float64[:,:](int64[:,:], int64[:,:], float64[:,:], int64[:,:],
                    uint8[:], int64[:], float64[:,:]),
       cache=True)
-def _d8_reverse_distance_numba(min_order, max_order, rdist, fdir,
-                               indegree, startnodes, weights):
+def _d8_reverse_distance_recur_numba(min_order, max_order, rdist, fdir,
+                                     indegree, startnodes, weights):
     n = startnodes.size
     for k in range(n):
         startnode = startnodes.flat[k]
@@ -546,6 +762,24 @@ def _d8_reverse_distance_numba(min_order, max_order, rdist, fdir,
                                        rdist, fdir, indegree, weights)
     return rdist
 
+@njit(float64[:,:](int64[:,:], int64[:,:], float64[:,:], int64[:,:],
+                   uint8[:], int64[:], float64[:,:]),
+      cache=True)
+def _d8_reverse_distance_iter_numba(min_order, max_order, rdist, fdir,
+                                    indegree, startnodes, weights):
+    n = startnodes.size
+    for k in range(n):
+        startnode = startnodes.flat[k]
+        endnode = fdir.flat[startnode]
+        while(indegree.flat[startnode] == 0):
+            min_order.flat[endnode] = min(min_order.flat[endnode], rdist.flat[startnode])
+            max_order.flat[endnode] = max(max_order.flat[endnode], rdist.flat[startnode])
+            indegree.flat[endnode] -= 1
+            rdist.flat[endnode] = max_order.flat[endnode] + weights.flat[endnode]
+            startnode = endnode
+            endnode = fdir.flat[startnode]
+    return rdist
+
 # Functions for 'resolve_flats'
 
 @njit(UniTuple(boolean[:,:], 3)(float64[:,:], int64[:]),
@@ -766,10 +1000,9 @@ def _d8_hand_recursion(child, parent, hand, offsets, r_dirmap, fdir):
             hand.flat[neighbor] = parent
             _d8_hand_recursion(neighbor, parent, hand, offsets, r_dirmap, fdir)
 
-@njit(int64[:,:](int64[:], int64[:,:], UniTuple(int64, 8)),
+@njit(int64[:,:](int64[:,:], boolean[:,:], UniTuple(int64, 8)),
       cache=True)
-def _d8_hand_recursive_numba(parents, fdir, dirmap):
-    n = parents.size
+def _d8_hand_recur_numba(fdir, mask, dirmap):
     offset = fdir.shape[1]
     offsets = np.array([-offset, 1 - offset, 1,
                         1 + offset, offset, - 1 + offset,
@@ -778,6 +1011,8 @@ def _d8_hand_recursive_numba(parents, fdir, dirmap):
                          dirmap[7], dirmap[0], dirmap[1],
                          dirmap[2], dirmap[3]])
     hand = -np.ones(fdir.shape, dtype=np.int64)
+    parents = np.flatnonzero(mask)
+    n = parents.size
     for i in range(n):
         parent = parents[i]
         hand.flat[parent] = parent
@@ -835,10 +1070,9 @@ def _dinf_hand_recursion(child, parent, hand, offsets, r_dirmap, fdir_0, fdir_1)
             hand.flat[neighbor] = parent
             _dinf_hand_recursion(neighbor, parent, hand, offsets, r_dirmap, fdir_0, fdir_1)
 
-@njit(int64[:,:](int64[:], int64[:,:], int64[:,:], UniTuple(int64, 8)),
+@njit(int64[:,:](int64[:,:], int64[:,:], boolean[:,:], UniTuple(int64, 8)),
       cache=True)
-def _dinf_hand_recursive_numba(parents, fdir_0, fdir_1, dirmap):
-    n = parents.size
+def _dinf_hand_recur_numba(fdir_0, fdir_1, mask, dirmap):
     offset = fdir_0.shape[1]
     offsets = np.array([-offset, 1 - offset, 1,
                         1 + offset, offset, - 1 + offset,
@@ -847,6 +1081,8 @@ def _dinf_hand_recursive_numba(parents, fdir_0, fdir_1, dirmap):
                          dirmap[7], dirmap[0], dirmap[1],
                          dirmap[2], dirmap[3]])
     hand = -np.ones(fdir_0.shape, dtype=np.int64)
+    parents = np.flatnonzero(mask)
+    n = parents.size
     for i in range(n):
         parent = parents[i]
         hand.flat[parent] = parent
@@ -890,8 +1126,8 @@ def _d8_streamorder_recursion(startnode, endnode, min_order, max_order,
 
 @njit(int64[:,:](int64[:,:], int64[:,:], int64[:,:], int64[:,:], uint8[:], uint8[:], int64[:]),
       cache=True)
-def _d8_streamorder_numba(min_order, max_order, order, fdir,
-                          indegree, orig_indegree, startnodes):
+def _d8_streamorder_recur_numba(min_order, max_order, order, fdir,
+                                indegree, orig_indegree, startnodes):
     n = startnodes.size
     for k in range(n):
         startnode = startnodes.flat[k]
@@ -900,6 +1136,27 @@ def _d8_streamorder_numba(min_order, max_order, order, fdir,
                                  fdir, indegree, orig_indegree)
     return order
 
+@njit(int64[:,:](int64[:,:], int64[:,:], int64[:,:], int64[:,:], uint8[:], uint8[:], int64[:]),
+      cache=True)
+def _d8_streamorder_iter_numba(min_order, max_order, order, fdir,
+                               indegree, orig_indegree, startnodes):
+    n = startnodes.size
+    for k in range(n):
+        startnode = startnodes.flat[k]
+        endnode = fdir.flat[startnode]
+        while (indegree.flat[startnode] == 0):
+            min_order.flat[endnode] = min(min_order.flat[endnode], order.flat[startnode])
+            max_order.flat[endnode] = max(max_order.flat[endnode], order.flat[startnode])
+            indegree.flat[endnode] -= 1
+            if ((min_order.flat[endnode] == max_order.flat[endnode]) and
+                (orig_indegree.flat[endnode] > 1)):
+                order.flat[endnode] = max_order.flat[endnode] + 1
+            else:
+                order.flat[endnode] = max_order.flat[endnode]
+            startnode = endnode
+            endnode = fdir.flat[startnode]
+    return order
+
 @njit(void(int64, int64, int64[:,:], uint8[:], uint8[:], List(List(int64)), List(int64)),
       cache=True)
 def _d8_stream_network_recursion(startnode, endnode, fdir, indegree,
@@ -918,7 +1175,7 @@ def _d8_stream_network_recursion(startnode, endnode, fdir, indegree,
 
 @njit(List(List(int64))(int64[:,:], uint8[:], uint8[:], int64[:]),
       cache=True)
-def _d8_stream_network_numba(fdir, indegree, orig_indegree, startnodes):
+def _d8_stream_network_recur_numba(fdir, indegree, orig_indegree, startnodes):
     n = startnodes.size
     profiles = [[0]]
     _ = profiles.pop()
@@ -930,6 +1187,26 @@ def _d8_stream_network_numba(fdir, indegree, orig_indegree, startnodes):
                                      orig_indegree, profiles, profile)
     return profiles
 
+@njit(List(List(int64))(int64[:,:], uint8[:], uint8[:], int64[:]),
+      cache=True)
+def _d8_stream_network_iter_numba(fdir, indegree, orig_indegree, startnodes):
+    n = startnodes.size
+    profiles = [[0]]
+    _ = profiles.pop()
+    for k in range(n):
+        startnode = startnodes.flat[k]
+        endnode = fdir.flat[startnode]
+        profile = [startnode]
+        while (indegree.flat[startnode] == 0):
+            profile.append(endnode)
+            indegree.flat[endnode] -= 1
+            if (orig_indegree[endnode] > 1):
+                profiles.append(profile)
+                profile = [endnode]
+            startnode = endnode
+            endnode = fdir.flat[startnode]
+    return profiles
+
 @njit(parallel=True)
 def _d8_cell_dh_numba(startnodes, endnodes, dem):
     n = startnodes.size
diff --git a/pysheds/io.py b/pysheds/io.py
index 05afe34..8d477ea 100644
--- a/pysheds/io.py
+++ b/pysheds/io.py
@@ -150,8 +150,9 @@ def read_raster(data, band=1, window=None, window_crs=None, mask_geometry=False,
 
 def to_ascii(data, file_name, target_view=None, delimiter=' ', fmt=None,
              interpolation='nearest', apply_input_mask=False,
-             apply_output_mask=True, affine=None, shape=None, crs=None,
-             mask=None, nodata=None, dtype=None, **kwargs):
+             apply_output_mask=True, inherit_nodata=True, affine=None,
+             shape=None, crs=None, mask=None, nodata=None, dtype=None,
+             **kwargs):
     """
     Writes a Raster object to a formatted ascii text file.
 
@@ -174,6 +175,9 @@ def to_ascii(data, file_name, target_view=None, delimiter=' ', fmt=None,
                         If True, mask the input Raster according to data.mask.
     apply_output_mask : bool
                         If True, mask the output Raster according to target_view.mask.
+    inherit_nodata : bool
+                     If True, output ascii inherits `nodata` value from `data`.
+                     If False, output ascii uses `nodata` value from `target_view`.
     affine : affine.Affine
                 Affine transformation matrix (overrides target_view.affine)
     shape : tuple of ints (length 2)
@@ -193,9 +197,9 @@ def to_ascii(data, file_name, target_view=None, delimiter=' ', fmt=None,
         target_view = data.viewfinder
     data = View.view(data, target_view, interpolation=interpolation,
                      apply_input_mask=apply_input_mask,
-                     apply_output_mask=apply_output_mask, affine=affine,
-                     shape=shape, crs=crs, mask=mask, nodata=nodata,
-                     dtype=dtype)
+                     apply_output_mask=apply_output_mask,
+                     inherit_nodata=inherit_nodata, affine=affine, shape=shape,
+                     crs=crs, mask=mask, nodata=nodata, dtype=dtype)
     try:
         assert (abs(data.affine.a) == abs(data.affine.e))
     except:
@@ -225,8 +229,9 @@ def to_ascii(data, file_name, target_view=None, delimiter=' ', fmt=None,
 
 def to_raster(data, file_name, target_view=None, profile=None, blockxsize=256,
               blockysize=256, interpolation='nearest', apply_input_mask=False,
-              apply_output_mask=True, affine=None, shape=None, crs=None,
-              mask=None, nodata=None, dtype=None, **kwargs):
+              apply_output_mask=True, inherit_nodata=True, affine=None,
+              shape=None, crs=None, mask=None, nodata=None, dtype=None,
+              **kwargs):
     """
     Writes gridded data to a raster.
 
@@ -251,6 +256,9 @@ def to_raster(data, file_name, target_view=None, profile=None, blockxsize=256,
                         If True, mask the input Raster according to data.mask.
     apply_output_mask : bool
                         If True, mask the output Raster according to target_view.mask.
+    inherit_nodata : bool
+                     If True, output Raster inherits `nodata` value from `data`.
+                     If False, output Raster uses `nodata` value from `target_view`.
     affine : affine.Affine
                 Affine transformation matrix (overrides target_view.affine)
     shape : tuple of ints (length 2)
@@ -268,9 +276,9 @@ def to_raster(data, file_name, target_view=None, profile=None, blockxsize=256,
         target_view = data.viewfinder
     data = View.view(data, target_view, interpolation=interpolation,
                      apply_input_mask=apply_input_mask,
-                     apply_output_mask=apply_output_mask, affine=affine,
-                     shape=shape, crs=crs, mask=mask, nodata=nodata,
-                     dtype=dtype)
+                     apply_output_mask=apply_output_mask,
+                     inherit_nodata=inherit_nodata, affine=affine, shape=shape,
+                     crs=crs, mask=mask, nodata=nodata, dtype=dtype)
     height, width = data.shape
     default_blockx = width
     default_profile = {
diff --git a/pysheds/sgrid.py b/pysheds/sgrid.py
index 9433439..d2c2861 100644
--- a/pysheds/sgrid.py
+++ b/pysheds/sgrid.py
@@ -262,10 +262,11 @@ def read_raster(self, data, band=1, window=None, window_crs=None,
                                       window_crs=window_crs, metadata=metadata,
                                       mask_geometry=mask_geometry, **kwargs)
 
-    def to_ascii(self, data, file_name, target_view=None, delimiter=' ', fmt=None,
-                interpolation='nearest', apply_input_mask=False,
-                apply_output_mask=True, affine=None, shape=None, crs=None,
-                mask=None, nodata=None, dtype=None, **kwargs):
+    def to_ascii(self, data, file_name, target_view=None, delimiter=' ',
+                 fmt=None, interpolation='nearest', apply_input_mask=False,
+                 apply_output_mask=True, inherit_nodata=True, affine=None,
+                 shape=None, crs=None, mask=None, nodata=None, dtype=None,
+                 **kwargs):
         """
         Writes a Raster object to a formatted ascii text file.
 
@@ -288,6 +289,9 @@ def to_ascii(self, data, file_name, target_view=None, delimiter=' ', fmt=None,
                             If True, mask the input Raster according to self.mask.
         apply_output_mask : bool
                             If True, mask the output Raster according to target_view.mask.
+        inherit_nodata : bool
+                         If True, output ascii inherits `nodata` value from `data`.
+                         If False, output ascii uses `nodata` value from grid's viewfinder.
         affine : affine.Affine
                     Affine transformation matrix (overrides target_view.affine)
         shape : tuple of ints (length 2)
@@ -309,15 +313,16 @@ def to_ascii(self, data, file_name, target_view=None, delimiter=' ', fmt=None,
                                    delimiter=delimiter, fmt=fmt, interpolation=interpolation,
                                    apply_input_mask=apply_input_mask,
                                    apply_output_mask=apply_output_mask,
+                                   inherit_nodata=inherit_nodata,
                                    affine=affine, shape=shape, crs=crs,
                                    mask=mask, nodata=nodata,
                                    dtype=dtype, **kwargs)
 
-    def to_raster(self, data, file_name, target_view=None, profile=None, view=True,
-                blockxsize=256, blockysize=256, interpolation='nearest',
-                apply_input_mask=False, apply_output_mask=True, affine=None,
-                shape=None, crs=None, mask=None, nodata=None, dtype=None,
-                **kwargs):
+    def to_raster(self, data, file_name, target_view=None, profile=None,
+                  blockxsize=256, blockysize=256, interpolation='nearest',
+                  apply_input_mask=False, apply_output_mask=True,
+                  inherit_nodata=True, affine=None, shape=None, crs=None,
+                  mask=None, nodata=None, dtype=None, **kwargs):
         """
         Writes gridded data to a raster.
 
@@ -342,6 +347,9 @@ def to_raster(self, data, file_name, target_view=None, profile=None, view=True,
                             If True, mask the input Raster according to self.mask.
         apply_output_mask : bool
                             If True, mask the output Raster according to target_view.mask.
+        inherit_nodata : bool
+                         If True, output Raster inherits `nodata` value from `data`.
+                         If False, output Raster uses `nodata` value from grid's viewfinder.
         affine : affine.Affine
                     Affine transformation matrix (overrides target_view.affine)
         shape : tuple of ints (length 2)
@@ -358,12 +366,12 @@ def to_raster(self, data, file_name, target_view=None, profile=None, view=True,
         if target_view is None:
             target_view = self.viewfinder
         return pysheds.io.to_raster(data, file_name, target_view=target_view,
-                                    profile=profile, view=view,
-                                    blockxsize=blockxsize,
+                                    profile=profile, blockxsize=blockxsize,
                                     blockysize=blockysize,
                                     interpolation=interpolation,
                                     apply_input_mask=apply_input_mask,
                                     apply_output_mask=apply_output_mask,
+                                    inherit_nodata=inherit_nodata,
                                     affine=affine, shape=shape, crs=crs,
                                     mask=mask, nodata=nodata, dtype=dtype,
                                     **kwargs)
@@ -420,7 +428,7 @@ def from_raster(cls, data, **kwargs):
             raise TypeError('`data` must be a Raster or str.')
 
     def view(self, data, data_view=None, target_view=None, interpolation='nearest',
-             apply_input_mask=False, apply_output_mask=True,
+             apply_input_mask=False, apply_output_mask=True, inherit_nodata=True,
              affine=None, shape=None, crs=None, mask=None, nodata=None,
              dtype=None, inherit_metadata=True, new_metadata={}, **kwargs):
         """
@@ -449,6 +457,9 @@ def view(self, data, data_view=None, target_view=None, interpolation='nearest',
                            If True, mask the input Raster according to data.mask.
         apply_output_mask : bool
                            If True, mask the output Raster according to grid.mask.
+        inherit_nodata : bool
+                         If True, output Raster inherits `nodata` value from `data` or `data_view`.
+                         If False, output Raster uses `nodata` value from grid's viewfinder.
         affine : affine.Affine
                  Affine transformation matrix (overrides target_view.affine)
         shape : tuple of ints (length 2)
@@ -462,7 +473,7 @@ def view(self, data, data_view=None, target_view=None, interpolation='nearest',
         dtype : numpy datatype
                 Desired datatype of the output array.
         inherit_metadata : bool
-                           If True, output Raster inherits metadata from input data.
+                           If True, output Raster inherits metadata from input Raster.
         new_metadata : dict
                        Optional metadata to add to output Raster.
 
@@ -490,6 +501,7 @@ def view(self, data, data_view=None, target_view=None, interpolation='nearest',
                         interpolation=interpolation,
                         apply_input_mask=apply_input_mask,
                         apply_output_mask=apply_output_mask,
+                        inherit_nodata=inherit_nodata,
                         affine=affine, shape=shape,
                         crs=crs, mask=mask, nodata=nodata,
                         dtype=dtype,
@@ -629,7 +641,8 @@ def _dinf_flowdir(self, dem, nodata_cells, nodata_out=np.nan, flats=-1, pits=-2,
                                     metadata=dem.metadata, nodata=nodata_out)
 
     def catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                  nodata_out=False, xytype='coordinate', routing='d8', snap='corner', **kwargs):
+                  nodata_out=False, xytype='coordinate', routing='d8', snap='corner',
+                  algorithm='iterative', **kwargs):
         """
         Delineates a watershed from a given pour point (x, y).
 
@@ -650,7 +663,7 @@ def catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4, 8, 16
                  [N, NE, E, SE, S, SW, W, NW]
         nodata_out : int or float
                      Value to indicate `no data` in output array.
-        xytype : 'coordinate' or 'index'
+        xytype : str
                  How to interpret parameters 'x' and 'y'.
                      'coordinate' : x and y represent geographic coordinates
                                     (will be passed to self.nearest_cell).
@@ -664,6 +677,10 @@ def catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4, 8, 16
                Function to use for self.nearest_cell:
                'corner' : numpy.around()
                'center' : numpy.floor()
+        algorithm : str
+                    Algorithm type to use:
+                    'iterative' : Use an iterative algorithm (recommended).
+                    'recursive' : Use a recursive algorithm.
 
         Additional keyword arguments (**kwargs) are passed to self.view.
 
@@ -693,14 +710,17 @@ def catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4, 8, 16
                                 .format(x, y, fdir.shape))
         if routing.lower() == 'd8':
             catch = self._d8_catchment(x, y, fdir=fdir, pour_value=pour_value, dirmap=dirmap,
-                                       nodata_out=nodata_out, xytype=xytype, snap=snap)
+                                       nodata_out=nodata_out, xytype=xytype, snap=snap,
+                                       algorithm=algorithm)
         elif routing.lower() == 'dinf':
             catch = self._dinf_catchment(x, y, fdir=fdir, pour_value=pour_value, dirmap=dirmap,
-                                         nodata_out=nodata_out, xytype=xytype, snap=snap)
+                                         nodata_out=nodata_out, xytype=xytype, snap=snap,
+                                         algorithm=algorithm)
         return catch
 
     def _d8_catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                      nodata_out=False, xytype='coordinate', snap='corner'):
+                      nodata_out=False, xytype='coordinate', snap='corner',
+                      algorithm='iterative'):
         # Pad the rim
         left, right, top, bottom = self._pop_rim(fdir, nodata=0)
         # If xytype is 'coordinate', delineate catchment based on cell nearest
@@ -708,7 +728,12 @@ def _d8_catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4, 8
         if xytype in {'label', 'coordinate'}:
             x, y = self.nearest_cell(x, y, fdir.affine, snap)
         # Delineate the catchment
-        catch = _self._d8_catchment_numba(fdir, (y, x), dirmap)
+        if algorithm.lower() == 'iterative':
+            catch = _self._d8_catchment_iter_numba(fdir, (y, x), dirmap)
+        elif algorithm.lower() == 'recursive':
+            catch = _self._d8_catchment_recur_numba(fdir, (y, x), dirmap)
+        else:
+            raise ValueError('Algorithm must be `iterative` or `recursive`.')
         if pour_value is not None:
             catch[y, x] = pour_value
         catch = self._output_handler(data=catch, viewfinder=fdir.viewfinder,
@@ -716,7 +741,8 @@ def _d8_catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4, 8
         return catch
 
     def _dinf_catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                        nodata_out=False, xytype='coordinate', snap='corner'):
+                        nodata_out=False, xytype='coordinate', snap='corner',
+                        algorithm='iterative'):
         # Find nodata cells
         nodata_cells = self._get_nodata_cells(fdir)
         # Split dinf flowdir
@@ -728,7 +754,12 @@ def _dinf_catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4,
         if xytype in {'label', 'coordinate'}:
             x, y = self.nearest_cell(x, y, fdir.affine, snap)
         # Delineate the catchment
-        catch = _self._dinf_catchment_numba(fdir_0, fdir_1, (y, x), dirmap)
+        if algorithm.lower() == 'iterative':
+            catch = _self._dinf_catchment_iter_numba(fdir_0, fdir_1, (y, x), dirmap)
+        elif algorithm.lower() == 'recursive':
+            catch = _self._dinf_catchment_recur_numba(fdir_0, fdir_1, (y, x), dirmap)
+        else:
+            raise ValueError('Algorithm must be `iterative` or `recursive`.')
         # if pour point needs to be a special value, set it
         if pour_value is not None:
             catch[y, x] = pour_value
@@ -737,7 +768,8 @@ def _dinf_catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4,
         return catch
 
     def accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                     nodata_out=0., efficiency=None, routing='d8', cycle_size=1, **kwargs):
+                     nodata_out=0., efficiency=None, routing='d8', cycle_size=1,
+                     algorithm='iterative', **kwargs):
         """
         Generates a flow accumulation raster. If no weights are provided, the value of each cell
         is equal to the number of upstream cells. If weights are provided, the value of each cell
@@ -768,6 +800,10 @@ def accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
                      that d-infinity routing can generate cycles that will cause
                      the accumulation algorithm to abort. These cycles are removed prior
                      to running the d-infinity accumulation algorithm.)
+        algorithm : str
+                    Algorithm type to use:
+                    'iterative' : Use an iterative algorithm (recommended).
+                    'recursive' : Use a recursive algorithm.
 
         Additional keyword arguments (**kwargs) are passed to self.view.
 
@@ -795,16 +831,16 @@ def accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         if routing.lower() == 'd8':
             acc = self._d8_accumulation(fdir, weights=weights, dirmap=dirmap,
                                         nodata_out=nodata_out,
-                                        efficiency=efficiency)
+                                        efficiency=efficiency, algorithm=algorithm)
         elif routing.lower() == 'dinf':
             acc = self._dinf_accumulation(fdir, weights=weights, dirmap=dirmap,
                                           nodata_out=nodata_out,
                                           efficiency=efficiency,
-                                          cycle_size=cycle_size)
+                                          cycle_size=cycle_size, algorithm=algorithm)
         return acc
 
     def _d8_accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                         nodata_out=0., efficiency=None, **kwargs):
+                         nodata_out=0., efficiency=None, algorithm='iterative', **kwargs):
         # Find nodata cells and invalid cells
         nodata_cells = self._get_nodata_cells(fdir)
         invalid_cells = ~np.in1d(fdir.ravel(), dirmap).reshape(fdir.shape)
@@ -829,17 +865,31 @@ def _d8_accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16,
         indegree = np.bincount(endnodes.ravel(), minlength=fdir.size).astype(np.uint8)
         # Set starting nodes to those with no predecessors
         startnodes = startnodes[(indegree == 0)]
-        # Compute accumulation
+        # Compute accumulation for no efficiency case
         if efficiency is None:
-            acc = _self._d8_accumulation_numba(acc, endnodes, indegree, startnodes)
+            if algorithm.lower() == 'iterative':
+                acc = _self._d8_accumulation_iter_numba(acc, endnodes, indegree, startnodes)
+            elif algorithm.lower() == 'recursive':
+                acc = _self._d8_accumulation_recur_numba(acc, endnodes, indegree, startnodes)
+            else:
+                raise ValueError('Algorithm must be `iterative` or `recursive`.')
+        # Compute accumulation for efficiency case
         else:
-            acc = _self._d8_accumulation_eff_numba(acc, endnodes, indegree, startnodes, eff)
+            if algorithm.lower() == 'iterative':
+                acc = _self._d8_accumulation_eff_iter_numba(acc, endnodes, indegree,
+                                                            startnodes, eff)
+            elif algorithm.lower() == 'recursive':
+                acc = _self._d8_accumulation_eff_recur_numba(acc, endnodes, indegree,
+                                                             startnodes, eff)
+            else:
+                raise ValueError('Algorithm must be `iterative` or `recursive`.')
         acc = self._output_handler(data=acc, viewfinder=fdir.viewfinder,
                                    metadata=fdir.metadata, nodata=nodata_out)
         return acc
 
     def _dinf_accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                           nodata_out=0., efficiency=None, cycle_size=1, **kwargs):
+                           nodata_out=0., efficiency=None, cycle_size=1, algorithm='iterative',
+                           **kwargs):
         # Find nodata cells and invalid cells
         nodata_cells = self._get_nodata_cells(fdir)
         # Split d-infinity grid
@@ -866,20 +916,37 @@ def _dinf_accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16
         indegree = (indegree_0 + indegree_1).astype(np.uint8)
         # Set starting nodes to those with no predecessors
         startnodes = startnodes[(indegree == 0)]
-        # Compute accumulation
+        # Compute accumulation for no efficiency case
         if efficiency is None:
-            acc = _self._dinf_accumulation_numba(acc, endnodes_0, endnodes_1, indegree,
-                                           startnodes, prop_0, prop_1)
+            if algorithm.lower() == 'iterative':
+                acc = _self._dinf_accumulation_iter_numba(acc, endnodes_0, endnodes_1,
+                                                          indegree, startnodes, prop_0,
+                                                          prop_1)
+            elif algorithm.lower() == 'recursive':
+                acc = _self._dinf_accumulation_recur_numba(acc, endnodes_0, endnodes_1,
+                                                           indegree, startnodes, prop_0,
+                                                           prop_1)
+            else:
+                raise ValueError('Algorithm must be `iterative` or `recursive`.')
+        # Compute accumulation for efficiency case
         else:
-            acc = _self._dinf_accumulation_eff_numba(acc, endnodes_0, endnodes_1, indegree,
-                                               startnodes, prop_0, prop_1, eff)
+            if algorithm.lower() == 'iterative':
+                acc = _self._dinf_accumulation_eff_iter_numba(acc, endnodes_0, endnodes_1,
+                                                            indegree, startnodes, prop_0,
+                                                            prop_1, eff)
+            elif algorithm.lower() == 'recursive':
+                acc = _self._dinf_accumulation_eff_recur_numba(acc, endnodes_0, endnodes_1,
+                                                               indegree, startnodes, prop_0,
+                                                               prop_1, eff)
+            else:
+                raise ValueError('Algorithm must be `iterative` or `recursive`.')
         acc = self._output_handler(data=acc, viewfinder=fdir.viewfinder,
                                    metadata=fdir.metadata, nodata=nodata_out)
         return acc
 
     def distance_to_outlet(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
                            nodata_out=np.nan, routing='d8', method='shortest',
-                           xytype='coordinate', snap='corner', **kwargs):
+                           xytype='coordinate', snap='corner', algorithm='iterative', **kwargs):
         """
         Generates a raster representing the (weighted) topological distance from each cell
         to the outlet, moving downstream.
@@ -918,6 +985,10 @@ def distance_to_outlet(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4,
                Function to use on array for indexing:
                'corner' : numpy.around()
                'center' : numpy.floor()
+        algorithm : str
+                    Algorithm type to use:
+                    'iterative' : Use an iterative algorithm (recommended).
+                    'recursive' : Use a recursive algorithm.
 
         Additional keyword arguments (**kwargs) are passed to self.view.
 
@@ -951,18 +1022,17 @@ def distance_to_outlet(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4,
             dist = self._d8_flow_distance(x=x, y=y, fdir=fdir, weights=weights,
                                           dirmap=dirmap, nodata_out=nodata_out,
                                           method=method, xytype=xytype,
-                                          snap=snap)
+                                          snap=snap, algorithm=algorithm)
         elif routing.lower() == 'dinf':
             dist = self._dinf_flow_distance(x=x, y=y, fdir=fdir, weights=weights,
                                             dirmap=dirmap, nodata_out=nodata_out,
                                             method=method, xytype=xytype,
-                                            snap=snap)
+                                            snap=snap, algorithm=algorithm)
         return dist
 
-    def _d8_flow_distance(self, x, y, fdir, weights=None,
-                          dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                          nodata_out=np.nan, method='shortest',
-                          xytype='coordinate', snap='corner', **kwargs):
+    def _d8_flow_distance(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
+                          nodata_out=np.nan, method='shortest', xytype='coordinate',
+                          snap='corner', algorithm='iterative', **kwargs):
         # Find nodata cells and invalid cells
         nodata_cells = self._get_nodata_cells(fdir)
         invalid_cells = ~np.in1d(fdir.ravel(), dirmap).reshape(fdir.shape)
@@ -973,15 +1043,19 @@ def _d8_flow_distance(self, x, y, fdir, weights=None,
             x, y = self.nearest_cell(x, y, fdir.affine, snap)
         if weights is None:
             weights = (~nodata_cells).reshape(fdir.shape).astype(np.float64)
-        dist = _self._d8_flow_distance_numba(fdir, weights, (y, x), dirmap)
+        if algorithm.lower() == 'iterative':
+            dist = _self._d8_flow_distance_iter_numba(fdir, weights, (y, x), dirmap)
+        elif algorithm.lower() == 'recursive':
+            dist = _self._d8_flow_distance_recur_numba(fdir, weights, (y, x), dirmap)
+        else:
+            raise ValueError('Algorithm must be `iterative` or `recursive`.')
         dist = self._output_handler(data=dist, viewfinder=fdir.viewfinder,
                                    metadata=fdir.metadata, nodata=nodata_out)
         return dist
 
-    def _dinf_flow_distance(self, x, y, fdir, weights=None,
-                          dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                          nodata_out=np.nan, method='shortest',
-                          xytype='coordinate', snap='corner', **kwargs):
+    def _dinf_flow_distance(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
+                            nodata_out=np.nan, method='shortest', xytype='coordinate',
+                            snap='corner', algorithm='iterative', **kwargs):
         # Find nodata cells
         nodata_cells = self._get_nodata_cells(fdir)
         # Split d-infinity grid
@@ -995,8 +1069,14 @@ def _dinf_flow_distance(self, x, y, fdir, weights=None,
             weights_0 = (~nodata_cells).reshape(fdir.shape).astype(np.float64)
             weights_1 = weights_0
         if method.lower() == 'shortest':
-            dist = _self._dinf_flow_distance_numba(fdir_0, fdir_1, weights_0,
-                                                   weights_1, (y, x), dirmap)
+            if algorithm.lower() == 'iterative':
+                dist = _self._dinf_flow_distance_iter_numba(fdir_0, fdir_1, weights_0,
+                                                            weights_1, (y, x), dirmap)
+            elif algorithm.lower() == 'recursive':
+                dist = _self._dinf_flow_distance_recur_numba(fdir_0, fdir_1, weights_0,
+                                                             weights_1, (y, x), dirmap)
+            else:
+                raise ValueError('Algorithm must be `iterative` or `recursive`.')
         else:
             raise NotImplementedError("Only implemented for shortest path distance.")
         # Prepare output
@@ -1005,7 +1085,8 @@ def _dinf_flow_distance(self, x, y, fdir, weights=None,
         return dist
 
     def compute_hand(self, fdir, dem, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                     nodata_out=None, routing='d8', return_index=False, **kwargs):
+                     nodata_out=None, routing='d8', return_index=False, algorithm='iterative',
+                     **kwargs):
         """
         Computes the height above nearest drainage (HAND), based on a flow direction grid,
         a digital elevation grid, and a grid containing the locations of drainage channels.
@@ -1035,6 +1116,10 @@ def compute_hand(self, fdir, dem, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
                          of the (topologically) nearest channel cell.
                        - If False, return a Raster where each cell indicates the elevation
                          above the (topologically) nearest channel cell.
+        algorithm : str
+                    Algorithm type to use:
+                    'iterative' : Use an iterative algorithm (recommended).
+                    'recursive' : Use a recursive algorithm.
 
         Additional keyword arguments (**kwargs) are passed to self.view.
 
@@ -1066,11 +1151,13 @@ def compute_hand(self, fdir, dem, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
                 nodata_out = np.nan
         # Compute height above nearest drainage
         if routing.lower() == 'd8':
-            hand = self._d8_compute_hand(fdir=fdir, mask=mask,
-                                         dirmap=dirmap, nodata_out=nodata_out)
+            hand = self._d8_compute_hand(fdir=fdir, mask=mask, dirmap=dirmap,
+                                         nodata_out=nodata_out,
+                                         algorithm=algorithm)
         elif routing.lower() == 'dinf':
             hand = self._dinf_compute_hand(fdir=fdir, mask=mask,
-                                           nodata_out=nodata_out)
+                                           nodata_out=nodata_out,
+                                           algorithm=algorithm)
         # If index is not desired, return heights
         if not return_index:
             hand = _self._assign_hand_heights_numba(hand, dem, nodata_out)
@@ -1079,7 +1166,7 @@ def compute_hand(self, fdir, dem, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         return hand
 
     def _d8_compute_hand(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                         nodata_out=-1):
+                         nodata_out=-1, algorithm='iterative'):
         # Find nodata cells and invalid cells
         nodata_cells = self._get_nodata_cells(fdir)
         invalid_cells = ~np.in1d(fdir.ravel(), dirmap).reshape(fdir.shape)
@@ -1088,13 +1175,18 @@ def _d8_compute_hand(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         fdir[invalid_cells] = 0
         dirleft, dirright, dirtop, dirbottom = self._pop_rim(fdir, nodata=0)
         maskleft, maskright, masktop, maskbottom = self._pop_rim(mask, nodata=False)
-        hand = _self._d8_hand_iter_numba(fdir, mask, dirmap)
+        if algorithm.lower() == 'iterative':
+            hand = _self._d8_hand_iter_numba(fdir, mask, dirmap)
+        elif algorithm.lower() == 'recursive':
+            hand = _self._d8_hand_recur_numba(fdir, mask, dirmap)
+        else:
+            raise ValueError('Algorithm must be `iterative` or `recursive`.')
         hand = self._output_handler(data=hand, viewfinder=fdir.viewfinder,
                                     metadata=fdir.metadata, nodata=-1)
         return hand
 
     def _dinf_compute_hand(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                         nodata_out=-1):
+                           nodata_out=-1, algorithm='iterative'):
         # Get nodata cells
         nodata_cells = self._get_nodata_cells(fdir)
         # Split dinf flowdir
@@ -1105,13 +1197,18 @@ def _dinf_compute_hand(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         dirleft_1, dirright_1, dirtop_1, dirbottom_1 = self._pop_rim(fdir_1,
                                                                      nodata=0)
         maskleft, maskright, masktop, maskbottom = self._pop_rim(mask, nodata=False)
-        hand = _self._dinf_hand_iter_numba(fdir_0, fdir_1, mask, dirmap)
+        if algorithm.lower() == 'iterative':
+            hand = _self._dinf_hand_iter_numba(fdir_0, fdir_1, mask, dirmap)
+        elif algorithm.lower() == 'recursive':
+            hand = _self._dinf_hand_recur_numba(fdir_0, fdir_1, mask, dirmap)
+        else:
+            raise ValueError('Algorithm must be `iterative` or `recursive`.')
         hand = self._output_handler(data=hand, viewfinder=fdir.viewfinder,
                                     metadata=fdir.metadata, nodata=-1)
         return hand
 
     def extract_river_network(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                              routing='d8', **kwargs):
+                              routing='d8', algorithm='iterative', **kwargs):
         """
         Generates river segments from accumulation and flow_direction arrays.
 
@@ -1128,6 +1225,10 @@ def extract_river_network(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32)
         routing : str
                   Routing algorithm to use:
                   'd8'   : D8 flow directions
+        algorithm : str
+                    Algorithm type to use:
+                    'iterative' : Use an iterative algorithm (recommended).
+                    'recursive' : Use a recursive algorithm.
 
         Additional keyword arguments (**kwargs) are passed to self.view.
 
@@ -1159,7 +1260,14 @@ def extract_river_network(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32)
         indegree = np.bincount(endnodes.ravel(), minlength=fdir.size).astype(np.uint8)
         orig_indegree = np.copy(indegree)
         startnodes = startnodes[(indegree == 0)]
-        profiles = _self._d8_stream_network_numba(endnodes, indegree, orig_indegree, startnodes)
+        if algorithm.lower() == 'iterative':
+            profiles = _self._d8_stream_network_iter_numba(endnodes, indegree,
+                                                           orig_indegree, startnodes)
+        elif algorithm.lower() == 'recursive':
+            profiles = _self._d8_stream_network_recur_numba(endnodes, indegree,
+                                                            orig_indegree, startnodes)
+        else:
+            raise ValueError('Algorithm must be `iterative` or `recursive`.')
         # Fill geojson dict with profiles
         featurelist = []
         for index, profile in enumerate(profiles):
@@ -1171,7 +1279,7 @@ def extract_river_network(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32)
         return geo
 
     def stream_order(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                     nodata_out=0, routing='d8', **kwargs):
+                     nodata_out=0, routing='d8', algorithm='iterative', **kwargs):
         """
         Computes the Strahler stream order.
 
@@ -1190,6 +1298,10 @@ def stream_order(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         routing : str
                   Routing algorithm to use:
                   'd8'   : D8 flow directions
+        algorithm : str
+                    Algorithm type to use:
+                    'iterative' : Use an iterative algorithm (recommended).
+                    'recursive' : Use a recursive algorithm.
 
         Additional keyword arguments (**kwargs) are passed to self.view.
 
@@ -1223,14 +1335,20 @@ def stream_order(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         min_order = np.full(fdir.shape, np.iinfo(np.int64).max, dtype=np.int64)
         max_order = np.ones(fdir.shape, dtype=np.int64)
         order = np.where(mask, 1, 0).astype(np.int64).reshape(fdir.shape)
-        order = _self._d8_streamorder_numba(min_order, max_order, order, endnodes,
-                                        indegree, orig_indegree, startnodes)
+        if algorithm.lower() == 'iterative':
+            order = _self._d8_streamorder_iter_numba(min_order, max_order, order, endnodes,
+                                                     indegree, orig_indegree, startnodes)
+        elif algorithm.lower() == 'recursive':
+            order = _self._d8_streamorder_recur_numba(min_order, max_order, order, endnodes,
+                                                      indegree, orig_indegree, startnodes)
+        else:
+            raise ValueError('Algorithm must be `iterative` or `recursive`.')
         order = self._output_handler(data=order, viewfinder=fdir.viewfinder,
                                      metadata=fdir.metadata, nodata=nodata_out)
         return order
 
     def distance_to_ridge(self, fdir, mask, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                          nodata_out=0, routing='d8', **kwargs):
+                          nodata_out=0, routing='d8', algorithm='iterative', **kwargs):
         """
         Generates a raster representing the (weighted) topological distance from each cell
         to its originating drainage divide, moving upstream.
@@ -1250,6 +1368,10 @@ def distance_to_ridge(self, fdir, mask, weights=None, dirmap=(64, 128, 1, 2, 4,
         routing : str
                   Routing algorithm to use:
                   'd8'   : D8 flow directions
+        algorithm : str
+                    Algorithm type to use:
+                    'iterative' : Use an iterative algorithm (recommended).
+                    'recursive' : Use a recursive algorithm.
 
         Additional keyword arguments (**kwargs) are passed to self.view.
 
@@ -1290,8 +1412,16 @@ def distance_to_ridge(self, fdir, mask, weights=None, dirmap=(64, 128, 1, 2, 4,
         min_order = np.full(fdir.shape, np.iinfo(np.int64).max, dtype=np.int64)
         max_order = np.ones(fdir.shape, dtype=np.int64)
         rdist = np.zeros(fdir.shape, dtype=np.float64)
-        rdist = _self._d8_reverse_distance_numba(min_order, max_order, rdist,
-                                            endnodes, indegree, startnodes, weights)
+        if algorithm.lower() == 'iterative':
+            rdist = _self._d8_reverse_distance_iter_numba(min_order, max_order, rdist,
+                                                        endnodes, indegree, startnodes,
+                                                        weights)
+        elif algorithm.lower() == 'recursive':
+            rdist = _self._d8_reverse_distance_recur_numba(min_order, max_order, rdist,
+                                                           endnodes, indegree, startnodes,
+                                                           weights)
+        else:
+            raise ValueError('Algorithm must be `iterative` or `recursive`.')
         rdist = self._output_handler(data=rdist, viewfinder=fdir.viewfinder,
                                      metadata=fdir.metadata, nodata=nodata_out)
         return rdist
diff --git a/pysheds/sview.py b/pysheds/sview.py
index b21a737..0822a6a 100644
--- a/pysheds/sview.py
+++ b/pysheds/sview.py
@@ -129,25 +129,25 @@ def viewfinder(self, new_viewfinder):
         self._viewfinder = new_viewfinder
 
     @property
-    def bbox(self):
-        return self.viewfinder.bbox
-
-    @property
-    def coords(self):
-        return self.viewfinder.coords
-
-    @property
-    def axes(self):
-        return self.viewfinder.axes
+    def affine(self):
+        return self.viewfinder.affine
 
-    @property
-    def view_shape(self):
-        return self.viewfinder.shape
+    @affine.setter
+    def affine(self, new_affine):
+        self.viewfinder.affine = new_affine
 
     @property
     def mask(self):
         return self.viewfinder.mask
 
+    @mask.setter
+    def mask(self, new_mask):
+        try:
+            assert (new_mask.shape == self.shape)
+        except:
+            raise ValueError('`mask` must have same dimensions as Raster.')
+        self.viewfinder.mask = new_mask
+
     @property
     def nodata(self):
         return self.viewfinder.nodata
@@ -160,9 +160,21 @@ def nodata(self, new_nodata):
     def crs(self):
         return self.viewfinder.crs
 
+    @crs.setter
+    def crs(self, new_crs):
+        self.viewfinder.crs = new_crs
+
+    @property
+    def bbox(self):
+        return self.viewfinder.bbox
+
     @property
-    def view_size(self):
-        return np.prod(self.viewfinder.shape)
+    def coords(self):
+        return self.viewfinder.coords
+
+    @property
+    def axes(self):
+        return self.viewfinder.axes
 
     @property
     def extent(self):
@@ -170,10 +182,6 @@ def extent(self):
         extent = (bbox[0], bbox[2], bbox[1], bbox[3])
         return extent
 
-    @property
-    def affine(self):
-        return self.viewfinder.affine
-
     @property
     def properties(self):
         property_dict = {
@@ -245,11 +253,13 @@ class ViewFinder():
     Attributes
     ==========
     affine : Affine transformation matrix (uses affine module).
-    shape : The shape of the raster (number of rows, number of columns).
+    shape : The shape of the raster (number of rows, number of columns). Note: if
+            the shape of the provided `mask` differs, the `mask` shape will be used.
     crs : The coordinate reference system.
     nodata : The value indicating `no data`.
     mask : A boolean array used to mask raster cells; may be used to indicate
-           which cells lie inside a catchment.
+           which cells lie inside a catchment. Note: if the shape of the provided
+           `mask` is different from `shape`, the shape of `mask` will be used.
     bbox : The bounding box of the raster (xmin, ymin, xmax, ymax).
     extent : The extent of the raster (xmin, xmax, ymin, ymax).
     size : The number of cells in the raster.
@@ -265,7 +275,6 @@ class ViewFinder():
     def __init__(self, affine=Affine(1., 0., 0., 0., 1., 0.), shape=(1,1),
                  nodata=0, mask=None, crs=pyproj.Proj(_pyproj_init)):
         self.affine = affine
-        self.shape = shape
         self.crs = crs
         self.nodata = nodata
         if mask is None:
@@ -308,18 +317,11 @@ def affine(self, new_affine):
 
     @property
     def shape(self):
-        return self._shape
+        return self.mask.shape
 
     @shape.setter
     def shape(self, new_shape):
-        try:
-            assert len(new_shape) == 2
-            assert isinstance(new_shape[0], int)
-            assert isinstance(new_shape[1], int)
-        except:
-            raise ValueError('`shape` must be an integer sequence of length 2.')
-        new_shape = tuple(new_shape)
-        self._shape = new_shape
+        raise AttributeError('Viewfinder shape cannot be changed after instantiation.')
 
     @property
     def mask(self):
@@ -327,15 +329,11 @@ def mask(self):
 
     @mask.setter
     def mask(self, new_mask):
-        try:
-            assert (new_mask.shape == self.shape)
-        except:
-            raise ValueError('`mask` shape must be the same as `self.shape`')
         try:
             assert (np.min_scalar_type(new_mask) <= np.dtype(np.bool8))
         except:
             raise TypeError('`mask` must be of boolean type')
-        new_mask = new_mask.astype(np.bool8)
+        new_mask = np.asarray(new_mask).astype(np.bool8)
         self._mask = new_mask
 
     @property
@@ -443,7 +441,7 @@ def __init__(self):
 
     @classmethod
     def view(cls, data, target_view, data_view=None, interpolation='nearest',
-             apply_input_mask=False, apply_output_mask=True,
+             apply_input_mask=False, apply_output_mask=True, inherit_nodata=True,
              affine=None, shape=None, crs=None, mask=None, nodata=None,
              dtype=None, inherit_metadata=True, new_metadata={}):
         """
@@ -467,6 +465,9 @@ def view(cls, data, target_view, data_view=None, interpolation='nearest',
                            If True, mask the input Raster according to data.mask.
         apply_output_mask : bool
                            If True, mask the output Raster according to grid.mask.
+        inherit_nodata : bool
+                         If True, output Raster inherits `nodata` value from `data_view`.
+                         If False, output Raster uses `nodata` value from `target_view`.
         affine : affine.Affine
                  Affine transformation matrix (overrides target_view.affine)
         shape : tuple of ints (length 2)
@@ -496,6 +497,10 @@ def view(cls, data, target_view, data_view=None, interpolation='nearest',
             except:
                 raise TypeError('`data` must be a Raster instance.')
             data_view = data.viewfinder
+        # By default, use dataset's `nodata` value
+        if nodata is None:
+            if inherit_nodata:
+                nodata = data_view.nodata
         # Override parameters of target view if desired
         target_view = cls._override_target_view(target_view,
                                                 affine=affine,
diff --git a/setup.py b/setup.py
index 4798396..98ce916 100644
--- a/setup.py
+++ b/setup.py
@@ -3,7 +3,7 @@
 from setuptools import setup
 
 setup(name='pysheds',
-      version='0.3',
+      version='0.3.1',
       description='🌎 Simple and fast watershed delineation in python.',
       author='Matt Bartos',
       author_email='mdbartos@umich.edu',
diff --git a/tests/test_grid.py b/tests/test_grid.py
index 8349c28..4fe1025 100644
--- a/tests/test_grid.py
+++ b/tests/test_grid.py
@@ -3,7 +3,7 @@
 import warnings
 import numpy as np
 from pysheds.grid import Grid
-from pysheds.sview import Raster
+from pysheds.view import Raster, ViewFinder
 from pysheds.rfsm import RFSM
 
 current_dir = os.path.dirname(os.path.realpath(__file__))
@@ -70,7 +70,7 @@ def test_constructors():
     newgrid = grid.from_ascii(dir_path, dtype=np.uint8, crs=crs)
     new_fdir = grid.read_ascii(dir_path, dtype=np.uint8, crs=crs)
     assert((fdir == new_fdir).all())
-    del newgrid
+    newgrid = Grid(viewfinder=fdir.viewfinder)
 
 def test_dtype():
     assert(fdir.dtype == np.uint8)
@@ -152,6 +152,10 @@ def test_computed_fdir_catch():
     catch_dinf = grid.catchment(x, y, fdir_dinf, dirmap=dirmap, routing='dinf',
                                 xytype='coordinate')
     assert(np.count_nonzero(catch_dinf) > 11300)
+    catch_d8_recur = grid.catchment(x, y, fdir_d8, dirmap=dirmap, routing='d8',
+                                    xytype='coordinate', algorithm='recursive')
+    catch_dinf_recur = grid.catchment(x, y, fdir_dinf, dirmap=dirmap, routing='dinf',
+                                      xytype='coordinate', algorithm='recursive')
 
 def test_accumulation():
     fdir = d.fdir
@@ -186,12 +190,16 @@ def test_accumulation():
     acc_dinf = grid.accumulation(fdir_dinf, dirmap=dirmap, routing='dinf')
     assert(acc_dinf.max() > 11300)
     # #set nodata to 1
-    eff = grid.view(dinf_eff)
-    eff[eff==dinf_eff.nodata] = 1
-    acc_dinf_eff = grid.accumulation(fdir_dinf, dirmap=dirmap,
-                                     routing='dinf', efficiency=eff)
+    # eff = grid.view(dinf_eff)
+    # eff[eff==dinf_eff.nodata] = 1
+    # acc_dinf_eff = grid.accumulation(fdir_dinf, dirmap=dirmap,
+    #                                  routing='dinf', efficiency=eff)
     # pos = np.where(grid.dinf_acc==grid.dinf_acc.max())
     # assert(np.round(grid.dinf_acc[pos] / grid.dinf_acc_eff[pos]) == 4.)
+    acc_d8_recur = grid.accumulation(fdir_d8, dirmap=dirmap, routing='d8',
+                                     algorithm='recursive')
+    acc_dinf_recur = grid.accumulation(fdir_dinf, dirmap=dirmap, routing='dinf',
+                                       algorithm='recursive')
     d.acc = acc
 
 def test_hand():
@@ -201,6 +209,10 @@ def test_hand():
     fdir_dinf = d.fdir_dinf
     hand_d8 = grid.compute_hand(fdir, dem, acc > 100, routing='d8')
     hand_dinf = grid.compute_hand(fdir_dinf, dem, acc > 100, routing='dinf')
+    hand_d8_recur = grid.compute_hand(fdir, dem, acc > 100, routing='d8',
+                                      algorithm='recursive')
+    hand_dinf_recur = grid.compute_hand(fdir_dinf, dem, acc > 100, routing='dinf',
+                                        algorithm='recursive')
 
 def test_distance_to_outlet():
     fdir = d.fdir
@@ -219,16 +231,23 @@ def test_distance_to_outlet():
                        dirmap=dirmap, xytype='label')
     grid.distance_to_outlet(x, y, fdir_dinf, dirmap=dirmap, weights=weights,
                        routing='dinf', xytype='label')
+    # Test recursive
+    grid.distance_to_outlet(x, y, fdir, dirmap=dirmap, xytype='coordinate',
+                            routing='d8', algorithm='recursive')
+    grid.distance_to_outlet(x, y, fdir_dinf, dirmap=dirmap, xytype='coordinate',
+                            routing='dinf', algorithm='recursive')
 
 def test_stream_order():
     fdir = d.fdir
     acc = d.acc
     order = grid.stream_order(fdir, acc > 100)
+    order = grid.stream_order(fdir, acc > 100, algorithm='recursive')
 
 def test_distance_to_ridge():
     fdir = d.fdir
     acc = d.acc
     order = grid.distance_to_ridge(fdir, acc > 100)
+    order = grid.distance_to_ridge(fdir, acc > 100, algorithm='recursive')
 
 def test_cell_dh():
     fdir = d.fdir
@@ -322,6 +341,7 @@ def test_extract_river_network():
     grid.clip_to(catch)
     rivers = grid.extract_river_network(catch, acc > 20)
     assert(isinstance(rivers, dict))
+    grid.extract_river_network(catch, acc > 20, algorithm='recursive')
     # TODO: Need more checks here. Check if endnodes equals next startnode
 
 def test_view_methods():
@@ -398,3 +418,8 @@ def test_polygonize_rasterize():
 #     waterlevel = rfsm.compute_waterlevel(input_vol)
 #     end_vol = (area*np.where(waterlevel, waterlevel - dem, 0)).sum()
 #     assert np.allclose(end_vol, input_vol.sum())
+
+def test_misc():
+    dem = d.dem
+    l, r, t, b = grid._pop_rim(dem, nodata=0)
+    grid._replace_rim(dem, l, r, t, b)