Preprocessing DEM with pysheds seems to produce incorrect accumulations

[groutr]

Using pysheds to fill depressions and resolve flats seems to produce incorrect flow directions and accumulations.
Preprocessing the same DEM with richdem seems to produce more usable results. The goal is to extract the stream from the bottom of the DEM.

The flow direction and accumulation steps in all cases have been computed by pysheds. I am comparing the accumulation results.

Original DEM:

Using pysheds' builtin fill depressions:

Using richdem's FillDepressions and then computing flowdir and accumulation with pysheds:

Using pysheds' builtin resolve flats:

Using richdem's ResolveFlats and then computing flowdir and accumulation with pysheds:

Any idea what could be happening?

[mdbartos]

Not sure. Did you try setting the epsilon value and number of iterations as suggested in #238?

Just making sure that you are doing fill_pits > fill_depressions > resolve_flats in that order.

I could take a shot at implementing the methods as they are programmed in richdem. Which methods specifically?

• https://github.com/r-barnes/richdem/blob/master/programs/tiled_flat_resolution/main.cpp
• https://github.com/r-barnes/richdem/blob/master/programs/parallel_priority_flood/Zhou2016pf.hpp
• https://github.com/r-barnes/richdem/blob/master/programs/parallel_priority_flood/Barnes2014pf.hpp

[mdbartos]

Another potential issue is that in this case a single depression/flat is touching multiple edges (3 in this case).

[groutr]

I should mention, that this issue arises on specific DEMs. I have observed that pysheds seems to struggle with DEMs that don't have dramatic changes in elevation. These set of DEMs that I'm working with are coastal areas.

The plots above are a 1000x1000 window of a much larger DEM (8000x8000 roughly). The smaller window was for plotting purposes only. The methods were run on the full DEM.

I am doing the operations in the correct order (except I'm not filling pits, I thought those were also handled by fill_depressions). So the order is fill_depressions > resolve_flats. However, in this specific scenario I was trying to isolate things, so the fill depressions plot is only fill_depressions and resolve flats is only resolve_flats. Though, running both fill_depressions and resolve_flats in sequence doesn't change the nature of the results. Each method in both packages were run with the default arguments. I didn't change the value of epsilon. In richdem, I am using the equivalent methods of FillDepressions > ResolveFlats.

[mdbartos]

See if the methods suggested in #238 solve the issue. The fact that it is struggling with flat areas makes me think that lowering the epsilon value could help.

[groutr]

@mdbartos I tried adjusting the epsilon value when resolving flats, however, it doesn't seem to fix the issue.
I then tried reimplementing fill_depressions using the Priority-Flood algorithm and I now get results very similar to richdem. I think the issue is using sk-image to fill depressions.

I can put together a PR later today or tomorrow.

[qpoterek]

Hello @groutr ! I had the same problem as described in this issue. Unfortunately, the tips suggested by @mdbartos didn't bring any improvement. I was wondering if you could kindly share your implementation of fill_depression mentioned in the previous comment 🥺

[groutr]

import heapq
import numpy as np

def priority_flood(dem):
    open = []
    closed = np.zeros_like(dem, dtype='bool', subok=False)
    
    y, x = dem.shape
    open.extend(zip(dem[0], range(x), repeat(0)))
    open.extend(zip(dem[y-1], range(x), repeat(y-1)))
    open.extend(zip(dem[:, 0], range(y), repeat(0)))
    open.extend(zip(dem[:, x-1], range(y), repeat(x-1)))
    heapq.heapify(open)
    closed[0] = True
    closed[y-1] = True
    closed[:, 0] = True
    closed[:, x-1] = True
    
    pq_pop = heapq.heappop
    pq_push = heapq.heappush
    row_offsets = np.array([-1, -1, 0, 1, 1, 1, 0, -1])
    col_offsets = np.array([0, 1, 1, 1, 0, -1, -1, -1])
    while open:
        c, y1, x1 = pq_pop(open)
        
        cols = x + col_offsets
        rows = y + row_offsets
        for r, c in zip(rows, cols):
            if 0 <= r < y and 0 <= c < x and not dem.mask[r, c]:
                if closed[r, c]:
                    continue
                dem[r, c] = max(dem[r, c], dem[y1, x1])
                closed[r, c] = True
                pq_push(open, (dem[r, c], r, c))

    return dem

It's a direct implementation of Algorithm 1 in the linked paper.

[mdbartos]

Hi @groutr, this is excellent. Would you be interested in creating a PR?

[groutr]

@mdbartos The code above was simply for testing agaist pysheds. For a PR, I wanted to implement the better, but slightly more complex algorithms in the paper, but unfortunately haven't had time yet.
Ultimately, I think Alg 2 (better priority-flood) and Alg 4 might be useful for pysheds.

[groutr]

I'm not familiar with optimizing numba code. My best attempt processes a DEM with shape of (12150, 9622) in about 40min. I suspect the heapq-based priority queue is really slowing things down.

Just for kicks, I wrote a cython version and it runs in about 10min on the same DEM.

[mdbartos]

I recently saw this tool for profiling numba code: https://pythonspeed.com/articles/numba-profiling/

Could be useful for diagnosing the slowdown.