accumulation issue

[rogerzzl]

Firstly many thanks to the contributors for the pysheds tool. I met an issue about generating accumulation and didnot figure it out yet.

my codes:

    grid = Grid.from_raster(inputRaster)
    dem = grid.read_raster(inputRaster)
    # Fill pits in DEM
    pit_filled_dem = grid.fill_pits(dem)
    # Fill depressions in DEM
    flooded_dem = grid.fill_depressions(pit_filled_dem)
    # Resolve flats in DEM
    inflated_dem = grid.resolve_flats(flooded_dem)
    # Specify directional mapping
    dirmap = (64, 128, 1, 2, 4, 8, 16, 32)
    fdir = grid.flowdir(inflated_dem, dirmap=dirmap)
    acc = grid.accumulation(fdir, dirmap=dirmap)
    # plot
    fig, ax = plt.subplots(figsize=(8, 6))
    fig.patch.set_alpha(0)
    plt.grid('on', zorder=0)
    acc_img = np.where(grid.mask, acc + 1, np.nan)
    im = ax.imshow(acc_img, extent=grid.extent, zorder=2,
                   cmap='cubehelix',
                   norm=colors.LogNorm(1, acc.max()))
    plt.colorbar(im, ax=ax, label='Upstream Cells')
    plt.title('Flow Accumulation')
    plt.xlabel('Longitude')
    plt.ylabel('Latitude')
    plt.show()

there was a fracture in the main channel as shown in the red box in my accumulation figure (

While I used ArcPy to test that the channels at this red box position should be linked together, and this is also the reality
.
).

[mdbartos]

Thanks @rogerzzl, my guess is that if you are using D8, the flats did not resolve properly. Try the fix listed here by using a smaller epsilon value: #224 (comment)

[rogerzzl]

I tried to use small epsilon value like resolve_flats(flooded_dem, eps=1e-7), but the acc result is still the same as the previous output

[rogerzzl]

oops, I tried to look back the codes in sgrid.py, the param max_iter might cause this issue, so once I tentatively added resolve_flats(flooded_dem,eps=1e-7, max_iter=5000), the right result is there!!

Thank you Dr. Matthew Bartos for your useful suggestions and efforts with the awesome tool Pysheds, I'll close this issue then.

All the best!