How to use landlab - 10 common examples

runtime = inputs.read_float("total_time")
    dt = inputs.read_float("dt")

    nt = int(runtime // dt)
    uplift_per_step = uplift_rate * dt

    mg = RasterModelGrid((nrows, ncols), xy_spacing=(dx, dx))

    mg.add_zeros("topographic__elevation", at="node")
    z = np.loadtxt(os.path.join(_THIS_DIR, "seddepinit.txt"))
    mg["node"]["topographic__elevation"] = z

    mg.set_closed_boundaries_at_grid_edges(True, False, True, False)

    fr = FlowAccumulator(mg, flow_director="D8")
    sde = SedDepEroder(mg, **inputs)

    for i in range(nt):
        mg.at_node["topographic__elevation"][mg.core_nodes] += uplift_per_step
        mg = fr.run_one_step()
        mg, _ = sde.erode(dt)

    z_tg = np.loadtxt(os.path.join(_THIS_DIR, "seddepz_tg.txt"))

    assert_array_almost_equal(
        mg.at_node["topographic__elevation"][mg.core_nodes], z_tg[mg.core_nodes]
    )

Returns
    -------
    ndarray of int
        The link ID as the nth vertical links.

    Examples
    --------
    >>> from landlab.grid.structured_quad.links import nth_vertical_link
    >>> shape = (3, 4)
    >>> nth_vertical_link(shape, 4)
    1
    >>> nth_vertical_link(shape, (3, 4, 11))
    array([0, 1, 5])
    """
    links = np.asarray(links, dtype=np.int)
    return as_id_array(
        (links // (2 * shape[1] - 1)) * shape[1]
        + links % (2 * shape[1] - 1)
        - (shape[1] - 1)
    )

Returns
    -------
    ndarray of int
        The link ID as the nth horizontal links.

    Examples
    --------
    >>> from landlab.components.overland_flow._links import nth_horizontal_link
    >>> shape = (3, 4)
    >>> nth_horizontal_link(shape, 16)
    8
    >>> nth_horizontal_link(shape, (1, 7, 8))
    array([1, 3, 4])
    """
    links = np.asarray(links, dtype=np.int)
    return as_id_array(
        (links // (2 * shape[1] - 1)) * (shape[1] - 1) + links % (2 * shape[1] - 1)
    )

profile_IDs,
            mg.at_node["flow__link_to_receiver_node"],
        )
        prf.plot_profiles(
            dists_upstr, profile_IDs, mg.at_node["topographic__elevation"]
        )

print("completed run to steady state...")
if show_figs_in_run:
    show()  # will cause a hang

# save a copy of the init conditions:
mg_init = deepcopy(mg)

# REinstantiate the components:
fr = FlowAccumulator(mg, flow_director="D8")
tle = TransportLimitedEroder(mg, input_file)
uplift_rate *= 10.  # accelerate tenfold
runtime = 200000.
nt = int(runtime // dt)
uplift_per_step = uplift_rate * dt
print("uplift per step: {0}".format(uplift_per_step))

x_profiles = []
z_profiles = []
S_profiles = []
A_profiles = []

# plot init conds
if make_output_plots:
    mg = fr.route_flow(grid=mg)
    pylab.figure("long_profile_anim")

im = pylab.plot(mg.dx * np.arange(nrows), elev_r[:, int(ncols // 2)])

    print('Completed loop ', i)

print('Completed the simulation. Plotting...')

# Finalize and plot:
# put a title on figure 4
pylab.figure(4)
pylab.title('N-S cross_section, linear diffusion')
pylab.xlabel('Distance')
pylab.ylabel('Elevation')

# figure 5 is the map of the final elevations
pylab.figure(5)
imshow_node_grid(mg, 'topographic__elevation')

# superpose this final form onto figure 3:
pylab.figure(3)
# turn the 1-D array of elevation values into a spatially accurate 2-D
# gridded format, for plotting
elev_r = mg.node_vector_to_raster(mg['node']['topographic__elevation'])
im = pylab.plot(mg.dx * np.arange(nrows), elev_r[:, int(ncols // 2)])
pylab.xlabel('Distance')
pylab.ylabel('Elevation')

pylab.show()  # this line displays all of the figures you've issued plot commands for, since you last called show()

def test_field_name_array_float_case6():
    """Topography as array, runoff rate as array"""
    mg = RasterModelGrid((5, 4), xy_spacing=(1, 1))
    topographic__elevation = np.array(
        [
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
            21.0,
            10.0,
            0.0,
            0.0,
            31.0,
            20.0,
            0.0,
            0.0,
            32.0,

def test_D_infinity_flat_closed_lower():
    mg = RasterModelGrid((5, 4), xy_spacing=(1, 1))
    z = mg.add_zeros("node", "topographic__elevation")
    z[mg.core_nodes] += 1
    mg.set_closed_boundaries_at_grid_edges(
        bottom_is_closed=True,
        left_is_closed=True,
        right_is_closed=True,
        top_is_closed=True,
    )

    fd = FlowDirectorDINF(mg)
    fd.run_one_step()

    node_ids = np.arange(mg.number_of_nodes)
    true_recievers = -1 * np.ones(fd.receivers.shape)
    true_recievers[:, 0] = node_ids

def test_netcdf_write_uint8_field_netcdf4(tmpdir):
    """Test write_netcdf with a grid that has an uint8 field."""
    field = RasterModelGrid((4, 3))
    field.add_field("node", "topographic__elevation", np.arange(12, dtype=np.uint8))

    with tmpdir.as_cwd():
        write_netcdf("test.nc", field, format="NETCDF4")

        root = nc.Dataset("test.nc", "r", format="NETCDF4")

        for name in ["topographic__elevation"]:
            assert name in root.variables
            assert_array_equal(root.variables[name][:].flatten(), field.at_node[name])
            assert root.variables[name][:].dtype == "uint8"

        root.close()

def test_providing_array_and_kwargs():
    mg = RasterModelGrid((10, 10))
    mg.add_zeros("node", "topographic__elevation")
    noise = np.random.rand(mg.size("node"))
    mg.at_node["topographic__elevation"] += noise
    fr = FlowAccumulator(mg, flow_director="D8")
    fr.run_one_step()
    mask = np.zeros(len(mg.at_node["topographic__elevation"]), dtype=np.uint8)
    mask[np.where(mg.at_node["drainage_area"] > 5)] = 1
    with pytest.warns(UserWarning):
        DrainageDensity(mg, channel__mask=mask, area_coefficient=1)
    with pytest.warns(UserWarning):
        DrainageDensity(mg, channel__mask=mask, slope_coefficient=1)
    with pytest.warns(UserWarning):
        DrainageDensity(mg, channel__mask=mask, area_exponent=1)
    with pytest.warns(UserWarning):
        DrainageDensity(mg, channel__mask=mask, slope_exponent=1)
    with pytest.warns(UserWarning):

def test_diagonal_list_with_scalar_arg():
    rmg = RasterModelGrid((5, 4))

    assert_array_equal(rmg.diagonal_adjacent_nodes_at_node[6], np.array([11, 9, 1, 3]))
    assert_array_equal(rmg.diagonal_adjacent_nodes_at_node[-1], np.array([X, X, 14, X]))
    assert_array_equal(
        rmg.diagonal_adjacent_nodes_at_node[-2], np.array([X, X, 13, 15])
    )