search

mdbartos / pysheds

:earth_americas: Simple and fast watershed delineation in python.
GNU General Public License v3.0
702 stars 188 forks source link

Role of nodata_cells in flowdir #209

Closed philippkraft closed 1 year ago

philippkraft commented 1 year ago

Currently nodata cells (as returned by sgrid._get_nodata_cells(..)) are intepreted by the three flow dir functions as obstacles for the flow, like here for D8:

https://github.com/mdbartos/pysheds/blob/9d960962bfe742c534da6577c3353169e6b6a501/pysheds/sgrid.py#L629-L630

However, for some applications nodata (or masked) cells should rather serve as outflow edges. I would be willing to contribute a pull request adding a parameter so users can indicate if they want their nodata cells to be understood either as obtacle (like now) or as an edge/sink. The ArcGIS flowdir-tool acts like the latter.

Use case: I want to write a fast USLE (or rather its German variant ABAG) calculator using free software without QGIS or GrassGIS dependencies. The erosive hillslope length is divided by landscape features like rivers, forests, streets with ditches etc. The USLE assumes surface runoff that enters these landscape features never reaches crop land again. I want to use nodata cells to indicate these sinks - my template ArcGIS model builder tool works that way.

Before starting, I would like to know:

1) are you interested? 2) do you forsee major complications on this journey? Currently I would only change the nodata values depending on the flag instead of dem.max() + 1 to dem.min() - 1.

Thank you for your software!

philippkraft commented 1 year ago

Found first a valid workaround and then a solution for me to avoid these problems. I leave it here in case others have a similar problem:

Work around

One can set the nodata values before calling the flowdir functions to values below outlet:

nodata_mask = dem == dem.nodata
dem[nodata_mask] = dem.min() - 1
flowdir = grid.flowdir(dem)
dem[nodata_mask] = dem.nodata
flowdir[nodata_mask] = flowdir.nodata

My final solution

Do not cut out the study area from the dem before processing, calculate flowdir the normal way and cut the flow intersections (barriers) into the flowdir raster. Calculate accumulation from that.

fd_cut = flowdir.copy()
fd_cut[nodata_mask | barrier] = fd_cut.nodata
acc_cut = grid.accumulation(fd_cut)