Question: How can I create shapefiles of the upstream area of >10,000 points?

[MauKruisheer]

Hi,

I would like to create polygons of the contributing area for > 10,000 points. I have all information that I need for this, including a 30m DEM, a river network (LineString geodataframe file) and derived flow direction raster. All compiled using pyflwdir.

I'm currently using the pyflwdir tool, but this takes around 5 minutes per point.. Has anyone any clues on how to approach this easier and quicker? Also, to insert the streams from a LineString geodataframe (instead of this flow.stream_order() > 4 command). This is my current script:

with rasterio.open(flwdir_fn, "r") as src:
flwdir = src.read(1)
crs = src.crs
flw = pyflwdir.from_array(
flwdir,
ftype="d8",
transform=src.transform,
latlon=crs.is_geographic,
cache=True,
)

for point in geo_df["geometry"].head():
x = point.x
y = point.y
print(x, y)
subbasins = flw.basins(xy=(x, y), streams=flw.stream_order() >= 4)

Many many thanks!

[DirkEilander]

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.

[MauKruisheer]

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