How to create a raster delineating all drainage basins using pysheds python

Question

I am going through the tutorial of "Pysheds" available in the following website
https://www.hatarilabs.com/ih-en/elevation-model-conditioning-and-stream-network-delimitation-with-python-and-pysheds-tutorial
Corresponding to all the steps mentioned in this tutorial, similar steps can be done in ArcMap (under Spatial Analyst --> Hydrology Tool).
Now the issue is I want to know is there any function available in Python Pysheds library to delineate all the drainage basins available in the DEM (after generation of flow direction map from DEM) as it is there in ArcMap's "Basin" tool under Hydrology.

Gevaert Joep · Answer

Ok, so my code while trying to solve the same issue. I got stuck at the error:
ValueError: Pour point (112, 269) is out of bounds for dataset with shape (176, 271).

where this is the code. Only a vector layer is needed, elevation is downloaded from SRTM using the elevation module:
#download DEM data from SRTM, with bounds as input
DEMLoc = '../data/Hadocha_SRTDEM.tif'
if (1):#not os.path.isfile(DEMLoc)):
    DEMLoc = os.getcwd() + '/' + DEMLoc
    bounds = [37.10384181194198, 9.545816522573169, 37.15290416533696, 9.57720460067058]
    elevation.clip(bounds=bounds,output=DEMLoc,product='SRTM1')
    gdal.Warp(destNameOrDestDS=DEMLoc, srcDSOrSrcDSTab=DEMLoc, dstSRS='EPSG:32637', resampleAlg='cubic', xRes=20, yRes=20)

grid = Grid.from_raster(DEMLoc, data_name='dem')
grid.fill_depressions(data='dem', out_name='filled')
grid.resolve_flats(data='filled',out_name='resolved')
grid.flowdir(data='resolved', out_name='dir')
grid.accumulation(data='dir', out_name='acc')

maxindex = np.where(grid.acc==np.max(grid.acc))
grid.catchment(x = maxindex[0][0], y = maxindex[1][0],data = 'dir', out_name='catch', xytype='index')

apologies for the long lines, very un-pythonic, bad me. Hope it helps on the way.

How to create a raster delineating all drainage basins using pysheds python

One Answer

Add your own answers!

Ask a Question