search

mdbartos / pysheds

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

catchment and flow_distance functions of pysheds give different areas for the same pour point #163

Closed yonaskd closed 2 years ago

yonaskd commented 2 years ago

Distance

catchment

Hi Everyone,

I am new for the pysheds. When i tried to use the example codes from mdbratos to delineate a watershed using the elevation data from my study area, it is giving different areas for catchment and flow_distance even if i used the same x, y (as both index and label). Did i missed something? Thanks for your insight.

Yonas

mdbartos commented 2 years ago

Hi Yonas,

I just pushed major changes yesterday. I’d recommend upgrading to v0.3 first. After that, can you post your code?

MDB

yonaskd commented 2 years ago

Hi Matt,

Thanks for the quick response. I will update the package. I have fixed this particular issue as described below. The problem happened when i used the entire watershed flow direction (grid.dir) for delineating the catchment and determining the flow distance as follow

grid = Grid.from_raster('dem.tif', data_name='dem') grid.fill_depressions(grid.dem, out_name='flooded_dem') grid.resolve_flats('flooded_dem', out_name='inflated_dem') dirmap = (64, 128, 1, 2, 4, 8, 16, 32) grid.flowdir(data='inflated_dem', out_name='dir', dirmap=dirmap)

c,r = 4288, 916 # pour point. This point was also used directly x = c, y = r and xytype ='index' and give the same catch and dist. x,y = grid.affine*(c,r) grid.catchment(x=x, y=y, data='dir', dirmap=dirmap, out_name='catch', recursionlimit=15000, xytype='label') grid.flow_distance(x=x, y=y, data='dir', dirmap=dirmap, out_name='dist', xytype='label', nodata_in=np.nan)

The above code is the one which provided different watershed areas for catchment and flow distance even though the x and y coordinate are the same. But, when i used the clipped flow direction (or 'catch' instead of 'dir') in the flow_distance function, it provides the same area as the catchment as expected.

grid.clip_to('catch') catch = grid.view('catch',nodata=np.nan) grid.flow_distance(x=x, y=y, data=catch, dirmap=dirmap, out_name='dist', xytype='label', nodata_in=np.nan)

The above code fixed the issue.

Side note, one thing i noted is that flow distance computation takes considerably longer time, which may become my next my main concern since the reason i am inclining to pysheds instead of arcpy is the potential to use the python multiprocessing functionality in HPC since i will perform thousands of watershed delineations and hydrological analyses. Thanks again.

Yonas

Distance2

mdbartos commented 2 years ago

Great, if you are interested in HPC I'd highly recommend upgrading. pysheds v0.3 is anywhere between 10-100x faster than previous versions due to numba integration.