Question: How can I create shapefiles of the upstream area of >10,000 points? about pyflwdir HOT 2 OPEN

hi @MauKruisheer,

You can pass multiple locations to the basins method as a tuple of x and y coordinates at once, instead of looping over your geodataframe as shown in this example. The streams argument is optional and used to snap outlet points the nearest downstream stream. You can also rasterize your own stream network using e.g. rasterio to get a stream network mask instead of using the a stream order based mask (or create the stream order mask outside the loop).

Note however, that if some points are within the basins of other points you would have to sum the contributing areas of each subbasin to get the total contributing area. You could use the snap method but provide a mask of the outlet points instead of the streams to snap points to the next downstream point and use this info to determine which areas you have to sum.

from pyflwdir.

hi @DirkEilander ,

Thanks for your answer. The first part of your answer seems to work for me: many basins are being created. The second part of your answer I tried to do, but seem to fail.. Can you explain/show how this snap method works? I am now doing this:

`
with rasterio.open(flwdir_fn, "r") as src:
flwdir, _ = mask(src, watershed_lvl0.geometry, invert=False)
flwdir = flwdir.reshape(flwdir.shape[0] * flwdir.shape[1], flwdir.shape[2])
crs = src.crs
flw = pyflwdir.from_array(
flwdir,
ftype="d8",
transform=src.transform,
latlon=crs.is_geographic,
cache=True,
)

streams = readFile(intermediate_dir / 'RiverBinary_R.tif')

subbasins = flw.basins(xy=(x, y))

gdf_bas = vectorize(subbasins.astype(np.int32), 0, flw.transform, src.crs, name="basin")
gdf_bas.to_file(intermediate_dir / 'test_basins.shp')

points = pd.merge(left=geo_df, right = gdf_bas, how='inner', left_index= True, right_on= "basin")
points.rename({'geometry_x': 'geometry'}, axis=1, inplace=True)

exampleR = rioxarray.open_rasterio(flwdir_fn, masked=True).squeeze()

points_R = shpToRaster(points, exampleR, "basin").to_array().squeeze()
points_Binary = points_R.values > 0

idxs1, dists = flw.snap(xy= (x, y), mask= points_Binary*1, unit="m")
`

from pyflwdir.