mdbartos / pysheds

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

Set the threshold of the catchment area? #225

Closed staceyzhang-cbcl closed 1 year ago

staceyzhang-cbcl commented 1 year ago

I wonder if there is an option to set the area threshold of the catchment?

I am working on delineating catchments/subcatchments of a large region. With a list of outlets along a river, the catchments delineated become larger as the outlets approach the main channel.

For example, the target catchments look like this. This was created using pcswmm's watershed delineation tool by setting a certain level of target discretization. image

But when I set the same pour point, the area is rather large. image

mdbartos commented 1 year ago

It is possible to limit the catchment size. One could add a counter to this function that checks for the number of total cells added:

https://github.com/mdbartos/pysheds/blob/32fb2efba94c5da4ee8e2556d484a218d352226b/pysheds/_sgrid.py#L316

Is that what you want to do though? How is the target level of discretization defined?

staceyzhang-cbcl commented 1 year ago

I took a look at the function, but not sure where to add the counter and when do use it. Is it before grid.catchment? But it does seem like what I want to do.

The catchment discretization seems to have something to do with the landscape heterogeneity, and therefore generates distributes output variables (Dehotin and Braud, 2008). What I was told is that this could limit the size of catchments to be delineated. So I think this function might be what I'm looking for. How do I use it?

staceyzhang-cbcl commented 1 year ago

Hi!

It turned out that I wasn't using the right way to put together the polygons into one shapefile. I should have sorted the polygons in a descending order, so that smaller polygons wouldn't be covered by larger ones. Thank you Matt for answering my question!

image

mdbartos commented 1 year ago

Thanks. If you want to create non-overlapping polygons, there is a code example in this issue.

(Note that this example uses the old convention of providing a data name instead of inputting a raster object to the function.)