mdbartos / pysheds

:earth_americas: Simple and fast watershed delineation in python.
GNU General Public License v3.0
708 stars 191 forks source link

grid.accumulation ValueError: 'list' argument must have no negative elements #175

Open fluviotect opened 2 years ago

fluviotect commented 2 years ago

ValueError Traceback (most recent call last) ~\AppData\Local\Temp/ipykernel_15104/190316472.py in 3 4 # Compute accumulation ----> 5 acc = grid.accumulation(fdir, apply_output_mask=False)

~\miniconda3\envs\CatchDelin\lib\site-packages\pysheds\sgrid.py in accumulation(self, fdir, weights, dirmap, nodata_out, efficiency, routing, cycle_size, kwargs) 794 efficiency = self._input_handler(efficiency, kwargs) 795 if routing.lower() == 'd8': --> 796 acc = self._d8_accumulation(fdir, weights=weights, dirmap=dirmap, 797 nodata_out=nodata_out, 798 efficiency=efficiency)

~\miniconda3\envs\CatchDelin\lib\site-packages\pysheds\sgrid.py in _d8_accumulation(self, fdir, weights, dirmap, nodata_out, efficiency, **kwargs) 827 eff = np.asarray(eff) 828 # Find indegree of all cells --> 829 indegree = np.bincount(endnodes.ravel(), minlength=fdir.size).astype(np.uint8) 830 # Set starting nodes to those with no predecessors 831 startnodes = startnodes[(indegree == 0)]

<__array_function__ internals> in bincount(*args, **kwargs) ValueError: 'list' argument must have no negative elements

I have tried several different approaches to generate the catchment:

catch = grid.catchment(x=x, y=y, fdir=fdir, xytype='coordinates')

# Clip the view to the catchment
grid.clip_to(catch)

# Compute accumulation
acc = grid.accumulation(fdir, apply_output_mask=False)
# Find the row and column index corresponding to the pour point
col, row = grid.nearest_cell(x, y)

# Delineate a catchment
catch_ind = grid.catchment(x=col, y=row, fdir=fdir, xytype='index')

# Clip the view to the catchment
grid.clip_to(catch_ind)

# Compute accumulation
acc = grid.accumulation(fdir, apply_output_mask=False)
# Compute accumulation
acc = grid.accumulation(fdir)

# Snap pour point to high accumulation cell
x_snap, y_snap = grid.snap_to_mask(acc > 10000, (x, y))

# Delineate the catchment
catch = grid.catchment(x=x_snap, y=y_snap, fdir=fdir, xytype='coordinate')

They all throw the same error. I have also tried:

fluviotect commented 2 years ago

this works:

fdir = grid.flowdir(proj_dem)
acc = grid.accumulation(fdir)

see Issue 174 for related code

fluviotect commented 2 years ago

I reckon the problem lies with pour-point snapping, but I'm not sure how to resolve in a manner that doesn't make more sense for me to do the delineation in QGIS or Arc in the first place.

The "catchment" is bizarre on a number of levels:

image

mdbartos commented 2 years ago

I don't have enough info to give an answer. Feel free to post the dataset.

fluviotect commented 2 years ago

Cool. What's your preference? upload externally and post a link?

mdbartos commented 2 years ago

Yes, external upload would be preferred.