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")