Using Pysheds on non-square Raster

Question

I am using Pysheds for watershed delineation and displaying potential flow networks. Everything works well (like in the GitHub Readme) as long as the raster (TIFF) is square:

But if I try to use a raster clipped to my desired extents, the code does not work anymore:

The occuring error is:

  File "...\Anaconda3\envs\geo_env\lib\site-packages\spyder_kernels\py3compat.py", line 356, in compat_exec
    exec(code, globals, locals)

  File "...\python\createshp.py", line 128, in <module>
    x_snap, y_snap = grid.snap_to_mask(acc > 1000, (x, y))

  File "...\Anaconda3\envs\geo_env\lib\site-packages\pysheds\sgrid.py", line 2394, in snap_to_mask
    return View.snap_to_mask(mask, xy, affine=affine,

  File "...\Anaconda3\envs\geo_env\lib\site-packages\pysheds\sview.py", line 791, in snap_to_mask
    return tree_xy[ix]

IndexError: index 0 is out of bounds for axis 0 with size 0

Why does it say index 0 is out of bounds, when my catchment point is inside my tiff-area. Am I missing something out?

In my code I am using the point read from a SHP-file for the "Delinate a catchment"-step, which works for the square tiff but not for the none-square tiff. The code I am trying to run:

# Loop through points file
wells = gpd.read_file(r'...\points.shp')
for index, row in wells.iterrows():
    # Extracting attriburtes from the shapefile
    num = row['Id']
    y = row['geometry'].y
    x = row['geometry'].x   
    
    
    # Read elevation raster
    grid = Grid.from_raster('...\none-square.tif')
    dem = grid.read_raster('...\none-square.tif')
    
    # Condition DEM
    # ----------------------
    # Fill pits in DEM
    pit_filled_dem = grid.fill_pits(dem)
    
    # Fill depressions in DEM
    flooded_dem = grid.fill_depressions(pit_filled_dem)
        
    # Resolve flats in DEM
    inflated_dem = grid.resolve_flats(flooded_dem)
    
    
    # Determine D8 flow directions from DEM
    # Specify directional mapping
             #N   NE   E  SE S  SW W   NW
    dirmap = (64, 128, 1, 2, 4, 8, 16, 32)
        
    # Compute flow directions
    fdir = grid.flowdir(inflated_dem, dirmap=dirmap)
    
    
    # Calculate flow accumulation
    acc = grid.accumulation(fdir, dirmap=dirmap)
    
    
    # Snap pour point to high accumulation cell
    x_snap, y_snap = grid.snap_to_mask(acc > 1000, (x, y))
    
    # Delineate the catchment
    catch = grid.catchment(x=x_snap, y=y_snap, fdir=fdir, dirmap=dirmap, 
                           xytype='coordinate')
    
    shapes = grid.polygonize()

Answer 1

I suspect it might have to do with the indexes pysheds creates in the background or maybe there are no points in the clipped raster that accumuluate over 1000 (have you checked the values for x_snap and y_snap?)

Could you try to define x_snap and y_snap with the following (This should work, it does on a non rectangualar DEM for me):

Specify pour point (greatest accumulation)

y_snap, x_snap = np.unravel_index(np.argmax(acc), acc.shape)

Delineate a catchment

catch = grid.catchment(x=x_snap , y=y_snap, fdir=fdir, xytype='index')