Error while doing Watershed analysis using Python

Question

I am trying to do watershed analysis using Python in jupyter notebook. Whenever I am passing commands regarding watershed delineation in Python I am getting output as raster matrix and when I am plotting it in graph using imshow() command I am getting squared images in graph and not the desired output. I am passing commands such as:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import geopandas as gpd
from pysheds.grid import Grid
import mplleaflet
%matplotlib inline
import tifffile as tiff

grid = Grid.from_raster(r'D:\...\ASTGTMV003_N19E073_dem.tif', data_name='dem')
grid.view('dem')

output :

Raster([[202., 204., 204., ..., 517., 524., 532.],
        [204., 202., 202., ..., 514., 519., 523.],
        [202., 201., 200., ..., 519., 516., 519.],
        ...,
        [  0.,   0.,   0., ..., 606., 606., 604.],
        [  0.,   0.,   0., ..., 607., 607., 603.],
        [  0.,   0.,   0., ..., 607., 612., 610.]], dtype=float32)

... And then while calculating catchment area:

Flow dir:

grid.flowdir(data='inflated_dem', out_name='dir', dirmap=dirmap)
grid.view('dir')

Output:

Raster([[  0.,   0.,   0., ...,   0.,   0.,   0.],
        [  0.,   2.,   2., ...,  -1.,   4.,   0.],
        [  0.,   4.,   1., ...,  64.,  -1.,   0.],
        ...,
        [  0.,   4.,   4., ...,   1., 128.,   0.],
        [  0.,  -1.,  -1., ..., 128.,   1.,   0.],
        [  0.,   0.,   0., ...,   0.,   0.,   0.]], dtype=float32)

But then I plotted it in image using:

plotFigure(grid.dir,'Flow Direction','viridis')

And I got image but while calculating catchment area By specifying x and y:

x, y = 73.6, 19.6
grid.catchment(data='dir', x=x, y=y, dirmap=dirmap, out_name='catch',
               recursionlimit=15000, xytype='label')
grid.clip_to('catch')
catch = grid.view('catch')
plt.imshow(catch)

I am getting output as

Why is this occurring?

Answer 1

So this is a fairly common error associated with pysheds that to my knowledge has not been fixed (I moved away from pysheds to TopoToolbox for exactly this error). There is an open issue on this still on github. It sounds like one solution is to add to the catchment function the following:

grid.catchment(data='dir', x=x, y=y, dirmap=dirmap, out_name='catch',
               recursionlimit=15000, xytype='label', snap = 'center')

The default is to snap the x,y pour point to the corner of a pixel, which is likely causing the odd watershed delineation. The default:

def nearest_cell(self, x, y, affine=None, snap='corner'): 

My experience with this has been that sometime it fixes it, sometimes it doesn't. I have not been able to figure out the exactly reasoning for it occurring, I think it has something to do with coordinate systems. Hopefully this helps.