Geographic Information Systems Asked on June 22, 2021
I am going through the tutorial of "Pysheds" available in the following website
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.
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.
Answered by Gevaert Joep on June 22, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP