Skip to content

Commit 7703f9a

Browse files
kbarnhartkbarnhart
committed
intermediate updates to tutorial
[skip ci]
1 parent 7201d98 commit 7703f9a

1 file changed

Lines changed: 160 additions & 42 deletions

File tree

‎lithology/introduction_to_lithology.ipynb‎

Lines changed: 160 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -52,28 +52,41 @@
5252
"source": [
5353
"## Part 1: Creating layered rock\n",
5454
"\n",
55-
"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",
55+
"First we will create an instance of a LithoLayers to learn how this component works. Both LithoLayers and Lithology work closely with a landlab ModelGrid and stores information about rock type at each node. \n",
5656
"\n",
5757
"To create LithoLayers you need the following information:\n",
5858
"\n",
59-
"1. A model grid that has the field `'topographic__elevation'` already created.\n",
60-
"2. A dictionary of rock property attributes that maps. \n",
61-
"3. A list of elevations, alled `z0` that your layers will go through at specified anchor point (default value for the anchor point is (0, 0). When `z0` is negative.... TODO.\n",
62-
"4. A functional form in two variables (x and y) that defines the shape of your surface. \n",
59+
"1. A model grid that has the field `'topographic__elevation'` already created. \n",
60+
"2. A list of elevations, called `'layer_elevations'` that the bottom of your layers will go through at specified plan-view anchor point (default value for the anchor point is (x, y) = (0, 0) and a list of IDs that indicate the rock type of that layer. When `'layer_elevations'` is negative that means that the layer goes through the anchor point above the topographic surface. These layers will be created where they extend below the topographic surface.\n",
61+
"3. A dictionary of rock property attributes that maps a rock ID type to property values. \n",
62+
"4. A functional form in x and y that defines the shape of your surface. \n",
6363
"\n",
6464
"The use of this function form makes it possible for any function of x and y to be passed to LithoLayers.\n",
65-
"In this example the option for an anticline is currently uncommented, though options for steep and shallow dipping layers are also provided. \n",
6665
"\n",
67-
"In this tutorial we will first make an example to help build intuition. If you are interested in reading more... [LithoLayers](http://landlab.readthedocs.io/en/release/landlab.layers.litholayers.html)\n",
66+
"In this tutorial we will first make an example to help build intuition and then do two more complex examples. Most of the functionality of Lithology and LithoLayers is shown in this tutorial, but if you want to read the full component documentation for LithoLayers, it can be found [here](http://landlab.readthedocs.io/en/release/landlab.layers.litholayers.html). Links to both components documentation can be found at the bottom of the tutorial.\n",
6867
"\n",
69-
"First, we create a RasterModelGrid with topography initialized with a small amount of random noise. "
68+
"First, we create a small RasterModelGrid with topography. "
69+
]
70+
},
71+
{
72+
"cell_type": "code",
73+
"execution_count": null,
74+
"metadata": {
75+
"collapsed": true
76+
},
77+
"outputs": [],
78+
"source": [
79+
"mg = RasterModelGrid((10, 10), 1)\n",
80+
"z = mg.add_zeros('node', 'topographic__elevation') "
7081
]
7182
},
7283
{
7384
"cell_type": "markdown",
7485
"metadata": {},
7586
"source": [
76-
"Make a very small grid, with a very easy example of layers varying, and illustrate that first so people get their confidence up. I know people can read the docstrings, but I think the point of the notebook is that people should not need to read the docstring with the notebook?"
87+
"Next we make our layer elevations. We will make 20 layers that are 5 meters thick. At the anchor point, half of the layers will be above the ground (`'layer_elevations'` will have negative values) and half will be below the ground (`'layer_elevations'` have positive values). \n",
88+
"\n",
89+
"We will make this with the [`np.arrange`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.arange.html) function. We will also make the bottom layer really really thick so that we won't be able to erode through through it. "
7790
]
7891
},
7992
{
@@ -84,21 +97,35 @@
8497
},
8598
"outputs": [],
8699
"source": [
87-
"mg = RasterModelGrid((100, 60), 200)\n",
88-
"z = mg.add_zeros('node', 'topographic__elevation') \n",
89-
"random_field = 0.01*np.random.randn(mg.size('node'))\n",
90-
"z += random_field - random_field.min()"
100+
"layer_elevations = 5. * np.arange(-10, 10)\n",
101+
"\n",
102+
"# we create a bottom layer that is very thick. \n",
103+
"layer_elevations[-1] = layer_elevations[-2] + 100"
91104
]
92105
},
93106
{
94107
"cell_type": "markdown",
95108
"metadata": {},
96109
"source": [
97-
"Next we.... \n",
98-
"\n",
99-
"Desribe attrs, z0s, bottom layer, and IDS. say what tile is. \n",
100-
"\n",
101-
"Remind people what it means when z0s are negative. "
110+
"Next we create alternating `0`s `1`s `2`s and `3`s with the [np.tile](https://docs.scipy.org/doc/numpy/reference/generated/numpy.tile.html) function. "
111+
]
112+
},
113+
{
114+
"cell_type": "code",
115+
"execution_count": null,
116+
"metadata": {
117+
"collapsed": true
118+
},
119+
"outputs": [],
120+
"source": [
121+
"layer_ids = np.tile([0,1,2,3], 5) "
122+
]
123+
},
124+
{
125+
"cell_type": "markdown",
126+
"metadata": {},
127+
"source": [
128+
"Our attributes dictionary has the following form. "
102129
]
103130
},
104131
{
@@ -110,26 +137,70 @@
110137
"outputs": [],
111138
"source": [
112139
"attrs = {'K_sp': {0: 0.0003,\n",
113-
" 1: 0.0001}}\n",
114-
"\n",
115-
"\n",
116-
"z0s = 50 * np.arange(-20, 20)\n",
117-
"\n",
118-
"# we create a bottom layer that is very thick. \n",
119-
"z0s[-1] = z0s[-2] + 10000\n",
120-
"\n",
121-
"\n",
122-
"ids = np.tile([0,1], 20) \n",
123-
"\n",
140+
" 1: 0.0001,\n",
141+
" 2: 0.0002,\n",
142+
" 3: 0.0004}}"
143+
]
144+
},
145+
{
146+
"cell_type": "markdown",
147+
"metadata": {},
148+
"source": [
149+
"`'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",
124150
"\n",
125-
"\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)`"
152+
]
153+
},
154+
{
155+
"cell_type": "code",
156+
"execution_count": null,
157+
"metadata": {
158+
"collapsed": true
159+
},
160+
"outputs": [],
161+
"source": [
162+
"func = lambda x, y : x + (2. * y)"
163+
]
164+
},
165+
{
166+
"cell_type": "markdown",
167+
"metadata": {},
168+
"source": [
169+
"Finally we construct our LithoLayers component."
170+
]
171+
},
172+
{
173+
"cell_type": "code",
174+
"execution_count": null,
175+
"metadata": {},
176+
"outputs": [],
177+
"source": [
178+
"lith = LithoLayers(mg, layer_elevations, layer_ids, function=func, attrs=attrs)"
126179
]
127180
},
128181
{
129182
"cell_type": "markdown",
130183
"metadata": {},
131184
"source": [
132-
"TODO, describe making functions"
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. "
186+
]
187+
},
188+
{
189+
"cell_type": "code",
190+
"execution_count": null,
191+
"metadata": {},
192+
"outputs": [],
193+
"source": [
194+
"imshow_grid(mg, 'K_sp')"
195+
]
196+
},
197+
{
198+
"cell_type": "markdown",
199+
"metadata": {},
200+
"source": [
201+
"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+
"\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. "
133204
]
134205
},
135206
{
@@ -140,15 +211,45 @@
140211
},
141212
"outputs": [],
142213
"source": [
143-
"func = lambda x, y : ((0.001*x)+(0.003*y))\n",
144-
"\n"
214+
"z -= 1.\n",
215+
"dz_ad = 0.\n",
216+
"lith.run_one_step(dz_advection = dz_ad)"
145217
]
146218
},
147219
{
148220
"cell_type": "markdown",
149221
"metadata": {},
150222
"source": [
151-
"Finally we construct our LithoLayers component. Make a plot and talk through it. "
223+
"We can re-plot the value of `'K_sp'`. We will see that the locatio of layers has changed. "
224+
]
225+
},
226+
{
227+
"cell_type": "code",
228+
"execution_count": null,
229+
"metadata": {},
230+
"outputs": [],
231+
"source": [
232+
"imshow_grid(mg, 'K_sp')"
233+
]
234+
},
235+
{
236+
"cell_type": "code",
237+
"execution_count": null,
238+
"metadata": {},
239+
"outputs": [],
240+
"source": [
241+
"ds = lith.to_xarray(np.arange(10))\n",
242+
"hvds_rock = hv.Dataset (ds.rock_type__id)"
243+
]
244+
},
245+
{
246+
"cell_type": "code",
247+
"execution_count": null,
248+
"metadata": {},
249+
"outputs": [],
250+
"source": [
251+
"%opts Image style(interpolation='bilinear', cmap='viridis') plot[colorbar=True]\n",
252+
"hvds_rock.to(hv.Image, ['x', 'y'])"
152253
]
153254
},
154255
{
@@ -244,7 +345,9 @@
244345
{
245346
"cell_type": "code",
246347
"execution_count": null,
247-
"metadata": {},
348+
"metadata": {
349+
"collapsed": true
350+
},
248351
"outputs": [],
249352
"source": [
250353
"imshow_grid(mg, 'K_sp')"
@@ -334,7 +437,9 @@
334437
{
335438
"cell_type": "code",
336439
"execution_count": null,
337-
"metadata": {},
440+
"metadata": {
441+
"collapsed": true
442+
},
338443
"outputs": [],
339444
"source": [
340445
"print(ds)"
@@ -352,7 +457,9 @@
352457
{
353458
"cell_type": "code",
354459
"execution_count": null,
355-
"metadata": {},
460+
"metadata": {
461+
"collapsed": true
462+
},
356463
"outputs": [],
357464
"source": [
358465
"ds.topographic__elevation"
@@ -402,7 +509,9 @@
402509
{
403510
"cell_type": "code",
404511
"execution_count": null,
405-
"metadata": {},
512+
"metadata": {
513+
"collapsed": true
514+
},
406515
"outputs": [],
407516
"source": [
408517
"imshow_grid(mg, 'topographic__elevation', cmap='viridis')"
@@ -422,7 +531,9 @@
422531
{
423532
"cell_type": "code",
424533
"execution_count": null,
425-
"metadata": {},
534+
"metadata": {
535+
"collapsed": true
536+
},
426537
"outputs": [],
427538
"source": [
428539
"hvds_topo = hv.Dataset(ds.topographic__elevation)\n",
@@ -442,7 +553,9 @@
442553
{
443554
"cell_type": "code",
444555
"execution_count": null,
445-
"metadata": {},
556+
"metadata": {
557+
"collapsed": true
558+
},
446559
"outputs": [],
447560
"source": [
448561
"%opts Image style(interpolation='bilinear', cmap='viridis') plot[colorbar=True]\n",
@@ -562,7 +675,9 @@
562675
{
563676
"cell_type": "code",
564677
"execution_count": null,
565-
"metadata": {},
678+
"metadata": {
679+
"collapsed": true
680+
},
566681
"outputs": [],
567682
"source": [
568683
"imshow_grid(mg2, 'topographic__elevation', cmap='viridis')"
@@ -580,7 +695,9 @@
580695
{
581696
"cell_type": "code",
582697
"execution_count": null,
583-
"metadata": {},
698+
"metadata": {
699+
"collapsed": true
700+
},
584701
"outputs": [],
585702
"source": [
586703
"volcanic_deposits = np.zeros(mg2.size('node'))\n",
@@ -636,6 +753,7 @@
636753
"cell_type": "code",
637754
"execution_count": null,
638755
"metadata": {
756+
"collapsed": true,
639757
"scrolled": true
640758
},
641759
"outputs": [],

0 commit comments

Comments
 (0)