|
11 | 11 | "cell_type": "markdown", |
12 | 12 | "metadata": {}, |
13 | 13 | "source": [ |
14 | | - "# Introduction to the RockBlock and LayeredRockBlock objects\n", |
| 14 | + "# Introduction to the Lithology and LithoLayers objects\n", |
15 | 15 | "\n", |
16 | | - "The RockBlock and LayeredRockBlock are two landlab objects meant to make it easier to work with spatially variable lithology that produces spatially variable parameter values (e.g. stream power erodability or diffusivity). \n", |
| 16 | + "Lithology and LithoLayers are two landlab components meant to make it easier to work with spatially variable lithology that produces spatially variable parameter values (e.g. stream power erodability or diffusivity). \n", |
17 | 17 | "\n", |
18 | | - "In this tutorial we will first use the LayeredRockBlock to erode either dipping layeres or an anticline. Then we will use RockBlock to create inverted topography. \n", |
| 18 | + "In this tutorial we will first use the LithoLayers to erode either dipping layeres or an anticline. Then we will use Lithology to create inverted topography. \n", |
19 | 19 | "\n", |
20 | 20 | "To start, we will import the necessary modules. A note: this tutorial uses the [HoloViews package](http://holoviews.org) for visualization. This package is a great tool for dealing with multidimentional annotated data (e.g. an xarray dataset). If you get an error on import, consider updating dask (this is what the author needed to do in April 2018). You will also need to have the [Bokeh](https://bokeh.pydata.org/en/latest/) and [Matplotlib](https://matplotlib.org) packages installed." |
21 | 21 | ] |
|
33 | 33 | "%matplotlib inline\n", |
34 | 34 | "\n", |
35 | 35 | "import holoviews as hv\n", |
36 | | - "hv.extension('bokeh', 'matplotlib')\n", |
| 36 | + "hv.notebook_extension('matplotlib')\n", |
37 | 37 | "\n", |
38 | 38 | "from landlab import RasterModelGrid\n", |
39 | | - "from landlab.layers import RockBlock, LayeredRockBlock\n", |
40 | | - "from landlab.components import FlowAccumulator, FastscapeEroder, LinearDiffuser\n", |
| 39 | + "from landlab.components import FlowAccumulator, FastscapeEroder, LinearDiffuser, Lithology, LithoLayers\n", |
41 | 40 | "from landlab.plot import imshow_grid" |
42 | 41 | ] |
43 | 42 | }, |
|
47 | 46 | "source": [ |
48 | 47 | "## Part 1: Eroding layered rock\n", |
49 | 48 | "\n", |
50 | | - "First we will create an instance of a LayedRockBlock. Both LayeredRockBlock and RockBlock work closely with a landlab ModelGrid and stores information about rock type at each node. \n", |
| 49 | + "First we will create an instance of a LithoLayers. Both LithoLayers and Lithology work closely with a landlab ModelGrid and stores information about rock type at each node. \n", |
51 | 50 | "\n", |
52 | | - "To create a LayeredRock you need the following information:\n", |
| 51 | + "To create LithoLayers you need the following information:\n", |
53 | 52 | "\n", |
54 | 53 | "1. A dictionary of rock property attributes. \n", |
55 | 54 | "2. A model grid that has the field `'topographic__elevation'` already created.\n", |
|
81 | 80 | "ids = np.tile([0,1], 20) \n", |
82 | 81 | "\n", |
83 | 82 | "# Anticline\n", |
84 | | - "rb = LayeredRockBlock(mg, z0s, ids, x0=7000, y0=7000, function = lambda x, y : ((0.002*x)**2+(0.001*y)**2), attrs=attrs)\n", |
| 83 | + "lith = LithoLayers(mg, z0s, ids, x0=7000, y0=7000, function = lambda x, y : ((0.002*x)**2+(0.001*y)**2), attrs=attrs)\n", |
85 | 84 | "\n", |
86 | 85 | "# Shallow dips\n", |
87 | | - "#rb = LayeredRockBlock(mg, z0s, ids, x0=5000, y0=5000, function = lambda x, y : ((0.001*x)+(0.003*y)), attrs=attrs)\n", |
| 86 | + "#lith = LithoLayers(mg, z0s, ids, x0=5000, y0=5000, function = lambda x, y : ((0.001*x)+(0.003*y)), attrs=attrs)\n", |
88 | 87 | "\n", |
89 | 88 | "# Steeper dips\n", |
90 | | - "#rb = LayeredRockBlock(mg, z0s, ids, x0=5000, y0=5000, function = lambda x, y : ((0.01*x)+(0.01*y)), attrs=attrs)" |
| 89 | + "#lith = LithoLayers(mg, z0s, ids, x0=5000, y0=5000, function = lambda x, y : ((0.01*x)+(0.01*y)), attrs=attrs)" |
91 | 90 | ] |
92 | 91 | }, |
93 | 92 | { |
94 | 93 | "cell_type": "markdown", |
95 | 94 | "metadata": {}, |
96 | 95 | "source": [ |
97 | | - "Now that we've created the LayeredRockBlock, model grid fields for each of the LayeredRockBlock attributes exist and have been set to the values of the rock exposed at the surface. \n", |
| 96 | + "Now that we've created LithoLayers, model grid fields for each of the LithoLayers attributes exist and have been set to the values of the rock exposed at the surface. \n", |
98 | 97 | "\n", |
99 | 98 | "Here we plot the value of `'K_sp'` as a function of the model grid. " |
100 | 99 | ] |
101 | 100 | }, |
102 | 101 | { |
103 | 102 | "cell_type": "code", |
104 | 103 | "execution_count": null, |
105 | | - "metadata": {}, |
| 104 | + "metadata": { |
| 105 | + "collapsed": true |
| 106 | + }, |
106 | 107 | "outputs": [], |
107 | 108 | "source": [ |
108 | 109 | "imshow_grid(mg, 'K_sp')" |
|
116 | 117 | "\n", |
117 | 118 | "Next, lets instantiate a FlowAccumulator and a FastscapeEroder to create a simple landscape evolution model. \n", |
118 | 119 | "\n", |
119 | | - "We will point the FastscapeEroder to the model grid field `'K_sp'` so that it will respond to the spatially variable erodabilities created by the LayeredRockBlock. \n", |
| 120 | + "We will point the FastscapeEroder to the model grid field `'K_sp'` so that it will respond to the spatially variable erodabilities created by LithoLayers. \n", |
120 | 121 | "\n", |
121 | 122 | "We will also instatiate an xarray dataset used to store the output of our model through time for visualization. " |
122 | 123 | ] |
|
150 | 151 | "cell_type": "markdown", |
151 | 152 | "metadata": {}, |
152 | 153 | "source": [ |
153 | | - "Here we run the model. In each time step we first run the FlowAccumulator to direct flow and accumulatate drainage area. Then the FastscapeEroder erodes the topography based on the stream power equation using the erodability value in the field`'K_sp'`. We create an uplift field that uplifts only the model grid's core nodes. After uplifting these core nodes, we update the LayeredRockBlock. Importantly, we must tell the LayeredRockBlock how it has been advected upward by uplift. \n", |
| 154 | + "Here we run the model. In each time step we first run the FlowAccumulator to direct flow and accumulatate drainage area. Then the FastscapeEroder erodes the topography based on the stream power equation using the erodability value in the field`'K_sp'`. We create an uplift field that uplifts only the model grid's core nodes. After uplifting these core nodes, we update LithoLayers. Importantly, we must tell the LithoLayers how it has been advected upward by uplift. \n", |
154 | 155 | "\n", |
155 | | - "`rb.run_one_step` has an optional argument `rock_id` to use when some material may be deposited. Since we are using the FastscapeEroder which is fully detachment limited, we don't need to set this. \n", |
| 156 | + "`lith.run_one_step` has an optional argument `rock_id` to use when some material may be deposited. Since we are using the FastscapeEroder which is fully detachment limited, we don't need to set this. \n", |
156 | 157 | "\n", |
157 | 158 | "Within each timestep we save information about the model for plotting. " |
158 | 159 | ] |
|
171 | 172 | " dz_ad = np.zeros(mg.size('node'))\n", |
172 | 173 | " dz_ad[mg.core_nodes] = U * dt\n", |
173 | 174 | " z += dz_ad\n", |
174 | | - " rb.run_one_step(dz_advection = dz_ad)\n", |
| 175 | + " lith.run_one_step(dz_advection = dz_ad)\n", |
175 | 176 | " \n", |
176 | 177 | " for of in out_fields:\n", |
177 | 178 | " ds[of][i,:,:] = mg['node'][of].reshape(mg.shape)" |
|
187 | 188 | { |
188 | 189 | "cell_type": "code", |
189 | 190 | "execution_count": null, |
190 | | - "metadata": {}, |
| 191 | + "metadata": { |
| 192 | + "collapsed": true |
| 193 | + }, |
191 | 194 | "outputs": [], |
192 | 195 | "source": [ |
193 | 196 | "imshow_grid(mg, 'topographic__elevation', cmap='viridis')" |
|
207 | 210 | { |
208 | 211 | "cell_type": "code", |
209 | 212 | "execution_count": null, |
210 | | - "metadata": {}, |
| 213 | + "metadata": { |
| 214 | + "collapsed": true |
| 215 | + }, |
211 | 216 | "outputs": [], |
212 | 217 | "source": [ |
213 | | - "hvds = hv.Dataset(ds)\n", |
214 | | - "hvds" |
| 218 | + "hvds_topo = hv.Dataset(ds.topographic__elevation)\n", |
| 219 | + "hvds_rock = hv.Dataset(ds.rock_type__id)\n", |
| 220 | + "hvds_topo" |
215 | 221 | ] |
216 | 222 | }, |
217 | 223 | { |
218 | 224 | "cell_type": "markdown", |
219 | 225 | "metadata": {}, |
220 | 226 | "source": [ |
221 | | - "Next we specify that we want two images, one showing rock type and one showing topographic elevation. A slider bar shows us model time in millions of years. " |
| 227 | + "Next we specify that we want two images, one showing rock type and one showing topographic elevation. A slider bar shows us model time in millions of years. \n", |
| 228 | + "\n", |
| 229 | + "Be patient. Running this next block may take a moment. HoloViews is rendering an image of all time slices so you can see an animated slider. This is pretty magical (but not instantaneous). " |
222 | 230 | ] |
223 | 231 | }, |
224 | 232 | { |
225 | 233 | "cell_type": "code", |
226 | 234 | "execution_count": null, |
227 | | - "metadata": {}, |
| 235 | + "metadata": { |
| 236 | + "collapsed": true |
| 237 | + }, |
228 | 238 | "outputs": [], |
229 | 239 | "source": [ |
230 | | - "%opts Image (cmap='viridis')\n", |
231 | | - "topo = hvds.to(hv.Image, ['x', 'y'], ['topographic__elevation'])\n", |
232 | | - "rock = hvds.to(hv.Image, ['x', 'y'], ['rock_type__id'])\n", |
| 240 | + "%opts Image style(interpolation='bilinear', cmap='viridis') plot[colorbar=True]\n", |
| 241 | + "topo = hvds_topo.to(hv.Image, ['x', 'y'])\n", |
| 242 | + "rock = hvds_rock.to(hv.Image, ['x', 'y'])\n", |
233 | 243 | "\n", |
234 | 244 | "topo + rock" |
235 | 245 | ] |
|
245 | 255 | "\n", |
246 | 256 | "## Part 2 Creation of Inverted Topography\n", |
247 | 257 | "\n", |
248 | | - "Here we will explore making inverted topography by eroding a RockBlock with constant properties for half of the model evaluation time, and then filling RockBlock in with resistant material only where the drainage area is large. This is meant as a simple example of filling in valleys with volcanic material. \n", |
| 258 | + "Here we will explore making inverted topography by eroding Lithology with constant properties for half of the model evaluation time, and then filling Lithology in with resistant material only where the drainage area is large. This is meant as a simple example of filling in valleys with volcanic material. \n", |
249 | 259 | "\n", |
250 | | - "All of the details of the options for creating a [RockBlock](http://landlab.readthedocs.io/en/release/landlab.layers.rockblock.html) can be found here. \n", |
| 260 | + "All of the details of the options for creating a [Lithology](http://landlab.readthedocs.io/en/release/landlab.layers.lithology.html) can be found here. \n", |
251 | 261 | "\n", |
252 | | - "In the next code block we make a new model and run it. Note that because we are using the LinearDiffuser which may result in depositing material we must specify a rock_id in the RockBlock `run_one_step` method.\n", |
| 262 | + "In the next code block we make a new model and run it. Note that because we are using the LinearDiffuser which may result in depositing material we must specify a rock_id in the Lithology `run_one_step` method.\n", |
253 | 263 | "\n", |
254 | 264 | "We also are handling the model grid boundary conditions differently than in the last example, setting the boundaries on the left and right sides to closed. " |
255 | 265 | ] |
|
276 | 286 | " 'D': {0: 0.4,\n", |
277 | 287 | " 1: 0.001}}\n", |
278 | 288 | "\n", |
279 | | - "rb2 = RockBlock(mg2, thicknesses2, ids2, attrs=attrs2)\n", |
| 289 | + "lith2 = Lithology(mg2, thicknesses2, ids2, attrs=attrs2)\n", |
280 | 290 | "\n", |
281 | 291 | "nts = 500\n", |
282 | 292 | "U = 0.005\n", |
|
305 | 315 | " dz_ad2 = np.zeros(mg2.size('node'))\n", |
306 | 316 | " dz_ad2[mg2.core_nodes] = U * dt\n", |
307 | 317 | " z2 += dz_ad2\n", |
308 | | - " rb2.run_one_step(dz_advection = dz_ad2, rock_id=0)\n", |
| 318 | + " lith2.run_one_step(dz_advection = dz_ad2, rock_id=0)\n", |
309 | 319 | " \n", |
310 | 320 | " for of in out_fields:\n", |
311 | 321 | " ds2[of][i,:,:] = mg2['node'][of].reshape(mg2.shape)" |
|
321 | 331 | { |
322 | 332 | "cell_type": "code", |
323 | 333 | "execution_count": null, |
324 | | - "metadata": {}, |
| 334 | + "metadata": { |
| 335 | + "collapsed": true |
| 336 | + }, |
325 | 337 | "outputs": [], |
326 | 338 | "source": [ |
327 | | - "imshow_grid(mg2, 'topographic__elevation')" |
| 339 | + "imshow_grid(mg2, 'topographic__elevation', cmap='viridis')" |
328 | 340 | ] |
329 | 341 | }, |
330 | 342 | { |
|
339 | 351 | { |
340 | 352 | "cell_type": "code", |
341 | 353 | "execution_count": null, |
342 | | - "metadata": {}, |
| 354 | + "metadata": { |
| 355 | + "collapsed": true |
| 356 | + }, |
343 | 357 | "outputs": [], |
344 | 358 | "source": [ |
345 | 359 | "volcanic_deposits = np.zeros(mg2.size('node'))\n", |
346 | 360 | "da_big_enough = mg2['node']['drainage_area']>5e4\n", |
347 | 361 | "\n", |
348 | 362 | "topo_difference_from_top = mg2['node']['topographic__elevation'].max() - mg2['node']['topographic__elevation']\n", |
349 | 363 | "\n", |
350 | | - "volcanic_deposits[da_big_enough] = 0.5 * topo_difference_from_top[da_big_enough]\n", |
| 364 | + "volcanic_deposits[da_big_enough] = 0.25 * topo_difference_from_top[da_big_enough]\n", |
351 | 365 | "volcanic_deposits[mg2.boundary_nodes] = 0.0\n", |
352 | 366 | "\n", |
353 | 367 | "z2 += volcanic_deposits\n", |
354 | | - "rb2.run_one_step(rock_id=1)\n", |
| 368 | + "lith2.run_one_step(rock_id=1)\n", |
355 | 369 | "\n", |
356 | 370 | "imshow_grid(mg2, volcanic_deposits)" |
357 | 371 | ] |
|
378 | 392 | " dz_ad2 = np.zeros(mg2.size('node'))\n", |
379 | 393 | " dz_ad2[mg2.core_nodes] = U * dt\n", |
380 | 394 | " z2 += dz_ad2\n", |
381 | | - " rb2.run_one_step(dz_advection = dz_ad2, rock_id=0)\n", |
| 395 | + " lith2.run_one_step(dz_advection = dz_ad2, rock_id=0)\n", |
382 | 396 | " \n", |
383 | 397 | " for of in out_fields:\n", |
384 | 398 | " ds2[of][i,:,:] = mg2['node'][of].reshape(mg2.shape)" |
|
394 | 408 | { |
395 | 409 | "cell_type": "code", |
396 | 410 | "execution_count": null, |
397 | | - "metadata": {}, |
| 411 | + "metadata": { |
| 412 | + "collapsed": true |
| 413 | + }, |
398 | 414 | "outputs": [], |
399 | 415 | "source": [ |
400 | | - "imshow_grid(mg2, 'topographic__elevation')" |
| 416 | + "imshow_grid(mg2, 'topographic__elevation', cmap='viridis')" |
401 | 417 | ] |
402 | 418 | }, |
403 | 419 | { |
|
410 | 426 | { |
411 | 427 | "cell_type": "code", |
412 | 428 | "execution_count": null, |
413 | | - "metadata": {}, |
| 429 | + "metadata": { |
| 430 | + "collapsed": true |
| 431 | + }, |
414 | 432 | "outputs": [], |
415 | 433 | "source": [ |
416 | | - "%opts Image (cmap='viridis')\n", |
| 434 | + "hvds_topo2 = hv.Dataset(ds2.topographic__elevation)\n", |
| 435 | + "hvds_rock2 = hv.Dataset(ds2.rock_type__id)\n", |
417 | 436 | "\n", |
418 | | - "hvds2 = hv.Dataset(ds2)\n", |
419 | | - "topo2 = hvds2.to(hv.Image, ['x', 'y'], ['topographic__elevation'])\n", |
420 | | - "rock2 = hvds2.to(hv.Image, ['x', 'y'], ['rock_type__id'])\n", |
| 437 | + "%opts Image style(interpolation='bilinear', cmap='viridis') plot[colorbar=True]\n", |
| 438 | + "topo2 = hvds_topo2.to(hv.Image, ['x', 'y'])\n", |
| 439 | + "rock2 = hvds_rock2.to(hv.Image, ['x', 'y'])\n", |
421 | 440 | "\n", |
422 | 441 | "topo2 + rock2" |
423 | 442 | ] |
|
437 | 456 | "\n", |
438 | 457 | "Nice work getting to the end of the tutorial!\n", |
439 | 458 | "\n", |
440 | | - "For more detailed information about the [RockBlock](http://landlab.readthedocs.io/en/release/landlab.layers.rockblock.html) and [LayeredRockBlock](http://landlab.readthedocs.io/en/release/landlab.layers.layeredrockblock.html) objects, check out their detailed documentation. \n", |
| 459 | + "For more detailed information about the [Lithology](http://landlab.readthedocs.io/en/release/landlab.layers.lithology.html) and [LithoLayers](http://landlab.readthedocs.io/en/release/landlab.layers.litholayers.html) objects, check out their detailed documentation. \n", |
441 | 460 | "\n", |
442 | 461 | "# **Click [here](https://github.com/landlab/landlab/wiki/Tutorials) for more Landlab tutorials**" |
443 | 462 | ] |
|
0 commit comments