<jupyter_start><jupyter_text>anuga to landlab: example 2In this example we use [Anuga](https://github.com/GeoscienceAustralia/anuga_core) to create a mesh with variable resolution. With Anuga, we are able to specify regions and specific grid cell resolutions within those regions. In this example we will demonstrate this by doing the following:1. Creation of an *Anuga* mesh with 2 regions of different resolutions2. Create simple ridge that is in the center of the domain and each slope has a different cell resolution3. Translate the *Anuga* mesh to a *landlab* grid4. Run the *landlab* flow accumulation routine5. Fill in elevation holes that were created in the mesh generation6. Re-run the flow accumulation routine  Import librariesThis notebook requires the *numpy*, *matplotlib*, *landlab* and *anuga* packages. See the [documentation](https://elbeejay.github.io/meshing-with-landlab/dependencies/index.html) for links to the source installation guides.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import anuga
import anuga.utilities.idwinterp as idt<jupyter_output><empty_output><jupyter_text>Generate the Anuga meshTo generate the Anuga mesh we first define a bounding polygon that defines the outer extents. Then a region within the domain is defined `c_region` by listing points in counter-clockwise order that form a closed shape. The resolution (cell area) is defined for a low resolution and higher resolution region. Then the mesh is visualized so this distinction in resolution can be observed.<jupyter_code># define domain extents
bounding_polygon = [[0.0, 0.0],
                    [100.0, 0.0],
                    [100.0, 100.0],
                    [0.0, 100.0]]
# define a region of the domain to model in higher-res
c_region = [[0.0, 50.0], [100.0, 50.0], [100.0, 100.0], [0.0, 100.0]]
# define boundaries
boundary_tags = {'bottom': [0,1],
                'right': [1,2],
                'top': [2,3],
                'left': [3,0]}
# define low and high res values
low_res = 50
high_res = 1
# create the domain
domain = anuga.create_domain_from_regions(bounding_polygon,
                                          boundary_tags,
                                          maximum_triangle_area=low_res,
                                          interior_regions=[[c_region, high_res]],
                                          mesh_filename='test.msh')
# visualize mesh
dplotter = anuga.Domain_plotter(domain)
plt.figure(figsize=(6,6))
plt.triplot(dplotter.triang, linewidth=0.2)
plt.axis('equal')
plt.show()<jupyter_output>Figure files for each frame will be stored in _plot<jupyter_text>Generate fake topography and assign to the meshHere we generate a simple sloped topography. This is defined using *numpy*, so it is on a regular grid. Built-in interpolation functionality from *anuga* is applied to interpolate from the regular grid to the *anuga* mesh. Lastly the elevation values are assigned to the mesh.<jupyter_code># create fake topography as a grid
topo_gridded = np.zeros((100, 100))
[a, b] = np.shape(topo_gridded)
for i in range(0, b):
    if i < b/2:
        topo_gridded[i, :] = i
    else:
        topo_gridded[i, :] = np.abs(i - 100)
        
grid_y, grid_x = np.mgrid[0:100, 0:100]

plt.figure()
plt.imshow(topo_gridded)
plt.colorbar()
plt.show()
# flatten
xval = grid_x.flatten()
yval = grid_y.flatten()
topoval = topo_gridded.flatten()
# interpolate gridded values to irregular grid
idwtree = idt.invdisttree(np.transpose((xval, yval)), topoval)
topo = idwtree(domain.centroid_coordinates, nnear=3, eps=0, p=1, weights=None)
# set values
domain.set_quantity('elevation', topo, location='centroids')<jupyter_output><empty_output><jupyter_text>Translate to landlabTo translate the mesh as well as the elevation values from *anuga* to *landlab*, first we need to identify the vertex coordinates and their corresponding elevation values. After that the `landlab.VoronoiDelaunayGrid` object is established using the coordinate values. This grid is plotted to ensure that the variable resolution is preserved. Then the elevation attribute is initialized for each node, and the translation of *anuga* elevation values to the *landlab* grid takes place. The resulting *landlab* grid with elevations is plotted next to the source elevation data to visually check them.<jupyter_code># get vertex coordinates and their elevation values
e = domain.get_quantity('elevation')
X, Y, A, V = e.get_vertex_values()
XY = np.column_stack((X,Y))
from landlab.grid import VoronoiDelaunayGrid
from landlab.plot.graph import plot_graph
worked = 0
while worked <= 25:
    try:
        grid = VoronoiDelaunayGrid(X,Y)
        worked = 100
    except Exception:
        worked += 1
        
print(worked)
plt.figure(figsize=(5, 5))
plt.title("Landlab's plot_graph(grid, at='link', with_id=False)")
plot_graph(grid, at="link", with_id=False)
z_vals = grid.add_zeros("topographic__elevation", at="node")
# translating elevation from anuga grid to landlab grid
for i in range(0, len(grid.node_x)):
    ind = np.where(XY==[grid.node_x[i],grid.node_y[i]])[0][0]
    grid.at_node['topographic__elevation'][i] = A[ind]
from landlab.plot.imshow import imshow_grid
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.title('Landlab Grid')
imshow_grid(grid, 'topographic__elevation', show_elements=True, cmap='viridis')
plt.subplot(1,2,2)
plt.title('Source Topo Data')
plt.imshow(topo_gridded)
plt.gca().invert_yaxis()
plt.colorbar()
plt.show()<jupyter_output>/home/jayh/miniconda3/envs/espin/lib/python3.8/site-packages/landlab/plot/imshow.py:267: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
  cmap.set_bad(color=color_for_closed)<jupyter_text>Running the flow accumulatorNow the `landlab.FlowAccumulator` module is initialized and run. Using this on the variable resolution grid shows the impact grid resolution can have on model results (even on a simple sloped topography). The final cell plots the source elevation data on the left, the generated landlab grid with elevations in the center, and the resulting flow accumulation map on the right.<jupyter_code>from landlab.components import FlowAccumulator
from landlab.components import FlowDirectorSteepest
fa = FlowAccumulator(grid, 'topographic__elevation',
                     flow_director=FlowDirectorSteepest)
fa.run_one_step()
imshow_grid(grid, 'drainage_area', show_elements=False, cmap='viridis')
from landlab.plot.imshow import imshow_grid
plt.figure(figsize=(15,5))
plt.subplot(1,3,1)
plt.title('Source Topo Data')
plt.imshow(topo_gridded)
plt.gca().invert_yaxis()
plt.colorbar()
plt.subplot(1,3,2)
plt.title('Landlab Grid')
imshow_grid(grid, 'topographic__elevation', show_elements=True, cmap='viridis')
plt.subplot(1,3,3)
plt.title('Landlab Flow Accumulation')
imshow_grid(grid, 'drainage_area', show_elements=False, cmap='viridis')
plt.show()<jupyter_output><empty_output><jupyter_text>Hole fillingIf we look closely at the center panel in the figure above, we might notice that some cells near the ridge have low elevation values. This is likely due to the conversion process between the mesh types, leaving some unintended 'holes' in the topography. While this is not ideal, it is not too different from the 'holes' that are commonly found in DEMs. As such, there are built-in *landlab* components for dealing with 'holes'. We will first visualize the domain to see the holes, then apply the `landlab.LakeMapperBarnes` component to fill in the holes, and take a look at the resulting grid.<jupyter_code>imshow_grid(grid, "topographic__elevation", show_elements=False, cmap='viridis')
plt.show()
from landlab.components import LakeMapperBarnes
z_init = z_vals.copy()
lmb = LakeMapperBarnes(grid, method='Steepest', surface=z_init, fill_flat=True,
                       redirect_flow_steepest_descent=False, track_lakes=False, ignore_overfill=True)
lmb.run_one_step()
imshow_grid(grid, "topographic__elevation", show_elements=False, cmap='viridis')
plt.show()<jupyter_output><empty_output><jupyter_text>Re-running the FlowAccumulatorIt is worth noting that the holes predominantly appeared in the 'coarse' region of the mesh, and therefore could be due to large irregular cells whose values were interpolated from a finer regular grid. Anyway, now that the holes are filled, we can re-run the `landlab.FlowAccumulator` component and take a look at how grid resolution impacts our results.<jupyter_code>fa = FlowAccumulator(grid, 'topographic__elevation',
                     flow_director=FlowDirectorSteepest)
fa.run_one_step()
imshow_grid(grid, 'drainage_area', show_elements=False, cmap='viridis')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09274284426414243, 'fraction_numerical': 0.015100292990759523, 'mean_word_length': 4.6781829814459375, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 27, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '204202', 'license': 'MIT', 'n_tokens_mistral': 2861, 'n_tokens_neox': 2673, 'n_words': 915, 'path': 'examples/anuga_landlab_02.ipynb', 'repo_name': 'elbeejay/meshing-with-landlab', 'size': '8874'}