D8 and Dinf flow directions look incorrect

[markwang0]

trafficstars

Describe the bug

D8 and Dinf flow directions calculated with pysheds do not appear as expected when run on a particular DEM. The flow direction rasters have diagonal streaks of values.

To Reproduce

Download the DEM that causes the problem:

$ curl https://utexas.box.com/shared/static/2ny8fbovuy2rhq1ddz7dqka86w5vcfzy.tif -Lo OC1mTest_PM_filtered.tif 

Use pysheds to calculate flow directions on the DEM:

from pysheds.grid import Grid
# read and process DEM
dem_path = "OC1mTest_PM_filtered.tif"
grid = Grid.from_raster(dem_path)
dem = grid.read_raster(dem_path)
flooded_dem = grid.fill_depressions(dem)
inflated_dem = grid.resolve_flats(flooded_dem)
# calculate and write flow direction rasters
fdir_inf = grid.flowdir(inflated_dem, routing="dinf")
fdir_d8 = grid.flowdir(inflated_dem, routing="d8")
grid.to_raster(fdir_inf, "pysheds_dinf_fdr.tif", blockxsize=16, blockysize=16)
grid.to_raster(fdir_d8, "pysheds_d8_fdr.tif", blockxsize=16, blockysize=16)

These are the resulting rasters:

Expected behavior

I would expect the results to look like these flow direction rasters calculated with TauDEM:

$ d8flowdir -fel OC1mTest_PM_filtered.tif -p taudem_d8_fdr.tif -sd8 taudem_d8_slp.tif        
$ dinfflowdir -fel OC1mTest_PM_filtered.tif -ang taudem_dinf_fdr.tif -slp taudem_dinf_slp.tif

Environment

M1 MacBook Pro Python 3.12.3 arm64

$ pip list
Package         Version
--------------- -----------
affine          2.4.0
attrs           23.2.0
certifi         2024.2.2
click           8.1.7
click-plugins   1.1.1
cligj           0.7.2
geojson         3.1.0
imageio         2.34.1
lazy_loader     0.4
llvmlite        0.42.0
networkx        3.3
numba           0.59.1
numpy           1.26.4
packaging       24.0
pandas          2.2.2
pillow          10.3.0
pip             24.0
pyparsing       3.1.2
pyproj          3.6.1
pysheds         0.3.5
python-dateutil 2.9.0.post0
pytz            2024.1
rasterio        1.3.10
scikit-image    0.23.2
scipy           1.13.0
setuptools      69.5.1
six             1.16.0
snuggs          1.4.7
tifffile        2024.4.24
tzdata          2024.1
wheel           0.43.0

[mdbartos]

Thanks Mark, this does look strange. Will look into this as soon as I get a chance. Note that there have been some updates to the depression filling algorithm on the main branch of the github repo (#243) that may not yet be reflected on the distributed packages. It may be worth cloning from github directly and seeing if this changes the results.

[markwang0]

Thank you @mdbartos. Here is the result of running the above code with a pysheds installation from GitHub:

[groutr]

@markwang0 I think the issue here might be filling depressions. Is TauDEM filling depressions? The depression filling step fills in a large portion of the DEM. On my machine, if I skip the filling depressions step, pysheds' routing produces better looking output.

Pysheds D8 with just resolve flats Pysheds dinf with just resolving flats

[markwang0]

That makes sense @groutr -- when I skip fill_depressions and just do resolve_flats, I get results that look like yours. I did not fill depressions when running TauDEM