How to use the landlab.CLOSED_BOUNDARY function in landlab

def test_fixed_gradient_and_value_boundary():
    """testing multiple boundary conditions with fixed links"""
    grid = RasterModelGrid((3, 4))
    grid["node"]["topographic__elevation"] = np.zeros(grid.number_of_nodes)
    grid["link"]["topographic__slope"] = np.zeros(grid.number_of_links)

    grid.set_closed_boundaries_at_grid_edges(False, True, False, True)
    grid.status_at_node[4] = FG
    grid.status_at_node[7] = FG
    grid.status_at_node[1] = 1

    assert_array_equal(
        grid.status_at_node, [CB, FV, CB, CB, FG, 0, 0, FG, CB, CB, CB, CB]
    )

    assert_array_equal(
        grid.status_at_link, [4, 4, 4, 4, 0, 4, 4, 2, 0, 2, 4, 4, 4, 4, 4, 4, 4]
    )

def test_set_status_with_array():
    """Test that active links are reset after changing the node status."""
    grid = RasterModelGrid((4, 5))

    assert_array_equal(
        grid.active_links,
        [5, 6, 7, 9, 10, 11, 12, 14, 15, 16, 18, 19, 20, 21, 23, 24, 25],
    )

    grid.status_at_node[:5] = CB
    assert_array_equal(
        grid.status_at_node,
        [CB, CB, CB, CB, CB, FV, 0, 0, 0, FV, FV, 0, 0, 0, FV, FV, FV, FV, FV, FV],
    )

    assert_array_equal(
        grid.active_links, [9, 10, 11, 12, 14, 15, 16, 18, 19, 20, 21, 23, 24, 25]
    )

def test_multiple_status_node(four_by_four_raster):
    four_by_four_raster.set_closed_boundaries_at_grid_edges(True, True, False, False)
    vals = constant(
        four_by_four_raster,
        "values",
        "node",
        where=[CORE_NODE, CLOSED_BOUNDARY],
        value=10.0,
    )
    true_array = np.array(
        [
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
            10.0,
            10.0,
            10.0,
            0.0,
            10.0,
            10.0,
            10.0,

from six.moves import range
import numpy as np
from matplotlib.pyplot import show, plot, figure
from landlab import RasterModelGrid, CLOSED_BOUNDARY, imshow_grid_at_node
from landlab.components import FlowRouter, SedDepEroder

mg = RasterModelGrid((30, 3), 200.)  # ((10, 3), 200.)
for edge in (mg.nodes_at_left_edge, mg.nodes_at_top_edge,
             mg.nodes_at_right_edge):
    mg.status_at_node[edge] = CLOSED_BOUNDARY

z = mg.add_zeros('node', 'topographic__elevation')
th = mg.add_zeros('node', 'channel_sediment__depth')
th += 0.0007

fr = FlowRouter(mg)
sde = SedDepEroder(mg, K_sp=1.e-5,
                   sed_dependency_type='almost_parabolic',
                   Qc='power_law', K_t=1.e-5)

z[:] = mg.node_y/10000.

initz = z.copy()

dt = 100.
up = 0.0005

dem_grid.set_nodata_nodes_to_inactive(z, 0)  # set nodata nodes to inactive bounds
outlet_node = dem_grid.grid_coords_to_node_id(outlet_row, outlet_column)

grid = RasterModelGrid(4, 5, 1.0)
grid.set_inactive_boundaries(False, True, True, True)
z = grid.add_zeros("node", "Land_Surface__Elevation")
z[6] = 4.5
z[7] = 3.
z[8] = 1.
z[11] = 4.
z[12] = 2.8
z[13] = 2.


# Get array of interior (active) node IDs
interior_nodes = np.where(grid.status_at_node != CLOSED_BOUNDARY)[0]

# Route flow
flow_router = FlowAccumulator(grid, flow_director="D8")
grid = flow_router.route_flow()

for i in range(grid.number_of_nodes):
    print(
        i,
        grid.node_x[i],
        grid.node_y[i],
        z[i],
        grid.status_at_node[i],
        r[i],
        a[i],
        q[i],
        ss[i],

zoutlet=z[outlet_loc]

    #some helper and parameter values
    total_run_time = 500000 #yr
    uplift_rate = 0.001 #m/yr
    rain_rate = 1 #m/yr
    storm_duration = 50 #years

    elapsed_time = 0 #years

    #set up fluvial incision component
    incisor = PowerLawIncision('input_file.txt',rg)

    #print "elevations before ", rg.node_vector_to_raster(z)
    #interior_nodes are the nodes on which you will be operating
    interior_nodes = np.where(rg.status_at_node != CLOSED_BOUNDARY)[0]

    while elapsed_time < total_run_time:
        #uplift the landscape
        z[interior_nodes] = z[interior_nodes]+uplift_rate * storm_duration
        z[outlet_loc]=zoutlet
        #erode the landscape
        z = incisor.run_one_storm(rg,z,rain_rate,storm_duration)
        #update the time
        elapsed_time = elapsed_time+storm_duration

    #below purely for plotting reasons
    #instantiate variable of type RouteFlowD8 Class
    flow_router = RouteFlowD8(len(z))
    #initial flow direction
    flowdirs, max_slopes = flow_router.calc_flowdirs(rg,z)
    #insantiate variable of type AccumFlow Class

#some helper and parameter values
    total_run_time = 10000 #yr
    one_twentieth_time = total_run_time/20
    uplift_rate = 0.001 #m/yr
    rain_rate = 1 #m/yr
    storm_duration = 50 #years

    elapsed_time = 0 #years

    #set up fluvial incision component
    incisor = PowerLawIncision('input_file.txt',rg)

    #print "elevations before ", rg.node_vector_to_raster(z)
    #interior_nodes are the nodes on which you will be operating
    interior_nodes = np.where(rg.status_at_node != CLOSED_BOUNDARY)[0]

    while elapsed_time < total_run_time:
        #erode the landscape
        z = incisor.run_one_storm(rg,z,rain_rate,storm_duration)
        #uplift the landscape
        z[interior_nodes] = z[interior_nodes]+uplift_rate * storm_duration
        #update the time
        elapsed_time = elapsed_time+storm_duration
        #print "elapsed_time", elapsed_time
        if elapsed_time%one_twentieth_time == 0:
            print("elapsed time",elapsed_time)
            elev_raster = rg.node_vector_to_raster(z,True)
            # Plot topography
            pylab.figure(22)
            im = pylab.imshow(elev_raster, cmap=pylab.cm.RdBu, extent=[0, nc*rg.dx, 0, nr*rg.dx])
            cb = pylab.colorbar(im)

shape = (mg.number_of_nodes, 1)

        plt.quiver(mg.x_of_node.reshape(shape), mg.y_of_node.reshape(shape),
               xdist.reshape(shape), ydist.reshape(shape),
               color=colors,
               angles='xy',
               scale_units='xy',
               scale=1,
               zorder=3)

    # Plot differen types of nodes:
    o, = plt.plot(mg.x_of_node[mg.status_at_node == CORE_NODE], mg.y_of_node[mg.status_at_node == CORE_NODE], 'b.', label='Core Nodes', zorder=4)
    fg, = plt.plot(mg.x_of_node[mg.status_at_node == FIXED_VALUE_BOUNDARY], mg.y_of_node[mg.status_at_node == FIXED_VALUE_BOUNDARY], 'c.', label='Fixed Gradient Nodes', zorder=5)
    fv, = plt.plot(mg.x_of_node[mg.status_at_node == FIXED_GRADIENT_BOUNDARY], mg.y_of_node[mg.status_at_node ==FIXED_GRADIENT_BOUNDARY], 'g.', label='Fixed Value Nodes', zorder=6)
    c, = plt.plot(mg.x_of_node[mg.status_at_node == CLOSED_BOUNDARY], mg.y_of_node[mg.status_at_node ==CLOSED_BOUNDARY], 'r.', label='Closed Nodes', zorder=7)

    node_id = np.arange(mg.number_of_nodes)
    flow_to_self = receivers[:,0]==node_id

    fts, = plt.plot(mg.x_of_node[flow_to_self], mg.y_of_node[flow_to_self], 'kx', markersize=6, label = 'Flows To Self', zorder=8)

    ax = plt.gca()

    ax.legend(labels = ['Core Nodes', 'Fixed Gradient Nodes', 'Fixed Value Nodes', 'Closed Nodes', 'Flows To Self'],
              handles = [o, fg, fv, c, fts], numpoints=1, loc='center left', bbox_to_anchor=(1.7, 0.5))
    sm = plt.cm.ScalarMappable(cmap=propColor, norm=plt.Normalize(vmin=0, vmax=1))
    sm._A = []
    cx = plt.colorbar(sm)
    cx.set_label('Proportion of Flow')
    plt.title(title)

h = self.h
        g = self.g
        n = self.n
        tau_crit = self.tau_crit
        q = self.q
        qs = self.qs
        tau = self.tau
        rhog = self.rhog
        alpha = self.alpha
        mpm = self.mpm
        zm = self.zm
        dtmax = self.dtmax
        erode_start_time = self.erode_start_time
        ten_thirds = 10./3.

        interior_cells = np.where(grid.status_at_node != CLOSED_BOUNDARY)[0]


        # Calculate the effective flow depth at active links. Bates et al. 2010
        # recommend using the difference between the highest water-surface
        # and the highest bed elevation between each pair of cells.
        zmax = grid.max_of_link_end_node_values(z)
        w = h+z   # water-surface height
        wmax = grid.max_of_link_end_node_values(w)
        hflow = wmax - zmax

        # Calculate water-surface slopes
        water_surface_slope = grid.calc_grad_of_active_link(w)

        # Calculate the unit discharges (Bates et al., eq 11)
        q = (q-g*hflow*dtmax*water_surface_slope)/ \
            (1.+g*hflow*dtmax*n*n*abs(q)/(hflow**ten_thirds))

elevs = return_array_at_node(grid, elevs)

    ### Step 1, some basic set-up, gathering information about the grid.

    # Calculate the number of nodes.
    num_nodes = len(elevs)

    # Set the number of receivers and facets.
    num_receivers = 2
    num_facets = 8

    # Create a node array
    node_id = np.arange(num_nodes)

    # find where there are closed nodes.
    closed_nodes = grid.status_at_node == CLOSED_BOUNDARY

    # create an array of the triangle numbers
    tri_numbers = np.arange(num_facets)

    ### Step 3, create some triangle datastructures because landlab (smartly)
    # makes it hard to deal with diagonals.

    # create list of triangle neighbors at node. Use orientation associated
    # with tarboton's 1997 algorithm, orthogonal link first, then diagonal.
    # has shape, (nnodes, 8 triangles, 2 neighbors)
    n_at_node = grid.adjacent_nodes_at_node
    dn_at_node = grid.diagonal_adjacent_nodes_at_node

    triangle_neighbors_at_node = np.stack(
        [
            np.vstack((n_at_node[:, 0], dn_at_node[:, 0])),