|
| 1 | +""" |
| 2 | + Plot a surface + points (clasts) on it |
| 3 | +
|
| 4 | + Use matplotlib Axes3D.plot_surface and Axes3D.scatter3D |
| 5 | +
|
| 6 | + Axes3D.plot_surface(X, Y, Z, *args, **kwargs) |
| 7 | + X, Y and Z are 2D arrays of similar shape. |
| 8 | + In our case, the plot is clearer if we don't plot boundary nodes |
| 9 | + (but that could bean option). So the shape is that of the LL grid |
| 10 | + core nodes. |
| 11 | +
|
| 12 | + X = dx * [[1, 2, 3, ..., nb_of_columns-1] |
| 13 | + [1, 2, 3, ..., nb_of_columns-1] |
| 14 | + ... |
| 15 | + [1, 2, 3, ..., nb_of_columns-1]] |
| 16 | +
|
| 17 | + Y = dy * [[1, 1, 1, ..., 1] |
| 18 | + [2, 2, 2, ..., 2] |
| 19 | + ... |
| 20 | + [nb_of_rows-1, ..., nb_of_rows-1]] |
| 21 | +
|
| 22 | + Z = [[z_core_node1, z_core_node2, ...] |
| 23 | + [..., ... ] |
| 24 | + ... |
| 25 | + [..., ] |
| 26 | +
|
| 27 | + The rstride and cstride kwargs set the stride used to sample the input data |
| 28 | + to generate the graph. If 1k by 1k arrays are passed in, the default values |
| 29 | + for the strides will result in a 100x100 grid being plotted. Defaults to 10. |
| 30 | +
|
| 31 | + The kwargs alpha sets the transparency factor (0 to 1). |
| 32 | +
|
| 33 | + The LL nodes are at the crossing of the white lines that form the surface |
| 34 | + (the surface is made of LL patches, not cells). |
| 35 | +
|
| 36 | +
|
| 37 | + Axes3D.scatter(xs, ys, zs=0, zdir='z', s=20, c=None, |
| 38 | + depthshade=True, *args, **kwargs) |
| 39 | + xs, ys, zs are 1D arrays defining the position of the points. |
| 40 | + s = size (in points, so relation to the space scales of the plot itself is |
| 41 | + not straightforward...) |
| 42 | + c = color |
| 43 | + depthshade = Whether or not to shade the scatter markers to give the |
| 44 | + appearance of depth. Default is True. |
| 45 | +""" |
| 46 | + |
| 47 | +from landlab import RasterModelGrid |
| 48 | +import numpy as np |
| 49 | + |
| 50 | +import matplotlib.pyplot as plt |
| 51 | +from mpl_toolkits.mplot3d import Axes3D |
| 52 | + |
| 53 | + |
| 54 | +### Create Raster Model Grid |
| 55 | +rows = 20 |
| 56 | +columns = 25 |
| 57 | +dx = 1 |
| 58 | +dy = 2 |
| 59 | + |
| 60 | +mg = RasterModelGrid((rows, columns), spacing=(dy, dx)) |
| 61 | + |
| 62 | +# Add elevation field |
| 63 | +z = mg.node_y*0.1 |
| 64 | +_ = mg.add_field('node', 'topographic__elevation', z) |
| 65 | + |
| 66 | +# Shape of core nodes: |
| 67 | +# to be modified for more complex grids (use a mask?) |
| 68 | +number_of_core_node_rows = mg.number_of_node_rows - 2 |
| 69 | +number_of_core_node_columns = mg.number_of_node_columns - 2 |
| 70 | + |
| 71 | +##################################################################### |
| 72 | + |
| 73 | +### Data for 3D plot of topographic surface |
| 74 | +xplot = mg.node_x[mg.core_nodes].reshape(( |
| 75 | + number_of_core_node_rows, number_of_core_node_columns)) |
| 76 | +yplot = mg.node_y[mg.core_nodes].reshape(( |
| 77 | + number_of_core_node_rows, number_of_core_node_columns)) |
| 78 | + |
| 79 | +# Elevation of core nodes, reshaped for 3D plot: |
| 80 | +zplot = mg.at_node['topographic__elevation'][mg.core_nodes].reshape( |
| 81 | + (number_of_core_node_rows, number_of_core_node_columns)) |
| 82 | + |
| 83 | +##################################################################### |
| 84 | + |
| 85 | +### 3D plot of elevation surface: |
| 86 | +# Figure and type of projection: |
| 87 | +fig = plt.figure(1) |
| 88 | +ax = plt.axes(projection='3d') |
| 89 | + |
| 90 | +# Plot surface: |
| 91 | +ax.plot_surface(xplot, yplot, zplot, cmap='binary', rstride=1, cstride=1, |
| 92 | + alpha=0.5) |
| 93 | + |
| 94 | +# Set initial view of the graph (elevation and azimuth): |
| 95 | +ax.view_init(elev=None, azim=-130) |
| 96 | + |
| 97 | +ax.set_xlabel('X axis') |
| 98 | +ax.set_ylabel('Y axis') |
| 99 | +ax.set_zlabel('Z axis') |
| 100 | + |
| 101 | +##################################################################### |
| 102 | + |
| 103 | +### Data for 3D plot of topographic surface |
| 104 | +# The Clast Set class (from Clast Tracker) provides |
| 105 | +# clast node ID, x and y coordinates, clast elevation, etc. |
| 106 | +# but for this example, we define them here: |
| 107 | + |
| 108 | +clast__node = np.array([382, 386, 386, 390, 392, 392, 392, 392]) |
| 109 | +clast__number = clast__node.size |
| 110 | + |
| 111 | +clast__x = mg.node_x[clast__node] |
| 112 | +clast__y= mg.node_y[clast__node] |
| 113 | + |
| 114 | +clast__elevation = np.array([3, 3, 3, 3, 3, 3, 3, 3]) |
| 115 | +clast__size = np.array([0.5, 0.5, 0.5, 1, 1, 1, 1, 1]) |
| 116 | + |
| 117 | +# For display purpose, clasts sizes are increased: |
| 118 | +# Marker size is in "points"... |
| 119 | +sizes = clast__size * 100 |
| 120 | + |
| 121 | +# Count the number of clast on each node (WILL PUT THIS IN THE CLAST TRACKER) |
| 122 | +clast__number_at_node = np.zeros(mg.number_of_nodes) |
| 123 | +for i in range(0, mg.number_of_nodes): |
| 124 | + clast__number_at_node[i] = list(clast__node).count(mg.nodes.reshape( |
| 125 | + clast__number_at_node.shape)[i]) |
| 126 | + |
| 127 | +# Assign colors to clast markers as a function of clast density on the node: |
| 128 | +clast__color = np.zeros(clast__number) |
| 129 | +for j in range(0, clast__number): |
| 130 | + clast__color[j] = clast__number_at_node[clast__node[j]] |
| 131 | + |
| 132 | +##################################################################### |
| 133 | + |
| 134 | +### 3D plot of points (scatter): |
| 135 | + |
| 136 | +plot = ax.scatter(clast__x, clast__y, clast__elevation, marker='H', |
| 137 | + s=sizes, c=clast__color, label=clast__color, cmap='cool') |
| 138 | +cbar = plt.colorbar(plot, ticks=[1, 2, 3, 4], shrink=0.7) |
| 139 | +cbar.set_label('# of clasts') |
0 commit comments