can you use grid_river_network to extract connection node coordinates?

[mnky9800n]

I am looking at grid_river_network and it seems to create a network of the river but I am actually interested in identifying a point where a tributary connects with the main river. Is there some function that identifies these points?

[mdbartos]

Thanks for raising this issue. Prior to v0.3, the function grid.extract_profiles allowed for this functionality. I can take a shot at re-implementing it this weekend.

[mdbartos]

Draft of functionality available at: https://github.com/mdbartos/pysheds/commit/959c79208070d20c6c1eaf6ef0a0c48a4464d23e

Will provide an example soon.

[mdbartos]

This functionality has been added via grid.extract_profiles in v0.3.3, now available on pypi (see PR #186).

Here's an example showing how it can be used to investigate the connectivity of a river network. Here, I am plotting the elevation profile of a contiguous chain of river segments:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors
import seaborn as sns
from pysheds.grid import Grid

grid = Grid.from_raster('elevation.tiff')
dem = grid.read_raster('elevation.tiff')

pit_filled_dem = grid.fill_pits(dem)
flooded_dem = grid.fill_depressions(pit_filled_dem)
inflated_dem = grid.resolve_flats(flooded_dem)
fdir = grid.flowdir(inflated_dem)
acc = grid.accumulation(fdir)

x, y = -97.294, 32.737
x_snap, y_snap = grid.snap_to_mask(acc > 1000, (x, y))
catch = grid.catchment(x=x_snap, y=y_snap, fdir=fdir, 
                       xytype='coordinate')
grid.clip_to(catch)
clipped_dem = grid.view(dem, apply_output_mask=False)

# Get profiles and connections
profiles, connections = grid.extract_profiles(fdir, acc > 50,
                                              include_endpoint=True)

# Iterate through all connections starting with index 14
index = 14
chain = []
while True:
    chain.append(profiles[index])
    next_index = connections[index]
    if index == next_index:
        break
    index = next_index

# View chain
chain

Output...

```python [[6485, 6484, 6650], [6650, 6816], [6816, 6815, 6981, 6980, 6979], [6979, 6811, 6810, 6642, 6474, 6473, 6472, 6304, 6136, 5969, 5802, 5635, 5467], [5467, 5301, 5302, 5135, 4968, 4801], [4801, 4634, 4467, 4300, 4134], [4134, 3968, 3801, 3635, 3469, 3470, 3304, 3137, 2970, 2803], [2803, 2635, 2467], [2467, 2299, 2131, 1963, 1796, 1629], [1629, 1462, 1296, 1129, 963, 796, 630, 464], [464, 296]] ```

Plotting the profiles:

# Plot segmented river profile
sns.set_palette('husl', len(chain))

x0 = 0
for profile in chain:
    y = clipped_dem.flat[profile]
    x = x0 + np.arange(len(profile))
    plt.plot(x, y)
    x0 += len(profile) - 1

plt.title('Segmented river profile')
plt.ylabel('Elevation [m]')
plt.xlabel('Distance [cells]')
plt.savefig('profile.png', bbox_inches='tight')