From a22592280da6f535fdf849577aa1c25f81ad7a9c Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Sat, 15 Jan 2022 00:39:44 -0600
Subject: [PATCH 01/27] Add MFD flowdir and catchment computation

---
 pysheds/_sgrid.py | 60 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 60 insertions(+)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index b60c5dc..93242b0 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -252,6 +252,35 @@ def _angle_to_d8_numba(angles, dirmap, nodata_cells):
         props_1.flat[i] = prop_1
     return fdirs_0, fdirs_1, props_0, props_1
 
+@njit(float64[:,:,:](float64[:,:], float64, float64, boolean[:,:], int64, int64),
+      parallel=True,
+      cache=True)
+def _mfd_flowdir_numba(dem, dx, dy, nodata_cells, nodata_out, p=1):
+    m, n = dem.shape
+    fdir = np.zeros((8, m, n), dtype=np.float64)
+    row_offsets = np.array([-1, -1, 0, 1, 1, 1, 0, -1])
+    col_offsets = np.array([0, 1, 1, 1, 0, -1, -1, -1])
+    dd = np.sqrt(dx**2 + dy**2)
+    distances = np.array([dy, dd, dx, dd, dy, dd, dx, dd])
+    for i in prange(1, m - 1):
+        for j in prange(1, n - 1):
+            if nodata_cells[i, j]:
+                fdir[:, i, j] = nodata_out
+            else:
+                elev = dem[i, j]
+                den = 0.
+                for k in range(8):
+                    row_offset = row_offsets[k]
+                    col_offset = col_offsets[k]
+                    distance = distances[k]
+                    num = (elev - dem[i + row_offset, j + col_offset])**p / distance
+                    if num > 0:
+                        fdir[k, i, j] = num
+                        den += num
+                if den > 0:
+                    fdir[:, i, j] /= den
+    return fdir
+
 # Functions for 'catchment'
 
 @njit(void(int64, boolean[:,:], int64[:,:], int64[:], int64[:]),
@@ -373,6 +402,37 @@ def _dinf_catchment_iter_numba(fdir_0, fdir_1, pour_point, dirmap):
                     queue.append(neighbor)
     return catch
 
+@njit(boolean[:,:](float64[:,:,:], UniTuple(int64, 2)),
+      cache=True)
+def _mfd_catchment_iter_numba(fdir, pour_point):
+    _, m, n = fdir.shape
+    mn = m * n
+    catch = np.zeros((m, n), dtype=np.bool8)
+    i, j = pour_point
+    ix = (i * n) + j
+    offsets = np.array([-n, 1 - n, 1,
+                        1 + n, n, - 1 + n,
+                        - 1, - 1 - n])
+    r_dirmap = np.array([4, 5, 6, 7,
+                         0, 1, 2, 3])
+    queue = [ix]
+    while queue:
+        parent = queue.pop()
+        catch.flat[parent] = True
+        neighbors = offsets + parent
+        for k in range(8):
+            neighbor = neighbors[k]
+            neighbor_dir = r_dirmap[k]
+            visited = catch.flat[neighbor]
+            if visited:
+                continue
+            else:
+                kix = neighbor + (neighbor_dir * mn)
+                points_to = fdir.flat[kix] > 0.
+                if points_to:
+                    queue.append(neighbor)
+    return catch
+
 # Functions for 'accumulation'
 
 @njit(void(int64, int64, float64[:,:], int64[:,:], uint8[:]),

From 4b117d189884d003266e569abd8c35948b96e55d Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Sat, 15 Jan 2022 01:33:00 -0600
Subject: [PATCH 02/27] First working version of mfd flow accumulation

---
 pysheds/_sgrid.py | 37 +++++++++++++++++++++++++++++++++++++
 1 file changed, 37 insertions(+)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index 93242b0..f197bbe 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -642,6 +642,43 @@ def _dinf_accumulation_eff_iter_numba(acc, fdir_0, fdir_1, indegree, startnodes,
                     queue.append(endnode_1)
     return acc
 
+@njit(parallel=True)
+def _mfd_bincount(fdir):
+    p, m, n = fdir.shape
+    mn = m * n
+    out = np.zeros(mn, dtype=np.uint8)
+    for i in range(p):
+        fdir_i = fdir[i]
+        for j in prange(mn):
+            endnode = fdir_i.flat[j]
+            if endnode != j:
+                out[endnode] += 1
+    return out
+
+@njit
+def _mfd_accumulation_iter_numba(acc, fdir, props, indegree, startnodes):
+    n = startnodes.size
+    queue = [0]
+    _ = queue.pop()
+    for k in range(n):
+        startnode = startnodes.flat[k]
+        queue.append(startnode)
+        while queue:
+            startnode = queue.pop()
+            for i in range(8):
+                fdir_i = fdir[i]
+                props_i = props[i]
+                endnode = fdir_i.flat[startnode]
+                if endnode == startnode:
+                    continue
+                else:
+                    prop = props_i.flat[startnode]
+                    acc.flat[endnode] += (prop * acc.flat[startnode])
+                    indegree.flat[endnode] -= 1
+                    if (indegree.flat[endnode] == 0):
+                        queue.append(endnode)
+    return acc
+
 # Functions for 'flow_distance'
 
 @njit(void(int64, int64[:,:], boolean[:,:], float64[:,:], float64[:,:],

From fa4d74a8ec8a655ef5b46ce8a0cb9e057bc2891d Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Sat, 15 Jan 2022 19:15:07 -0600
Subject: [PATCH 03/27] Add mfd distance_to_outlet

---
 pysheds/_sgrid.py | 36 ++++++++++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index f197bbe..a4dd121 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -831,6 +831,42 @@ def _dinf_flow_distance_iter_numba(fdir_0, fdir_1, weights_0, weights_1,
                     queue.append(neighbor)
     return dist
 
+# TODO: Weights should actually by (8, m, n)
+# neighbor_dist = parent_dist + weights.flat[kix]
+@njit(float64[:,:](float64[:,:,:], UniTuple(int64, 2), float64[:,:]),
+      cache=True)
+def _mfd_flow_distance_iter_numba(fdir, pour_point, weights):
+    _, m, n = fdir.shape
+    mn = m * n
+    dist = np.full((m, n), np.inf, dtype=np.float64)
+    i, j = pour_point
+    ix = (i * n) + j
+    offsets = np.array([-n, 1 - n, 1,
+                        1 + n, n, - 1 + n,
+                        - 1, - 1 - n])
+    r_dirmap = np.array([4, 5, 6, 7,
+                         0, 1, 2, 3])
+    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]
+            neighbor_dir = r_dirmap[k]
+            current_neighbor_dist = dist.flat[neighbor]
+            kix = neighbor + (neighbor_dir * mn)
+            points_to = fdir.flat[kix] > 0.
+            if points_to:
+                neighbor_dist = parent_dist + weights.flat[neighbor]
+                if (neighbor_dist < current_neighbor_dist):
+                    dist.flat[neighbor] = neighbor_dist
+                    queue.append(neighbor)
+    return dist
+
+# Functions for 'reverse_flow_distance'
+
 @njit(void(int64, int64, int64[:,:], int64[:,:], float64[:,:],
            int64[:,:], uint8[:], float64[:,:]),
       cache=True)

From 9e8f34eb22251f0b799ce7c38a10eaaf9a61cbed Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Sat, 15 Jan 2022 20:51:23 -0600
Subject: [PATCH 04/27] Add heap to distance_to_outlet

---
 pysheds/_sgrid.py | 43 +++++++++++++++++++++++++++++++++++++++++++
 pysheds/sgrid.py  |  2 +-
 2 files changed, 44 insertions(+), 1 deletion(-)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index a4dd121..15c3ba7 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -1,3 +1,4 @@
+from heapq import heappop, heappush
 import numpy as np
 from numba import njit, prange
 from numba.types import float64, int64, uint32, uint16, uint8, boolean, UniTuple, Tuple, List, void
@@ -831,6 +832,48 @@ def _dinf_flow_distance_iter_numba(fdir_0, fdir_1, weights_0, weights_1,
                     queue.append(neighbor)
     return dist
 
+@njit(float64[:,:](int64[:,:], int64[:,:], float64[:,:], float64[:,:],
+                   UniTuple(int64, 2), UniTuple(int64, 8)),
+      cache=True)
+def _dinf_flow_distance_heap_numba(fdir_0, fdir_1, weights_0, weights_1,
+                                   pour_point, dirmap):
+    dist = np.full(fdir_0.shape, np.inf, dtype=np.float64)
+    visited = np.zeros(fdir_0.shape, dtype=np.bool8)
+    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 = [(0., ix)]
+    while queue:
+        parent_dist, parent = heappop(queue)
+        visited.flat[parent] = True
+        neighbors = offsets + parent
+        for k in range(8):
+            neighbor = neighbors[k]
+            if visited.flat[neighbor]:
+                continue
+            else:
+                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
+                        heappush(queue, (neighbor_dist_0, 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
+                        heappush(queue, (neighbor_dist_1, neighbor))
+    return dist
+
 # TODO: Weights should actually by (8, m, n)
 # neighbor_dist = parent_dist + weights.flat[kix]
 @njit(float64[:,:](float64[:,:,:], UniTuple(int64, 2), float64[:,:]),
diff --git a/pysheds/sgrid.py b/pysheds/sgrid.py
index d2c2861..651ce97 100644
--- a/pysheds/sgrid.py
+++ b/pysheds/sgrid.py
@@ -1070,7 +1070,7 @@ def _dinf_flow_distance(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4
             weights_1 = weights_0
         if method.lower() == 'shortest':
             if algorithm.lower() == 'iterative':
-                dist = _self._dinf_flow_distance_iter_numba(fdir_0, fdir_1, weights_0,
+                dist = _self._dinf_flow_distance_heap_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,

From ac196e23504b3e3dbe0a20d951d0f86a13032c06 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Sun, 16 Jan 2022 00:12:47 -0600
Subject: [PATCH 05/27] Add heap to mfd distance_to_outlet

---
 pysheds/_sgrid.py | 38 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index 15c3ba7..77d84b3 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -908,6 +908,44 @@ def _mfd_flow_distance_iter_numba(fdir, pour_point, weights):
                     queue.append(neighbor)
     return dist
 
+# TODO: Weights should actually by (8, m, n)
+# neighbor_dist = parent_dist + weights.flat[kix]
+@njit(float64[:,:](float64[:,:,:], UniTuple(int64, 2), float64[:,:]),
+      cache=True)
+def _mfd_flow_distance_heap_numba(fdir, pour_point, weights):
+    _, m, n = fdir.shape
+    mn = m * n
+    dist = np.full((m, n), np.inf, dtype=np.float64)
+    visited = np.zeros((m, n), dtype=np.bool8)
+    i, j = pour_point
+    ix = (i * n) + j
+    offsets = np.array([-n, 1 - n, 1,
+                        1 + n, n, - 1 + n,
+                        - 1, - 1 - n])
+    r_dirmap = np.array([4, 5, 6, 7,
+                         0, 1, 2, 3])
+    dist.flat[ix] = 0.
+    queue = [(0., ix)]
+    while queue:
+        parent_dist, parent = heappop(queue)
+        visited.flat[parent] = True
+        neighbors = offsets + parent
+        for k in range(8):
+            neighbor = neighbors[k]
+            if visited.flat[neighbor]:
+                continue
+            else:
+                neighbor_dir = r_dirmap[k]
+                current_neighbor_dist = dist.flat[neighbor]
+                kix = neighbor + (neighbor_dir * mn)
+                points_to = fdir.flat[kix] > 0.
+                if points_to:
+                    neighbor_dist = parent_dist + weights.flat[neighbor]
+                    if (neighbor_dist < current_neighbor_dist):
+                        dist.flat[neighbor] = neighbor_dist
+                        heappush(queue, (neighbor_dist, neighbor))
+    return dist
+
 # Functions for 'reverse_flow_distance'
 
 @njit(void(int64, int64, int64[:,:], int64[:,:], float64[:,:],

From 844c938c2f571ee015f8331f6ad60a862e126286 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Sun, 16 Jan 2022 01:31:10 -0600
Subject: [PATCH 06/27] Add MFD HAND

---
 pysheds/_sgrid.py | 40 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index 77d84b3..008444b 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -1305,6 +1305,46 @@ def _dinf_hand_recur_numba(fdir_0, fdir_1, mask, dirmap):
         _dinf_hand_recursion(parent, parent, hand, offsets, r_dirmap, fdir_0, fdir_1)
     return hand
 
+@njit(int64[:,:](float64[:,:,:], boolean[:,:]),
+      cache=True)
+def _mfd_hand_iter_numba(fdir, mask):
+    _, m, n = fdir.shape
+    mn = m * n
+    offsets = np.array([-n, 1 - n, 1,
+                        1 + n, n, - 1 + n,
+                        - 1, - 1 - n])
+    r_dirmap = np.array([4, 5, 6, 7,
+                         0, 1, 2, 3])
+    hand = -np.ones((m, n), dtype=np.int64)
+    cur_queue = []
+    next_queue = []
+    for i in range(hand.size):
+        if mask.flat[i]:
+            hand.flat[i] = i
+            cur_queue.append(i)
+    while True:
+        if not cur_queue:
+            break
+        while cur_queue:
+            k = cur_queue.pop()
+            neighbors = offsets + k
+            for j in range(8):
+                neighbor = neighbors[j]
+                visited = (hand.flat[neighbor] >= 0)
+                if visited:
+                    continue
+                else:
+                    neighbor_dir = r_dirmap[j]
+                    kix = neighbor + (neighbor_dir * mn)
+                    points_to = fdir.flat[kix] > 0.
+                    if points_to:
+                        hand.flat[neighbor] = hand.flat[k]
+                        next_queue.append(neighbor)
+        while next_queue:
+            next_cell = next_queue.pop()
+            cur_queue.append(next_cell)
+    return hand
+
 @njit(float64[:,:](int64[:,:], float64[:,:], float64),
       parallel=True,
       cache=True)

From e641ec5f9d9f7f04865361f4b93147427b3afa10 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Sun, 16 Jan 2022 15:22:57 -0600
Subject: [PATCH 07/27] Add d-inf support for distance_to_ridge

---
 pysheds/_sgrid.py | 50 +++++++++++++++++++++-------
 pysheds/sgrid.py  | 83 +++++++++++++++++++++++++++++++++++++----------
 2 files changed, 105 insertions(+), 28 deletions(-)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index 008444b..85c55af 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -948,45 +948,43 @@ def _mfd_flow_distance_heap_numba(fdir, pour_point, weights):
 
 # Functions for 'reverse_flow_distance'
 
-@njit(void(int64, int64, int64[:,:], int64[:,:], float64[:,:],
+@njit(void(int64, int64, int64[:,:], float64[:,:],
            int64[:,:], uint8[:], float64[:,:]),
       cache=True)
-def _d8_reverse_distance_recursion(startnode, endnode, min_order, max_order,
+def _d8_reverse_distance_recursion(startnode, endnode, max_order,
                                    rdist, fdir, indegree, weights):
-    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
     if indegree.flat[endnode] == 0:
         rdist.flat[endnode] = max_order.flat[endnode] + weights.flat[endnode]
         new_startnode = endnode
         new_endnode = fdir.flat[new_startnode]
-        _d8_reverse_distance_recursion(new_startnode, new_endnode, min_order,
-                                       max_order, rdist, fdir, indegree, weights)
+        _d8_reverse_distance_recursion(new_startnode, new_endnode, max_order,
+                                       rdist, fdir, indegree, weights)
 
-@njit(float64[:,:](int64[:,:], int64[:,:], float64[:,:], int64[:,:],
+@njit(float64[:,:](int64[:,:], float64[:,:], int64[:,:],
                    uint8[:], int64[:], float64[:,:]),
       cache=True)
-def _d8_reverse_distance_recur_numba(min_order, max_order, rdist, fdir,
+def _d8_reverse_distance_recur_numba(max_order, rdist, fdir,
                                      indegree, startnodes, weights):
     n = startnodes.size
     for k in range(n):
         startnode = startnodes.flat[k]
         endnode = fdir.flat[startnode]
-        _d8_reverse_distance_recursion(startnode, endnode, min_order, max_order,
+        _d8_reverse_distance_recursion(startnode, endnode, max_order,
                                        rdist, fdir, indegree, weights)
     return rdist
 
-@njit(float64[:,:](int64[:,:], int64[:,:], float64[:,:], int64[:,:],
+@njit(float64[:,:](int64[:,:], float64[:,:], int64[:,:],
                    uint8[:], int64[:], float64[:,:]),
       cache=True)
-def _d8_reverse_distance_iter_numba(min_order, max_order, rdist, fdir,
+def _d8_reverse_distance_iter_numba(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]
@@ -994,6 +992,36 @@ def _d8_reverse_distance_iter_numba(min_order, max_order, rdist, fdir,
             endnode = fdir.flat[startnode]
     return rdist
 
+# TODO: This should probably have two weights vectors
+@njit(float64[:,:](float64[:,:], int64[:,:], int64[:,:],
+                   uint8[:], int64[:], float64[:,:]),
+      cache=True)
+def _dinf_reverse_distance_iter_numba(rdist, fdir_0, fdir_1,
+                                      indegree, startnodes, weights):
+    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]
+            rdist.flat[endnode_0] = max(rdist.flat[endnode_0],
+                                        rdist.flat[startnode] + weights.flat[endnode_0])
+            rdist.flat[endnode_1] = max(rdist.flat[endnode_1],
+                                        rdist.flat[startnode] + weights.flat[endnode_1])
+            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 rdist
+
 # Functions for 'resolve_flats'
 
 @njit(UniTuple(boolean[:,:], 3)(float64[:,:], int64[:]),
diff --git a/pysheds/sgrid.py b/pysheds/sgrid.py
index 651ce97..4241859 100644
--- a/pysheds/sgrid.py
+++ b/pysheds/sgrid.py
@@ -1041,6 +1041,7 @@ def _d8_flow_distance(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4,
         fdir[invalid_cells] = 0
         if xytype in {'label', 'coordinate'}:
             x, y = self.nearest_cell(x, y, fdir.affine, snap)
+        # TODO: Should this be ones for all cells?
         if weights is None:
             weights = (~nodata_cells).reshape(fdir.shape).astype(np.float64)
         if algorithm.lower() == 'iterative':
@@ -1066,6 +1067,7 @@ def _dinf_flow_distance(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4
             weights_0 = weights
             weights_1 = weights
         else:
+            # TODO: Should this be ones for all cells?
             weights_0 = (~nodata_cells).reshape(fdir.shape).astype(np.float64)
             weights_1 = weights_0
         if method.lower() == 'shortest':
@@ -1347,8 +1349,9 @@ def stream_order(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
                                      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', algorithm='iterative', **kwargs):
+    def distance_to_ridge(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
+                          nodata_out=np.nan, routing='d8', algorithm='iterative', cycle_size=1,
+                          **kwargs):
         """
         Generates a raster representing the (weighted) topological distance from each cell
         to its originating drainage divide, moving upstream.
@@ -1383,43 +1386,89 @@ def distance_to_ridge(self, fdir, mask, weights=None, dirmap=(64, 128, 1, 2, 4,
         """
         if routing.lower() == 'd8':
             fdir_overrides = {'dtype' : np.int64, 'nodata' : fdir.nodata}
+        elif routing.lower() == 'dinf':
+            fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
         else:
-            raise NotImplementedError('Only implemented for D8 routing.')
-        mask_overrides = {'dtype' : np.bool8, 'nodata' : False}
+            raise NotImplementedError('Routing method must be one of: `d8`, `dinf`')
         kwargs.update(fdir_overrides)
         fdir = self._input_handler(fdir, **kwargs)
-        kwargs.update(mask_overrides)
-        mask = self._input_handler(mask, **kwargs)
         if weights is not None:
             weights_overrides = {'dtype' : np.float64, 'nodata' : weights.nodata}
             kwargs.update(weights_overrides)
             weights = self._input_handler(weights, **kwargs)
+        if routing.lower() == 'd8':
+            rdist = self._d8_distance_to_ridge(fdir=fdir, weights=weights,
+                                               dirmap=dirmap, algorithm=algorithm,
+                                               nodata_out=nodata_out)
+        elif routing.lower() == 'dinf':
+            rdist = self._dinf_distance_to_ridge(fdir=fdir, weights=weights,
+                                                 dirmap=dirmap, algorithm=algorithm,
+                                                 nodata_out=nodata_out, cycle_size=cycle_size)
+        return rdist
+
+    def _d8_distance_to_ridge(self, fdir, weights, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
+                              algorithm='iterative', nodata_out=np.nan, **kwargs):
         # Find nodata cells and invalid cells
         nodata_cells = self._get_nodata_cells(fdir)
         invalid_cells = ~np.in1d(fdir.ravel(), dirmap).reshape(fdir.shape)
         # Set nodata cells to zero
         fdir[nodata_cells] = 0
         fdir[invalid_cells] = 0
+        # TODO: Should this be ones for all cells?
         if weights is None:
             weights = (~nodata_cells).reshape(fdir.shape).astype(np.float64)
-        maskleft, maskright, masktop, maskbottom = self._pop_rim(mask, nodata=0)
-        masked_fdir = np.where(mask, fdir, 0).astype(np.int64)
         startnodes = np.arange(fdir.size, dtype=np.int64)
-        endnodes = _self._flatten_fdir_numba(masked_fdir, dirmap).reshape(fdir.shape)
+        endnodes = _self._flatten_fdir_numba(fdir, dirmap).reshape(fdir.shape)
         indegree = np.bincount(endnodes.ravel()).astype(np.uint8)
-        orig_indegree = np.copy(indegree)
         startnodes = startnodes[(indegree == 0)]
-        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)
         if algorithm.lower() == 'iterative':
-            rdist = _self._d8_reverse_distance_iter_numba(min_order, max_order, rdist,
-                                                        endnodes, indegree, startnodes,
-                                                        weights)
+            rdist = _self._d8_reverse_distance_iter_numba(max_order, rdist,
+                                                          endnodes, indegree,
+                                                          startnodes, weights)
+        elif algorithm.lower() == 'recursive':
+            rdist = _self._d8_reverse_distance_recur_numba(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
+
+    def _dinf_distance_to_ridge(self, fdir, weights, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
+                                algorithm='iterative', nodata_out=np.nan, cycle_size=1,
+                                **kwargs):
+        # Find nodata cells and invalid cells
+        nodata_cells = self._get_nodata_cells(fdir)
+        # Split d-infinity grid
+        fdir_0, fdir_1, prop_0, prop_1 = _self._angle_to_d8_numba(fdir, dirmap, nodata_cells)
+        # Get matching of start and end nodes
+        startnodes = np.arange(fdir.size, dtype=np.int64)
+        endnodes_0 = _self._flatten_fdir_numba(fdir_0, dirmap).reshape(fdir.shape)
+        endnodes_1 = _self._flatten_fdir_numba(fdir_1, dirmap).reshape(fdir.shape)
+        # Remove cycles
+        _self._dinf_fix_cycles_numba(endnodes_0, endnodes_1, cycle_size)
+        # Find indegree of all cells
+        indegree_0 = np.bincount(endnodes_0.ravel(), minlength=fdir.size)
+        indegree_1 = np.bincount(endnodes_1.ravel(), minlength=fdir.size)
+        indegree = (indegree_0 + indegree_1).astype(np.uint8)
+        # Set starting nodes to those with no predecessors
+        startnodes = startnodes[(indegree == 0)]
+        # TODO: Should this be ones for all cells?
+        if weights is None:
+            weights = np.ones(fdir_0.shape, dtype=np.float64)
+        rdist = np.zeros(fdir_0.shape, dtype=np.float64)
+        if algorithm.lower() == 'iterative':
+            rdist = _self._dinf_reverse_distance_iter_numba(rdist,
+                                                            endnodes_0,
+                                                            endnodes_1,
+                                                            indegree,
+                                                            startnodes,
+                                                            weights)
         elif algorithm.lower() == 'recursive':
-            rdist = _self._d8_reverse_distance_recur_numba(min_order, max_order, rdist,
-                                                           endnodes, indegree, startnodes,
-                                                           weights)
+            raise NotImplementedError('Recursive algorithm not implemented for distance_to_ridge.')
         else:
             raise ValueError('Algorithm must be `iterative` or `recursive`.')
         rdist = self._output_handler(data=rdist, viewfinder=fdir.viewfinder,

From 5df327722dd726d1130593a368b09f47d7b34fb4 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Sun, 16 Jan 2022 15:36:23 -0600
Subject: [PATCH 08/27] Simplify D8 distance_to_ridge

---
 pysheds/_sgrid.py | 32 +++++++++++++++-----------------
 pysheds/sgrid.py  | 13 ++++++-------
 2 files changed, 21 insertions(+), 24 deletions(-)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index 85c55af..012ec86 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -948,46 +948,44 @@ def _mfd_flow_distance_heap_numba(fdir, pour_point, weights):
 
 # Functions for 'reverse_flow_distance'
 
-@njit(void(int64, int64, int64[:,:], float64[:,:],
-           int64[:,:], uint8[:], float64[:,:]),
+@njit(void(int64, int64, float64[:,:], int64[:,:], uint8[:], float64[:,:]),
       cache=True)
-def _d8_reverse_distance_recursion(startnode, endnode, max_order,
-                                   rdist, fdir, indegree, weights):
-    max_order.flat[endnode] = max(max_order.flat[endnode], rdist.flat[startnode])
+def _d8_reverse_distance_recursion(startnode, endnode, rdist, fdir, indegree,
+                                   weights):
+    rdist.flat[endnode] = max(rdist.flat[endnode],
+                              rdist.flat[startnode] + weights.flat[endnode])
     indegree.flat[endnode] -= 1
     if indegree.flat[endnode] == 0:
-        rdist.flat[endnode] = max_order.flat[endnode] + weights.flat[endnode]
         new_startnode = endnode
         new_endnode = fdir.flat[new_startnode]
-        _d8_reverse_distance_recursion(new_startnode, new_endnode, max_order,
-                                       rdist, fdir, indegree, weights)
+        _d8_reverse_distance_recursion(new_startnode, new_endnode, rdist, fdir,
+                                       indegree, weights)
 
-@njit(float64[:,:](int64[:,:], float64[:,:], int64[:,:],
+@njit(float64[:,:](float64[:,:], int64[:,:],
                    uint8[:], int64[:], float64[:,:]),
       cache=True)
-def _d8_reverse_distance_recur_numba(max_order, rdist, fdir,
+def _d8_reverse_distance_recur_numba(rdist, fdir,
                                      indegree, startnodes, weights):
     n = startnodes.size
     for k in range(n):
         startnode = startnodes.flat[k]
         endnode = fdir.flat[startnode]
-        _d8_reverse_distance_recursion(startnode, endnode, max_order,
-                                       rdist, fdir, indegree, weights)
+        _d8_reverse_distance_recursion(startnode, endnode, rdist, fdir,
+                                       indegree, weights)
     return rdist
 
-@njit(float64[:,:](int64[:,:], float64[:,:], int64[:,:],
+@njit(float64[:,:](float64[:,:], int64[:,:],
                    uint8[:], int64[:], float64[:,:]),
       cache=True)
-def _d8_reverse_distance_iter_numba(max_order, rdist, fdir,
-                                    indegree, startnodes, weights):
+def _d8_reverse_distance_iter_numba(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):
-            max_order.flat[endnode] = max(max_order.flat[endnode], rdist.flat[startnode])
+            rdist.flat[endnode] = max(rdist.flat[endnode],
+                                      rdist.flat[startnode] + weights.flat[endnode])
             indegree.flat[endnode] -= 1
-            rdist.flat[endnode] = max_order.flat[endnode] + weights.flat[endnode]
             startnode = endnode
             endnode = fdir.flat[startnode]
     return rdist
diff --git a/pysheds/sgrid.py b/pysheds/sgrid.py
index 4241859..71edc48 100644
--- a/pysheds/sgrid.py
+++ b/pysheds/sgrid.py
@@ -1421,16 +1421,15 @@ def _d8_distance_to_ridge(self, fdir, weights, dirmap=(64, 128, 1, 2, 4, 8, 16,
         endnodes = _self._flatten_fdir_numba(fdir, dirmap).reshape(fdir.shape)
         indegree = np.bincount(endnodes.ravel()).astype(np.uint8)
         startnodes = startnodes[(indegree == 0)]
-        max_order = np.ones(fdir.shape, dtype=np.int64)
         rdist = np.zeros(fdir.shape, dtype=np.float64)
         if algorithm.lower() == 'iterative':
-            rdist = _self._d8_reverse_distance_iter_numba(max_order, rdist,
-                                                          endnodes, indegree,
-                                                          startnodes, weights)
+            rdist = _self._d8_reverse_distance_iter_numba(rdist, endnodes,
+                                                          indegree, startnodes,
+                                                          weights)
         elif algorithm.lower() == 'recursive':
-            rdist = _self._d8_reverse_distance_recur_numba(max_order, rdist,
-                                                           endnodes, indegree,
-                                                           startnodes, weights)
+            rdist = _self._d8_reverse_distance_recur_numba(rdist, endnodes,
+                                                           indegree, startnodes,
+                                                           weights)
         else:
             raise ValueError('Algorithm must be `iterative` or `recursive`.')
         rdist = self._output_handler(data=rdist, viewfinder=fdir.viewfinder,

From ce5c1703db599b6e3bb57a0886da2fce9d0a62a9 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Mon, 17 Jan 2022 20:07:07 -0600
Subject: [PATCH 09/27] Add mfd distance_to_ridge

---
 pysheds/_sgrid.py | 30 ++++++++++++++++++++++++++++--
 1 file changed, 28 insertions(+), 2 deletions(-)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index 012ec86..3b108e9 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -643,7 +643,7 @@ def _dinf_accumulation_eff_iter_numba(acc, fdir_0, fdir_1, indegree, startnodes,
                     queue.append(endnode_1)
     return acc
 
-@njit(parallel=True)
+@njit(uint8[:](int64[:,:,:]), parallel=True)
 def _mfd_bincount(fdir):
     p, m, n = fdir.shape
     mn = m * n
@@ -656,7 +656,8 @@ def _mfd_bincount(fdir):
                 out[endnode] += 1
     return out
 
-@njit
+@njit(float64[:,:](float64[:,:], int64[:,:,:], float64[:,:,:], uint8[:], int64[:]),
+      cache=True)
 def _mfd_accumulation_iter_numba(acc, fdir, props, indegree, startnodes):
     n = startnodes.size
     queue = [0]
@@ -1020,6 +1021,31 @@ def _dinf_reverse_distance_iter_numba(rdist, fdir_0, fdir_1,
                     queue.append(endnode_1)
     return rdist
 
+@njit(float64[:,:](float64[:,:], int64[:,:,:], uint8[:], int64[:], float64[:,:]),
+      cache=True)
+def _mfd_reverse_distance_iter_numba(rdist, fdir, indegree, startnodes, weights):
+    n = startnodes.size
+    queue = [0]
+    _ = queue.pop()
+    for k in range(n):
+        startnode = startnodes.flat[k]
+        queue.append(startnode)
+        while queue:
+            startnode = queue.pop()
+            for i in range(8):
+                fdir_i = fdir[i]
+                endnode = fdir_i.flat[startnode]
+                if endnode == startnode:
+                    continue
+                else:
+                    weight = weights.flat[startnode]
+                    rdist.flat[endnode] = max(rdist.flat[endnode],
+                                              rdist.flat[startnode] + weights.flat[endnode])
+                    indegree.flat[endnode] -= 1
+                    if (indegree.flat[endnode] == 0):
+                        queue.append(endnode)
+    return rdist
+
 # Functions for 'resolve_flats'
 
 @njit(UniTuple(boolean[:,:], 3)(float64[:,:], int64[:]),

From 442c20a1ea067a4019c1b9dc7526b7e4337ef566 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Mon, 17 Jan 2022 21:26:53 -0600
Subject: [PATCH 10/27] Add new flatten_fdir for mfd

---
 pysheds/_sgrid.py | 51 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 51 insertions(+)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index 3b108e9..7f4f8ce 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -1709,6 +1709,57 @@ def _flatten_fdir_no_boundary(fdir, dirmap):
         flat_fdir.flat[k] = k + offset
     return flat_fdir
 
+# TODO: Assumes pits and flats are removed
+@njit(int64[:,:,:](float64[:,:,:]),
+      parallel=True,
+      cache=True)
+def _flatten_mfd_fdir_numba(fdir):
+    p, r, c = fdir.shape
+    n = r * c
+    flat_fdir = np.zeros((p, r, c), dtype=np.int64)
+    offsets = np.array([0 - c, 1 - c, 1 + 0, 1 + c,
+                        0 + c, -1 + c, -1 + 0, -1 - c],
+                       dtype=np.int64)
+    left_map = np.array([offsets[0], offsets[1], offsets[2], offsets[3],
+                         offsets[4], 0, 0, 0],
+                        dtype=np.int64)
+    right_map = np.array([offsets[0], 0, 0, 0,
+                          offsets[4], offsets[5], offsets[6], offsets[7]],
+                         dtype=np.int64)
+    top_map = np.array([0, 0, offsets[2], offsets[3],
+                        offsets[4], offsets[5], offsets[6], 0],
+                       dtype=np.int64)
+    bottom_map = np.array([offsets[0], offsets[1], offsets[2], 0,
+                           0, 0, offsets[6], offsets[7]],
+                          dtype=np.int64)
+    for i in prange(8):
+        for k in prange(n):
+            kix = k + (i * n)
+            cell_value = fdir.flat[kix]
+            if cell_value == 0:
+                offset = 0
+            else:
+                on_left = ((k % c) == 0)
+                on_right = (((k + 1) % c) == 0)
+                on_top = (k < c)
+                on_bottom = (k > (n - c - 1))
+                on_boundary = (on_left | on_right | on_top | on_bottom)
+                if on_boundary:
+                    # TODO: This seems like it could cause errors at corner points
+                    # TODO: Check if offset is already zero
+                    if on_left:
+                        offset = left_map[i]
+                    if on_right and (offset != 0):
+                        offset = right_map[i]
+                    if on_top and (offset != 0):
+                        offset = top_map[i]
+                    if on_bottom and (offset != 0):
+                        offset = bottom_map[i]
+                else:
+                    offset = offsets[i]
+            flat_fdir.flat[kix] = k + offset
+    return flat_fdir
+
 @njit
 def _construct_matching(fdir, dirmap):
     n = fdir.size

From aafa193a6060401fe066cea95f95daebd65082a6 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Tue, 18 Jan 2022 19:57:22 -0600
Subject: [PATCH 11/27] Add view support for 3d raster

---
 pysheds/sview.py | 78 +++++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 68 insertions(+), 10 deletions(-)

diff --git a/pysheds/sview.py b/pysheds/sview.py
index 0822a6a..e59d5a8 100644
--- a/pysheds/sview.py
+++ b/pysheds/sview.py
@@ -243,6 +243,16 @@ def to_crs(self, new_crs, **kwargs):
                             data_view=self.viewfinder, **kwargs)
         return new_raster
 
+class MultiRaster(Raster):
+    def __new__(cls, input_array, viewfinder=None, metadata={}):
+        return super().__new__(cls, input_array,
+                               viewfinder=viewfinder,
+                               metadata=metadata)
+
+    # TODO: Anything relying on shape[0] and shape[1] will be wrong
+    def to_crs(self, new_crs, **kwargs):
+        raise NotImplementedError('Not Implemented for MultiRaster.')
+
 class ViewFinder():
     """
     Class that defines a spatial reference system for a Raster or Grid instance.
@@ -286,6 +296,7 @@ def __eq__(self, other):
         if isinstance(other, ViewFinder):
             is_eq = True
             is_eq &= (self.affine == other.affine)
+            # TODO: This shape will not work for multiraster
             is_eq &= (self.shape[0] == other.shape[0])
             is_eq &= (self.shape[1] == other.shape[1])
             is_eq &= (self.crs == other.crs)
@@ -406,6 +417,7 @@ def is_congruent_with(self, other):
         if isinstance(other, ViewFinder):
             is_congruent = True
             is_congruent &= (self.affine == other.affine)
+            # TODO: This won't work for ndim > 2
             is_congruent &= (self.shape[0] == other.shape[0])
             is_congruent &= (self.shape[1] == other.shape[1])
             is_congruent &= (self.crs == other.crs)
@@ -456,7 +468,7 @@ def view(cls, data, target_view, data_view=None, interpolation='nearest',
         target_view : ViewFinder
                       The desired spatial reference system.
         data_view : ViewFinder
-                    The spatial reference system of the data. Defaults to the Raster dataset's
+                    The spatial reference system of the data. Defaults to the Raster dataset
                     `viewfinder` attribute.
         interpolation : 'nearest', 'linear'
                         Interpolation method to be used if spatial reference systems
@@ -490,12 +502,16 @@ def view(cls, data, target_view, data_view=None, interpolation='nearest',
         out : Raster
               View of the input Raster at the provided target view.
         """
-        # If no data view given, use data's view
+        # If no data view given, use data view
+        is_raster = isinstance(data, Raster)
+        is_multiraster = isinstance(data, MultiRaster)
+        is_2d_array = isinstance(data, np.ndarray) and (data.ndim == 2)
+        is_3d_array = isinstance(data, np.ndarray) and (data.ndim == 3)
         if data_view is None:
             try:
-                assert(isinstance(data, Raster))
+                assert (is_raster or is_multiraster)
             except:
-                raise TypeError('`data` must be a Raster instance.')
+                raise TypeError('`data` must be a Raster or MultiRaster instance.')
             data_view = data.viewfinder
         # By default, use dataset's `nodata` value
         if nodata is None:
@@ -513,9 +529,32 @@ def view(cls, data, target_view, data_view=None, interpolation='nearest',
                                     dtype=dtype,
                                     interpolation=interpolation)
         # Mask input data if desired
-        if apply_input_mask:
+        has_input_mask = not data_view.mask.all()
+        if apply_input_mask and has_input_mask:
+            # TODO: Is this necessary?
             arr = np.where(data_view.mask, data, target_view.nodata).astype(dtype)
             data = Raster(arr, data.viewfinder, metadata=data.metadata)
+        if is_multiraster or is_3d_array:
+            out = cls._view_multiraster(data, target_view, data_view=data_view,
+                                        interpolation=interpolation,
+                                        apply_output_mask=apply_output_mask,
+                                        dtype=dtype)
+        elif is_raster or is_2d_array:
+            out = cls._view_raster(data, target_view, data_view=data_view,
+                                   interpolation=interpolation,
+                                   apply_output_mask=apply_output_mask,
+                                   dtype=dtype)
+        else:
+            raise TypeError('`data` must be a Raster, MultiRaster, 2D array or 3D array.')
+        # Write metadata
+        if inherit_metadata:
+            out.metadata.update(data.metadata)
+        out.metadata.update(new_metadata)
+        return out
+
+    @classmethod
+    def _view_raster(cls, data, target_view, data_view=None, interpolation='nearest',
+                     apply_output_mask=True, dtype=None):
         # If data view and target view are the same, return a copy of the data
         if data_view.is_congruent_with(target_view):
             out = cls._view_same_viewfinder(data, data_view, target_view, dtype,
@@ -525,10 +564,27 @@ def view(cls, data, target_view, data_view=None, interpolation='nearest',
             out = cls._view_different_viewfinder(data, data_view, target_view, dtype,
                                                  apply_output_mask=apply_output_mask,
                                                  interpolation=interpolation)
-        # Write metadata
-        if inherit_metadata:
-            out.metadata.update(data.metadata)
-        out.metadata.update(new_metadata)
+        return out
+
+    @classmethod
+    def _view_multiraster(cls, data, target_view, data_view=None, interpolation='nearest',
+                          apply_output_mask=True, dtype=None):
+        k, m, n = data.shape
+        out = np.full((k, *target_view.shape), target_view.nodata, dtype=dtype)
+        out_mask = np.ones((k, *target_view.shape), dtype=np.bool8)
+        for i in range(k):
+            slice_viewfinder = ViewFinder(affine=data_view.affine, mask=data_view.mask[i],
+                                          nodata=data_view.nodata, crs=data_view.crs)
+            data_i = Raster(np.asarray(data[i]), viewfinder=slice_viewfinder)
+            view_i = cls._view_raster(data_i, target_view, slice_viewfinder,
+                                      interpolation=interpolation,
+                                      apply_output_mask=apply_output_mask,
+                                      dtype=dtype)
+            out[i, :, :] = view_i
+            out_mask[i, :, :] = view_i.mask
+        target_viewfinder = ViewFinder(affine=target_view.affine, mask=out_mask,
+                                       nodata=target_view.nodata, crs=target_view.crs)
+        out = MultiRaster(out, viewfinder=target_viewfinder)
         return out
 
     @classmethod
@@ -807,7 +863,7 @@ def _override_dtype(cls, data, target_view, dtype=None, interpolation='nearest')
     @classmethod
     def _view_same_viewfinder(cls, data, data_view, target_view, dtype,
                               apply_output_mask=True):
-        if apply_output_mask:
+        if (apply_output_mask) and (not target_view.mask.all()):
             out = np.where(target_view.mask, data, target_view.nodata).astype(dtype)
         else:
             out = np.asarray(data.copy(), dtype=dtype)
@@ -824,6 +880,7 @@ def _view_different_viewfinder(cls, data, data_view, target_view, dtype,
         else:
             out = cls._view_different_crs(out, data, data_view,
                                           target_view, interpolation)
+        # TODO: Skip if mask is True everywhere
         if apply_output_mask:
             np.place(out, ~target_view.mask, target_view.nodata)
         out = Raster(out, target_view)
@@ -839,6 +896,7 @@ def _view_same_crs(cls, view, data, data_view, target_view, interpolation='neare
         x_ix, _ = cls.affine_transform(inv_affine, x,
                                        np.zeros(target_view.shape[1],
                                                 dtype=np.float64))
+        # TODO: Does this work for rotated data?
         if interpolation == 'nearest':
             view = _self._view_fill_by_axes_nearest_numba(data, view, y_ix, x_ix)
         elif interpolation == 'linear':

From 29dc19d09a5de158f971cf7e4e5cecdebef69025 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Tue, 18 Jan 2022 20:35:29 -0600
Subject: [PATCH 12/27] Increase support for multiraster and ensure correct
 dimensions

---
 pysheds/sview.py | 88 ++++++++++++++++++++++++++++++++++++------------
 1 file changed, 67 insertions(+), 21 deletions(-)

diff --git a/pysheds/sview.py b/pysheds/sview.py
index e59d5a8..ee89cd2 100644
--- a/pysheds/sview.py
+++ b/pysheds/sview.py
@@ -49,6 +49,11 @@ class Raster(np.ndarray):
     """
 
     def __new__(cls, input_array, viewfinder=None, metadata={}):
+        try:
+            assert isinstance(input_array, np.ndarray)
+            assert input_array.ndim == 2
+        except:
+            raise TypeError('Input must be a 2-dimensional array-like object.')
         # Handle case where input is a Raster itself
         if isinstance(input_array, Raster):
             input_array, viewfinder, metadata = cls._handle_raster_input(input_array,
@@ -218,7 +223,7 @@ def to_crs(self, new_crs, **kwargs):
         old_crs = self.crs
         dx = self.affine.a
         dy = self.affine.e
-        m, n = self.shape
+        m, n = self.shape[-2], self.shape[-1]
         Y, X = np.mgrid[0:m, 0:n]
         top = np.column_stack([X[0, :], Y[0, :]])
         bottom = np.column_stack([X[-1, :], Y[-1, :]])
@@ -228,7 +233,7 @@ def to_crs(self, new_crs, **kwargs):
         xi, yi = boundary[:,0], boundary[:,1]
         xb, yb = View.affine_transform(self.affine, xi, yi)
         xb_p, yb_p = pyproj.transform(old_crs, new_crs, xb, yb,
-                                    errcheck=True, always_xy=True)
+                                      errcheck=True, always_xy=True)
         x0_p = xb_p.min() if (dx > 0) else xb_p.max()
         y0_p = yb_p.min() if (dy > 0) else yb_p.max()
         xn_p = xb_p.max() if (dx > 0) else xb_p.min()
@@ -240,18 +245,56 @@ def to_crs(self, new_crs, **kwargs):
                                     nodata=self.nodata, mask=self.mask,
                                     crs=new_crs)
         new_raster = View.view(self, target_view=new_viewfinder,
-                            data_view=self.viewfinder, **kwargs)
+                               data_view=self.viewfinder, **kwargs)
         return new_raster
 
 class MultiRaster(Raster):
     def __new__(cls, input_array, viewfinder=None, metadata={}):
-        return super().__new__(cls, input_array,
-                               viewfinder=viewfinder,
-                               metadata=metadata)
-
-    # TODO: Anything relying on shape[0] and shape[1] will be wrong
-    def to_crs(self, new_crs, **kwargs):
-        raise NotImplementedError('Not Implemented for MultiRaster.')
+        try:
+            assert isinstance(input_array, np.ndarray)
+            if (input_array.ndim == 2):
+                input_array = input_array.reshape(1, *input_array.shape)
+            assert input_array.ndim == 3
+        except:
+            raise TypeError('Input must be a 2d or 3d array-like object.')
+        # Handle case where input is a Raster itself
+        if isinstance(input_array, Raster):
+            input_array, viewfinder, metadata = cls._handle_raster_input(input_array,
+                                                                         viewfinder,
+                                                                         metadata)
+        # Create a numpy array from the input
+        obj = np.asarray(input_array).view(cls)
+        # If no viewfinder provided, construct one congruent with the array shape
+        if viewfinder is None:
+            viewfinder = ViewFinder(shape=obj.shape)
+        # If a viewfinder is provided, ensure that it is a viewfinder...
+        else:
+            try:
+                assert(isinstance(viewfinder, ViewFinder))
+            except:
+                raise ValueError("Must initialize with a ViewFinder.")
+            # Ensure that viewfinder shape is correct...
+            try:
+                assert viewfinder.shape == obj.shape
+            except:
+                raise ValueError('Viewfinder and array shape must be the same.')
+        # Test typing of array
+        try:
+            assert not np.issubdtype(obj.dtype, np.object_)
+            assert not np.issubdtype(obj.dtype, np.flexible)
+        except:
+            raise TypeError('`object` and `flexible` dtypes not allowed.')
+        try:
+            assert np.min_scalar_type(viewfinder.nodata) <= obj.dtype
+        except:
+            raise TypeError('`nodata` value not representable in dtype of array.')
+        # Don't allow original viewfinder and metadata to be modified
+        viewfinder = viewfinder.copy()
+        metadata = metadata.copy()
+        # Set attributes of array
+        obj._viewfinder = viewfinder
+        obj.metadata = metadata
+        return obj
 
 class ViewFinder():
     """
@@ -296,9 +339,7 @@ def __eq__(self, other):
         if isinstance(other, ViewFinder):
             is_eq = True
             is_eq &= (self.affine == other.affine)
-            # TODO: This shape will not work for multiraster
-            is_eq &= (self.shape[0] == other.shape[0])
-            is_eq &= (self.shape[1] == other.shape[1])
+            is_eq &= (self.shape == other.shape)
             is_eq &= (self.crs == other.crs)
             is_eq &= (self.mask == other.mask).all()
             if np.isnan(self.nodata):
@@ -379,7 +420,7 @@ def size(self):
     def bbox(self):
         shape = self.shape
         xmin, ymax = View.affine_transform(self.affine, 0, 0)
-        xmax, ymin = View.affine_transform(self.affine, shape[1], shape[0])
+        xmax, ymin = View.affine_transform(self.affine, shape[-1], shape[-2])
         _bbox = (xmin, ymin, xmax, ymax)
         return _bbox
 
@@ -411,15 +452,15 @@ def properties(self):
 
     @property
     def axes(self):
-        return View.axes(self.affine, self.shape)
+        shape = self.shape
+        return View.axes(self.affine, (shape[-2], shape[-1]))
 
     def is_congruent_with(self, other):
         if isinstance(other, ViewFinder):
             is_congruent = True
             is_congruent &= (self.affine == other.affine)
-            # TODO: This won't work for ndim > 2
-            is_congruent &= (self.shape[0] == other.shape[0])
-            is_congruent &= (self.shape[1] == other.shape[1])
+            is_congruent &= (self.shape[-2] == other.shape[-2])
+            is_congruent &= (self.shape[-1] == other.shape[-1])
             is_congruent &= (self.crs == other.crs)
             return is_congruent
         else:
@@ -790,6 +831,10 @@ def clip_to_mask(cls, data, mask=None, pad=(0,0,0,0)):
               A Raster dataset clipped to the bounding box of the non-null entries
               in the given mask.
         """
+        try:
+            assert (data.ndim == 2)
+        except:
+            raise ValueError('Data must be 2-dimensional')
         try:
             for value in pad:
                 assert (isinstance(value, int))
@@ -863,7 +908,8 @@ def _override_dtype(cls, data, target_view, dtype=None, interpolation='nearest')
     @classmethod
     def _view_same_viewfinder(cls, data, data_view, target_view, dtype,
                               apply_output_mask=True):
-        if (apply_output_mask) and (not target_view.mask.all()):
+        has_output_mask = not target_view.mask.all()
+        if (apply_output_mask) and (has_output_mask):
             out = np.where(target_view.mask, data, target_view.nodata).astype(dtype)
         else:
             out = np.asarray(data.copy(), dtype=dtype)
@@ -880,8 +926,8 @@ def _view_different_viewfinder(cls, data, data_view, target_view, dtype,
         else:
             out = cls._view_different_crs(out, data, data_view,
                                           target_view, interpolation)
-        # TODO: Skip if mask is True everywhere
-        if apply_output_mask:
+        has_output_mask = not target_view.mask.all()
+        if (apply_output_mask) and (has_output_mask):
             np.place(out, ~target_view.mask, target_view.nodata)
         out = Raster(out, target_view)
         return out

From 7f7a2e294339eedc2dc43aca197ff108c74bb1fb Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Tue, 18 Jan 2022 23:01:22 -0600
Subject: [PATCH 13/27] Add more checks on multiraster

---
 pysheds/sview.py | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/pysheds/sview.py b/pysheds/sview.py
index ee89cd2..d48f8ef 100644
--- a/pysheds/sview.py
+++ b/pysheds/sview.py
@@ -50,7 +50,9 @@ class Raster(np.ndarray):
 
     def __new__(cls, input_array, viewfinder=None, metadata={}):
         try:
+            # MultiRaster must be subclass of ndarray
             assert isinstance(input_array, np.ndarray)
+            # Ensure MultiRaster is 2D
             assert input_array.ndim == 2
         except:
             raise TypeError('Input must be a 2-dimensional array-like object.')
@@ -251,9 +253,12 @@ def to_crs(self, new_crs, **kwargs):
 class MultiRaster(Raster):
     def __new__(cls, input_array, viewfinder=None, metadata={}):
         try:
+            # MultiRaster must be subclass of ndarray
             assert isinstance(input_array, np.ndarray)
+            # If 2D, upcast to 3D
             if (input_array.ndim == 2):
                 input_array = input_array.reshape(1, *input_array.shape)
+            # Ensure MultiRaster is 3D
             assert input_array.ndim == 3
         except:
             raise TypeError('Input must be a 2d or 3d array-like object.')
@@ -475,6 +480,8 @@ def view(self, raster, **kwargs):
         target_view = self
         return View.view(raster, data_view, target_view, **kwargs)
 
+# TODO: Ensure that target_view cannot be 3-dimensional
+# TODO: Or use only last two entries of viewfinder.shape
 class View():
     """
     Class containing methods for manipulating views of gridded datasets.
@@ -543,11 +550,17 @@ def view(cls, data, target_view, data_view=None, interpolation='nearest',
         out : Raster
               View of the input Raster at the provided target view.
         """
-        # If no data view given, use data view
+        # TODO: Temporary check on target_view
+        try:
+            assert (len(target_view.shape) == 2)
+        except:
+            raise ValueError('Target view must be 2-dimensional.')
+        # TODO: Use type instead of isinstance
         is_raster = isinstance(data, Raster)
         is_multiraster = isinstance(data, MultiRaster)
         is_2d_array = isinstance(data, np.ndarray) and (data.ndim == 2)
         is_3d_array = isinstance(data, np.ndarray) and (data.ndim == 3)
+        # If no data view given, use data view
         if data_view is None:
             try:
                 assert (is_raster or is_multiraster)

From c65f8d1dbbb91aee9a8b006cb75aed54bbf0bc4d Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Wed, 19 Jan 2022 00:27:13 -0600
Subject: [PATCH 14/27] Add mfd flowdir, catchment, and accumulation to grid
 class

---
 pysheds/_sgrid.py |   2 +-
 pysheds/sgrid.py  | 130 ++++++++++++++++++++++++++++++++++++++++++----
 2 files changed, 120 insertions(+), 12 deletions(-)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index 7f4f8ce..a101b5d 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -253,7 +253,7 @@ def _angle_to_d8_numba(angles, dirmap, nodata_cells):
         props_1.flat[i] = prop_1
     return fdirs_0, fdirs_1, props_0, props_1
 
-@njit(float64[:,:,:](float64[:,:], float64, float64, boolean[:,:], int64, int64),
+@njit(float64[:,:,:](float64[:,:], float64, float64, boolean[:,:], float64, int64),
       parallel=True,
       cache=True)
 def _mfd_flowdir_numba(dem, dx, dy, nodata_cells, nodata_out, p=1):
diff --git a/pysheds/sgrid.py b/pysheds/sgrid.py
index 71edc48..caab2b2 100644
--- a/pysheds/sgrid.py
+++ b/pysheds/sgrid.py
@@ -29,7 +29,7 @@
 import pysheds.io
 
 # Import viewing functions
-from pysheds.sview import Raster
+from pysheds.sview import Raster, MultiRaster
 from pysheds.sview import View, ViewFinder
 
 # Import numba functions
@@ -611,6 +611,12 @@ def flowdir(self, dem, routing='d8', flats=-1, pits=-2, nodata_out=None,
             fdir = self._dinf_flowdir(dem=dem, nodata_cells=nodata_cells,
                                       nodata_out=nodata_out, flats=flats,
                                       pits=pits, dirmap=dirmap)
+        elif routing.lower() == 'mfd':
+            if nodata_out is None:
+                nodata_out = np.nan
+            fdir = self._mfd_flowdir(dem=dem, nodata_cells=nodata_cells,
+                                     nodata_out=nodata_out, flats=flats,
+                                     pits=pits, dirmap=dirmap)
         else:
             raise ValueError('Routing method must be one of: `d8`, `dinf`')
         fdir.metadata.update(default_metadata)
@@ -640,6 +646,21 @@ def _dinf_flowdir(self, dem, nodata_cells, nodata_out=np.nan, flats=-1, pits=-2,
         return self._output_handler(data=fdir, viewfinder=dem.viewfinder,
                                     metadata=dem.metadata, nodata=nodata_out)
 
+    def _mfd_flowdir(self, dem, nodata_cells, nodata_out=np.nan, flats=-1, pits=-2,
+                     dirmap=(64, 128, 1, 2, 4, 8, 16, 32)):
+        # Make sure nothing flows to the nodata cells
+        dem[nodata_cells] = dem.max() + 1
+        dx = abs(dem.affine.a)
+        dy = abs(dem.affine.e)
+        # TODO: Allow p value to be changed
+        fdir = _self._mfd_flowdir_numba(dem, dx, dy, nodata_cells,
+                                        nodata_out, p=1)
+        # Create new mask to match new shape
+        new_mask = np.tile(dem.viewfinder.mask, (8, 1, 1))
+        return self._output_handler(data=fdir, viewfinder=dem.viewfinder,
+                                    metadata=dem.metadata, nodata=nodata_out,
+                                    mask=new_mask)
+
     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',
                   algorithm='iterative', **kwargs):
@@ -695,8 +716,10 @@ def catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4, 8, 16
             input_overrides = {'dtype' : np.int64, 'nodata' : fdir.nodata}
         elif routing.lower() == 'dinf':
             input_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
+        elif routing.lower() == 'mfd':
+            input_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
         else:
-            raise ValueError('Routing method must be one of: `d8`, `dinf`')
+            raise ValueError('Routing method must be one of: `d8`, `dinf`, `mfd`')
         kwargs.update(input_overrides)
         fdir = self._input_handler(fdir, **kwargs)
         xmin, ymin, xmax, ymax = fdir.bbox
@@ -716,6 +739,10 @@ def catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4, 8, 16
             catch = self._dinf_catchment(x, y, fdir=fdir, pour_value=pour_value, dirmap=dirmap,
                                          nodata_out=nodata_out, xytype=xytype, snap=snap,
                                          algorithm=algorithm)
+        elif routing.lower() == 'mfd':
+            catch = self._mfd_catchment(x, y, fdir=fdir, pour_value=pour_value, dirmap=dirmap,
+                                         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),
@@ -767,6 +794,31 @@ def _dinf_catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4,
                                      metadata=fdir.metadata, nodata=nodata_out)
         return catch
 
+    def _mfd_catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
+                       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
+        # to given geographic coordinate
+        if xytype in {'label', 'coordinate'}:
+            x, y = self.nearest_cell(x, y, fdir.affine, snap)
+        # Delineate the catchment
+        if algorithm.lower() == 'iterative':
+            catch = _self._mfd_catchment_iter_numba(fdir, (y, x))
+        elif algorithm.lower() == 'recursive':
+            raise NotImplementedError('Recursive algorithm not implemented.')
+        else:
+            raise ValueError('Algorithm must be `iterative` or `recursive`.')
+        if pour_value is not None:
+            catch[y, x] = pour_value
+        # Create new mask because dimension reduced
+        new_mask = fdir.mask[0]
+        catch = self._output_handler(data=catch, viewfinder=fdir.viewfinder,
+                                     metadata=fdir.metadata, nodata=nodata_out,
+                                     mask=new_mask)
+        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,
                      algorithm='iterative', **kwargs):
@@ -816,6 +868,8 @@ def accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
             fdir_overrides = {'dtype' : np.int64, 'nodata' : fdir.nodata}
         elif routing.lower() == 'dinf':
             fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
+        elif routing.lower() == 'mfd':
+            fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
         else:
             raise ValueError('Routing method must be one of: `d8`, `dinf`')
         kwargs.update(fdir_overrides)
@@ -837,6 +891,10 @@ def accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
                                           nodata_out=nodata_out,
                                           efficiency=efficiency,
                                           cycle_size=cycle_size, algorithm=algorithm)
+        elif routing.lower() == 'mfd':
+            acc = self._mfd_accumulation(fdir, weights=weights, dirmap=dirmap,
+                                         nodata_out=nodata_out,
+                                         efficiency=efficiency, algorithm=algorithm)
         return acc
 
     def _d8_accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
@@ -944,6 +1002,53 @@ def _dinf_accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16
                                    metadata=fdir.metadata, nodata=nodata_out)
         return acc
 
+    def _mfd_accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
+                          nodata_out=0., efficiency=None, algorithm='iterative', **kwargs):
+        # Find nodata cells and invalid cells
+        nodata_cells = self._get_nodata_cells(fdir)
+        # Set nodata cells to zero
+        fdir[nodata_cells] = 0.
+        # Start and end nodes
+        startnodes = np.arange(fdir[0].size, dtype=np.int64)
+        props = fdir
+        endnodes = _self._flatten_mfd_fdir_numba(props)
+        if weights is not None:
+            acc = weights.astype(np.float64).reshape(fdir[0].shape)
+        # Otherwise, initialize accumulation array to ones where valid cells exist
+        else:
+            acc = np.ones(fdir[0].shape, dtype=np.float64)
+        acc = np.asarray(acc)
+        # If using efficiency, initialize array
+        if efficiency is not None:
+            eff = efficiency.astype(np.float64).reshape(fdir[0].shape)
+            eff = np.asarray(eff)
+        # Find indegree of all cells
+        indegree = _self._mfd_bincount(endnodes)
+        # Set starting nodes to those with no predecessors
+        startnodes = startnodes[(indegree == 0)]
+        # Compute accumulation for no efficiency case
+        if efficiency is None:
+            if algorithm.lower() == 'iterative':
+                acc = _self._mfd_accumulation_iter_numba(acc, endnodes, props,
+                                                         indegree, startnodes)
+            elif algorithm.lower() == 'recursive':
+                raise NotImplementedError('Not Implemented.')
+            else:
+                raise ValueError('Algorithm must be `iterative` or `recursive`.')
+        # Compute accumulation for efficiency case
+        else:
+            if algorithm.lower() == 'iterative':
+                raise NotImplementedError('Not Implemented.')
+            elif algorithm.lower() == 'recursive':
+                raise NotImplementedError('Not Implemented.')
+            else:
+                raise ValueError('Algorithm must be `iterative` or `recursive`.')
+        new_mask = fdir.mask[0]
+        acc = self._output_handler(data=acc, viewfinder=fdir.viewfinder,
+                                   metadata=fdir.metadata, nodata=nodata_out,
+                                   mask=new_mask)
+        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', algorithm='iterative', **kwargs):
@@ -2069,7 +2174,10 @@ def _output_handler(self, data, viewfinder, metadata={}, **kwargs):
         for param, value in kwargs.items():
             if (value is not None) and (hasattr(new_view, param)):
                 setattr(new_view, param, value)
-        dataset = Raster(data, new_view, metadata=metadata)
+        if (data.ndim == 2):
+            dataset = Raster(data, new_view, metadata=metadata)
+        elif (data.ndim == 3):
+            dataset = MultiRaster(data, new_view, metadata=metadata)
         return dataset
 
     def _get_nodata_cells(self, data):
@@ -2087,15 +2195,15 @@ def _get_nodata_cells(self, data):
     def _pop_rim(self, data, nodata=0):
         left, right, top, bottom = (data[:,0].copy(), data[:,-1].copy(),
                                     data[0,:].copy(), data[-1,:].copy())
-        data[:,0] = nodata
-        data[:,-1] = nodata
-        data[0,:] = nodata
-        data[-1,:] = nodata
+        data[...,  :,  0] = nodata
+        data[...,  :, -1] = nodata
+        data[...,  0,  :] = nodata
+        data[..., -1,  :] = nodata
         return left, right, top, bottom
 
     def _replace_rim(self, data, left, right, top, bottom):
-        data[:,0] = left
-        data[:,-1] = right
-        data[0,:] = top
-        data[-1,:] = bottom
+        data[...,  :,  0] = left
+        data[...,  :, -1] = right
+        data[...,  0,  :] = top
+        data[..., -1,  :] = bottom
         return None

From d610df7283d93c6a31dcb0181b099aacbc8253c8 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Sun, 23 Jan 2022 13:57:21 -0600
Subject: [PATCH 15/27] Add mfd to distance_to_outlet, compute_hand, and
 distance_to_ridge

---
 pysheds/_sgrid.py |  25 +++++++++++
 pysheds/sgrid.py  | 109 ++++++++++++++++++++++++++++++++++++++++++++--
 pysheds/sview.py  |   5 ++-
 3 files changed, 134 insertions(+), 5 deletions(-)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index a101b5d..8f2a13f 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -681,6 +681,31 @@ def _mfd_accumulation_iter_numba(acc, fdir, props, indegree, startnodes):
                         queue.append(endnode)
     return acc
 
+@njit(float64[:,:](float64[:,:], int64[:,:,:], float64[:,:,:], uint8[:], int64[:], float64[:,:]),
+      cache=True)
+def _mfd_accumulation_eff_iter_numba(acc, fdir, props, indegree, startnodes, 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()
+            for i in range(8):
+                fdir_i = fdir[i]
+                props_i = props[i]
+                endnode = fdir_i.flat[startnode]
+                if endnode == startnode:
+                    continue
+                else:
+                    prop = props_i.flat[startnode]
+                    acc.flat[endnode] += (prop * acc.flat[startnode] * eff.flat[startnode])
+                    indegree.flat[endnode] -= 1
+                    if (indegree.flat[endnode] == 0):
+                        queue.append(endnode)
+    return acc
+
 # Functions for 'flow_distance'
 
 @njit(void(int64, int64[:,:], boolean[:,:], float64[:,:], float64[:,:],
diff --git a/pysheds/sgrid.py b/pysheds/sgrid.py
index caab2b2..63418c5 100644
--- a/pysheds/sgrid.py
+++ b/pysheds/sgrid.py
@@ -1038,7 +1038,8 @@ def _mfd_accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16,
         # Compute accumulation for efficiency case
         else:
             if algorithm.lower() == 'iterative':
-                raise NotImplementedError('Not Implemented.')
+                acc = _self._mfd_accumulation_eff_iter_numba(acc, endnodes, props,
+                                                             indegree, startnodes, eff)
             elif algorithm.lower() == 'recursive':
                 raise NotImplementedError('Not Implemented.')
             else:
@@ -1106,6 +1107,8 @@ def distance_to_outlet(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4,
             input_overrides = {'dtype' : np.int64, 'nodata' : fdir.nodata}
         elif routing.lower() == 'dinf':
             input_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
+        elif routing.lower() == 'mfd':
+            input_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
         else:
             raise ValueError('Routing method must be one of: `d8`, `dinf`')
         kwargs.update(input_overrides)
@@ -1133,6 +1136,11 @@ def distance_to_outlet(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4,
                                             dirmap=dirmap, nodata_out=nodata_out,
                                             method=method, xytype=xytype,
                                             snap=snap, algorithm=algorithm)
+        elif routing.lower() == 'mfd':
+            dist = self._mfd_flow_distance(x=x, y=y, fdir=fdir, weights=weights,
+                                           dirmap=dirmap, nodata_out=nodata_out,
+                                           method=method, xytype=xytype,
+                                           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),
@@ -1191,6 +1199,35 @@ def _dinf_flow_distance(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4
                                    metadata=fdir.metadata, nodata=nodata_out)
         return dist
 
+    def _mfd_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):
+        # Pad the rim
+        left, right, top, bottom = self._pop_rim(fdir, nodata=0)
+        # Populate weights
+        if weights is None:
+            weights = np.ones(fdir[0].shape, dtype=np.float64)
+        # If xytype is 'coordinate', delineate catchment based on cell nearest
+        # to given geographic coordinate
+        if xytype in {'label', 'coordinate'}:
+            x, y = self.nearest_cell(x, y, fdir.affine, snap)
+        # Compute flow distances
+        if method.lower() == 'shortest':
+            if algorithm.lower() == 'iterative':
+                dist = _self._mfd_flow_distance_heap_numba(fdir, (y, x), weights)
+            elif algorithm.lower() == 'recursive':
+                raise NotImplementedError('Recursive algorithm not implemented.')
+            else:
+                raise ValueError('Algorithm must be `iterative` or `recursive`.')
+        else:
+            raise NotImplementedError("Only implemented for shortest path distance.")
+        # Create new mask because dimension reduced
+        new_mask = fdir.mask[0]
+        dist = self._output_handler(data=dist, viewfinder=fdir.viewfinder,
+                                    metadata=fdir.metadata, nodata=nodata_out,
+                                    mask=new_mask)
+        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, algorithm='iterative',
                      **kwargs):
@@ -1240,6 +1277,8 @@ def compute_hand(self, fdir, dem, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
             fdir_overrides = {'dtype' : np.int64, 'nodata' : fdir.nodata}
         elif routing.lower() == 'dinf':
             fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
+        elif routing.lower() == 'mfd':
+            fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
         else:
             raise ValueError('Routing method must be one of: `d8`, `dinf`')
         dem_overrides = {'dtype' : np.float64, 'nodata' : dem.nodata}
@@ -1265,11 +1304,16 @@ def compute_hand(self, fdir, dem, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
             hand = self._dinf_compute_hand(fdir=fdir, mask=mask,
                                            nodata_out=nodata_out,
                                            algorithm=algorithm)
+        elif routing.lower() == 'mfd':
+            hand = self._mfd_compute_hand(fdir=fdir, mask=mask,
+                                           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)
-            hand = self._output_handler(data=hand, viewfinder=fdir.viewfinder,
-                                        metadata=fdir.metadata, nodata=nodata_out)
+            hand_idx = hand
+            hand = _self._assign_hand_heights_numba(hand_idx, dem, nodata_out)
+            hand = self._output_handler(data=hand, viewfinder=hand_idx.viewfinder,
+                                        metadata=hand_idx.metadata, nodata=nodata_out)
         return hand
 
     def _d8_compute_hand(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
@@ -1314,6 +1358,25 @@ def _dinf_compute_hand(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
                                     metadata=fdir.metadata, nodata=-1)
         return hand
 
+    def _mfd_compute_hand(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
+                          nodata_out=-1, algorithm='iterative'):
+        # Get nodata cells
+        nodata_cells = self._get_nodata_cells(fdir)
+        # Pad the rim
+        dirleft, dirright, dirtop, dirbottom = self._pop_rim(fdir, nodata=0.)
+        maskleft, maskright, masktop, maskbottom = self._pop_rim(mask, nodata=False)
+        if algorithm.lower() == 'iterative':
+            hand = _self._mfd_hand_iter_numba(fdir, mask)
+        elif algorithm.lower() == 'recursive':
+            raise NotImplementedError('Recursive algorithm not implemented.')
+        else:
+            raise ValueError('Algorithm must be `iterative` or `recursive`.')
+        new_mask = fdir.mask[0]
+        hand = self._output_handler(data=hand, viewfinder=fdir.viewfinder,
+                                    metadata=fdir.metadata, nodata=-1,
+                                    mask=new_mask)
+        return hand
+
     def extract_river_network(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
                               routing='d8', algorithm='iterative', **kwargs):
         """
@@ -1493,6 +1556,8 @@ def distance_to_ridge(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16,
             fdir_overrides = {'dtype' : np.int64, 'nodata' : fdir.nodata}
         elif routing.lower() == 'dinf':
             fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
+        elif routing.lower() == 'mfd':
+            fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
         else:
             raise NotImplementedError('Routing method must be one of: `d8`, `dinf`')
         kwargs.update(fdir_overrides)
@@ -1509,6 +1574,10 @@ def distance_to_ridge(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16,
             rdist = self._dinf_distance_to_ridge(fdir=fdir, weights=weights,
                                                  dirmap=dirmap, algorithm=algorithm,
                                                  nodata_out=nodata_out, cycle_size=cycle_size)
+        elif routing.lower() == 'mfd':
+            rdist = self._mfd_distance_to_ridge(fdir=fdir, weights=weights,
+                                                dirmap=dirmap, algorithm=algorithm,
+                                                nodata_out=nodata_out, cycle_size=cycle_size)
         return rdist
 
     def _d8_distance_to_ridge(self, fdir, weights, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
@@ -1579,6 +1648,38 @@ def _dinf_distance_to_ridge(self, fdir, weights, dirmap=(64, 128, 1, 2, 4, 8, 16
                                      metadata=fdir.metadata, nodata=nodata_out)
         return rdist
 
+    def _mfd_distance_to_ridge(self, fdir, weights, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
+                               algorithm='iterative', nodata_out=np.nan, cycle_size=1,
+                               **kwargs):
+        # Find nodata cells and invalid cells
+        nodata_cells = self._get_nodata_cells(fdir)
+        # Set nodata cells to zero
+        fdir[nodata_cells] = 0.
+        # Start and end nodes
+        startnodes = np.arange(fdir[0].size, dtype=np.int64)
+        props = fdir
+        endnodes = _self._flatten_mfd_fdir_numba(props)
+        if weights is None:
+            weights = np.ones(fdir[0].shape, dtype=np.float64)
+        # Find indegree of all cells
+        indegree = _self._mfd_bincount(endnodes)
+        # Set starting nodes to those with no predecessors
+        startnodes = startnodes[(indegree == 0)]
+        # Compute distance to ridge
+        rdist = np.zeros(fdir[0].shape, dtype=np.float64)
+        if algorithm.lower() == 'iterative':
+            rdist = _self._mfd_reverse_distance_iter_numba(rdist, endnodes, indegree,
+                                                           startnodes, weights)
+        elif algorithm.lower() == 'recursive':
+            raise NotImplementedError('Not Implemented.')
+        else:
+            raise ValueError('Algorithm must be `iterative` or `recursive`.')
+        new_mask = fdir.mask[0]
+        rdist = self._output_handler(data=rdist, viewfinder=fdir.viewfinder,
+                                     metadata=fdir.metadata, nodata=nodata_out,
+                                     mask=new_mask)
+        return rdist
+
     def cell_dh(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
                 nodata_out=np.nan, routing='d8', **kwargs):
         """
diff --git a/pysheds/sview.py b/pysheds/sview.py
index d48f8ef..4132450 100644
--- a/pysheds/sview.py
+++ b/pysheds/sview.py
@@ -626,6 +626,7 @@ def _view_multiraster(cls, data, target_view, data_view=None, interpolation='nea
         k, m, n = data.shape
         out = np.full((k, *target_view.shape), target_view.nodata, dtype=dtype)
         out_mask = np.ones((k, *target_view.shape), dtype=np.bool8)
+        # TODO: Is it faster to concatenate, or write to empty array?
         for i in range(k):
             slice_viewfinder = ViewFinder(affine=data_view.affine, mask=data_view.mask[i],
                                           nodata=data_view.nodata, crs=data_view.crs)
@@ -918,10 +919,12 @@ def _override_dtype(cls, data, target_view, dtype=None, interpolation='nearest')
             raise TypeError('`object` and `flexible` dtypes not allowed.')
         return dtype
 
+    # TODO: Can speed this up by giving option to not copy
     @classmethod
     def _view_same_viewfinder(cls, data, data_view, target_view, dtype,
                               apply_output_mask=True):
         has_output_mask = not target_view.mask.all()
+        # TODO: pass in an empty array and fill to avoid copying
         if (apply_output_mask) and (has_output_mask):
             out = np.where(target_view.mask, data, target_view.nodata).astype(dtype)
         else:
@@ -955,7 +958,7 @@ def _view_same_crs(cls, view, data, data_view, target_view, interpolation='neare
         x_ix, _ = cls.affine_transform(inv_affine, x,
                                        np.zeros(target_view.shape[1],
                                                 dtype=np.float64))
-        # TODO: Does this work for rotated data?
+        # TODO: Does this work for rotated data? Or is it for axis-aligned data only?
         if interpolation == 'nearest':
             view = _self._view_fill_by_axes_nearest_numba(data, view, y_ix, x_ix)
         elif interpolation == 'linear':

From c0139058727e665934b7970bb1b3838a5be73d29 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Sun, 23 Jan 2022 15:33:13 -0600
Subject: [PATCH 16/27] Update docstrings and error messages for mfd

---
 pysheds/sgrid.py | 23 +++++++++++++++--------
 1 file changed, 15 insertions(+), 8 deletions(-)

diff --git a/pysheds/sgrid.py b/pysheds/sgrid.py
index 63418c5..6f2c8de 100644
--- a/pysheds/sgrid.py
+++ b/pysheds/sgrid.py
@@ -582,6 +582,7 @@ def flowdir(self, dem, routing='d8', flats=-1, pits=-2, nodata_out=None,
                   Routing algorithm to use:
                   'd8'   : D8 flow directions
                   'dinf' : D-infinity flow directions
+                  'mfd'  : Multiple flow directions
 
         Additional keyword arguments (**kwargs) are passed to self.view.
 
@@ -618,7 +619,7 @@ def flowdir(self, dem, routing='d8', flats=-1, pits=-2, nodata_out=None,
                                      nodata_out=nodata_out, flats=flats,
                                      pits=pits, dirmap=dirmap)
         else:
-            raise ValueError('Routing method must be one of: `d8`, `dinf`')
+            raise ValueError('Routing method must be one of: `d8`, `dinf`, `mfd`')
         fdir.metadata.update(default_metadata)
         return fdir
 
@@ -694,6 +695,7 @@ def catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4, 8, 16
                   Routing algorithm to use:
                   'd8'   : D8 flow directions
                   'dinf' : D-infinity flow directions
+                  'mfd'  : Multiple flow directions
         snap : str
                Function to use for self.nearest_cell:
                'corner' : numpy.around()
@@ -847,6 +849,7 @@ def accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
                   Routing algorithm to use:
                   'd8'   : D8 flow directions
                   'dinf' : D-infinity flow directions
+                  'mfd'  : Multiple flow directions
         cycle_size : int
                      Maximum length of cycles to check for in d-infinity grids. (Note
                      that d-infinity routing can generate cycles that will cause
@@ -871,7 +874,7 @@ def accumulation(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         elif routing.lower() == 'mfd':
             fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
         else:
-            raise ValueError('Routing method must be one of: `d8`, `dinf`')
+            raise ValueError('Routing method must be one of: `d8`, `dinf`, `mfd`')
         kwargs.update(fdir_overrides)
         fdir = self._input_handler(fdir, **kwargs)
         if weights is not None:
@@ -1078,6 +1081,7 @@ def distance_to_outlet(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4,
                   Routing algorithm to use:
                   'd8'   : D8 flow directions
                   'dinf' : D-infinity flow directions
+                  'mfd'  : Multiple flow directions
         xytype : 'coordinate' or 'index'
                  How to interpret parameters 'x' and 'y'.
                      'coordinate' : x and y represent geographic coordinates
@@ -1110,7 +1114,7 @@ def distance_to_outlet(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4,
         elif routing.lower() == 'mfd':
             input_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
         else:
-            raise ValueError('Routing method must be one of: `d8`, `dinf`')
+            raise ValueError('Routing method must be one of: `d8`, `dinf`, `mfd`')
         kwargs.update(input_overrides)
         fdir = self._input_handler(fdir, **kwargs)
         if weights is not None:
@@ -1253,7 +1257,8 @@ def compute_hand(self, fdir, dem, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         routing : str
                   Routing algorithm to use:
                   'd8'   : D8 flow directions
-                  'dinf' : D-infinity flow directions (not implemented)
+                  'dinf' : D-infinity flow directions
+                  'mfd'  : Multiple flow directions
         return_index : bool
                        Boolean value indicating desired output.
                        - If True, return a Raster where each cell indicates the index
@@ -1280,7 +1285,7 @@ def compute_hand(self, fdir, dem, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         elif routing.lower() == 'mfd':
             fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
         else:
-            raise ValueError('Routing method must be one of: `d8`, `dinf`')
+            raise ValueError('Routing method must be one of: `d8`, `dinf`, `mfd`')
         dem_overrides = {'dtype' : np.float64, 'nodata' : dem.nodata}
         mask_overrides = {'dtype' : np.bool8, 'nodata' : False}
         kwargs.update(fdir_overrides)
@@ -1411,7 +1416,7 @@ def extract_river_network(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32)
         if routing.lower() == 'd8':
             fdir_overrides = {'dtype' : np.int64, 'nodata' : fdir.nodata}
         else:
-            raise NotImplementedError('Only implemented for D8 routing.')
+            raise NotImplementedError('Only implemented for `d8` routing.')
         mask_overrides = {'dtype' : np.bool8, 'nodata' : False}
         kwargs.update(fdir_overrides)
         fdir = self._input_handler(fdir, **kwargs)
@@ -1483,7 +1488,7 @@ def stream_order(self, fdir, mask, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         if routing.lower() == 'd8':
             fdir_overrides = {'dtype' : np.int64, 'nodata' : fdir.nodata}
         else:
-            raise NotImplementedError('Only implemented for D8 routing.')
+            raise NotImplementedError('Only implemented for `d8` routing.')
         mask_overrides = {'dtype' : np.bool8, 'nodata' : False}
         kwargs.update(fdir_overrides)
         fdir = self._input_handler(fdir, **kwargs)
@@ -1539,6 +1544,8 @@ def distance_to_ridge(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16,
         routing : str
                   Routing algorithm to use:
                   'd8'   : D8 flow directions
+                  'dinf' : D-infinity flow directions
+                  'mfd'  : Multiple flow directions
         algorithm : str
                     Algorithm type to use:
                     'iterative' : Use an iterative algorithm (recommended).
@@ -1559,7 +1566,7 @@ def distance_to_ridge(self, fdir, weights=None, dirmap=(64, 128, 1, 2, 4, 8, 16,
         elif routing.lower() == 'mfd':
             fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
         else:
-            raise NotImplementedError('Routing method must be one of: `d8`, `dinf`')
+            raise NotImplementedError('Routing method must be one of: `d8`, `dinf`, `mfd`')
         kwargs.update(fdir_overrides)
         fdir = self._input_handler(fdir, **kwargs)
         if weights is not None:

From 6fcddb499b82b62821fd1b4b9cbd2eccf81591bf Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Sun, 23 Jan 2022 16:32:43 -0600
Subject: [PATCH 17/27] Rewrite view functions to reduce copying

---
 pysheds/_sview.py | 20 ++++++++++++++-----
 pysheds/sview.py  | 51 +++++++++++++++++++++++++----------------------
 2 files changed, 42 insertions(+), 29 deletions(-)

diff --git a/pysheds/_sview.py b/pysheds/_sview.py
index 0ff1a5f..c667ad3 100644
--- a/pysheds/_sview.py
+++ b/pysheds/_sview.py
@@ -3,17 +3,19 @@
 from numba.types import float64, UniTuple
 
 @njit(parallel=True)
-def _view_fill_numba(data, out, y_ix, x_ix, y_passed, x_passed):
+def _view_fill_numba(data, out, y_ix, x_ix, y_passed, x_passed, nodata):
     n = x_ix.size
     m = y_ix.size
     for i in prange(m):
         for j in prange(n):
             if (y_passed[i]) & (x_passed[j]):
                 out[i, j] = data[y_ix[i], x_ix[j]]
+            else:
+                out[i, j] = nodata
     return out
 
 @njit(parallel=True)
-def _view_fill_by_axes_nearest_numba(data, out, y_ix, x_ix):
+def _view_fill_by_axes_nearest_numba(data, out, y_ix, x_ix, nodata):
     m, n = y_ix.size, x_ix.size
     M, N = data.shape
     # Currently need to use inplace form of round
@@ -27,10 +29,12 @@ def _view_fill_by_axes_nearest_numba(data, out, y_ix, x_ix):
         for j in prange(n):
             if (y_in_bounds[i]) & (x_in_bounds[j]):
                 out[i, j] = data[y_near[i], x_near[j]]
+            else:
+                out[i, j] = nodata
     return out
 
 @njit(parallel=True)
-def _view_fill_by_axes_linear_numba(data, out, y_ix, x_ix):
+def _view_fill_by_axes_linear_numba(data, out, y_ix, x_ix, nodata):
     m, n = y_ix.size, x_ix.size
     M, N = data.shape
     # Find which cells are in bounds
@@ -63,10 +67,12 @@ def _view_fill_by_axes_linear_numba(data, out, y_ix, x_ix):
                          + ( ( 1 - tx[j] ) * ty[i] * ll )
                          + ( tx[j] * ty[i] * lr ) )
                 out[i, j] = value
+            else:
+                out[i, j] = nodata
     return out
 
 @njit(parallel=True)
-def _view_fill_by_entries_nearest_numba(data, out, y_ix, x_ix):
+def _view_fill_by_entries_nearest_numba(data, out, y_ix, x_ix, nodata):
     m, n = y_ix.size, x_ix.size
     M, N = data.shape
     # Currently need to use inplace form of round
@@ -81,10 +87,12 @@ def _view_fill_by_entries_nearest_numba(data, out, y_ix, x_ix):
     for i in prange(n):
         if (y_in_bounds[i]) & (x_in_bounds[i]):
             out.flat[i] = data[y_near[i], x_near[i]]
+        else:
+            out.flat[i] = nodata
     return out
 
 @njit(parallel=True)
-def _view_fill_by_entries_linear_numba(data, out, y_ix, x_ix):
+def _view_fill_by_entries_linear_numba(data, out, y_ix, x_ix, nodata):
     m, n = y_ix.size, x_ix.size
     M, N = data.shape
     # Find which cells are in bounds
@@ -118,6 +126,8 @@ def _view_fill_by_entries_linear_numba(data, out, y_ix, x_ix):
                     + ( ( 1 - tx[i] ) * ty[i] * ll )
                     + ( tx[i] * ty[i] * lr ) )
             out.flat[i] = value
+        else:
+            out.flat[i] = nodata
     return out
 
 @njit(UniTuple(float64[:], 2)(UniTuple(float64, 9), float64[:], float64[:]), parallel=True)
diff --git a/pysheds/sview.py b/pysheds/sview.py
index 4132450..9dfec5c 100644
--- a/pysheds/sview.py
+++ b/pysheds/sview.py
@@ -608,34 +608,37 @@ def view(cls, data, target_view, data_view=None, interpolation='nearest',
 
     @classmethod
     def _view_raster(cls, data, target_view, data_view=None, interpolation='nearest',
-                     apply_output_mask=True, dtype=None):
+                     apply_output_mask=True, dtype=None, out=None):
+        # Create an output array to fill
+        if out is None:
+            out = np.empty(target_view.shape, dtype=dtype)
         # If data view and target view are the same, return a copy of the data
         if data_view.is_congruent_with(target_view):
-            out = cls._view_same_viewfinder(data, data_view, target_view, dtype,
+            out = cls._view_same_viewfinder(data, data_view, target_view, out, dtype,
                                             apply_output_mask=apply_output_mask)
         # If data view and target view are different...
         else:
-            out = cls._view_different_viewfinder(data, data_view, target_view, dtype,
+            out = cls._view_different_viewfinder(data, data_view, target_view, out, dtype,
                                                  apply_output_mask=apply_output_mask,
                                                  interpolation=interpolation)
         return out
 
     @classmethod
     def _view_multiraster(cls, data, target_view, data_view=None, interpolation='nearest',
-                          apply_output_mask=True, dtype=None):
+                          apply_output_mask=True, dtype=None, out=None):
         k, m, n = data.shape
-        out = np.full((k, *target_view.shape), target_view.nodata, dtype=dtype)
+        if out is None:
+            out = np.empty((k, *target_view.shape), dtype=dtype)
         out_mask = np.ones((k, *target_view.shape), dtype=np.bool8)
-        # TODO: Is it faster to concatenate, or write to empty array?
         for i in range(k):
             slice_viewfinder = ViewFinder(affine=data_view.affine, mask=data_view.mask[i],
                                           nodata=data_view.nodata, crs=data_view.crs)
             data_i = Raster(np.asarray(data[i]), viewfinder=slice_viewfinder)
+            # Write to out array in-place
             view_i = cls._view_raster(data_i, target_view, slice_viewfinder,
                                       interpolation=interpolation,
                                       apply_output_mask=apply_output_mask,
-                                      dtype=dtype)
-            out[i, :, :] = view_i
+                                      dtype=dtype, out=out[i, :, :])
             out_mask[i, :, :] = view_i.mask
         target_viewfinder = ViewFinder(affine=target_view.affine, mask=out_mask,
                                        nodata=target_view.nodata, crs=target_view.crs)
@@ -921,21 +924,19 @@ def _override_dtype(cls, data, target_view, dtype=None, interpolation='nearest')
 
     # TODO: Can speed this up by giving option to not copy
     @classmethod
-    def _view_same_viewfinder(cls, data, data_view, target_view, dtype,
+    def _view_same_viewfinder(cls, data, data_view, target_view, out, dtype,
                               apply_output_mask=True):
+        out[:] = data
         has_output_mask = not target_view.mask.all()
-        # TODO: pass in an empty array and fill to avoid copying
         if (apply_output_mask) and (has_output_mask):
-            out = np.where(target_view.mask, data, target_view.nodata).astype(dtype)
-        else:
-            out = np.asarray(data.copy(), dtype=dtype)
-        out = Raster(out, target_view)
-        return out
+            out[~target_view.mask] = target_view.nodata
+        out_raster = Raster(out, target_view)
+        return out_raster
 
     @classmethod
-    def _view_different_viewfinder(cls, data, data_view, target_view, dtype,
+    def _view_different_viewfinder(cls, data, data_view, target_view, out, dtype,
                                    apply_output_mask=True, interpolation='nearest'):
-        out = np.full(target_view.shape, target_view.nodata, dtype=dtype)
+        # TODO: May need to fill with nodata here
         if (data_view.crs == target_view.crs):
             out = cls._view_same_crs(out, data, data_view,
                                      target_view, interpolation)
@@ -944,9 +945,9 @@ def _view_different_viewfinder(cls, data, data_view, target_view, dtype,
                                           target_view, interpolation)
         has_output_mask = not target_view.mask.all()
         if (apply_output_mask) and (has_output_mask):
-            np.place(out, ~target_view.mask, target_view.nodata)
-        out = Raster(out, target_view)
-        return out
+            out[~target_view.mask] = target_view.nodata
+        out_raster = Raster(out, target_view)
+        return out_raster
 
     @classmethod
     def _view_same_crs(cls, view, data, data_view, target_view, interpolation='nearest'):
@@ -958,11 +959,12 @@ def _view_same_crs(cls, view, data, data_view, target_view, interpolation='neare
         x_ix, _ = cls.affine_transform(inv_affine, x,
                                        np.zeros(target_view.shape[1],
                                                 dtype=np.float64))
+        nodata = target_view.nodata
         # TODO: Does this work for rotated data? Or is it for axis-aligned data only?
         if interpolation == 'nearest':
-            view = _self._view_fill_by_axes_nearest_numba(data, view, y_ix, x_ix)
+            view = _self._view_fill_by_axes_nearest_numba(data, view, y_ix, x_ix, nodata)
         elif interpolation == 'linear':
-            view = _self._view_fill_by_axes_linear_numba(data, view, y_ix, x_ix)
+            view = _self._view_fill_by_axes_linear_numba(data, view, y_ix, x_ix, nodata)
         else:
             raise ValueError('Interpolation method must be one of: `nearest`, `linear`')
         return view
@@ -974,10 +976,11 @@ def _view_different_crs(cls, view, data, data_view, target_view, interpolation='
                                   errcheck=True, always_xy=True)
         inv_affine = ~data_view.affine
         x_ix, y_ix = cls.affine_transform(inv_affine, xt, yt)
+        nodata = target_view.nodata
         if interpolation == 'nearest':
-            view = _self._view_fill_by_entries_nearest_numba(data, view, y_ix, x_ix)
+            view = _self._view_fill_by_entries_nearest_numba(data, view, y_ix, x_ix, nodata)
         elif interpolation == 'linear':
-            view = _self._view_fill_by_entries_linear_numba(data, view, y_ix, x_ix)
+            view = _self._view_fill_by_entries_linear_numba(data, view, y_ix, x_ix, nodata)
         else:
             raise ValueError('Interpolation method must be one of: `nearest`, `linear`')
         return view

From 89ae01d01adc15878a3373f079943f1b766dabe2 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Sun, 23 Jan 2022 19:55:51 -0600
Subject: [PATCH 18/27] Add tests for mfd methods

---
 tests/test_grid.py | 25 +++++++++++++++++++++++--
 1 file changed, 23 insertions(+), 2 deletions(-)

diff --git a/tests/test_grid.py b/tests/test_grid.py
index 4fe1025..1d25f90 100644
--- a/tests/test_grid.py
+++ b/tests/test_grid.py
@@ -131,6 +131,11 @@ def test_dinf_flowdir():
     fdir_dinf = grid.flowdir(inflated_dem, dirmap=dirmap, routing='dinf')
     d.fdir_dinf = fdir_dinf
 
+def test_mfd_flowdir():
+    inflated_dem = d.inflated_dem
+    fdir_mfd = grid.flowdir(inflated_dem, dirmap=dirmap, routing='mfd')
+    d.fdir_mfd = fdir_mfd
+
 def test_clip_pad():
     catch = d.catch
     grid.clip_to(catch)
@@ -143,15 +148,17 @@ def test_clip_pad():
 def test_computed_fdir_catch():
     fdir_d8 = d.fdir_d8
     fdir_dinf = d.fdir_dinf
+    fdir_mfd = d.fdir_mfd
     catch_d8 = grid.catchment(x, y, fdir_d8, dirmap=dirmap, routing='d8',
                               xytype='coordinate')
     assert(np.count_nonzero(catch_d8) > 11300)
     # Reference routing
-    catch_d8 = grid.catchment(x, y, fdir_d8, dirmap=dirmap, routing='d8',
-                              xytype='coordinate')
     catch_dinf = grid.catchment(x, y, fdir_dinf, dirmap=dirmap, routing='dinf',
                                 xytype='coordinate')
     assert(np.count_nonzero(catch_dinf) > 11300)
+    catch_mfd = grid.catchment(x, y, fdir_mfd, dirmap=dirmap, routing='mfd',
+                               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',
@@ -163,6 +170,7 @@ def test_accumulation():
     catch = d.catch
     fdir_d8 = d.fdir_d8
     fdir_dinf = d.fdir_dinf
+    fdir_mfd = d.fdir_mfd
     # TODO: This breaks if clip_to's padding of dir is nonzero
     grid.clip_to(fdir)
     acc = grid.accumulation(fdir, dirmap=dirmap, routing='d8')
@@ -189,6 +197,8 @@ def test_accumulation():
     assert(acc_d8.max() > 11300)
     acc_dinf = grid.accumulation(fdir_dinf, dirmap=dirmap, routing='dinf')
     assert(acc_dinf.max() > 11300)
+    acc_mfd = grid.accumulation(fdir_mfd, dirmap=dirmap, routing='mfd')
+    assert(acc_mfd.max() > 11200)
     # #set nodata to 1
     # eff = grid.view(dinf_eff)
     # eff[eff==dinf_eff.nodata] = 1
@@ -207,8 +217,10 @@ def test_hand():
     dem = d.dem
     acc = d.acc
     fdir_dinf = d.fdir_dinf
+    fdir_mfd = d.fdir_mfd
     hand_d8 = grid.compute_hand(fdir, dem, acc > 100, routing='d8')
     hand_dinf = grid.compute_hand(fdir_dinf, dem, acc > 100, routing='dinf')
+    hand_mfd = grid.compute_hand(fdir_mfd, dem, acc > 100, routing='mfd')
     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',
@@ -218,6 +230,7 @@ def test_distance_to_outlet():
     fdir = d.fdir
     catch = d.catch
     fdir_dinf = d.fdir_dinf
+    fdir_mfd = d.fdir_mfd
     grid.clip_to(catch)
     dist = grid.distance_to_outlet(x, y, fdir, dirmap=dirmap, xytype='coordinate')
     assert(dist[np.isfinite(dist)].max() == max_distance_d8)
@@ -227,10 +240,14 @@ def test_distance_to_outlet():
     weights = Raster(2 * np.ones(grid.shape), grid.viewfinder)
     grid.distance_to_outlet(x, y, fdir_dinf, dirmap=dirmap, routing='dinf',
                        xytype='coordinate')
+    grid.distance_to_outlet(x, y, fdir_mfd, dirmap=dirmap, routing='mfd',
+                       xytype='coordinate')
     grid.distance_to_outlet(x, y, fdir, weights=weights,
                        dirmap=dirmap, xytype='label')
     grid.distance_to_outlet(x, y, fdir_dinf, dirmap=dirmap, weights=weights,
                        routing='dinf', xytype='label')
+    grid.distance_to_outlet(x, y, fdir_mfd, dirmap=dirmap, weights=weights,
+                       routing='mfd', xytype='label')
     # Test recursive
     grid.distance_to_outlet(x, y, fdir, dirmap=dirmap, xytype='coordinate',
                             routing='d8', algorithm='recursive')
@@ -246,8 +263,12 @@ def test_stream_order():
 def test_distance_to_ridge():
     fdir = d.fdir
     acc = d.acc
+    fdir_dinf = d.fdir_dinf
+    fdir_mfd = d.fdir_mfd
     order = grid.distance_to_ridge(fdir, acc > 100)
     order = grid.distance_to_ridge(fdir, acc > 100, algorithm='recursive')
+    order = grid.distance_to_ridge(fdir_dinf, acc > 100, routing='dinf')
+    order = grid.distance_to_ridge(fdir_mfd, acc > 100, routing='mfd')
 
 def test_cell_dh():
     fdir = d.fdir

From f81f1d60440249300e3b13870badf691020dc7b4 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Tue, 25 Jan 2022 00:34:51 -0600
Subject: [PATCH 19/27] Fix issue with efficiency test

---
 tests/test_grid.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/tests/test_grid.py b/tests/test_grid.py
index 834bb7e..7b3a85c 100644
--- a/tests/test_grid.py
+++ b/tests/test_grid.py
@@ -155,7 +155,6 @@ def test_accumulation():
     # D8 flow accumulation without efficiency
     # external flow direction
     fdir = d.fdir
-    eff = d.eff
     catch = d.catch
     fdir_d8 = d.fdir_d8
     fdir_dinf = d.fdir_dinf

From d5276552e83315de5bd597a775a7b84508863a6e Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Tue, 25 Jan 2022 18:08:35 -0600
Subject: [PATCH 20/27] Update flow direction docs with MFD

---
 docs/flow-directions.md | 116 +++++++++++++++++++++++++++++++++++++++-
 1 file changed, 115 insertions(+), 1 deletion(-)

diff --git a/docs/flow-directions.md b/docs/flow-directions.md
index 0aca37f..66afeeb 100644
--- a/docs/flow-directions.md
+++ b/docs/flow-directions.md
@@ -135,7 +135,121 @@ Raster([[  nan,   nan,   nan, ...,   nan,   nan,   nan],
 
 Note that each entry takes a value between 0 and 2π, with `np.nan` representing unknown flow directions.
 
-Note that you must also specify `routing=dinf` when using `grid.catchment` or `grid.accumulation` with a D-infinity output grid.
+Note that you must also specify `routing='dinf'` when using other functions that use the flow direction grid such as `grid.catchment` or `grid.accumulation`.
+
+## Multiple flow directions (MFD)
+
+The multiple flow direction (MFD) routing scheme partitions flow from each cell to among as many as eight of its neighbors. The proportion of the cell that flows to each of its neighbors is proportional to the height gradient between the neighboring cells.
+
+MFD routing can be selected by using the keyword argument `routing='mfd'`.
+
+```python
+fdir = grid.flowdir(inflated_dem, routing='mfd')
+```
+
+<details>
+<summary>Output...</summary>
+<p>
+
+<pre>
+fdir
+
+MultiRaster([[[0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.37428797, 0.41595555, ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.35360402, 0.42297009, ..., 0.06924557,
+               0.        , 0.        ],
+              ...,
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ]],
+
+             [[0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.36288963, 0.33088875, ..., 0.06863035,
+               0.        , 0.        ],
+              [0.        , 0.40169546, 0.36123674, ..., 0.23938736,
+               0.17013502, 0.        ],
+              ...,
+              [0.        , 0.00850506, 0.10102002, ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.        , 0.04147018, ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ]],
+
+             [[0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.14276847, 0.06932945, ..., 0.48528137,
+               0.39806072, 0.        ],
+              [0.        , 0.1217316 , 0.06042334, ..., 0.4193337 ,
+               0.48612365, 0.        ],
+              ...,
+              [0.        , 0.1683329 , 0.28176027, ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.13663963, 0.24437534, ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ]],
+
+             ...,
+
+             [[0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              ...,
+              [0.        , 0.20829874, 0.04770285, ..., 0.29010027,
+               0.31952507, 0.        ],
+              [0.        , 0.20128372, 0.11750307, ..., 0.23404662,
+               0.28716789, 0.        ],
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ]],
+
+             [[0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              ...,
+              [0.        , 0.        , 0.        , ..., 0.03460397,
+               0.06389793, 0.        ],
+              [0.        , 0.0151827 , 0.        , ..., 0.        ,
+               0.06675575, 0.        ],
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ]],
+
+             [[0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.12005394, 0.18382625, ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.12296892, 0.15536983, ..., 0.        ,
+               0.        , 0.        ],
+              ...,
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ],
+              [0.        , 0.        , 0.        , ..., 0.        ,
+               0.        , 0.        ]]])
+</pre>
+
+</p>
+</details>
+
+<br>
+
+Note that if the original DEM is of shape `(m, n)`, then the returned flow direction grid will be of shape `(8, m, n)`. The first dimension corresponds to the eight flow directions, with index 0 corresponding to North, index 1 corresponding to Northeast, index 2 corresponding to East, and so on until index 7 corresponding to Northwest. The value of a given array element (k, i, j) represents the proportion of flow that is transferred from cell (i, j) to the neighboring cell k (where k ∈ [0, 7]).
+
+Note that you must also specify `routing='mfd'` when using other functions that use the flow direction grid such as `grid.catchment` or `grid.accumulation`.
+
 
 ## Effect of map projections on routing
 

From 3f142f9cd2e59c4dd25d0aebed89037f9e209ca2 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Tue, 25 Jan 2022 19:25:49 -0600
Subject: [PATCH 21/27] Fix pour point out of bounds issue for mfd

---
 pysheds/sgrid.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pysheds/sgrid.py b/pysheds/sgrid.py
index 6f2c8de..8e721b7 100644
--- a/pysheds/sgrid.py
+++ b/pysheds/sgrid.py
@@ -730,7 +730,7 @@ def catchment(self, x, y, fdir, pour_value=None, dirmap=(64, 128, 1, 2, 4, 8, 16
                 raise ValueError('Pour point ({}, {}) is out of bounds for dataset with bbox {}.'
                                 .format(x, y, (xmin, ymin, xmax, ymax)))
         elif xytype == 'index':
-            if (x < 0) or (y < 0) or (x >= fdir.shape[1]) or (y >= fdir.shape[0]):
+            if (x < 0) or (y < 0) or (x >= fdir.shape[-1]) or (y >= fdir.shape[-2]):
                 raise ValueError('Pour point ({}, {}) is out of bounds for dataset with shape {}.'
                                 .format(x, y, fdir.shape))
         if routing.lower() == 'd8':
@@ -1127,7 +1127,7 @@ def distance_to_outlet(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4,
                 raise ValueError('Pour point ({}, {}) is out of bounds for dataset with bbox {}.'
                                 .format(x, y, (xmin, ymin, xmax, ymax)))
         elif xytype == 'index':
-            if (x < 0) or (y < 0) or (x >= fdir.shape[1]) or (y >= fdir.shape[0]):
+            if (x < 0) or (y < 0) or (x >= fdir.shape[-1]) or (y >= fdir.shape[-2]):
                 raise ValueError('Pour point ({}, {}) is out of bounds for dataset with shape {}.'
                                 .format(x, y, fdir.shape))
         if routing.lower() == 'd8':

From 4ae44b3b793866093b9698bfd3c2b3deaeea3b12 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Tue, 25 Jan 2022 19:26:07 -0600
Subject: [PATCH 22/27] Update performance tests with mfd

---
 README.md     | 29 ++++++++++++++++++-----------
 docs/index.md | 29 ++++++++++++++++++-----------
 2 files changed, 36 insertions(+), 22 deletions(-)

diff --git a/README.md b/README.md
index c0cc2a7..f7b35fa 100644
--- a/README.md
+++ b/README.md
@@ -451,17 +451,24 @@ Performance benchmarks on a 2015 MacBook Pro (M: million, K: thousand):
 
 | Function                | Routing | Number of cells          | Run time |
 | ----------------------- | ------- | ------------------------ | -------- |
-| `flowdir`               | D8      |  36M                     | 1.09 [s] |
-| `flowdir`               | DINF    |  36M                     | 6.64 [s] |
-| `accumulation`          | D8      |  36M                     | 3.65 [s] |
-| `accumulation`          | DINF    |  36M                     | 16.2 [s] |
-| `catchment`             | D8      |  9.76M                   | 3.43 [s] |
-| `catchment`             | DINF    |  9.76M                   | 5.41 [s] |
-| `distance_to_outlet`    | D8      |  9.76M                   | 4.74 [s] |
-| `distance_to_outlet`    | DINF    |  9.76M                   | 1 [m] 13 [s] |
-| `distance_to_ridge`     | D8      |  36M                     | 6.83 [s] |
-| `hand`                  | D8      |  36M total, 730K channel | 12.9 [s] |
-| `hand`                  | DINF    |  36M total, 770K channel | 18.7 [s] |
+| `flowdir`               | D8      |  36M                     | 1.14 [s] |
+| `flowdir`               | DINF    |  36M                     | 7.01 [s] |
+| `flowdir`               | MFD     |  36M                     | 4.21 [s] |
+| `accumulation`          | D8      |  36M                     | 3.44 [s] |
+| `accumulation`          | DINF    |  36M                     | 14.9 [s] |
+| `accumulation`          | MFD     |  36M                     | 32.5 [s] |
+| `catchment`             | D8      |  9.76M                   | 2.19 [s] |
+| `catchment`             | DINF    |  9.76M                   | 3.5  [s] |
+| `catchment`             | MFD     |  9.76M                   | 17.1 [s] |
+| `distance_to_outlet`    | D8      |  9.76M                   | 2.98 [s] |
+| `distance_to_outlet`    | DINF    |  9.76M                   | 5.49 [s] |
+| `distance_to_outlet`    | MFD     |  9.76M                   | 13.1 [s] |
+| `distance_to_ridge`     | D8      |  36M                     | 4.53 [s] |
+| `distance_to_ridge`     | DINF    |  36M                     | 14.5 [s] |
+| `distance_to_ridge`     | MFD     |  36M                     | 31.3 [s] |
+| `hand`                  | D8      |  36M total, 730K channel | 12.3 [s] |
+| `hand`                  | DINF    |  36M total, 770K channel | 15.8 [s] |
+| `hand`                  | MFD     |  36M total, 770K channel | 29.8 [s] |
 | `stream_order`          | D8      |  36M total, 1M channel   | 3.99 [s] |
 | `extract_river_network` | D8      |  36M total, 345K channel | 4.07 [s] |
 | `detect_pits`           | N/A     |  36M                     | 1.80 [s] |
diff --git a/docs/index.md b/docs/index.md
index 7b8d14a..8fc03bc 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -455,17 +455,24 @@ Performance benchmarks on a 2015 MacBook Pro (M: million, K: thousand):
 
 | Function                | Routing | Number of cells          | Run time |
 | ----------------------- | ------- | ------------------------ | -------- |
-| `flowdir`               | D8      |  36M                     | 1.09 [s] |
-| `flowdir`               | DINF    |  36M                     | 6.64 [s] |
-| `accumulation`          | D8      |  36M                     | 3.65 [s] |
-| `accumulation`          | DINF    |  36M                     | 16.2 [s] |
-| `catchment`             | D8      |  9.76M                   | 3.43 [s] |
-| `catchment`             | DINF    |  9.76M                   | 5.41 [s] |
-| `distance_to_outlet`    | D8      |  9.76M                   | 4.74 [s] |
-| `distance_to_outlet`    | DINF    |  9.76M                   | 1 [m] 13 [s] |
-| `distance_to_ridge`     | D8      |  36M                     | 6.83 [s] |
-| `hand`                  | D8      |  36M total, 730K channel | 12.9 [s] |
-| `hand`                  | DINF    |  36M total, 770K channel | 18.7 [s] |
+| `flowdir`               | D8      |  36M                     | 1.14 [s] |
+| `flowdir`               | DINF    |  36M                     | 7.01 [s] |
+| `flowdir`               | MFD     |  36M                     | 4.21 [s] |
+| `accumulation`          | D8      |  36M                     | 3.44 [s] |
+| `accumulation`          | DINF    |  36M                     | 14.9 [s] |
+| `accumulation`          | MFD     |  36M                     | 32.5 [s] |
+| `catchment`             | D8      |  9.76M                   | 2.19 [s] |
+| `catchment`             | DINF    |  9.76M                   | 3.5  [s] |
+| `catchment`             | MFD     |  9.76M                   | 17.1 [s] |
+| `distance_to_outlet`    | D8      |  9.76M                   | 2.98 [s] |
+| `distance_to_outlet`    | DINF    |  9.76M                   | 5.49 [s] |
+| `distance_to_outlet`    | MFD     |  9.76M                   | 13.1 [s] |
+| `distance_to_ridge`     | D8      |  36M                     | 4.53 [s] |
+| `distance_to_ridge`     | DINF    |  36M                     | 14.5 [s] |
+| `distance_to_ridge`     | MFD     |  36M                     | 31.3 [s] |
+| `hand`                  | D8      |  36M total, 730K channel | 12.3 [s] |
+| `hand`                  | DINF    |  36M total, 770K channel | 15.8 [s] |
+| `hand`                  | MFD     |  36M total, 770K channel | 29.8 [s] |
 | `stream_order`          | D8      |  36M total, 1M channel   | 3.99 [s] |
 | `extract_river_network` | D8      |  36M total, 345K channel | 4.07 [s] |
 | `detect_pits`           | N/A     |  36M                     | 1.80 [s] |

From 19fe7af0996981c8d145366faa1f4d3b1a78dc3c Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Tue, 25 Jan 2022 22:57:18 -0600
Subject: [PATCH 23/27] Add MFD to dh, cell distances, and slopes

---
 pysheds/_sgrid.py | 34 +++++++++++++++++++++++++
 pysheds/sgrid.py  | 65 ++++++++++++++++++++++++++++++++++++++++++-----
 2 files changed, 92 insertions(+), 7 deletions(-)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index 8f2a13f..d1321b1 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -1562,6 +1562,23 @@ def _dinf_cell_dh_numba(startnodes, endnodes_0, endnodes_1, props_0, props_1, de
                       prop_1 * (dem.flat[startnode] - dem.flat[endnode_1]))
     return dh
 
+@njit(parallel=True)
+def _mfd_cell_dh_numba(startnodes, endnodes, props, dem):
+    k, m, n = props.shape
+    mn = m * n
+    N = startnodes.size
+    dh = np.zeros((m, n), dtype=np.float64)
+    for i in prange(N):
+        startnode = startnodes.flat[i]
+        elev = dem.flat[startnode]
+        for j in prange(k):
+            kix = startnode + (j * mn)
+            endnode = endnodes.flat[kix]
+            prop = props.flat[kix]
+            neighbor_elev = dem.flat[endnode]
+            dh.flat[startnode] += prop * (elev - neighbor_elev)
+    return dh
+
 @njit(parallel=True)
 def _d8_cell_distances_numba(fdir, dirmap, dx, dy):
     n = fdir.size
@@ -1596,6 +1613,23 @@ def _dinf_cell_distances_numba(fdir_0, fdir_1, prop_0, prop_1, dirmap, dx, dy):
         cdist.flat[k] = dist_k
     return cdist
 
+@njit(parallel=True)
+def _mfd_cell_distances_numba(startnodes, endnodes, props, dx, dy):
+    k, m, n = props.shape
+    mn = m * n
+    N = startnodes.size
+    dd = np.sqrt(dx**2 + dy**2)
+    distances = (dy, dd, dx, dd, dy, dd, dx, dd)
+    cdist = np.zeros((m, n), dtype=np.float64)
+    for i in prange(N):
+        startnode = startnodes.flat[i]
+        for j in prange(k):
+            kix = startnode + (j * mn)
+            endnode = endnodes.flat[kix]
+            prop = props.flat[kix]
+            cdist.flat[startnode] += prop * distances[j]
+    return cdist
+
 @njit(parallel=True)
 def _cell_slopes_numba(dh, cdist):
     n = dh.size
diff --git a/pysheds/sgrid.py b/pysheds/sgrid.py
index 8e721b7..66bad12 100644
--- a/pysheds/sgrid.py
+++ b/pysheds/sgrid.py
@@ -1708,7 +1708,8 @@ def cell_dh(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         routing : str
                   Routing algorithm to use:
                   'd8'   : D8 flow directions
-                  'dinf' : D-infinity flow directions (not implemented)
+                  'dinf' : D-infinity flow directions
+                  'mfd'  : Multiple flow directions
 
         Additional keyword arguments (**kwargs) are passed to self.view.
 
@@ -1721,8 +1722,10 @@ def cell_dh(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
             fdir_overrides = {'dtype' : np.int64, 'nodata' : fdir.nodata}
         elif routing.lower() == 'dinf':
             fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
+        elif routing.lower() == 'mfd':
+            fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
         else:
-            raise ValueError('Routing method must be one of: `d8`, `dinf`')
+            raise ValueError('Routing method must be one of: `d8`, `dinf`, `mfd`')
         dem_overrides = {'dtype' : np.float64, 'nodata' : dem.nodata}
         kwargs.update(fdir_overrides)
         fdir = self._input_handler(fdir, **kwargs)
@@ -1734,6 +1737,9 @@ def cell_dh(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         elif routing.lower() == 'dinf':
             dh = self._dinf_cell_dh(dem=dem, fdir=fdir, dirmap=dirmap,
                                     nodata_out=nodata_out)
+        elif routing.lower() == 'mfd':
+            dh = self._mfd_cell_dh(dem=dem, fdir=fdir, dirmap=dirmap,
+                                    nodata_out=nodata_out)
         return dh
 
     def _d8_cell_dh(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
@@ -1753,7 +1759,7 @@ def _d8_cell_dh(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         return dh
 
     def _dinf_cell_dh(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
-                         nodata_out=np.nan):
+                      nodata_out=np.nan):
         # Get nodata cells
         nodata_cells = self._get_nodata_cells(fdir)
         # Split dinf flowdir
@@ -1771,6 +1777,23 @@ def _dinf_cell_dh(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
                                   metadata=fdir.metadata, nodata=nodata_out)
         return dh
 
+    def _mfd_cell_dh(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
+                     nodata_out=np.nan):
+        # Find nodata cells and invalid cells
+        nodata_cells = self._get_nodata_cells(fdir)
+        # Set nodata cells to zero
+        fdir[nodata_cells] = 0.
+        # Start and end nodes
+        startnodes = np.arange(fdir[0].size, dtype=np.int64)
+        props = fdir
+        endnodes = _self._flatten_mfd_fdir_numba(props)
+        dh = _self._mfd_cell_dh_numba(startnodes, endnodes, props, dem)
+        new_mask = fdir.mask[0]
+        dh = self._output_handler(data=dh, viewfinder=fdir.viewfinder,
+                                  metadata=fdir.metadata, nodata=nodata_out,
+                                  mask=new_mask)
+        return dh
+
     def cell_distances(self, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32), nodata_out=np.nan,
                        routing='d8', **kwargs):
         """
@@ -1789,7 +1812,8 @@ def cell_distances(self, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32), nodata_out=
         routing : str
                   Routing algorithm to use:
                   'd8'   : D8 flow directions
-                  'dinf' : D-infinity flow directions (not implemented)
+                  'dinf' : D-infinity flow directions
+                  'mfd'  : Multiple flow directions
 
         Additional keyword arguments (**kwargs) are passed to self.view.
 
@@ -1803,8 +1827,10 @@ def cell_distances(self, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32), nodata_out=
             fdir_overrides = {'dtype' : np.int64, 'nodata' : fdir.nodata}
         elif routing.lower() == 'dinf':
             fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
+        elif routing.lower() == 'mfd':
+            fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
         else:
-            raise ValueError('Routing method must be one of: `d8`, `dinf`')
+            raise ValueError('Routing method must be one of: `d8`, `dinf`, `mfd`')
         kwargs.update(fdir_overrides)
         fdir = self._input_handler(fdir, **kwargs)
         if routing.lower() == 'd8':
@@ -1813,6 +1839,9 @@ def cell_distances(self, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32), nodata_out=
         elif routing.lower() == 'dinf':
             cdist = self._dinf_cell_distances(fdir=fdir, dirmap=dirmap,
                                            nodata_out=nodata_out)
+        elif routing.lower() == 'mfd':
+            cdist = self._mfd_cell_distances(fdir=fdir, dirmap=dirmap,
+                                             nodata_out=nodata_out)
         return cdist
 
     def _d8_cell_distances(self, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
@@ -1849,6 +1878,25 @@ def _dinf_cell_distances(self, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
                                      metadata=fdir.metadata, nodata=nodata_out)
         return cdist
 
+    def _mfd_cell_distances(self, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
+                            nodata_out=np.nan):
+        # Find nodata cells and invalid cells
+        nodata_cells = self._get_nodata_cells(fdir)
+        # Set nodata cells to zero
+        fdir[nodata_cells] = 0.
+        # Start and end nodes
+        startnodes = np.arange(fdir[0].size, dtype=np.int64)
+        props = fdir
+        endnodes = _self._flatten_mfd_fdir_numba(props)
+        dx = abs(fdir.affine.a)
+        dy = abs(fdir.affine.e)
+        cdist = _self._mfd_cell_distances_numba(startnodes, endnodes, props, dx, dy)
+        new_mask = fdir.mask[0]
+        cdist = self._output_handler(data=cdist, viewfinder=fdir.viewfinder,
+                                     metadata=fdir.metadata, nodata=nodata_out,
+                                     mask=new_mask)
+        return cdist
+
     def cell_slopes(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32), nodata_out=np.nan,
                     routing='d8', **kwargs):
         """
@@ -1870,7 +1918,8 @@ def cell_slopes(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32), nodata_ou
         routing : str
                   Routing algorithm to use:
                   'd8'   : D8 flow directions
-                  'dinf' : D-infinity flow directions (not implemented)
+                  'dinf' : D-infinity flow directions
+                  'mfd'  : Multiple flow directions
 
         Additional keyword arguments (**kwargs) are passed to self.view.
 
@@ -1884,8 +1933,10 @@ def cell_slopes(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32), nodata_ou
             fdir_overrides = {'dtype' : np.int64, 'nodata' : fdir.nodata}
         elif routing.lower() == 'dinf':
             fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
+        elif routing.lower() == 'mfd':
+            fdir_overrides = {'dtype' : np.float64, 'nodata' : fdir.nodata}
         else:
-            raise ValueError('Routing method must be one of: `d8`, `dinf`')
+            raise ValueError('Routing method must be one of: `d8`, `dinf`, `mfd`')
         dem_overrides = {'dtype' : np.float64, 'nodata' : dem.nodata}
         kwargs.update(fdir_overrides)
         fdir = self._input_handler(fdir, **kwargs)

From 349b6060bffdb0d9376cb98f2a5b10cf04716919 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Tue, 25 Jan 2022 23:23:23 -0600
Subject: [PATCH 24/27] Add type signatures to dh, cell distances, and slopes

---
 pysheds/_sgrid.py | 29 ++++++++++++++++++++++-------
 pysheds/sgrid.py  |  8 ++++----
 2 files changed, 26 insertions(+), 11 deletions(-)

diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index d1321b1..bc491ba 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -1538,7 +1538,9 @@ def _d8_stream_network_iter_numba(fdir, indegree, orig_indegree, startnodes):
             endnode = fdir.flat[startnode]
     return profiles
 
-@njit(parallel=True)
+@njit(float64[:,:](int64[:,:], int64[:,:], float64[:,:]),
+      parallel=True,
+      cache=True)
 def _d8_cell_dh_numba(startnodes, endnodes, dem):
     n = startnodes.size
     dh = np.zeros_like(dem)
@@ -1548,7 +1550,9 @@ def _d8_cell_dh_numba(startnodes, endnodes, dem):
         dh.flat[k] = dem.flat[startnode] - dem.flat[endnode]
     return dh
 
-@njit(parallel=True)
+@njit(float64[:,:](int64[:,:], int64[:,:], int64[:,:], float64[:,:], float64[:,:], float64[:,:]),
+      parallel=True,
+      cache=True)
 def _dinf_cell_dh_numba(startnodes, endnodes_0, endnodes_1, props_0, props_1, dem):
     n = startnodes.size
     dh = np.zeros(dem.shape, dtype=np.float64)
@@ -1562,7 +1566,9 @@ def _dinf_cell_dh_numba(startnodes, endnodes_0, endnodes_1, props_0, props_1, de
                       prop_1 * (dem.flat[startnode] - dem.flat[endnode_1]))
     return dh
 
-@njit(parallel=True)
+@njit(float64[:,:](int64[:,:], int64[:,:,:], float64[:,:,:], float64[:,:]),
+      parallel=True,
+      cache=True)
 def _mfd_cell_dh_numba(startnodes, endnodes, props, dem):
     k, m, n = props.shape
     mn = m * n
@@ -1579,7 +1585,9 @@ def _mfd_cell_dh_numba(startnodes, endnodes, props, dem):
             dh.flat[startnode] += prop * (elev - neighbor_elev)
     return dh
 
-@njit(parallel=True)
+@njit(float64[:,:](int64[:,:], UniTuple(int64, 8), float64, float64),
+      parallel=True,
+      cache=True)
 def _d8_cell_distances_numba(fdir, dirmap, dx, dy):
     n = fdir.size
     cdist = np.zeros(fdir.shape, dtype=np.float64)
@@ -1593,7 +1601,10 @@ def _d8_cell_distances_numba(fdir, dirmap, dx, dy):
         cdist.flat[k] = dist_map[fdir_k]
     return cdist
 
-@njit(parallel=True)
+@njit(float64[:,:](int64[:,:], int64[:,:], float64[:,:], float64[:,:], UniTuple(int64, 8),
+                   float64, float64),
+      parallel=True,
+      cache=True)
 def _dinf_cell_distances_numba(fdir_0, fdir_1, prop_0, prop_1, dirmap, dx, dy):
     n = fdir_0.size
     cdist = np.zeros(fdir_0.shape, dtype=np.float64)
@@ -1613,7 +1624,9 @@ def _dinf_cell_distances_numba(fdir_0, fdir_1, prop_0, prop_1, dirmap, dx, dy):
         cdist.flat[k] = dist_k
     return cdist
 
-@njit(parallel=True)
+@njit(float64[:,:](int64[:,:], int64[:,:,:], float64[:,:,:], float64, float64),
+      parallel=True,
+      cache=True)
 def _mfd_cell_distances_numba(startnodes, endnodes, props, dx, dy):
     k, m, n = props.shape
     mn = m * n
@@ -1630,7 +1643,9 @@ def _mfd_cell_distances_numba(startnodes, endnodes, props, dx, dy):
             cdist.flat[startnode] += prop * distances[j]
     return cdist
 
-@njit(parallel=True)
+@njit(float64[:,:](float64[:,:], float64[:,:]),
+      parallel=True,
+      cache=True)
 def _cell_slopes_numba(dh, cdist):
     n = dh.size
     slopes = np.zeros(dh.shape, dtype=np.float64)
diff --git a/pysheds/sgrid.py b/pysheds/sgrid.py
index 66bad12..3bbfa25 100644
--- a/pysheds/sgrid.py
+++ b/pysheds/sgrid.py
@@ -1751,7 +1751,7 @@ def _d8_cell_dh(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         fdir[nodata_cells] = 0
         fdir[invalid_cells] = 0
         dirleft, dirright, dirtop, dirbottom = self._pop_rim(fdir, nodata=0)
-        startnodes = np.arange(fdir.size, dtype=np.int64)
+        startnodes = np.arange(fdir.size, dtype=np.int64).reshape(fdir.shape)
         endnodes = _self._flatten_fdir_numba(fdir, dirmap).reshape(fdir.shape)
         dh = _self._d8_cell_dh_numba(startnodes, endnodes, dem)
         dh = self._output_handler(data=dh, viewfinder=fdir.viewfinder,
@@ -1769,7 +1769,7 @@ def _dinf_cell_dh(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
                                                                      nodata=0)
         dirleft_1, dirright_1, dirtop_1, dirbottom_1 = self._pop_rim(fdir_1,
                                                                      nodata=0)
-        startnodes = np.arange(fdir.size, dtype=np.int64)
+        startnodes = np.arange(fdir.size, dtype=np.int64).reshape(fdir.shape)
         endnodes_0 = _self._flatten_fdir_numba(fdir_0, dirmap).reshape(fdir.shape)
         endnodes_1 = _self._flatten_fdir_numba(fdir_1, dirmap).reshape(fdir.shape)
         dh = _self._dinf_cell_dh_numba(startnodes, endnodes_0, endnodes_1, prop_0, prop_1, dem)
@@ -1784,7 +1784,7 @@ def _mfd_cell_dh(self, dem, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         # Set nodata cells to zero
         fdir[nodata_cells] = 0.
         # Start and end nodes
-        startnodes = np.arange(fdir[0].size, dtype=np.int64)
+        startnodes = np.arange(fdir[0].size, dtype=np.int64).reshape(fdir[0].shape)
         props = fdir
         endnodes = _self._flatten_mfd_fdir_numba(props)
         dh = _self._mfd_cell_dh_numba(startnodes, endnodes, props, dem)
@@ -1885,7 +1885,7 @@ def _mfd_cell_distances(self, fdir, dirmap=(64, 128, 1, 2, 4, 8, 16, 32),
         # Set nodata cells to zero
         fdir[nodata_cells] = 0.
         # Start and end nodes
-        startnodes = np.arange(fdir[0].size, dtype=np.int64)
+        startnodes = np.arange(fdir[0].size, dtype=np.int64).reshape(fdir[0].shape)
         props = fdir
         endnodes = _self._flatten_mfd_fdir_numba(props)
         dx = abs(fdir.affine.a)

From 2964c8d75a5cbea5cc4fd802989dc4bb25f7d08d Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Tue, 25 Jan 2022 23:34:41 -0600
Subject: [PATCH 25/27] Add tests for MFD dh, cell distances, and slopes

---
 tests/test_grid.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/tests/test_grid.py b/tests/test_grid.py
index 7b3a85c..28a3e6d 100644
--- a/tests/test_grid.py
+++ b/tests/test_grid.py
@@ -307,23 +307,29 @@ def test_distance_to_ridge():
 def test_cell_dh():
     fdir = d.fdir
     fdir_dinf = d.fdir_dinf
+    fdir_mfd = d.fdir_mfd
     dem = d.dem
     dh_d8 = grid.cell_dh(dem, fdir, routing='d8')
     dh_dinf = grid.cell_dh(dem, fdir_dinf, routing='dinf')
+    dh_mfd = grid.cell_dh(dem, fdir_mfd, routing='mfd')
 
 def test_cell_distances():
     fdir = d.fdir
     fdir_dinf = d.fdir_dinf
+    fdir_mfd = d.fdir_mfd
     dem = d.dem
     cdist_d8 = grid.cell_distances(fdir, routing='d8')
     cdist_dinf = grid.cell_distances(fdir_dinf, routing='dinf')
+    cdist_mfd = grid.cell_distances(fdir_mfd, routing='mfd')
 
 def test_cell_slopes():
     fdir = d.fdir
     fdir_dinf = d.fdir_dinf
+    fdir_mfd = d.fdir_mfd
     dem = d.dem
     slopes_d8 = grid.cell_slopes(dem, fdir, routing='d8')
     slopes_dinf = grid.cell_slopes(dem, fdir_dinf, routing='dinf')
+    slopes_mfd = grid.cell_slopes(dem, fdir_mfd, routing='mfd')
 
 # def test_set_nodata():
 #     grid.set_nodata('dir', 0)

From b0f8bf5c6dffc528f5487448628284ba16e07a1d Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Wed, 26 Jan 2022 22:05:22 -0600
Subject: [PATCH 26/27] Minor cleanup; bump version number

---
 pysheds/__init__.py |  2 +-
 pysheds/_sgrid.py   |  2 ++
 pysheds/sview.py    | 20 ++++++++------------
 setup.py            |  2 +-
 4 files changed, 12 insertions(+), 14 deletions(-)

diff --git a/pysheds/__init__.py b/pysheds/__init__.py
index 260c070..f9aa3e1 100644
--- a/pysheds/__init__.py
+++ b/pysheds/__init__.py
@@ -1 +1 @@
-__version__ = "0.3.1"
+__version__ = "0.3.2"
diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index bc491ba..dd4242f 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -1745,6 +1745,8 @@ def _flatten_fdir_numba(fdir, dirmap):
         on_bottom = (k > (n - c - 1))
         on_boundary = (on_left | on_right | on_top | on_bottom)
         if on_boundary:
+            # TODO: This seems like it could cause errors at corner points
+            # TODO: Check if offset is already zero
             if on_left:
                 offset = left_map[cell_dir]
             if on_right:
diff --git a/pysheds/sview.py b/pysheds/sview.py
index 9dfec5c..a13d13a 100644
--- a/pysheds/sview.py
+++ b/pysheds/sview.py
@@ -480,8 +480,6 @@ def view(self, raster, **kwargs):
         target_view = self
         return View.view(raster, data_view, target_view, **kwargs)
 
-# TODO: Ensure that target_view cannot be 3-dimensional
-# TODO: Or use only last two entries of viewfinder.shape
 class View():
     """
     Class containing methods for manipulating views of gridded datasets.
@@ -585,7 +583,7 @@ def view(cls, data, target_view, data_view=None, interpolation='nearest',
         # Mask input data if desired
         has_input_mask = not data_view.mask.all()
         if apply_input_mask and has_input_mask:
-            # TODO: Is this necessary?
+            # TODO: Rewrite to avoid copying.
             arr = np.where(data_view.mask, data, target_view.nodata).astype(dtype)
             data = Raster(arr, data.viewfinder, metadata=data.metadata)
         if is_multiraster or is_3d_array:
@@ -738,10 +736,10 @@ def axes(cls, affine, shape):
         y, x : tuple
                y- and x-coordinates of axes
         """
-        y_ix = np.arange(shape[0])
-        x_ix = np.arange(shape[1])
-        x_null = np.zeros(shape[0])
-        y_null = np.zeros(shape[1])
+        y_ix = np.arange(shape[-2])
+        x_ix = np.arange(shape[-1])
+        x_null = np.zeros(shape[-2])
+        y_null = np.zeros(shape[-1])
         x, _ = cls.affine_transform(affine, x_ix, y_null)
         _, y = cls.affine_transform(affine, x_null, y_ix)
         return y, x
@@ -922,7 +920,6 @@ def _override_dtype(cls, data, target_view, dtype=None, interpolation='nearest')
             raise TypeError('`object` and `flexible` dtypes not allowed.')
         return dtype
 
-    # TODO: Can speed this up by giving option to not copy
     @classmethod
     def _view_same_viewfinder(cls, data, data_view, target_view, out, dtype,
                               apply_output_mask=True):
@@ -936,7 +933,6 @@ def _view_same_viewfinder(cls, data, data_view, target_view, out, dtype,
     @classmethod
     def _view_different_viewfinder(cls, data, data_view, target_view, out, dtype,
                                    apply_output_mask=True, interpolation='nearest'):
-        # TODO: May need to fill with nodata here
         if (data_view.crs == target_view.crs):
             out = cls._view_same_crs(out, data, data_view,
                                      target_view, interpolation)
@@ -954,13 +950,13 @@ def _view_same_crs(cls, view, data, data_view, target_view, interpolation='neare
         y, x = target_view.axes
         inv_affine = ~data_view.affine
         _, y_ix = cls.affine_transform(inv_affine,
-                                       np.zeros(target_view.shape[0],
+                                       np.zeros(target_view.shape[-2],
                                                 dtype=np.float64), y)
         x_ix, _ = cls.affine_transform(inv_affine, x,
-                                       np.zeros(target_view.shape[1],
+                                       np.zeros(target_view.shape[-1],
                                                 dtype=np.float64))
         nodata = target_view.nodata
-        # TODO: Does this work for rotated data? Or is it for axis-aligned data only?
+        # TODO: Check that this works for rotated data.
         if interpolation == 'nearest':
             view = _self._view_fill_by_axes_nearest_numba(data, view, y_ix, x_ix, nodata)
         elif interpolation == 'linear':
diff --git a/setup.py b/setup.py
index 98ce916..5f62d67 100644
--- a/setup.py
+++ b/setup.py
@@ -3,7 +3,7 @@
 from setuptools import setup
 
 setup(name='pysheds',
-      version='0.3.1',
+      version='0.3.2',
       description='🌎 Simple and fast watershed delineation in python.',
       author='Matt Bartos',
       author_email='mdbartos@umich.edu',

From c683077c2f87dae10623ac7f87d023e2e4d40457 Mon Sep 17 00:00:00 2001
From: Matt Bartos <mdbartos@umich.edu>
Date: Wed, 26 Jan 2022 22:24:28 -0600
Subject: [PATCH 27/27] Rename heap functions; add more speed tests to readme

---
 README.md         |  3 ++
 docs/index.md     |  3 ++
 pysheds/_sgrid.py | 72 -----------------------------------------------
 pysheds/sgrid.py  |  4 +--
 4 files changed, 8 insertions(+), 74 deletions(-)

diff --git a/README.md b/README.md
index f7b35fa..a495edd 100644
--- a/README.md
+++ b/README.md
@@ -478,10 +478,13 @@ Performance benchmarks on a 2015 MacBook Pro (M: million, K: thousand):
 | `resolve_flats`         | N/A     |  36M                     | 9.56 [s] |
 | `cell_dh`               | D8      |  36M                     | 2.34 [s] |
 | `cell_dh`               | DINF    |  36M                     | 4.92 [s] |
+| `cell_dh`               | MFD     |  36M                     | 30.1 [s] |
 | `cell_distances`        | D8      |  36M                     | 1.11 [s] |
 | `cell_distances`        | DINF    |  36M                     | 2.16 [s] |
+| `cell_distances`        | MFD     |  36M                     | 26.8 [s] |
 | `cell_slopes`           | D8      |  36M                     | 4.01 [s] |
 | `cell_slopes`           | DINF    |  36M                     | 10.2 [s] |
+| `cell_slopes`           | MFD     |  36M                     | 58.7 [s] |
 
 Speed tests were run on a conditioned DEM from the HYDROSHEDS DEM repository
 (linked above as `elevation.tiff`).
diff --git a/docs/index.md b/docs/index.md
index 8fc03bc..5db33ca 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -482,10 +482,13 @@ Performance benchmarks on a 2015 MacBook Pro (M: million, K: thousand):
 | `resolve_flats`         | N/A     |  36M                     | 9.56 [s] |
 | `cell_dh`               | D8      |  36M                     | 2.34 [s] |
 | `cell_dh`               | DINF    |  36M                     | 4.92 [s] |
+| `cell_dh`               | MFD     |  36M                     | 30.1 [s] |
 | `cell_distances`        | D8      |  36M                     | 1.11 [s] |
 | `cell_distances`        | DINF    |  36M                     | 2.16 [s] |
+| `cell_distances`        | MFD     |  36M                     | 26.8 [s] |
 | `cell_slopes`           | D8      |  36M                     | 4.01 [s] |
 | `cell_slopes`           | DINF    |  36M                     | 10.2 [s] |
+| `cell_slopes`           | MFD     |  36M                     | 58.7 [s] |
 
 Speed tests were run on a conditioned DEM from the HYDROSHEDS DEM repository
 (linked above as `elevation.tiff`).
diff --git a/pysheds/_sgrid.py b/pysheds/_sgrid.py
index dd4242f..1fe424f 100644
--- a/pysheds/_sgrid.py
+++ b/pysheds/_sgrid.py
@@ -826,44 +826,6 @@ def _dinf_flow_distance_recur_numba(fdir_0, fdir_1, weights_0, weights_1,
 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(float64[:,:](int64[:,:], int64[:,:], float64[:,:], float64[:,:],
-                   UniTuple(int64, 2), UniTuple(int64, 8)),
-      cache=True)
-def _dinf_flow_distance_heap_numba(fdir_0, fdir_1, weights_0, weights_1,
-                                   pour_point, dirmap):
-    dist = np.full(fdir_0.shape, np.inf, dtype=np.float64)
     visited = np.zeros(fdir_0.shape, dtype=np.bool8)
     r_dirmap = np.array([dirmap[4], dirmap[5], dirmap[6],
                          dirmap[7], dirmap[0], dirmap[1],
@@ -905,40 +867,6 @@ def _dinf_flow_distance_heap_numba(fdir_0, fdir_1, weights_0, weights_1,
 @njit(float64[:,:](float64[:,:,:], UniTuple(int64, 2), float64[:,:]),
       cache=True)
 def _mfd_flow_distance_iter_numba(fdir, pour_point, weights):
-    _, m, n = fdir.shape
-    mn = m * n
-    dist = np.full((m, n), np.inf, dtype=np.float64)
-    i, j = pour_point
-    ix = (i * n) + j
-    offsets = np.array([-n, 1 - n, 1,
-                        1 + n, n, - 1 + n,
-                        - 1, - 1 - n])
-    r_dirmap = np.array([4, 5, 6, 7,
-                         0, 1, 2, 3])
-    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]
-            neighbor_dir = r_dirmap[k]
-            current_neighbor_dist = dist.flat[neighbor]
-            kix = neighbor + (neighbor_dir * mn)
-            points_to = fdir.flat[kix] > 0.
-            if points_to:
-                neighbor_dist = parent_dist + weights.flat[neighbor]
-                if (neighbor_dist < current_neighbor_dist):
-                    dist.flat[neighbor] = neighbor_dist
-                    queue.append(neighbor)
-    return dist
-
-# TODO: Weights should actually by (8, m, n)
-# neighbor_dist = parent_dist + weights.flat[kix]
-@njit(float64[:,:](float64[:,:,:], UniTuple(int64, 2), float64[:,:]),
-      cache=True)
-def _mfd_flow_distance_heap_numba(fdir, pour_point, weights):
     _, m, n = fdir.shape
     mn = m * n
     dist = np.full((m, n), np.inf, dtype=np.float64)
diff --git a/pysheds/sgrid.py b/pysheds/sgrid.py
index 3bbfa25..cf1ffc4 100644
--- a/pysheds/sgrid.py
+++ b/pysheds/sgrid.py
@@ -1189,7 +1189,7 @@ def _dinf_flow_distance(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4
             weights_1 = weights_0
         if method.lower() == 'shortest':
             if algorithm.lower() == 'iterative':
-                dist = _self._dinf_flow_distance_heap_numba(fdir_0, fdir_1, weights_0,
+                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,
@@ -1218,7 +1218,7 @@ def _mfd_flow_distance(self, x, y, fdir, weights=None, dirmap=(64, 128, 1, 2, 4,
         # Compute flow distances
         if method.lower() == 'shortest':
             if algorithm.lower() == 'iterative':
-                dist = _self._mfd_flow_distance_heap_numba(fdir, (y, x), weights)
+                dist = _self._mfd_flow_distance_iter_numba(fdir, (y, x), weights)
             elif algorithm.lower() == 'recursive':
                 raise NotImplementedError('Recursive algorithm not implemented.')
             else: