def find_contributing_basins(chainids, HB_gdf):
    """
    Given an input GeoDataFrame of HydroBasins shapefiles and a list of chainids
    denoting which basins are part of the chain, this function walks upstream
    from the upstream-most basin by following the "NEXT_DOWN" attribute until
    all possible basins are included. This process is repeated, but stops when
    the previous basin is encountered. The result is a list of sets, where each
    set contains the INCREMENTAL basin indices for each subbasin. i.e. the
    most-downstream subbasin would be found by unioning all the sets.

    IMPORTANT: chainids must be arranged in US->DS direction.

    Parameters
    ----------
    chainids : list
        Denotes which basins are part of the chain
    HB_gdf : GeoDataFrame
        HydroBasins shapefiles

    Returns
    -------
    list of sets
        each set contains the incremental basin indices for each subbasin
    """
    subbasin_idcs = []
    visited_subbasins = set()

    for idx in chainids:
        sb_idcs = set([idx])
        sb_check = set(HB_gdf[HB_gdf.NEXT_DOWN == HB_gdf.HYBAS_ID[idx]].index)
        while sb_check:
            idx_check = sb_check.pop()

            if idx_check in visited_subbasins:
                continue

            sb_idcs.add(idx_check)

            basin_id_check = HB_gdf.HYBAS_ID[idx_check]
            sb_check = (
                sb_check | set(HB_gdf[HB_gdf.NEXT_DOWN == basin_id_check].index) - visited_subbasins
            )

        # Store the incremental indices
        subbasin_idcs.append(sb_idcs)

        # Update the visited subbasins (so we don't redo them)
        visited_subbasins = visited_subbasins | sb_idcs

    return subbasin_idcs