Isolate watershedΒΆ
In this script the watershed starting from the padma river, is isolated from the HydroSheds dataset.
[ ]:
import geopandas as gp
import shapely as sh
Load data
[ ]:
#water_orig = gp.read_file('as_riv_15s/as_riv_15s.shp')
[ ]:
water = water_orig.copy()
[ ]:
water = water.set_index('ARCID')
[ ]:
water.head()
[ ]:
spatial_index = water.sindex
Function to find the inflowing arc:
[ ]:
def get_inflow(spatial_index, gdf, arcid):
line = gdf.loc[arcid]['geometry']
possible_matches_index = list(spatial_index.intersection(line.bounds))
possible_matches = gdf.iloc[possible_matches_index]
precise_matches = possible_matches[possible_matches.intersects(line.boundary[0])]
return precise_matches
Function to get the corresponding id
[ ]:
def get_inflow_arcid(spatial_index, gdf, arcid):
matches = get_inflow(spatial_index, gdf, arcid)
ids = set(matches.index)
ids.remove(arcid)
return ids
[ ]:
get_inflow(spatial_index,water,770777)
[ ]:
get_inflow_arcid(spatial_index,water,770777)
Code to traverse dataset starting from arc 770777, which is the arc ending in the ocean:
[ ]:
i = 0
search_set = {770777}
network_set = search_set.copy()
import time
start = time.time()
while search_set:
new_search_set = set()
for arcid in search_set.copy():
new_search_set.update(get_inflow_arcid(spatial_index,water,arcid))
if new_search_set:
network_set.update(new_search_set)
search_set.remove(arcid)
i=i+1
print(i, end="\r")
if new_search_set:
search_set.update(new_search_set)
end = time.time()
print(end - start)
[ ]:
i
[ ]:
selection = water.loc[network_set]
[ ]:
selection = selection.reset_index()
Save watershed to shapefile
[ ]:
selection.to_file("out/padma.shp")