4

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')

enter image description here

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 enter image description here

Why is this occurring?

2
  • If I'm not mistaken, the extent of your ASTER DTM tile 'ASTGTMV003_N19E073_dem.tif' is from 19 N to 18 N. The coordinates in the file name are the NW corner. But you chose a y value of 19.6. That would be outside (north) of the tile. – Micha Jan 19 at 9:08
  • Actually I am getting 'y' extent from 19.0-20.0(19.0,19.2,19.4,19.6,19.8,20.0)and ' x' extent from 73.0-74.0 (73.0,73.2,73.4,73.6,73.8,74.0) – Omkar Kadlag Jan 19 at 9:42
0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.