giswqs
commented
2 years ago @Joeperdefloep Nice work! This would be a great addition to the package.
Open feefladder opened 2 years ago
giswqs
commented
2 years ago @Joeperdefloep Nice work! This would be a great addition to the package.
feefladder
commented
2 years ago ooh,quick reply!
some other important notes I forgot to mention previously:
d8_methods.hpp. output with x0=0,y0=0,x1=4,y0=4:
2 2 0 0 0
0 2 2 0 0
0 0 2 2 0
0 0 0 2 2
0 0 0 0 2
103 100 105 46 104 <-- where do these come from?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...
d8_methods.hpp included in catchment_delineation_generic.hpp only for sgn() function. maybe move to some common file?catch_delineation, because I had a no such file or directory: "catch_delineation" error and then just wrote a simple one by handTODO notes where I had some questions inside the code.feefladder
commented
2 years ago reference for better line-drawing, including slopes >1: https://www.geeksforgeeks.org/bresenhams-line-generation-algorithm/
feefladder
commented
2 years ago 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! sample code:
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:
Warning! No geotransform defined. Choosing a standard one! (Top left cell's top let corner at <0,0>; cells are 1x1.)
Warning! No geotransform defined. Choosing a standard one! (Top left cell's top let corner at <0,0>; cells are 1x1.)
Warning! No geotransform defined. Choosing a standard one! (Top left cell's top let corner at <0,0>; cells are 1x1.)
Warning! No geotransform defined. Choosing a standard one! (Top left cell's top let corner at <0,0>; cells are 1x1.)
Warning! No geotransform defined. Choosing a standard one! (Top left cell's top let corner at <0,0>; cells are 1x1.)
[[0 0 0 0 0]
[0 1 1 1 0]
[2 1 1 1 0]
[0 1 1 1 0]
[0 0 0 0 0]]
Warning! No geotransform defined. Choosing a standard one! (Top left cell's top let corner at <0,0>; cells are 1x1.)
Warning! No geotransform defined. Choosing a standard one! (Top left cell's top let corner at <0,0>; cells are 1x1.)
Warning! No geotransform defined. Choosing a standard one! (Top left cell's top let corner at <0,0>; cells are 1x1.)
[[0 0 1 1 1]
[0 1 1 1 1]
[2 1 1 1 1]
[0 1 1 1 1]
[0 0 1 1 1]]but also a lot of these when I start richdem:
HDF5-DIAG: Error detected in HDF5 (1.10.7) thread 1:
#000: H5T.c line 1929 in H5Tcopy(): not a datatype or dataset
major: Invalid arguments to routine
minor: Inappropriate typeand a segfault at the end.
I think this has something to do with my installation, rather than the PR
feefladder
commented
2 years ago @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] += 1hope 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.hppfunctions. 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 cases0cells, to give more user flexibility,flowProportionsmethod 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!