|
15 | 15 | "\n", |
16 | 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 LithoLayers to erode either dipping layeres or an anticline. Then we will use Lithology to create inverted topography. \n", |
| 18 | + "This tutorial is meant for users who have some experience using Landlab components.\n", |
19 | 19 | "\n", |
20 | | - "We will also use [xarray](https://xarray.pydata.org/en/stable/) to store and annotate our model output. \n", |
| 20 | + "In this tutorial we will explore the creation of layered topography. After an example that will let you see how LithoLayers works, we will work through two more complicated examples. In the first example, we use the LithoLayers to erode either dipping layeres or an anticline. Then we will use Lithology to create inverted topography. \n", |
21 | 21 | "\n", |
22 | | - "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." |
| 22 | + "We will use [xarray](https://xarray.pydata.org/en/stable/) to store and annotate our model output. While we won't extensively discuss the use of xarray, some background will be provided. \n", |
| 23 | + "\n", |
| 24 | + "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.\n", |
| 25 | + "\n", |
| 26 | + "In testing we've seen some users have an error raised related to the matplotlib backend. " |
23 | 27 | ] |
24 | 28 | }, |
25 | 29 | { |
|
76 | 80 | }, |
77 | 81 | "outputs": [], |
78 | 82 | "source": [ |
79 | | - "mg = RasterModelGrid((10, 10), 1)\n", |
| 83 | + "mg = RasterModelGrid((10, 15), 1)\n", |
80 | 84 | "z = mg.add_zeros('node', 'topographic__elevation') " |
81 | 85 | ] |
82 | 86 | }, |
|
148 | 152 | "source": [ |
149 | 153 | "`'K_sp'` is the attribute that we want to track through the layered rock, `0`, `1`, `2`, `3` are the rock type IDs, and `0.0003` and `0.0001` are the values for `'K_sp'` for the rock types `0` and `1`. \n", |
150 | 154 | "\n", |
151 | | - "Finally, we define our function. Here we will use [lambda expression](https://docs.python.org/3/tutorial/controlflow.html#lambda-expressions) to create a small anonymous function. In this case we define a function of `x` and `y` that returns the value `x + (2. * y)`" |
| 155 | + "Finally, we define our function. Here we will use [lambda expression](https://docs.python.org/3/tutorial/controlflow.html#lambda-expressions) to create a small anonymous function. In this case we define a function of `x` and `y` that returns the value `x + (2. * y)`.\n", |
| 156 | + "\n", |
| 157 | + "This means that planar rock layers will dip into the ground to the North-North-East. By changing this functional form, we can make more complicated rock layers." |
152 | 158 | ] |
153 | 159 | }, |
154 | 160 | { |
|
166 | 172 | "cell_type": "markdown", |
167 | 173 | "metadata": {}, |
168 | 174 | "source": [ |
169 | | - "Finally we construct our LithoLayers component." |
| 175 | + "Finally we construct our LithoLayers component by passing the correct arguments." |
170 | 176 | ] |
171 | 177 | }, |
172 | 178 | { |
173 | 179 | "cell_type": "code", |
174 | 180 | "execution_count": null, |
175 | | - "metadata": {}, |
| 181 | + "metadata": { |
| 182 | + "collapsed": true |
| 183 | + }, |
176 | 184 | "outputs": [], |
177 | 185 | "source": [ |
178 | 186 | "lith = LithoLayers(mg, layer_elevations, layer_ids, function=func, attrs=attrs)" |
|
182 | 190 | "cell_type": "markdown", |
183 | 191 | "metadata": {}, |
184 | 192 | "source": [ |
185 | | - "LithoLayers will make sure that the model grid has at-node grid fields with the layer attribute names. In this case, this means that the model grid will now include a grid field called `'K_sp'`. We can plot this with the Landlab [imshow_grid](http://landlab.readthedocs.io/en/release/landlab.plot.html#landlab.plot.imshow.imshow_grid) function. " |
| 193 | + "LithoLayers will make sure that the model grid has at-node grid fields with the layer attribute names. In this case, this means that the model grid will now include a grid field called `'K_sp'` and a field called `'rock_type__id'`. We can plot these with the Landlab [imshow_grid](http://landlab.readthedocs.io/en/release/landlab.plot.html#landlab.plot.imshow.imshow_grid) function. " |
186 | 194 | ] |
187 | 195 | }, |
188 | 196 | { |
|
191 | 199 | "metadata": {}, |
192 | 200 | "outputs": [], |
193 | 201 | "source": [ |
194 | | - "imshow_grid(mg, 'K_sp')" |
| 202 | + "imshow_grid(mg, 'rock_type__id', cmap='viridis')" |
195 | 203 | ] |
196 | 204 | }, |
197 | 205 | { |
198 | 206 | "cell_type": "markdown", |
199 | 207 | "metadata": {}, |
200 | 208 | "source": [ |
| 209 | + "As you can see, we have layers that strike East-South-East. Since we can only see the surface expression of the lyaers, we can't infer the dip direction or magnitude from the plot. \n", |
| 210 | + "\n", |
201 | 211 | "If the topographic surface erodes, then you will want to update LithoLayers. Like most Landlab components, LithoLayers uses a `run_one_step` method to update. \n", |
202 | 212 | "\n", |
203 | | - "Next we will erode the topography by decrementing the variable `z` which point to the topographic elevation of our model grid by 1. If the rock mass is being advected up or down by an external force (e.g. tectonic rock uplift), then then advection must be specified. The `dz_advection` argument can be a single value or an array of size number-of-nodes. " |
| 213 | + "Next we will erode the topography by decrementing the variable `z` which point to the topographic elevation of our model grid by an amount 1. If the rock mass is being advected up or down by an external force (e.g. tectonic rock uplift), then then advection must be specified. The `dz_advection` argument can be a single value or an array of size number-of-nodes. " |
204 | 214 | ] |
205 | 215 | }, |
206 | 216 | { |
|
220 | 230 | "cell_type": "markdown", |
221 | 231 | "metadata": {}, |
222 | 232 | "source": [ |
223 | | - "We can re-plot the value of `'K_sp'`. We will see that the locatio of layers has changed. " |
| 233 | + "We can re-plot the value of `'K_sp'`. We will see that the location of the surface expression of the rock layers has changed. As we expect, the location has changed in a way that is consistent with dipping to the NNE. " |
224 | 234 | ] |
225 | 235 | }, |
226 | 236 | { |
|
229 | 239 | "metadata": {}, |
230 | 240 | "outputs": [], |
231 | 241 | "source": [ |
232 | | - "imshow_grid(mg, 'K_sp')" |
| 242 | + "imshow_grid(mg, 'rock_type__id', cmap='viridis')" |
233 | 243 | ] |
234 | 244 | }, |
235 | 245 | { |
236 | | - "cell_type": "code", |
237 | | - "execution_count": null, |
| 246 | + "cell_type": "markdown", |
238 | 247 | "metadata": {}, |
239 | | - "outputs": [], |
240 | 248 | "source": [ |
241 | | - "ds = lith.to_xarray(np.arange(10))\n", |
242 | | - "hvds_rock = hv.Dataset (ds.rock_type__id)" |
| 249 | + "Next we will create a xarray dataset that has 3D information about our Lithology to help visualize the layers in space. We will use the `rock_cube_to_xarray` method of the LithoLayers component. \n", |
| 250 | + "\n", |
| 251 | + "We will then convert this xarray dataset into a HoloViews dataset so we can visualize the result. " |
243 | 252 | ] |
244 | 253 | }, |
245 | 254 | { |
|
248 | 257 | "metadata": {}, |
249 | 258 | "outputs": [], |
250 | 259 | "source": [ |
251 | | - "%opts Image style(interpolation='bilinear', cmap='viridis') plot[colorbar=True]\n", |
| 260 | + "ds = lith.rock_cube_to_xarray(np.arange(30))\n", |
| 261 | + "hvds_rock = hv.Dataset (ds.rock_type__id)\n", |
| 262 | + "\n", |
| 263 | + "%opts Image style(cmap='viridis') plot[colorbar=True]\n", |
252 | 264 | "hvds_rock.to(hv.Image, ['x', 'y'])" |
253 | 265 | ] |
254 | 266 | }, |
| 267 | + { |
| 268 | + "cell_type": "markdown", |
| 269 | + "metadata": {}, |
| 270 | + "source": [ |
| 271 | + "The slider allows us to change the depth below the topographic surface.\n", |
| 272 | + "\n", |
| 273 | + "We can also plot the cube of rock created with LithoLayers as a cross section. " |
| 274 | + ] |
| 275 | + }, |
255 | 276 | { |
256 | 277 | "cell_type": "code", |
257 | 278 | "execution_count": null, |
| 279 | + "metadata": {}, |
| 280 | + "outputs": [], |
| 281 | + "source": [ |
| 282 | + "%opts Image style(cmap='viridis') plot[colorbar=True, invert_yaxis=True]\n", |
| 283 | + "hvds_rock.to(hv.Image, ['x', 'z'])" |
| 284 | + ] |
| 285 | + }, |
| 286 | + { |
| 287 | + "cell_type": "markdown", |
258 | 288 | "metadata": { |
259 | 289 | "collapsed": true |
260 | 290 | }, |
| 291 | + "source": [ |
| 292 | + "Hopefuly this gives you a sense of how LithoLayers works. The next two blocks of code have all the steps we just worked through in one place. \n", |
| 293 | + "\n", |
| 294 | + "Try modifying the layer thicknesses, the size of the grid, the function used to create the form of the layers, and the location of the anchor point to gain intuition for how you can use LithoLayers to create different types of layered rock. " |
| 295 | + ] |
| 296 | + }, |
| 297 | + { |
| 298 | + "cell_type": "code", |
| 299 | + "execution_count": null, |
| 300 | + "metadata": {}, |
| 301 | + "outputs": [], |
| 302 | + "source": [ |
| 303 | + "# Parameters that control the size and shape of the model grid\n", |
| 304 | + "number_of_rows = 50\n", |
| 305 | + "number_of_columns = 50\n", |
| 306 | + "dx = 1\n", |
| 307 | + "\n", |
| 308 | + "# Parameters that control the LithoLayers\n", |
| 309 | + "\n", |
| 310 | + "# the layer shape function\n", |
| 311 | + "func = lambda x, y : (0.5 * x)**2 + (0.5 * y)**2\n", |
| 312 | + "\n", |
| 313 | + "# the layer thicknesses\n", |
| 314 | + "layer_thickness = 50.\n", |
| 315 | + "\n", |
| 316 | + "# the location of the anchor point\n", |
| 317 | + "x0 = 25\n", |
| 318 | + "y0 = 25\n", |
| 319 | + "\n", |
| 320 | + "# the resolution at which you sample to create the plan view and cros-section view figures.\n", |
| 321 | + "sample_depths = np.arange(0, 30, 1)\n", |
| 322 | + "\n", |
| 323 | + "\n", |
| 324 | + "# create the model grid\n", |
| 325 | + "mg = RasterModelGrid((number_of_rows, number_of_columns), dx)\n", |
| 326 | + "z = mg.add_zeros('node', 'topographic__elevation') \n", |
| 327 | + "\n", |
| 328 | + "# set up LithoLayers inputs\n", |
| 329 | + "layer_ids = np.tile([0,1,2,3], 5) \n", |
| 330 | + "layer_elevations = layer_thickness * np.arange(-10, 10)\n", |
| 331 | + "layer_elevations[-1] = layer_elevations[-2] + 100\n", |
| 332 | + "attrs = {'K_sp': {0: 0.0003,\n", |
| 333 | + " 1: 0.0001,\n", |
| 334 | + " 2: 0.0002,\n", |
| 335 | + " 3: 0.0004}}\n", |
| 336 | + "\n", |
| 337 | + "# create LithoLayers\n", |
| 338 | + "lith = LithoLayers(mg, layer_elevations, layer_ids, x0=x0, y0=y0, function=func, attrs=attrs)\n", |
| 339 | + "\n", |
| 340 | + "# get the rock-cube data structure and plot\n", |
| 341 | + "ds = lith.rock_cube_to_xarray(sample_depths)\n", |
| 342 | + "hvds_rock = hv.Dataset (ds.rock_type__id)\n", |
| 343 | + "\n", |
| 344 | + "# make a plan view image\n", |
| 345 | + "%opts Image style(cmap='viridis') plot[colorbar=True]\n", |
| 346 | + "hvds_rock.to(hv.Image, ['x', 'y'])" |
| 347 | + ] |
| 348 | + }, |
| 349 | + { |
| 350 | + "cell_type": "markdown", |
| 351 | + "metadata": {}, |
| 352 | + "source": [ |
| 353 | + "You can also make a cross section of this new LithoLayers component. " |
| 354 | + ] |
| 355 | + }, |
| 356 | + { |
| 357 | + "cell_type": "code", |
| 358 | + "execution_count": null, |
| 359 | + "metadata": {}, |
261 | 360 | "outputs": [], |
262 | | - "source": [] |
| 361 | + "source": [ |
| 362 | + "%opts Image style(cmap='viridis') plot[colorbar=True, invert_yaxis=True]\n", |
| 363 | + "hvds_rock.to(hv.Image, ['x', 'z'])" |
| 364 | + ] |
263 | 365 | }, |
264 | 366 | { |
265 | 367 | "cell_type": "markdown", |
|
268 | 370 | "## Part 2: Creation of a landscape evolution model with LithoLayers\n", |
269 | 371 | "\n", |
270 | 372 | "\n", |
271 | | - "Introductory text.\n" |
| 373 | + "In this next section, we will run LithoLayers with components used for a simple Landscape Evolution Model. \n", |
| 374 | + "\n", |
| 375 | + "We will start by creating the grid." |
272 | 376 | ] |
273 | 377 | }, |
274 | 378 | { |
|
282 | 386 | "mg = RasterModelGrid((100, 60), 200)\n", |
283 | 387 | "z = mg.add_zeros('node', 'topographic__elevation') \n", |
284 | 388 | "random_field = 0.01*np.random.randn(mg.size('node'))\n", |
285 | | - "z += random_field - random_field.min()\n", |
286 | | - "\n", |
| 389 | + "z += random_field - random_field.min()" |
| 390 | + ] |
| 391 | + }, |
| 392 | + { |
| 393 | + "cell_type": "markdown", |
| 394 | + "metadata": {}, |
| 395 | + "source": [ |
| 396 | + "Next we set all the parameters for LithoLayers. Here we have two types of rock with different erodabilities. " |
| 397 | + ] |
| 398 | + }, |
| 399 | + { |
| 400 | + "cell_type": "code", |
| 401 | + "execution_count": null, |
| 402 | + "metadata": { |
| 403 | + "collapsed": true |
| 404 | + }, |
| 405 | + "outputs": [], |
| 406 | + "source": [ |
287 | 407 | "attrs = {'K_sp': {0: 0.0003,\n", |
288 | 408 | " 1: 0.0001}}\n", |
289 | 409 | "\n", |
290 | | - "\n", |
291 | 410 | "z0s = 50 * np.arange(-20, 20)\n", |
292 | | - "\n", |
293 | | - "# we create a bottom layer that is very thick. \n", |
294 | 411 | "z0s[-1] = z0s[-2] + 10000\n", |
295 | 412 | "\n", |
296 | | - "\n", |
297 | | - "ids = np.tile([0,1], 20) \n", |
298 | | - "\n", |
299 | | - "\n", |
300 | | - "\n", |
| 413 | + "ids = np.tile([0,1], 20) " |
| 414 | + ] |
| 415 | + }, |
| 416 | + { |
| 417 | + "cell_type": "markdown", |
| 418 | + "metadata": {}, |
| 419 | + "source": [ |
| 420 | + "There are three functional forms that you can choose between. Here we define each of them. " |
| 421 | + ] |
| 422 | + }, |
| 423 | + { |
| 424 | + "cell_type": "code", |
| 425 | + "execution_count": null, |
| 426 | + "metadata": { |
| 427 | + "collapsed": true |
| 428 | + }, |
| 429 | + "outputs": [], |
| 430 | + "source": [ |
301 | 431 | "# Anticline\n", |
302 | 432 | "anticline_func = lambda x, y : ((0.002*x)**2+(0.001*y)**2)\n", |
303 | 433 | "\n", |
|
312 | 442 | "cell_type": "markdown", |
313 | 443 | "metadata": {}, |
314 | 444 | "source": [ |
315 | | - "Intervening text." |
| 445 | + "The default option is to make an anticline, but you can comment/uncomment lines to choose a different functional form. " |
316 | 446 | ] |
317 | 447 | }, |
318 | 448 | { |
|
0 commit comments