-
Notifications
You must be signed in to change notification settings - Fork 71
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Added catchment delineation functions and exposed to Python #53
base: master
Are you sure you want to change the base?
Conversation
@joeperdefloep Nice work! This would be a great addition to the package. |
ooh,quick reply! some other important notes I forgot to mention previously:
EDIT: so this is caused by the line extending to the edge of the grid... If I make a grid of 6 and then draw a line until y1=4, there is no problem, but an only 0 line. I think this may need some more work...
|
reference for better line-drawing, including slopes >1: https://www.geeksforgeeks.org/bresenhams-line-generation-algorithm/ |
c45532c
to
e96a122
Compare
8289e25
to
c51ebaf
Compare
Ok... So I found out that my commit was 45+ commits behind the main branch, which made some things obsolete that I wrote before, and gave some headaches while fixing. It all works now! import richdem as rd
import numpy as np
dem = rd.rdarray(np.array([
[1,2,3,4,5],
[1,2,3,4,5],
[1,2,3,4,5],
[1,2,3,4,5],
[1,2,3,4,5],]),no_data=0)
mask = rd.rdarray(np.array([
[0,0,0,0,0],
[0,0,0,0,0],
[1,0,0,0,0],
[0,0,0,0,0],
[0,0,0,0,0]]),no_data=0)
print(rd.UpslopeCells(dem,mask=mask,method="Quinn"))
print(rd.UpslopeCells(dem,mask=mask,method="mflow")) output:
but also a lot of these when I start richdem:
and a segfault at the end.
full output:
I think this has something to do with my installation, rather than the PR
|
@giswqs Since it didn't get replies or merged quickly, I also implemented it in Cython, albeit less generic: %%cython -+
from libcpp.queue cimport queue
from libcpp.vector cimport vector
from libc.math cimport isnan
cdef struct xypoint:
long x
long y
# 3 2 1
# 4 X 0
# 5 6 7
cdef int[8] dx = [1, 1, 0, -1, -1, -1, 0, 1]
cdef int[8] dy = [0, 1, 1, 1, 0, -1, -1, -1]
def upslope_cells(double[:,:] mask, double[:,:] dem) -> void:
cdef queue[xypoint] qu
for i in range(mask.shape[0]):
for j in range(mask.shape[1]):
if mask[i,j] != 0:
qu.push(xypoint(i,j))
while not qu.empty():
point = qu.front()
qu.pop()
for i in range(8):
x = point.x + dx[i]
y = point.y + dy[i]
if x < 0 or x > mask.shape[0]:
continue
if y < 0 or y > mask.shape[1]:
continue
if mask[x,y] == 0:
if dem[x,y] > dem[point.x,point.y]:
qu.push(xypoint(x,y))
mask[x,y] += 1 hope it helps! |
Added catchment delineation, because I needed it.
edit: Actually, this PR finds all upslope cells, given an input, using multiple possible methods. This is different than catchment delineation, because two catchments can actually overlap if a divergent flow metric is used. That functionality is actually already present in a modified algorithm for depression-filling
I saw that it was a todo, so took the liberty to implement it.
I needed the functionality myself to also use a mask in stead of a line.
some notes, in no particular order:
catchment_delineation.hpp
functions. This is actually important, because the raster should be filled with nodata values in order for the algorithm to work properly. I also do this there, but not super nice yet, because I wrote a 2-liner helper function that only works for 3/4 cases0
cells, to give more user flexibility,flowProportions
method was fast enough for what I needed. Can add on request. (Also I do not see an easy way to implement Dinf without copying a lot of code)FlowProportions
, which includes a redundantcopyFromWrapped
, but this seemed easier/more understandable than re-writing.It could be that some documentation for c++ is not fully optimal yet, but I can work on that. For now I just wanted to make this pull request. EDIT: documentation should be rather good.
hope you like it!