Testing the pysheds with new data set MERIT-Hydro; second 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/).

The following code does not work for an outlet (the first outlet location); I think this might be related to clip_to functionality and it fails in the extract_river_network.

However it works for an outlet a bit more downstream (which the first outlet that fails is nested and inside)

Thank you in advance.

# 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 = -90.7775, 31.50306 # outlet that doesnt work
x, y = -90.7999, 31.48755 # outlet that works! the outlet that doesnt work is nest inside this outlet that works

# read the bil or geotiff of direction from MERIT-Hydro Direction (Yamazaki et al., 2019)
grid = Grid.from_raster('n30w095_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_2A.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_2B.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_2C.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_2D.png')

and here is the underlying direction map: n30w095_dir.tif.zip

Error is:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-6-7820009a4787> in <module>
     35 # extracting river netwrok
     36 branches = None
---> 37 branches = grid.extract_river_network(fdir='catch', acc='acc',
     38                                       threshold=100, dirmap=dirmap)
     39 

~/opt/anaconda3/envs/myenv38/lib/python3.8/site-packages/pysheds/grid.py in extract_river_network(self, fdir, acc, threshold, dirmap, nodata_in, routing, apply_mask, ignore_metadata, **kwargs)
   2765             A = scipy.sparse.lil_matrix((fdir.size, fdir.size))
   2766             for i,j in zip(start[no_fork], end[no_fork]):
-> 2767                 A[i,j] = 1
   2768             n_components, labels = csgraph.connected_components(A)
   2769             u, inverse, c = np.unique(labels, return_inverse=True, return_counts=True)

~/opt/anaconda3/envs/myenv38/lib/python3.8/site-packages/scipy/sparse/lil.py in __setitem__(self, key, x)
    329             if x.size > 1:
    330                 raise ValueError("Trying to assign a sequence to an item")
--> 331             return self._set_intXint(key[0], key[1], x)
    332         # Everything else takes the normal path.
    333         IndexMixin.__setitem__(self, key, x)

~/opt/anaconda3/envs/myenv38/lib/python3.8/site-packages/scipy/sparse/lil.py in _set_intXint(self, row, col, x)
    298 
    299     def _set_intXint(self, row, col, x):
--> 300         _csparsetools.lil_insert(self.shape[0], self.shape[1], self.rows,
    301                                  self.data, row, col, x)
    302 

_csparsetools.pyx in scipy.sparse._csparsetools.lil_insert()

_csparsetools.pyx in scipy.sparse._csparsetools.lil_insert()

IndexError: column index (118580) out of bounds