Testing the pysheds with new data set MERIT-Hydro; first issue

[ShervanGharari]

Hello,

I am trying to replicate the basin delineation. I have face some issue using the newly released data set, MERIT-HYDRO from Yamazaki et al, 2019 (http://hydro.iis.u-tokyo.ac.jp/~yamadai/MERIT_Hydro/).

There are few issues with the following code that tries to delineate a basin from the direction of MERIT-Hydro data set.

# import necessary packages
from pysheds.grid import Grid
import matplotlib.pyplot as plt
import numpy as np
import json
import geopandas as gpd
from shapely.geometry import Point
from geopandas import GeoSeries
# set the font and font size for plots
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams.update({'font.size': 16})

# input
dirmap = (64,  128,  1,   2,    4,   8,    16,  32) # there data include 0 and 247 as no data apart from dir values
x, y = -71.8342, 41.4748 # outlet

# read the bil or geotiff of direction from MERIT-Hydro Direction (Yamazaki et al., 2019)
grid = Grid.from_raster('n40w075_dir.tif', data_name='dir',  nodata=0)

# Compute perliminary accumulation for snapping the gauges on the river network
grid.accumulation(data='dir', out_name='acc', dirmap=dirmap)

# putting the outlet point exactly on the river network
xy = np.column_stack([x, y])
new_xy = grid.snap_to_mask(grid.acc > 1000, xy, return_dist=False)
new_xs, new_ys = new_xy[:,0], new_xy[:,1]
print('corrected gauges on the river network :',new_xs, new_ys)

# Delineate the catchment
grid.catchment(data='dir', x=new_xs, y=new_ys, dirmap=dirmap, out_name='catch',
               recursionlimit=150000, xytype='label')
grid.clip_to('catch')

# extracting river netwrok
branches = None
branches = grid.extract_river_network(fdir='catch', acc='acc',
                                      threshold=100, dirmap=dirmap)

# show the deliniated basin
plt.figure(figsize = (7,7))
plt.imshow(grid.catch.astype(int), extent=grid.extent, cmap='jet')#, zorder=1)
plt.colorbar(label='Masekd flow direction (-)')
plt.grid(zorder=0)
plt.title('Masekd flow direction map')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.tight_layout()
plt.savefig('Figure_1A.png')

# show the direction for the clip_to area
plt.figure(figsize = (7,7))
plt.imshow(grid.dir.astype(int), extent=grid.extent, cmap='jet')#, zorder=1)
plt.colorbar(label='Masekd flow direction (-)')
plt.grid(zorder=0)
plt.title('Masekd flow direction map')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.tight_layout()
plt.savefig('Figure_1B.png')

# show the flow accumulation
plt.figure(figsize = (7,7))
plt.imshow(grid.acc.astype(int), extent=grid.extent, cmap='jet')#, zorder=1)
plt.colorbar(label='Masekd flow accumulation (-)')
plt.grid(zorder=0)
plt.title('Masekd flow accumulation map')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.tight_layout()
plt.savefig('Figure_1C.png')

# plotting the river segments
fig, ax = plt.subplots(figsize=(7,7))
plt.grid('on', zorder=0)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('River network (>1000 accumulation)')
plt.xlim(grid.bbox[0], grid.bbox[2]) # the bounding box from grid
plt.ylim(grid.bbox[1], grid.bbox[3]) # the bounding box from grid
ax.set_aspect('equal')

for branch in branches['features']:
    line = np.asarray(branch['geometry']['coordinates'])
    plt.plot(line[:, 0], line[:, 1])
plt.savefig('Figure_1D.png')

The visualization result in few issues; first of all the clip_to seems not working as the direction keeps it extend. Also the clipped dir is not shown properly. and thirdly the river segment has a weird loop at the end of the most downstream segment.

I believe I am doing everything correctly. What can go wrong here. Thank you in advance.

The underlying flow direction is uploaded here.