diff --git a/setup.py b/setup.py
index 7fa77d19a..7c406d9ba 100644
--- a/setup.py
+++ b/setup.py
@@ -24,11 +24,11 @@
 if platform.system() == 'Darwin':
     compiler_args = []
     compiler_and_linker_args = ['-stdlib=libc++', '-std=c++20']
-    library_dirs = []
+    library_dirs = [f'{os.environ['CONDA_PREFIX']}/lib']
 else:
     compiler_args = ['/std:c++20']
-    library_dirs = [r"C:\Users\runneradmin\micromamba\envs\env\Library\lib"]
-    include_dirs.append(r"C:\Users\runneradmin\micromamba\envs\env\Library\include")
+    library_dirs = [f'{os.environ['CONDA_PREFIX']}/Library/lib']
+    include_dirs.append(f'{os.environ['CONDA_PREFIX']}/Library/include')
 
 
 class build_py(_build_py):
@@ -62,7 +62,7 @@ def run(self):
             extra_link_args=compiler_and_linker_args,
             language='c++',
             libraries=['gdal'],
-            library_dirs=library_dirs,
+            library_dirs=["/Users/emily/mambaforge/envs/plugin/envs/main/lib"],
             define_macros=[("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION")]
         ) for package, module, package_compiler_args in [
             ('delineateit', 'delineateit_core', []),
@@ -79,6 +79,7 @@ def run(self):
             ('scenic_quality', 'viewshed', ['-ffp-contract=off']),
             ('ndr', 'ndr_core', []),
             ('sdr', 'sdr_core', []),
+            ('sdr', 'sdr_core_main', []),
             ('seasonal_water_yield', 'seasonal_water_yield_core', [])
         ]
     ], compiler_directives={'language_level': '3'}),
diff --git a/src/natcap/invest/sdr/sdr_core_main.pyx b/src/natcap/invest/sdr/sdr_core_main.pyx
new file mode 100644
index 000000000..91a16ba3b
--- /dev/null
+++ b/src/natcap/invest/sdr/sdr_core_main.pyx
@@ -0,0 +1,672 @@
+import logging
+import os
+
+import numpy
+import pygeoprocessing
+cimport numpy
+cimport cython
+from osgeo import gdal
+
+from cpython.mem cimport PyMem_Malloc, PyMem_Free
+from cython.operator cimport dereference as deref
+from cython.operator cimport preincrement as inc
+from libc.time cimport time as ctime
+from libcpp.pair cimport pair
+from libcpp.set cimport set as cset
+from libcpp.list cimport list as clist
+from libcpp.stack cimport stack
+cimport libc.math as cmath
+
+cdef extern from "time.h" nogil:
+    ctypedef int time_t
+    time_t time(time_t*)
+
+LOGGER = logging.getLogger(__name__)
+
+
+# cmath is supposed to have M_SQRT2, but tests have been failing recently
+# due to a missing symbol.
+cdef double SQRT2 = cmath.sqrt(2)
+cdef double PI = 3.141592653589793238462643383279502884
+# This module creates rasters with a memory xy block size of 2**BLOCK_BITS
+cdef int BLOCK_BITS = 8
+# Number of raster blocks to hold in memory at once per Managed Raster
+cdef int MANAGED_RASTER_N_BLOCKS = 2**6
+
+# These offsets are for the neighbor rows and columns according to the
+# ordering: 3 2 1
+#           4 x 0
+#           5 6 7
+cdef int *ROW_OFFSETS = [0, -1, -1, -1,  0,  1, 1, 1]
+cdef int *COL_OFFSETS = [1,  1,  0, -1, -1, -1, 0, 1]
+
+
+# this is a least recently used cache written in C++ in an external file,
+# exposing here so _ManagedRaster can use it
+cdef extern from "LRUCache.h" nogil:
+    cdef cppclass LRUCache[KEY_T, VAL_T]:
+        LRUCache(int)
+        void put(KEY_T&, VAL_T&, clist[pair[KEY_T,VAL_T]]&)
+        clist[pair[KEY_T,VAL_T]].iterator begin()
+        clist[pair[KEY_T,VAL_T]].iterator end()
+        bint exist(KEY_T &)
+        VAL_T get(KEY_T &)
+
+# this ctype is used to store the block ID and the block buffer as one object
+# inside Managed Raster
+ctypedef pair[int, double*] BlockBufferPair
+
+# a class to allow fast random per-pixel access to a raster for both setting
+# and reading pixels.  Copied from src/pygeoprocessing/routing/routing.pyx,
+# revision 891288683889237cfd3a3d0a1f09483c23489fca.
+cdef class _ManagedRaster:
+    cdef LRUCache[int, double*]* lru_cache
+    cdef cset[int] dirty_blocks
+    cdef int block_xsize
+    cdef int block_ysize
+    cdef int block_xmod
+    cdef int block_ymod
+    cdef int block_xbits
+    cdef int block_ybits
+    cdef long raster_x_size
+    cdef long raster_y_size
+    cdef int block_nx
+    cdef int block_ny
+    cdef int write_mode
+    cdef bytes raster_path
+    cdef int band_id
+    cdef int closed
+
+    def __cinit__(self, raster_path, band_id, write_mode):
+        """Create new instance of Managed Raster.
+
+        Args:
+            raster_path (char*): path to raster that has block sizes that are
+                powers of 2. If not, an exception is raised.
+            band_id (int): which band in `raster_path` to index. Uses GDAL
+                notation that starts at 1.
+            write_mode (boolean): if true, this raster is writable and dirty
+                memory blocks will be written back to the raster as blocks
+                are swapped out of the cache or when the object deconstructs.
+
+        Returns:
+            None.
+        """
+        raster_info = pygeoprocessing.get_raster_info(raster_path)
+        self.raster_x_size, self.raster_y_size = raster_info['raster_size']
+        self.block_xsize, self.block_ysize = raster_info['block_size']
+        self.block_xmod = self.block_xsize-1
+        self.block_ymod = self.block_ysize-1
+
+        if not (1 <= band_id <= raster_info['n_bands']):
+            err_msg = (
+                "Error: band ID (%s) is not a valid band number. "
+                "This exception is happening in Cython, so it will cause a "
+                "hard seg-fault, but it's otherwise meant to be a "
+                "ValueError." % (band_id))
+            print(err_msg)
+            raise ValueError(err_msg)
+        self.band_id = band_id
+
+        if (self.block_xsize & (self.block_xsize - 1) != 0) or (
+                self.block_ysize & (self.block_ysize - 1) != 0):
+            # If inputs are not a power of two, this will at least print
+            # an error message. Unfortunately with Cython, the exception will
+            # present itself as a hard seg-fault, but I'm leaving the
+            # ValueError in here at least for readability.
+            err_msg = (
+                "Error: Block size is not a power of two: "
+                "block_xsize: %d, %d, %s. This exception is happening"
+                "in Cython, so it will cause a hard seg-fault, but it's"
+                "otherwise meant to be a ValueError." % (
+                    self.block_xsize, self.block_ysize, raster_path))
+            print(err_msg)
+            raise ValueError(err_msg)
+
+        self.block_xbits = numpy.log2(self.block_xsize)
+        self.block_ybits = numpy.log2(self.block_ysize)
+        self.block_nx = (
+            self.raster_x_size + (self.block_xsize) - 1) // self.block_xsize
+        self.block_ny = (
+            self.raster_y_size + (self.block_ysize) - 1) // self.block_ysize
+
+        self.lru_cache = new LRUCache[int, double*](MANAGED_RASTER_N_BLOCKS)
+        self.raster_path = <bytes> raster_path
+        self.write_mode = write_mode
+        self.closed = 0
+
+    def __dealloc__(self):
+        """Deallocate _ManagedRaster.
+
+        This operation manually frees memory from the LRUCache and writes any
+        dirty memory blocks back to the raster if `self.write_mode` is True.
+        """
+        self.close()
+
+    def close(self):
+        """Close the _ManagedRaster and free up resources.
+
+            This call writes any dirty blocks to disk, frees up the memory
+            allocated as part of the cache, and frees all GDAL references.
+
+            Any subsequent calls to any other functions in _ManagedRaster will
+            have undefined behavior.
+        """
+        if self.closed:
+            return
+        self.closed = 1
+        cdef int xi_copy, yi_copy
+        cdef numpy.ndarray[double, ndim=2] block_array = numpy.empty(
+            (self.block_ysize, self.block_xsize))
+        cdef double *double_buffer
+        cdef int block_xi
+        cdef int block_yi
+        # initially the win size is the same as the block size unless
+        # we're at the edge of a raster
+        cdef int win_xsize
+        cdef int win_ysize
+
+        # we need the offsets to subtract from global indexes for cached array
+        cdef int xoff
+        cdef int yoff
+
+        cdef clist[BlockBufferPair].iterator it = self.lru_cache.begin()
+        cdef clist[BlockBufferPair].iterator end = self.lru_cache.end()
+        if not self.write_mode:
+            while it != end:
+                # write the changed value back if desired
+                PyMem_Free(deref(it).second)
+                inc(it)
+            return
+
+        raster = gdal.OpenEx(
+            self.raster_path, gdal.GA_Update | gdal.OF_RASTER)
+        raster_band = raster.GetRasterBand(self.band_id)
+
+        # if we get here, we're in write_mode
+        cdef cset[int].iterator dirty_itr
+        while it != end:
+            double_buffer = deref(it).second
+            block_index = deref(it).first
+
+            # write to disk if block is dirty
+            dirty_itr = self.dirty_blocks.find(block_index)
+            if dirty_itr != self.dirty_blocks.end():
+                self.dirty_blocks.erase(dirty_itr)
+                block_xi = block_index % self.block_nx
+                block_yi = block_index // self.block_nx
+
+                # we need the offsets to subtract from global indexes for
+                # cached array
+                xoff = block_xi << self.block_xbits
+                yoff = block_yi << self.block_ybits
+
+                win_xsize = self.block_xsize
+                win_ysize = self.block_ysize
+
+                # clip window sizes if necessary
+                if xoff+win_xsize > self.raster_x_size:
+                    win_xsize = win_xsize - (
+                        xoff+win_xsize - self.raster_x_size)
+                if yoff+win_ysize > self.raster_y_size:
+                    win_ysize = win_ysize - (
+                        yoff+win_ysize - self.raster_y_size)
+
+                for xi_copy in xrange(win_xsize):
+                    for yi_copy in xrange(win_ysize):
+                        block_array[yi_copy, xi_copy] = (
+                            double_buffer[
+                                (yi_copy << self.block_xbits) + xi_copy])
+                raster_band.WriteArray(
+                    block_array[0:win_ysize, 0:win_xsize],
+                    xoff=xoff, yoff=yoff)
+            PyMem_Free(double_buffer)
+            inc(it)
+        raster_band.FlushCache()
+        raster_band = None
+        raster = None
+
+    cdef inline void set(self, long xi, long yi, double value):
+        """Set the pixel at `xi,yi` to `value`."""
+        cdef int block_xi = xi >> self.block_xbits
+        cdef int block_yi = yi >> self.block_ybits
+        # this is the flat index for the block
+        cdef int block_index = block_yi * self.block_nx + block_xi
+        if not self.lru_cache.exist(block_index):
+            self._load_block(block_index)
+        self.lru_cache.get(
+            block_index)[
+                ((yi & (self.block_ymod))<<self.block_xbits) +
+                (xi & (self.block_xmod))] = value
+        if self.write_mode:
+            dirty_itr = self.dirty_blocks.find(block_index)
+            if dirty_itr == self.dirty_blocks.end():
+                self.dirty_blocks.insert(block_index)
+
+    cdef inline double get(self, long xi, long yi):
+        """Return the value of the pixel at `xi,yi`."""
+        cdef int block_xi = xi >> self.block_xbits
+        cdef int block_yi = yi >> self.block_ybits
+        # this is the flat index for the block
+        cdef int block_index = block_yi * self.block_nx + block_xi
+        if not self.lru_cache.exist(block_index):
+            self._load_block(block_index)
+        return self.lru_cache.get(
+            block_index)[
+                ((yi & (self.block_ymod))<<self.block_xbits) +
+                (xi & (self.block_xmod))]
+
+    cdef void _load_block(self, int block_index) except *:
+        cdef int block_xi = block_index % self.block_nx
+        cdef int block_yi = block_index // self.block_nx
+
+        # we need the offsets to subtract from global indexes for cached array
+        cdef int xoff = block_xi << self.block_xbits
+        cdef int yoff = block_yi << self.block_ybits
+
+        cdef int xi_copy, yi_copy
+        cdef numpy.ndarray[double, ndim=2] block_array
+        cdef double *double_buffer
+        cdef clist[BlockBufferPair] removed_value_list
+
+        # determine the block aligned xoffset for read as array
+
+        # initially the win size is the same as the block size unless
+        # we're at the edge of a raster
+        cdef int win_xsize = self.block_xsize
+        cdef int win_ysize = self.block_ysize
+
+        # load a new block
+        if xoff+win_xsize > self.raster_x_size:
+            win_xsize = win_xsize - (xoff+win_xsize - self.raster_x_size)
+        if yoff+win_ysize > self.raster_y_size:
+            win_ysize = win_ysize - (yoff+win_ysize - self.raster_y_size)
+
+        raster = gdal.OpenEx(self.raster_path, gdal.OF_RASTER)
+        raster_band = raster.GetRasterBand(self.band_id)
+        block_array = raster_band.ReadAsArray(
+            xoff=xoff, yoff=yoff, win_xsize=win_xsize,
+            win_ysize=win_ysize).astype(
+            numpy.float64)
+        raster_band = None
+        raster = None
+        double_buffer = <double*>PyMem_Malloc(
+            (sizeof(double) << self.block_xbits) * win_ysize)
+        for xi_copy in xrange(win_xsize):
+            for yi_copy in xrange(win_ysize):
+                double_buffer[(yi_copy<<self.block_xbits)+xi_copy] = (
+                    block_array[yi_copy, xi_copy])
+        self.lru_cache.put(
+            <int>block_index, <double*>double_buffer, removed_value_list)
+
+        if self.write_mode:
+            raster = gdal.OpenEx(
+                self.raster_path, gdal.GA_Update | gdal.OF_RASTER)
+            raster_band = raster.GetRasterBand(self.band_id)
+
+        block_array = numpy.empty(
+            (self.block_ysize, self.block_xsize), dtype=numpy.double)
+        while not removed_value_list.empty():
+            # write the changed value back if desired
+            double_buffer = removed_value_list.front().second
+
+            if self.write_mode:
+                block_index = removed_value_list.front().first
+
+                # write back the block if it's dirty
+                dirty_itr = self.dirty_blocks.find(block_index)
+                if dirty_itr != self.dirty_blocks.end():
+                    self.dirty_blocks.erase(dirty_itr)
+
+                    block_xi = block_index % self.block_nx
+                    block_yi = block_index // self.block_nx
+
+                    xoff = block_xi << self.block_xbits
+                    yoff = block_yi << self.block_ybits
+
+                    win_xsize = self.block_xsize
+                    win_ysize = self.block_ysize
+
+                    if xoff+win_xsize > self.raster_x_size:
+                        win_xsize = win_xsize - (
+                            xoff+win_xsize - self.raster_x_size)
+                    if yoff+win_ysize > self.raster_y_size:
+                        win_ysize = win_ysize - (
+                            yoff+win_ysize - self.raster_y_size)
+
+                    for xi_copy in xrange(win_xsize):
+                        for yi_copy in xrange(win_ysize):
+                            block_array[yi_copy, xi_copy] = double_buffer[
+                                (yi_copy << self.block_xbits) + xi_copy]
+                    raster_band.WriteArray(
+                        block_array[0:win_ysize, 0:win_xsize],
+                        xoff=xoff, yoff=yoff)
+            PyMem_Free(double_buffer)
+            removed_value_list.pop_front()
+
+        if self.write_mode:
+            raster_band = None
+            raster = None
+
+
+def calculate_sediment_deposition(
+        mfd_flow_direction_path, e_prime_path, f_path, sdr_path,
+        target_sediment_deposition_path):
+    """Calculate sediment deposition layer.
+
+    This algorithm outputs both sediment deposition (t_i) and flux (f_i)::
+
+        t_i  =      dr_i  * (sum over j ∈ J of f_j * p(i,j)) + E'_i
+
+        f_i  = (1 - dr_i) * (sum over j ∈ J of f_j * p(i,j)) + E'_i
+
+
+                (sum over k ∈ K of SDR_k * p(i,k)) - SDR_i
+        dr_i = --------------------------------------------
+                              (1 - SDR_i)
+
+    where:
+
+    - ``p(i,j)`` is the proportion of flow from pixel ``i`` into pixel ``j``
+    - ``J`` is the set of pixels that are immediate upslope neighbors of
+      pixel ``i``
+    - ``K`` is the set of pixels that are immediate downslope neighbors of
+      pixel ``i``
+    - ``E'`` is ``USLE * (1 - SDR)``, the amount of sediment loss from pixel
+      ``i`` that doesn't reach a stream (``e_prime_path``)
+    - ``SDR`` is the sediment delivery ratio (``sdr_path``)
+
+    ``f_i`` is recursively defined in terms of ``i``'s upslope neighbors.
+    The algorithm begins from seed pixels that are local high points and so
+    have no upslope neighbors. It works downslope from each seed pixel,
+    only adding a pixel to the stack when all its upslope neighbors are
+    already calculated.
+
+    Note that this function is designed to be used in the context of the SDR
+    model. Because the algorithm is recursive upslope and downslope of each
+    pixel, nodata values in the SDR input would propagate along the flow path.
+    This case is not handled because we assume the SDR and flow dir inputs
+    will come from the SDR model and have nodata in the same places.
+
+    Args:
+        mfd_flow_direction_path (string): a path to a raster with
+            pygeoprocessing.routing MFD flow direction values.
+        e_prime_path (string): path to a raster that shows sources of
+            sediment that wash off a pixel but do not reach the stream.
+        f_path (string): path to a raster that shows the sediment flux
+            on a pixel for sediment that does not reach the stream.
+        sdr_path (string): path to Sediment Delivery Ratio raster.
+        target_sediment_deposition_path (string): path to created that
+            shows where the E' sources end up across the landscape.
+
+    Returns:
+        None.
+
+    """
+    LOGGER.info('Calculate sediment deposition')
+    cdef float target_nodata = -1
+    pygeoprocessing.new_raster_from_base(
+        mfd_flow_direction_path, target_sediment_deposition_path,
+        gdal.GDT_Float32, [target_nodata])
+    pygeoprocessing.new_raster_from_base(
+        mfd_flow_direction_path, f_path,
+        gdal.GDT_Float32, [target_nodata])
+
+    cdef _ManagedRaster mfd_flow_direction_raster = _ManagedRaster(
+        mfd_flow_direction_path, 1, False)
+    cdef _ManagedRaster e_prime_raster = _ManagedRaster(
+        e_prime_path, 1, False)
+    cdef _ManagedRaster sdr_raster = _ManagedRaster(sdr_path, 1, False)
+    cdef _ManagedRaster f_raster = _ManagedRaster(f_path, 1, True)
+    cdef _ManagedRaster sediment_deposition_raster = _ManagedRaster(
+        target_sediment_deposition_path, 1, True)
+
+    # given the pixel neighbor numbering system
+    #  3 2 1
+    #  4 x 0
+    #  5 6 7
+    # if a pixel `x` has a neighbor `n` in position `i`,
+    # then `n`'s neighbor in position `inflow_offsets[i]`
+    # is the original pixel `x`
+    cdef int *inflow_offsets = [4, 5, 6, 7, 0, 1, 2, 3]
+
+    cdef long n_cols, n_rows
+    flow_dir_info = pygeoprocessing.get_raster_info(mfd_flow_direction_path)
+    n_cols, n_rows = flow_dir_info['raster_size']
+    cdef int mfd_nodata = 0
+    cdef stack[long] processing_stack
+    cdef float sdr_nodata = pygeoprocessing.get_raster_info(
+        sdr_path)['nodata'][0]
+    cdef float e_prime_nodata = pygeoprocessing.get_raster_info(
+        e_prime_path)['nodata'][0]
+    cdef long col_index, row_index, win_xsize, win_ysize, xoff, yoff
+    cdef long global_col, global_row, j, k
+    cdef long flat_index
+    cdef long seed_col = 0
+    cdef long seed_row = 0
+    cdef long neighbor_row, neighbor_col, ds_neighbor_row, ds_neighbor_col
+    cdef int flow_val, neighbor_flow_val, ds_neighbor_flow_val
+    cdef int flow_weight, neighbor_flow_weight
+    cdef float flow_sum, neighbor_flow_sum
+    cdef float downslope_sdr_weighted_sum, sdr_i, sdr_j
+    cdef float p_j, p_val
+    cdef unsigned long n_pixels_processed = 0
+    cdef time_t last_log_time = ctime(NULL)
+
+    for offset_dict in pygeoprocessing.iterblocks(
+            (mfd_flow_direction_path, 1), offset_only=True, largest_block=0):
+        win_xsize = offset_dict['win_xsize']
+        win_ysize = offset_dict['win_ysize']
+        xoff = offset_dict['xoff']
+        yoff = offset_dict['yoff']
+
+        if ctime(NULL) - last_log_time > 5.0:
+            last_log_time = ctime(NULL)
+            LOGGER.info('Sediment deposition %.2f%% complete', 100 * (
+                n_pixels_processed / float(n_cols * n_rows)))
+
+        for row_index in range(win_ysize):
+            seed_row = yoff + row_index
+            for col_index in range(win_xsize):
+                seed_col = xoff + col_index
+                # check if this is a good seed pixel ( a local high point)
+                if mfd_flow_direction_raster.get(seed_col, seed_row) == mfd_nodata:
+                    continue
+                seed_pixel = 1
+                # iterate over each of the pixel's neighbors
+                for j in range(8):
+                    # skip if the neighbor is outside the raster bounds
+                    neighbor_row = seed_row + ROW_OFFSETS[j]
+                    if neighbor_row < 0 or neighbor_row >= n_rows:
+                        continue
+                    neighbor_col = seed_col + COL_OFFSETS[j]
+                    if neighbor_col < 0 or neighbor_col >= n_cols:
+                        continue
+                    # skip if the neighbor's flow direction is undefined
+                    neighbor_flow_val = <int>mfd_flow_direction_raster.get(
+                        neighbor_col, neighbor_row)
+                    if neighbor_flow_val == mfd_nodata:
+                        continue
+                    # if the neighbor flows into it, it's not a local high
+                    # point and so can't be a seed pixel
+                    neighbor_flow_weight = (
+                        neighbor_flow_val >> (inflow_offsets[j]*4)) & 0xF
+                    if neighbor_flow_weight > 0:
+                        seed_pixel = 0  # neighbor flows in, not a seed
+                        break
+
+                # if this can be a seed pixel and hasn't already been
+                # calculated, put it on the stack
+                if seed_pixel and sediment_deposition_raster.get(seed_col, seed_row) == target_nodata:
+                    processing_stack.push(seed_row * n_cols + seed_col)
+
+                while processing_stack.size() > 0:
+                    # loop invariant: cell has all upslope neighbors
+                    # processed. this is true for seed pixels because they
+                    # have no upslope neighbors.
+                    flat_index = processing_stack.top()
+                    processing_stack.pop()
+                    global_row = flat_index // n_cols
+                    global_col = flat_index % n_cols
+
+                    # (sum over j ∈ J of f_j * p(i,j) in the equation for t_i)
+                    # calculate the upslope f_j contribution to this pixel,
+                    # the weighted sum of flux flowing onto this pixel from
+                    # all neighbors
+                    f_j_weighted_sum = 0
+                    for j in range(8):
+                        neighbor_row = global_row + ROW_OFFSETS[j]
+                        if neighbor_row < 0 or neighbor_row >= n_rows:
+                            continue
+                        neighbor_col = global_col + COL_OFFSETS[j]
+                        if neighbor_col < 0 or neighbor_col >= n_cols:
+                            continue
+
+                        # see if there's an inflow from the neighbor to the
+                        # pixel
+                        neighbor_flow_val = (
+                            <int>mfd_flow_direction_raster.get(
+                                neighbor_col, neighbor_row))
+                        neighbor_flow_weight = (
+                            neighbor_flow_val >> (inflow_offsets[j]*4)) & 0xF
+                        if neighbor_flow_weight > 0:
+                            f_j = f_raster.get(neighbor_col, neighbor_row)
+                            if f_j == target_nodata:
+                                continue
+                            # sum up the neighbor's flow dir values in each
+                            # direction.
+                            # flow dir values are relative to the total
+                            neighbor_flow_sum = 0
+                            for k in range(8):
+                                neighbor_flow_sum += (
+                                    neighbor_flow_val >> (k*4)) & 0xF
+                            # get the proportion of the neighbor's flow that
+                            # flows into the original pixel
+                            p_val = neighbor_flow_weight / neighbor_flow_sum
+                            # add the neighbor's flux value, weighted by the
+                            # flow proportion
+                            f_j_weighted_sum += p_val * f_j
+
+                    # calculate sum of SDR values of immediate downslope
+                    # neighbors, weighted by proportion of flow into each
+                    # neighbor
+                    # (sum over k ∈ K of SDR_k * p(i,k) in the equation above)
+                    downslope_sdr_weighted_sum = 0
+                    flow_val = <int>mfd_flow_direction_raster.get(
+                        global_col, global_row)
+                    flow_sum = 0
+                    for k in range(8):
+                        flow_sum += (flow_val >> (k*4)) & 0xF
+
+                    # iterate over the neighbors again
+                    for j in range(8):
+                        # skip if neighbor is outside the raster boundaries
+                        neighbor_row = global_row + ROW_OFFSETS[j]
+                        if neighbor_row < 0 or neighbor_row >= n_rows:
+                            continue
+                        neighbor_col = global_col + COL_OFFSETS[j]
+                        if neighbor_col < 0 or neighbor_col >= n_cols:
+                            continue
+                        # if it is a downslope neighbor, add to the sum and
+                        # check if it can be pushed onto the stack yet
+                        flow_weight = (flow_val >> (j*4)) & 0xF
+                        if flow_weight > 0:
+                            sdr_j = sdr_raster.get(neighbor_col, neighbor_row)
+                            if sdr_j == sdr_nodata:
+                                continue
+                            if sdr_j == 0:
+                                # this means it's a stream, for SDR deposition
+                                # purposes, we set sdr to 1 to indicate this
+                                # is the last step on which to retain sediment
+                                sdr_j = 1
+                            p_j = flow_weight / flow_sum
+                            downslope_sdr_weighted_sum += sdr_j * p_j
+
+                            # check if we can add neighbor j to the stack yet
+                            #
+                            # if there is a downslope neighbor it
+                            # couldn't have been pushed on the processing
+                            # stack yet, because the upslope was just
+                            # completed
+                            upslope_neighbors_processed = 1
+                            # iterate over each neighbor-of-neighbor
+                            for k in range(8):
+                                # no need to push the one we're currently
+                                # calculating back onto the stack
+                                if inflow_offsets[k] == j:
+                                    continue
+                                # skip if neighbor-of-neighbor is outside
+                                # raster bounds
+                                ds_neighbor_row = (
+                                    neighbor_row + ROW_OFFSETS[k])
+                                if ds_neighbor_row < 0 or ds_neighbor_row >= n_rows:
+                                    continue
+                                ds_neighbor_col = (
+                                    neighbor_col + COL_OFFSETS[k])
+                                if ds_neighbor_col < 0 or ds_neighbor_col >= n_cols:
+                                    continue
+                                # if any upslope neighbor of j hasn't been
+                                # calculated, we can't push j onto the stack
+                                # yet
+                                ds_neighbor_flow_val = (
+                                    <int>mfd_flow_direction_raster.get(
+                                        ds_neighbor_col, ds_neighbor_row))
+                                if (ds_neighbor_flow_val >> (
+                                        inflow_offsets[k]*4)) & 0xF > 0:
+                                    if (sediment_deposition_raster.get(
+                                            ds_neighbor_col, ds_neighbor_row) ==
+                                            target_nodata):
+                                        upslope_neighbors_processed = 0
+                                        break
+                            # if all upslope neighbors of neighbor j are
+                            # processed, we can push j onto the stack.
+                            if upslope_neighbors_processed:
+                                processing_stack.push(
+                                    neighbor_row * n_cols +
+                                    neighbor_col)
+
+                    # nodata pixels should propagate to the results
+                    sdr_i = sdr_raster.get(global_col, global_row)
+                    if sdr_i == sdr_nodata:
+                        continue
+                    e_prime_i = e_prime_raster.get(global_col, global_row)
+                    if e_prime_i == e_prime_nodata:
+                        continue
+
+                    # This condition reflects property A in the user's guide.
+                    if downslope_sdr_weighted_sum < sdr_i:
+                        # i think this happens because of our low resolution
+                        # flow direction, it's okay to zero out.
+                        downslope_sdr_weighted_sum = sdr_i
+
+                    # these correspond to the full equations for
+                    # dr_i, t_i, and f_i given in the docstring
+                    if sdr_i == 1:
+                        # This reflects property B in the user's guide and is
+                        # an edge case to avoid division-by-zero.
+                        dr_i = 1
+                    else:
+                        dr_i = (downslope_sdr_weighted_sum - sdr_i) / (1 - sdr_i)
+
+                    # Lisa's modified equations
+                    t_i = dr_i * f_j_weighted_sum  # deposition, a.k.a trapped sediment
+                    f_i = (1 - dr_i) * f_j_weighted_sum + e_prime_i  # flux
+
+                    # On large flow paths, it's possible for dr_i, f_i and t_i
+                    # to have very small negative values that are numerically
+                    # equivalent to 0. These negative values were raising
+                    # questions on the forums and it's easier to clamp the
+                    # values here than to explain IEEE 754.
+                    if dr_i < 0:
+                        dr_i = 0
+                    if t_i < 0:
+                        t_i = 0
+                    if f_i < 0:
+                        f_i = 0
+
+                    sediment_deposition_raster.set(global_col, global_row, t_i)
+                    f_raster.set(global_col, global_row, f_i)
+        n_pixels_processed += (win_xsize * win_ysize)
+
+    LOGGER.info('Sediment deposition 100% complete')
+    sediment_deposition_raster.close()