Source code for landlab.components.lateral_erosion.node_finder

import numpy as np


[docs] def angle_finder(grid, dn, cn, rn): """Find the interior angle between two vectors on a grid. Parameters ---------- grid : ModelGrid A landlab grid. dn : int or array of int Node or nodes at the end of the first vector. cn : int or array of int Node or nodes at the vertex between vectors. rn : int or array of int Node or nodes at the end of the second vector. Returns ------- float or array of float Angle between vectors (in radians). Examples -------- >>> import numpy as np >>> from landlab import RasterModelGrid >>> from landlab.components.lateral_erosion.node_finder import angle_finder >>> grid = RasterModelGrid((3, 4)) >>> np.rad2deg(angle_finder(grid, 8, 5, 0)) 90.0 >>> np.rad2deg(angle_finder(grid, (8, 9, 10, 6), 5, 6)) array([ 135., 90., 45., 0.]) """ vertex = np.take(grid.x_of_node, cn), np.take(grid.y_of_node, cn) vec_1 = [ np.take(grid.x_of_node, dn) - vertex[0], np.take(grid.y_of_node, dn) - vertex[1], ] vec_2 = [ np.take(grid.x_of_node, rn) - vertex[0], np.take(grid.y_of_node, rn) - vertex[1], ] return np.arccos( (vec_1[0] * vec_2[0] + vec_1[1] * vec_2[1]) / np.sqrt((vec_1[0] ** 2 + vec_1[1] ** 2) * (vec_2[0] ** 2 + vec_2[1] ** 2)) )
[docs] def forty_five_node(donor, i, receiver, neighbors, diag_neigh): radcurv_angle = 0.67 lat_node = 0 # In Landlab 2019: diagonal list goes [NE, NW, SW, SE]. Node list are ordered as [E,N,W,S] # if water flows SE-N OR if flow NE-S or E-NW or E-SW, erode west node if ( donor == diag_neigh[0] and receiver == neighbors[3] or donor == diag_neigh[3] and receiver == neighbors[1] or donor == neighbors[0] and receiver == diag_neigh[2] or donor == neighbors[0] and receiver == diag_neigh[1] ): lat_node = neighbors[2] # if flow is from SW-N or NW-S or W-NE or W-SE, erode east node elif ( donor == diag_neigh[1] and receiver == neighbors[3] or donor == diag_neigh[2] and receiver == neighbors[1] or donor == neighbors[2] and receiver == diag_neigh[3] or donor == neighbors[2] and receiver == diag_neigh[0] ): lat_node = neighbors[0] # if flow is from SE-W or SW-E or S-NE or S-NW, erode north node elif ( donor == diag_neigh[3] and receiver == neighbors[2] or donor == diag_neigh[2] and receiver == neighbors[0] or donor == neighbors[3] and receiver == diag_neigh[0] or donor == neighbors[3] and receiver == diag_neigh[1] ): lat_node = neighbors[1] # if flow is from NE-W OR NW-E or N-SE or N-SW, erode south node elif ( donor == diag_neigh[0] and receiver == neighbors[2] or donor == diag_neigh[1] and receiver == neighbors[0] or donor == neighbors[1] and receiver == diag_neigh[3] or donor == neighbors[1] and receiver == diag_neigh[2] ): lat_node = neighbors[3] return lat_node, radcurv_angle
[docs] def ninety_node(donor, i, receiver, link_list, neighbors, diag_neigh): # if flow is 90 degrees if donor in diag_neigh and receiver in diag_neigh: radcurv_angle = 1.37 # if flow is NE-SE or NW-SW, erode south node if ( donor == diag_neigh[0] and receiver == diag_neigh[3] or donor == diag_neigh[1] and receiver == diag_neigh[2] ): lat_node = neighbors[3] # if flow is SW-NW or SE-NE, erode north node elif ( donor == diag_neigh[2] and receiver == diag_neigh[1] or donor == diag_neigh[3] and receiver == diag_neigh[0] ): lat_node = neighbors[1] # if flow is SW-SE or NW-NE, erode east node elif ( donor == diag_neigh[2] and receiver == diag_neigh[3] or donor == diag_neigh[1] and receiver == diag_neigh[0] ): lat_node = neighbors[0] # if flow is SE-SW or NE-NW, erode west node elif ( donor == diag_neigh[3] and receiver == diag_neigh[2] or donor == diag_neigh[0] and receiver == diag_neigh[1] ): lat_node = neighbors[2] elif donor not in diag_neigh and receiver not in diag_neigh: radcurv_angle = 1.37 # if flow is from east, erode west node if donor == neighbors[0]: lat_node = neighbors[2] # if flow is from north, erode south node elif donor == neighbors[1]: lat_node = neighbors[3] # if flow is from west, erode east node elif donor == neighbors[2]: lat_node = neighbors[0] # if flow is from south, erode north node elif donor == neighbors[3]: lat_node = neighbors[1] return lat_node, radcurv_angle
[docs] def straight_node(donor, i, receiver, neighbors, diag_neigh): # ***FLOW LINK IS STRAIGHT, NORTH TO SOUTH***# if donor == neighbors[1] or donor == neighbors[3]: # print "flow is stright, N-S from ", donor, " to ", flowdirs[i] radcurv_angle = 0.23 # neighbors are ordered E,N,W, S # if the west cell is boundary (neighbors=-1), erode from east node if neighbors[2] == -1: lat_node = neighbors[0] elif neighbors[0] == -1: lat_node = neighbors[2] else: # if could go either way, choose randomly. 0 goes East, 1 goes west ran_num = np.random.randint(0, 2) if ran_num == 0: lat_node = neighbors[0] if ran_num == 1: lat_node = neighbors[2] # ***FLOW LINK IS STRAIGHT, EAST-WEST**# elif donor == neighbors[0] or donor == neighbors[2]: radcurv_angle = 0.23 # Node list are ordered as [E,N,W,S] # if the north cell is boundary (neighbors=-1), erode from south node if neighbors[1] == -1: lat_node = neighbors[3] elif neighbors[3] == -1: lat_node = neighbors[1] else: # if could go either way, choose randomly. 0 goes south, 1 goes north ran_num = np.random.randint(0, 2) if ran_num == 0: lat_node = neighbors[1] if ran_num == 1: lat_node = neighbors[3] # if flow is straight across diagonal, choose node to erode at random elif donor in diag_neigh and receiver in diag_neigh: radcurv_angle = 0.23 if receiver == diag_neigh[0]: poss_diag_nodes = neighbors[0 : 1 + 1] elif receiver == diag_neigh[1]: poss_diag_nodes = neighbors[1 : 2 + 1] elif receiver == diag_neigh[2]: poss_diag_nodes = neighbors[2 : 3 + 1] elif receiver == diag_neigh[3]: poss_diag_nodes = [neighbors[3], neighbors[0]] ran_num = np.random.randint(0, 2) if ran_num == 0: lat_node = poss_diag_nodes[0] if ran_num == 1: lat_node = poss_diag_nodes[1] return lat_node, radcurv_angle
[docs] def node_finder(grid, i, flowdirs, drain_area): """Find lateral neighbor node of the primary node for straight, 45 degree, and 90 degree channel segments. Parameters ---------- grid : ModelGrid A Landlab grid object i : int node ID of primary node flowdirs : array Flow direction array drain_area : array drainage area array Returns ------- lat_node : int node ID of lateral node radcurv_angle : float inverse radius of curvature of channel at lateral node """ # receiver node of flow is flowdirs[i] receiver = flowdirs[i] # find indicies of where flowdirs=i to find donor nodes. # will donor nodes always equal the index of flowdir list? inflow = np.where(flowdirs == i) # if there are more than 1 donors, find the one with largest drainage area if len(inflow[0]) > 1: drin = drain_area[inflow] drmax = max(drin) maxinfl = inflow[0][np.where(drin == drmax)] # if donor nodes have same drainage area, choose one randomly if len(maxinfl) > 1: ran_num = np.random.randint(0, len(maxinfl)) maxinfln = maxinfl[ran_num] donor = [maxinfln] else: donor = maxinfl # if inflow is empty, no donor elif len(inflow[0]) == 0: donor = i # else donor is the only inflow else: donor = inflow[0] # now we have chosen donor cell, next figure out if inflow/outflow lines are # straight, 45, or 90 degree angle. and figure out which node to erode link_list = grid.links_at_node[i] # this gives list of active neighbors for specified node # the order of this list is: [E,N,W,S] neighbors = grid.active_adjacent_nodes_at_node[i] # this gives list of all diagonal neighbors for specified node # the order of this list is: [NE,NW,SW,SE] diag_neigh = grid.diagonal_adjacent_nodes_at_node[i] angle_diff = np.rad2deg(angle_finder(grid, donor, i, receiver)) if (donor == flowdirs[i]) or (donor == i): # this is a sink. no lateral ero radcurv_angle = 0.0 lat_node = 0 elif np.isclose(angle_diff, 0.0) or np.isclose(angle_diff, 180.0): [lat_node, radcurv_angle] = straight_node( donor, i, receiver, neighbors, diag_neigh ) elif np.isclose(angle_diff, 45.0) or np.isclose(angle_diff, 135.0): [lat_node, radcurv_angle] = forty_five_node( donor, i, receiver, neighbors, diag_neigh ) elif np.isclose(angle_diff, 90.0): [lat_node, radcurv_angle] = ninety_node( donor, i, receiver, link_list, neighbors, diag_neigh ) else: lat_node = 0 radcurv_angle = 0.0 dx = grid.dx # INVERSE radius of curvature. radcurv_angle = radcurv_angle / dx return int(lat_node), radcurv_angle