accept raster / array of pour points as input to catchment

[mdbartos]

accept raster / array of pour points as input to catchment

[mdbartos]

Hi @cfandel, thanks for your question. You can do this using a for loop. If your points are not topologically ordered, you can do something like this:

# Import modules
import numpy as np
import matplotlib.pyplot as plt
from pysheds.grid import Grid
import seaborn as sns
np.random.seed(0)

# Load flow direction grid
grid = Grid.from_raster('../data/n30w100_dir', data_name='dir')

# Generate the small sample catchment
x, y = -97.294167, 32.73750
grid.catchment(data='dir', x=x, y=y, out_name='catch',
               recursionlimit=15000, xytype='label', nodata_out=0)
grid.clip_to('catch')

# Get some example pour points
grid.accumulation(data='dir', out_name='acc', apply_mask=False)
branches = grid.extract_river_network('catch', 'acc')
points = np.vstack([geom['geometry']['coordinates'][-1] for geom in branches['features']])
points = points[np.random.choice(len(points), size=10, replace=False), :]
# Add outlet point
points = np.vstack([points, np.asarray([x, y])])

# Delineate catchments
labels = np.zeros(grid.shape, dtype=int)
catchments = []
for point in points:
    x, y = point
    catchment = grid.catchment(x=x, y=y, data='dir', xytype='label', inplace=False)
    catchments.append(catchment)
    labels += (catchment != 0)

# Find non-overlapping catchments
for index, catchment in enumerate(catchments):
    mask = (catchment > 0)
    label = labels[mask].min()
    catchments[index] = np.where(mask & (labels == label), catchment, 0)

labeled_catchments = sum([np.where(catchment, index + 1, 0) 
                          for index, catchment in enumerate(catchments)]).astype(float)
labeled_catchments[labeled_catchments == 0] = np.nan

# Plot results
fig, ax = plt.subplots(figsize=(8,6))
plt.imshow(labeled_catchments, cmap='cubehelix', zorder=1, extent=grid.extent)
plt.colorbar()
plt.scatter(points[:, 0], points[:, 1], zorder=2, c='r')

[mdbartos]

Do you mean that you want the output of the catchment function to be the contributing area of site A + the contributing area of site B + contributing area of site C, etc.?

[smnorris]

Pretty much. Ideally I want catchment output to be contributing area of all non-zeros or no-nodata in an input 'pour point' raster. This reduces need for cleaning DEM if I have good stream data already.

See in_pour_point_data parameter here

For a raster, this represents cells above which the contributing area, or catchment, will be determined. 

[smnorris]

Turns out that just adding one line fixes my issue for all test sites.

# Set mask to blurred mask
mask = blurred_mask
# Create a view onto the DEM array
dem = grid.view('dem', dtype=np.float64, nodata=np.nan)
# Subtract dz where mask is nonzero
dem[(mask > 0)] -= dz*mask[(mask > 0)]
# Smooth the DEM in masked area - trying to force the DEM to go downhill
# along the stream
dem[(mask > 0)] = ndimage.filters.minimum_filter(dem, 2)[(mask > 0)]

[mdbartos]

You've tried grid.fill_depressions followed by grid.resolve_flats right? That should force everything to route to the boundaries of the dataset.

Edit: nevermind, I read the original post again.

[smnorris]

This is my code, I can provide the DEM & streams if you want. The minimum filter seemed to help a bit.
https://gist.github.com/ca7e6e9fc9511b1096b7d105a69d365b

[mdbartos]

Hi @smnorris

If you want to send the data to me by email (see profile), I can take a look at it. Just let me know what problem(s) you are still having and need to resolve.

Thanks,
MDB

[smnorris]

Here are a handful of DEMs where I'm having the problem, plus the stream segments I'm burning in (and would have used as pour point rasters with the ArcGIS method)
https://hillcrestgeo.ca/outgoing/public

This shows what I'm doing, using the first DEM:
https://gist.github.com/39efa810e510acf418268c2438abeb68

Here is a view of the inputs and outputs in QGIS (point is the pour point, input rasterized stream segment in red). The streams shown on the bottom half of the image all flow to the pour point.

[cfandel]

Hi,

I actually would like the opposite functionality. I have a bunch of points, and I'd like to be able to pass lists of x and y coordinates for those points to grid.catchment, but get the catchment for each point, not for all the points together - similar to what the ArcMap Spatial Analyst Watershed tool does if you give it multiple pour points. The lower points should have catchments that exclude areas already captured by the higher points. Is this possible?

Should look something like this:

(notebook where I am trying this out can be found here: https://github.com/cfandel/gottesacker/blob/master/subwatersheds.ipynb)

[mdbartos]

Because there are several possible interpretations for what result an array of pour points should generate, I would prefer not to overload the inputs to the catchment function. However, I am interested in creating another function that is analogous to the watershed functions from scikit-image or scipy.ndimage:

http://scikit-image.org/docs/dev/auto_examples/segmentation/plot_watershed.html
http://scipy-lectures.org/advanced/image_processing/auto_examples/plot_watershed_segmentation.html

In other words, supply an array of labels and delineate the contributing area for each labeled region.

[cfandel]

This is perfect thank you! I was getting stuck on how to get the subcatchments to not overlap. Do you mind if I take the code chunk above and modify/incorporate it into what I'm working on?

[mdbartos]

Sure! Go ahead. There is probably a more memory-efficient method if you can guarantee that the points are sorted (i.e. a point occurs later in the list if it is downstream of another point).

[cfandel]

They're definitely not going to be sorted, unfortunately (they're karst conduit inlet nodes), but this should be efficient enough. Thanks so much!