HydroSHEDS¶
HydroSHEDS stands for Hydrological data and maps based on SHuttle Elevation Derivatives at multiple Scales and provides consistent hydrographic information on a global scale. It is mainly developed by the WWF and contains various geodata sets including river networks, watershed boundaries, drainage directions, and flow accumulations.
HydroRIVERS¶
HydroRIVERS provides a vectorized line network of all global rivers.
HydroBASINS¶
HydroBASINS is a series of polygon layers of watershed boundaries on a global scale. The consistent covarage, size and hirarchy of the nested sub-basins at different scales is supported with a coding scheme that allows analysis (Lehner and Grill, 2013).
hybas = gpd.read_file('data/shp/hydrosheds/hybas_eu_lev00_v1c/hybas_eu_lev00_v1c.shp')
hybas = hybas.set_index('HYBAS_ID')
rivers = gpd.read_file('data/shp/hydrosheds/eu_riv_15s/eu_riv_15s.shp')
Let us define a point along the Rhine river and look at the corresponding main basin:
lon,lat = (7.60942,47.49567)
point = Point(lon,lat)
ID = hybas.loc[hybas.intersects(point)].index[0]
mainbas = hybas.loc[hybas['MAIN_BAS']==hybas.loc[ID,'MAIN_BAS']]
mainrivers = rivers.clip(mainbas)
fig,axs = plt.subplots(1,2,figsize=(15,5))
ax=axs[0]
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.plot(ax=ax,color='lightgrey')
ax.set_xlim(-12,30)
ax.set_ylim(40,60)
mainbas.dissolve().plot(ax=ax)
mainrivers[mainrivers.UP_CELLS>50000].plot(ax=ax,color='#FF7032')
ax.scatter(x=lon,y=lat,color='red');
ax2=axs[1]
mainbas.boundary.plot(ax=ax2)
ax2.set_xlim(6.8,8.5)
ax2.set_ylim(47,47.8)
mainrivers[mainrivers.UP_CELLS>10000].plot(ax=ax2,color='#FF7032',lw=3)
ax2.scatter(x=lon,y=lat,color='red');
Paffstetter code¶
Is a crucial coding system for nesting and topological concepts. Basically, a basin is subdivided in 9 smaller basins and labeled from 1 to 9. Each sub-basin is divided again by 9 and labeled by adding the number to the previous code. Like this, the more levels of subdividing the more digits the code will have (Verdin and Verdin, 1999). The HydroBasins dataset come with 12 levels.
The pfaf_decode script provides a function to check if a basin is upstream. Unfortunately, the script does not provide the tweek for the modified Pfafstetter code used by hydroBasins. I adapted the script following the R HydroCode package. Basically, after matching the first equal digits the test for odd numbers in the trailing part leads to a false isUpstream result when the number contains a 0. This condition has to be changed in "pfafstatter.py" line 69:
for dd in range(nth,len(pfaf_b)):
if int(pfaf_b[dd]) % 2 == 0:
to:
for dd in range(nth,len(pfaf_b)):
if (int(pfaf_b[dd]) % 2 == 0) and (int(pfaf_b[dd]) != 0):
Now we can apply the check_upstream() function to all equal or higher levels of our point of interest (POI).
levstart = int(len(str(ID).rstrip('0')))
res = mainbas[mainbas.index == ID]
for level in range(levstart,13):
strlev = str(level)
pfaf_12 = str(hybas.loc[ID,'PFAF_'+strlev])
tmp = mainbas[[pfaf.check_upstream(pfaf_a,pfaf_12,includeClose=True,verbose=False) for pfaf_a in mainbas['PFAF_'+strlev].astype(str)]]
res = pd.concat([res,tmp]).drop_duplicates()
fig2,axs = plt.subplots(1,2,figsize=(15,5))
ax2=axs[0]
mainbas.boundary.plot(ax=ax2)
ax2.set_xlim(6.8,8.5)
ax2.set_ylim(47,47.8)
mainrivers[mainrivers.UP_CELLS>10000].plot(ax=ax2,color='#FF7032',lw=3)
res.plot(ax=ax2,color='k')
ax2.scatter(x=lon,y=lat,color='red')
ax3=axs[1]
res.boundary.plot(ax=ax3)
res.plot(ax=ax3,color='k')
mainrivers.clip(res).plot(ax=ax3,color='#FF7032')
ax3.scatter(x=lon,y=lat,color='red');
Test¶
Finally we write a function and test a set of spots.
def get_watershed(point):
#get hybas ID
ID = hybas.loc[hybas.intersects(point)].index[0]
#main basin subset
mainbas = hybas.loc[hybas['MAIN_BAS']==hybas.loc[ID,'MAIN_BAS']]
#get the level
levstart = int(len(str(ID).rstrip('0')))
#loop through higher levels
res = mainbas[mainbas.index == ID]
for level in range(levstart,13):
strlev = str(level)
pfaf_12 = str(hybas.loc[ID,'PFAF_'+strlev])
#make a upstream mask with check_upstream
mask = [pfaf.check_upstream(pfaf_a,pfaf_12,includeClose=True) for pfaf_a in mainbas['PFAF_'+strlev].astype(str)]
tmp = mainbas[mask]
res = pd.concat([res,tmp]).drop_duplicates()
return res
#test some spots
spots = gpd.read_file('data/shp/spots/test_spots.shp')
sheds = gpd.GeoDataFrame()
for i,spot in spots.iterrows():
tmp = get_watershed(spot.geometry).dissolve()
tmp['name'] = spot['Name']
sheds = pd.concat([sheds,tmp.loc[:,['name','geometry']]])
sheds = sheds.reset_index(drop=True)
m = folium.Map(location=(lat,lon), zoom_start=5, tiles='CartoDB positron')
folium.GeoJson(spots,name='Spots').add_to(m)
folium.GeoJson(sheds,style_function = lambda f: {
'color':'black',
'weight':2,
'fill':True,
'opacity':0.9,
},name='Watershed').add_to(m)
folium.LayerControl().add_to(m)
m
References¶
Lehner, B., Grill G. (2013): Global river hydrography and network routing: baseline data and new approaches to study the world’s large river systems. Hydrological Processes, 27(15): 2171–2186. Data is available at www.hydrosheds.org.
Verdin, K. L., Verdin J. P. (1999): A topological system for delineation and codification of the Earth’s river basins. Journal of Hydrology, 218: 1–12. doi: 10.1016/S00221694(99)000116.