Commit 19c8a630 authored by Guillaume Cordonnier's avatar Guillaume Cordonnier

fix bug in boundary conditions for receivers

parent 5e4ae44b
......@@ -15,6 +15,7 @@ def reset_receivers(receivers, nnodes):
@numba.njit
def compute_receivers_d8(receivers, dist2receivers, elevation,
active_nodes,
nx, ny, dx, dy):
# queen (D8) neighbor lookup
dr = np.array([-1, -1, -1, 0, 0, 1, 1, 1], dtype=np.intp)
......@@ -24,19 +25,24 @@ def compute_receivers_d8(receivers, dist2receivers, elevation,
tiny = np.finfo(elevation.dtype).tiny
for r in range(1, ny - 1):
for c in range(1, nx - 1):
for r in range(0, ny):
for c in range(0, nx):
inode = r * nx + c
if not active_nodes[inode]:
continue
slope_max = tiny
for k in range(8):
ineighbor = (r + dr[k]) * nx + (c + dc[k])
slope = (elevation[inode] - elevation[ineighbor]) / length[k]
if r + dr[k] >=0 and r + dr[k]<ny and c + dc[k] >=0 and c + dc[k] <nx:
ineighbor = (r + dr[k]) * nx + (c + dc[k])
slope = (elevation[inode] - elevation[ineighbor]) / length[k]
if slope > slope_max:
slope_max = slope
receivers[inode] = ineighbor
dist2receivers[inode] = length[k]
if slope > slope_max:
slope_max = slope
receivers[inode] = ineighbor
dist2receivers[inode] = length[k]
@numba.njit
......@@ -455,7 +461,7 @@ def _sediment_lakes(sills, basin_stack, basin_parent, basins, water_height, elev
continue
#find sill
i_sill = sills[b, 0] if elevation[sills[b, 0]] >= elevation[sills[b, 1]] else sills[b, 1]
i_sill = sills[b, 0] if water_height[sills[b, 0]] >= water_height[sills[b, 1]] else sills[b, 1]
parse[parse_begin] = i_sill
parse_end += 1
......@@ -463,9 +469,11 @@ def _sediment_lakes(sills, basin_stack, basin_parent, basins, water_height, elev
parsed[i_sill] = -1
while (parse_begin != parse_end):
i_node = parse[parse_begin]
parse_begin +=1
for nb_dir_i in range (8):
......@@ -551,7 +559,7 @@ def correct_flowrouting(receivers, dist2receivers, ndonors, donors,
basin_stack_size = _orient_basin_tree(conn_basins, conn_nodes, nbasins, basin0, mstree,
sills, basin_stack,sediment_lakes, elevation, basin_parent)
if sediment_lakes:
......@@ -569,7 +577,7 @@ def correct_flowrouting(receivers, dist2receivers, ndonors, donors,
if reorder:
reset_receivers(receivers, nnodes)
compute_receivers_d8(receivers, dist2receivers,
water_height, nx, ny, dx, dy)
water_height, active_nodes, nx, ny, dx, dy)
else:
_update_pits_receivers(receivers, dist2receivers, outlets,
......@@ -579,6 +587,7 @@ def correct_flowrouting(receivers, dist2receivers, ndonors, donors,
if reorder or not sediment_lakes:
compute_donors(ndonors, donors, receivers, nnodes)
compute_stack(stack, ndonors, donors, receivers, nnodes)
def compute_water_level(elevation, nx, ny, dx, dy, water_height, active_nodes = None, method='mst_linear', sediment_slope = 1e-10):
......@@ -612,7 +621,7 @@ def compute_water_level(elevation, nx, ny, dx, dy, water_height, active_nodes =
reset_receivers(receivers, nnodes)
compute_receivers_d8(receivers, dist2receivers,
elevation, nx, ny, dx, dy)
elevation, active_nodes, nx, ny, dx, dy)
compute_donors(ndonors, donors, receivers, nnodes)
compute_stack(stack, ndonors, donors, receivers, nnodes)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment