// buffer nodes by tolerance then generate n-sized chunks of nodes
n = int(len(G) ** 0.5)
chunks = get_chunks(gs_nodes.loc[idx].buffer(tolerance), n)
// unary_union each chunk then unary_union those unary unions
merged_nodes = gpd.GeoSeries(chunk.unary_union for chunk in chunks).unary_union
After Change
// finally, unary_union each chunk then unary_union those unary unions
n = int(len(gs_nodes) ** 0.5)
idx = pd.DataFrame((gs_nodes.x, gs_nodes.y)).T.sort_values([0, 1]).index
buffers = gs_nodes.loc[idx].buffer(tolerance)
chunks = (buffers.iloc[i : i + n] for i in range(0, len(buffers), n))
merged_nodes = gpd.GeoSeries(chunk.unary_union for chunk in chunks).unary_union
// if only a single node results, make it iterable to convert to GeoSeries
if isinstance(merged_nodes, Polygon):