high_edge_cells is returning empty array

Hi,

While testing pysheds HAND with an NED DEM, we are getting an unexpected error using the code that ran successfully with SRTM DEM. The method we are using is the same as the one in an earlier ticket: https://github.com/mdbartos/pysheds/issues/113

The DEM information and link to the numpy array is available as indicated in the code lines below.

for l  in labels_of_labeled_flood_mask:
  #some code to clip dem to rectangular area given the labeled flood mask
  #this gives me variables: dem, dem_proj4, dem_gT
  #dem = https://www.dropbox.com/s/tvenso6dk0sopp4/dem.npy?dl=0
  #dem_proj4 = Proj('+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs', preserve_units=True)
  #dem_gT = (-96.630708, 0.0002893451223967838, 0.0, 30.0527557, 0.0, -0.00028934512285926984)
  # clip_gT()=(-96.52307161446839, 0.0002893451223967838, 0.0, 29.607742901042442, 0.0, -0.00028934512285926984)
  #cip_gT() calculates clipped geoTransform (gT) based on image coordinates, output above.
  dirmap = (64, 128, 1, 2, 4, 8, 16, 32)
  affine=Affine.from_gdal(*clip_gT(dem_gT, min0, max0, min1, max1, method='image'))
  mask=np.ones(dem.shape, dtype=np.bool)
  grid=Grid(shape=dem.shape,affine=affine, crs=dem_proj4, mask=mask)
  grid.add_gridded_data(dem, data_name='dem',affine=affine, crs=dem_proj4, mask=mask)
  grid.fill_depressions('dem', out_name='flooded_dem')
  grid.resolve_flats('flooded_dem', out_name='inflated_dem')
  grid.flowdir(data='inflated_dem', out_name='dir', dirmap=dirmap)
  grid.accumulation(data='dir', dirmap=dirmap, out_name='acc')
  grid.compute_hand('dir', 'inflated_dem', grid.acc >100, out_name='hand')
  hand_output=grid.view('hand')

The DEM looks normal, and the error happens in the grid.resolve_flats step. Below is the grid.flooded_dem.

grid_flooded_dem

Here is the stack trace:

> <ipython-input-14-150d2d14908b>(105)<module>()
    103             grid.add_gridded_data(dem, data_name='dem',affine=affine, crs=dem_proj4, mask=mask)
    104             grid.fill_depressions('dem', out_name='flooded_dem')
--> 105             grid.resolve_flats('flooded_dem', out_name='inflated_dem')
    106             grid.flowdir(data='inflated_dem', out_name='dir', dirmap=dirmap)
    107             grid.accumulation(data='dir', dirmap=dirmap, out_name='acc')

  /opt/conda/lib/python3.7/site-packages/pysheds/grid.py(3300)resolve_flats()
   3298                 dem_mask = np.where(dem.ravel() == nodata_in)[0]
   3299         inside = self._inside_indices(dem, mask=dem_mask)
-> 3300         drainage_result = self._drainage_gradient(dem, inside)
   3301         drainage_grad, flats, high_edge_cells, low_edge_cells, labels, diff = drainage_result
   3302         drainage_grad = drainage_grad.astype(np.float)

  /opt/conda/lib/python3.7/site-packages/pysheds/grid.py(3240)_drainage_gradient()
   3238         grad_from_higher = self._grad_from_higher(high_edge_cells, inner_neighbors, diff,
   3239                                                   fdir_defined, in_bounds, labels, numlabels,
-> 3240                                                   crosswalk, inside)
   3241         grad_towards_lower = self._grad_towards_lower(low_edge_cells, inner_neighbors, diff,
   3242                                                       fdir_defined, in_bounds, labels, numlabels,

  /opt/conda/lib/python3.7/site-packages/pysheds/grid.py(3129)_grad_from_higher()
   3127                           fdir_defined, in_bounds, labels, numlabels, crosswalk, inside):
   3128         z = np.zeros_like(labels)
-> 3129         max_iter = np.bincount(labels.ravel())[1:].max()
   3130         u = high_edge_cells.copy()
   3131         z.flat[inside[u]] = 1

  /opt/conda/lib/python3.7/site-packages/numpy/core/_methods.py(30)_amax()
     28 def _amax(a, axis=None, out=None, keepdims=False,
     29           initial=_NoValue, where=True):
---> 30     return umr_maximum(a, axis, None, out, keepdims, initial, where)
     31 
     32 def _amin(a, axis=None, out=None, keepdims=False,

Poking around a bit, it looks like the high_edge_cells is an empty array.

ipdb> !high_edge_cells
array([], dtype=int64)

That is generated by self._get_high_edge_cells(diff, fdir_defined). Looks like the high_edge_cells_bool is all False:

ipdb> !(diff < 0).any(axis=0)
Raster([ True,  True,  True, ...,  True,  True,  True])
ipdb> (~fdir_defined & (diff < 0).any(axis=0))
Raster([False, False, False, ..., False, False, False])
ipdb> (~fdir_defined & (diff < 0).any(axis=0)).any()
Raster(False)
ipdb> !fdir_defined.shape
(2072,)
ipdb> !fdir_defined.any()
Raster(True)

Any ideas what might be the issue and how to proceed in such a case? Let me know if there is anything else that may help you find the issue.

Best, batu.

mdbartos / pysheds

high_edge_cells is returning empty array #118