Skip to content

Commit 413555a

Browse files
kbarnhartkbarnhart
committed
initial commit of rockblock tutorial
first attempt of this accidentaly had a run tutorial that was 24 MB. [skip ci]
1 parent 958abce commit 413555a

1 file changed

Lines changed: 386 additions & 0 deletions

File tree

‎rockblock/introduction_to_rockblock.ipynb‎

Lines changed: 386 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,386 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Introduction to the RockBlock and LayeredRockBlock objects\n"
8+
]
9+
},
10+
{
11+
"cell_type": "code",
12+
"execution_count": null,
13+
"metadata": {},
14+
"outputs": [],
15+
"source": [
16+
"import os\n",
17+
"import numpy as np\n",
18+
"import xarray as xr\n",
19+
"import dask\n",
20+
"%matplotlib inline\n",
21+
"\n",
22+
"import holoviews as hv\n",
23+
"hv.extension('bokeh', 'matplotlib')\n",
24+
"\n",
25+
"from landlab import RasterModelGrid\n",
26+
"from landlab.layers import RockBlock, LayeredRockBlock\n",
27+
"from landlab.components import FlowAccumulator, FastscapeEroder\n",
28+
"from landlab.plot import imshow_grid"
29+
]
30+
},
31+
{
32+
"cell_type": "markdown",
33+
"metadata": {},
34+
"source": []
35+
},
36+
{
37+
"cell_type": "code",
38+
"execution_count": null,
39+
"metadata": {
40+
"collapsed": true
41+
},
42+
"outputs": [],
43+
"source": [
44+
"attrs = {'K_sp': {0: 0.0003,\n",
45+
" 1: 0.0001}}\n",
46+
"\n",
47+
"mg = RasterModelGrid((50, 50), 100)\n",
48+
"z = mg.add_zeros('node', 'topographic__elevation') \n",
49+
"random_field = 0.01*np.random.randn(mg.size('node'))\n",
50+
"z += random_field - random_field.min()\n",
51+
"\n",
52+
"z0s = 50 * np.arange(-20, 20)\n",
53+
"z0s[-1] = z0s[-2] + 10000\n",
54+
"ids = np.tile([0,1], 20) \n",
55+
"\n",
56+
"# Anticline\n",
57+
"rb = LayeredRockBlock(mg, z0s, ids, x0=5000, y0=5000, function = lambda x, y : ((0.001*x)**2+(0.003*y)**2), attrs=attrs)\n",
58+
"\n",
59+
"# Shallow dips\n",
60+
"#rb = LayeredRockBlock(mg, z0s, ids, x0=5000, y0=5000, function = lambda x, y : ((0.001*x)+(0.003*y)), attrs=attrs)\n",
61+
"\n",
62+
"# Steeper dips\n",
63+
"#rb = LayeredRockBlock(mg, z0s, ids, x0=5000, y0=5000, function = lambda x, y : ((0.01*x)+(0.01*y)), attrs=attrs)"
64+
]
65+
},
66+
{
67+
"cell_type": "markdown",
68+
"metadata": {},
69+
"source": []
70+
},
71+
{
72+
"cell_type": "code",
73+
"execution_count": null,
74+
"metadata": {
75+
"collapsed": true
76+
},
77+
"outputs": [],
78+
"source": [
79+
"imshow_grid(mg, 'K_sp')"
80+
]
81+
},
82+
{
83+
"cell_type": "markdown",
84+
"metadata": {},
85+
"source": []
86+
},
87+
{
88+
"cell_type": "code",
89+
"execution_count": null,
90+
"metadata": {
91+
"collapsed": true
92+
},
93+
"outputs": [],
94+
"source": [
95+
"nts = 400\n",
96+
"U = 0.001\n",
97+
"dt = 1000\n",
98+
"\n",
99+
"fa = FlowAccumulator(mg)\n",
100+
"sp = FastscapeEroder(mg, K_sp='K_sp')\n",
101+
"\n",
102+
"out_fields = ['topographic__elevation',\n",
103+
" 'K_sp', \n",
104+
" 'rock_type__id']\n",
105+
"\n",
106+
"ds = xr.Dataset(data_vars={'topographic__elevation' : (('y', 'x', 't'), np.empty((mg.shape[0], mg.shape[1], nts))),\n",
107+
" 'K_sp' : (('y', 'x', 't'), np.empty((mg.shape[0], mg.shape[1], nts))),\n",
108+
" 'rock_type__id': (('y', 'x', 't'), np.empty((mg.shape[0], mg.shape[1], nts)))},\n",
109+
" coords={'x': mg.x_of_node.reshape(mg.shape)[0,:],\n",
110+
" 'y': mg.y_of_node.reshape(mg.shape)[:, 1],\n",
111+
" 't': dt*np.arange(nts)})"
112+
]
113+
},
114+
{
115+
"cell_type": "markdown",
116+
"metadata": {},
117+
"source": []
118+
},
119+
{
120+
"cell_type": "code",
121+
"execution_count": null,
122+
"metadata": {
123+
"collapsed": true
124+
},
125+
"outputs": [],
126+
"source": [
127+
"for i in range(nts):\n",
128+
" fa.run_one_step()\n",
129+
" sp.run_one_step(dt = dt)\n",
130+
" dz_ad = np.zeros(mg.size('node'))\n",
131+
" dz_ad[mg.core_nodes] = U * dt\n",
132+
" z += dz_ad\n",
133+
" rb.run_one_step(dz_advection = dz_ad)\n",
134+
" \n",
135+
" for of in out_fields:\n",
136+
" ds[of][:,:,i] = mg['node'][of].reshape(mg.shape)\n"
137+
]
138+
},
139+
{
140+
"cell_type": "markdown",
141+
"metadata": {},
142+
"source": []
143+
},
144+
{
145+
"cell_type": "code",
146+
"execution_count": null,
147+
"metadata": {
148+
"collapsed": true
149+
},
150+
"outputs": [],
151+
"source": [
152+
"imshow_grid(mg, 'topographic__elevation')"
153+
]
154+
},
155+
{
156+
"cell_type": "markdown",
157+
"metadata": {},
158+
"source": []
159+
},
160+
{
161+
"cell_type": "code",
162+
"execution_count": null,
163+
"metadata": {
164+
"collapsed": true
165+
},
166+
"outputs": [],
167+
"source": [
168+
"hvds = hv.Dataset(ds)\n",
169+
"hvds"
170+
]
171+
},
172+
{
173+
"cell_type": "markdown",
174+
"metadata": {},
175+
"source": []
176+
},
177+
{
178+
"cell_type": "code",
179+
"execution_count": null,
180+
"metadata": {
181+
"collapsed": true
182+
},
183+
"outputs": [],
184+
"source": [
185+
"%opts Image (cmap='viridis')\n",
186+
"\n",
187+
"topo = hvds.to(hv.Image, ['x', 'y'], ['topographic__elevation'])\n",
188+
"rock = hvds.to(hv.Image, ['x', 'y'], ['rock_type__id'])\n",
189+
"\n",
190+
"topo + rock"
191+
]
192+
},
193+
{
194+
"cell_type": "markdown",
195+
"metadata": {
196+
"collapsed": true
197+
},
198+
"source": [
199+
"Make inverted topography by filling RockBlock in with resistant material only where the DA is large. "
200+
]
201+
},
202+
{
203+
"cell_type": "code",
204+
"execution_count": null,
205+
"metadata": {
206+
"collapsed": true
207+
},
208+
"outputs": [],
209+
"source": [
210+
"mg2 = RasterModelGrid((50, 50), 100)\n",
211+
"z2 = mg2.add_zeros('node', 'topographic__elevation') \n",
212+
"random_field = 0.01*np.random.randn(mg2.size('node'))\n",
213+
"z2 += random_field - random_field.min()\n",
214+
"\n",
215+
"thicknesses2 = [1000]\n",
216+
"ids2 =[0]\n",
217+
"\n",
218+
"attrs2 = {'K_sp': {0: 0.0001,\n",
219+
" 1: 0.00001}}\n",
220+
"\n",
221+
"rb2 = RockBlock(mg2, thicknesses2, ids2, attrs=attrs2)\n",
222+
"\n",
223+
"nts = 400\n",
224+
"U = 0.001\n",
225+
"dt = 1000\n",
226+
"\n",
227+
"fa2 = FlowAccumulator(mg2)\n",
228+
"sp2 = FastscapeEroder(mg2, K_sp='K_sp')\n",
229+
"\n",
230+
"out_fields = ['topographic__elevation',\n",
231+
" 'K_sp', \n",
232+
" 'rock_type__id']\n",
233+
"\n",
234+
"ds2 = xr.Dataset(data_vars={'topographic__elevation' : (('y', 'x', 't'), np.empty((mg2.shape[0], mg2.shape[1], nts))),\n",
235+
" 'K_sp' : (('y', 'x', 't'), np.empty((mg2.shape[0], mg2.shape[1], nts))),\n",
236+
" 'rock_type__id': (('y', 'x', 't'), np.empty((mg2.shape[0], mg2.shape[1], nts)))},\n",
237+
" coords={'x': mg2.x_of_node.reshape(mg2.shape)[0,:],\n",
238+
" 'y': mg2.y_of_node.reshape(mg2.shape)[:, 1],\n",
239+
" 't': dt*np.arange(nts)})\n",
240+
"half_nts = int(nts/2)\n",
241+
"for i in range(half_nts):\n",
242+
" fa2.run_one_step()\n",
243+
" sp2.run_one_step(dt = dt)\n",
244+
" dz_ad2 = np.zeros(mg2.size('node'))\n",
245+
" dz_ad2[mg2.core_nodes] = U * dt\n",
246+
" z2 += dz_ad2\n",
247+
" rb2.run_one_step(dz_advection = dz_ad2)\n",
248+
" \n",
249+
" for of in out_fields:\n",
250+
" ds2[of][:,:,i] = mg2['node'][of].reshape(mg2.shape)\n"
251+
]
252+
},
253+
{
254+
"cell_type": "markdown",
255+
"metadata": {},
256+
"source": []
257+
},
258+
{
259+
"cell_type": "code",
260+
"execution_count": null,
261+
"metadata": {},
262+
"outputs": [],
263+
"source": [
264+
"imshow_grid(mg2, 'topographic__elevation')"
265+
]
266+
},
267+
{
268+
"cell_type": "markdown",
269+
"metadata": {},
270+
"source": []
271+
},
272+
{
273+
"cell_type": "code",
274+
"execution_count": null,
275+
"metadata": {},
276+
"outputs": [],
277+
"source": [
278+
"volcano_deposits = np.zeros(mg2.size('node'))\n",
279+
"da_big_enough = mg2['node']['drainage_area']>1e4\n",
280+
"volcano_deposits[da_big_enough] = 30.0\n",
281+
"volcano_deposits[mg2.boundary_nodes] = 0.0\n",
282+
"rb2.add_layer(volcano_deposits, rock_id=1)\n",
283+
"\n",
284+
"imshow_grid(mg2, da_big_enough)"
285+
]
286+
},
287+
{
288+
"cell_type": "markdown",
289+
"metadata": {},
290+
"source": []
291+
},
292+
{
293+
"cell_type": "code",
294+
"execution_count": null,
295+
"metadata": {
296+
"collapsed": true
297+
},
298+
"outputs": [],
299+
"source": [
300+
"for i in range(half_nts, nts):\n",
301+
" fa2.run_one_step()\n",
302+
" sp2.run_one_step(dt = dt)\n",
303+
" dz_ad = np.zeros(mg2.size('node'))\n",
304+
" dz_ad[mg2.core_nodes] = U * dt\n",
305+
" z2 += dz_ad\n",
306+
" rb2.run_one_step(dz_advection = dz_ad)\n",
307+
" \n",
308+
" for of in out_fields:\n",
309+
" ds2[of][:,:,i] = mg2['node'][of].reshape(mg2.shape)"
310+
]
311+
},
312+
{
313+
"cell_type": "markdown",
314+
"metadata": {},
315+
"source": []
316+
},
317+
{
318+
"cell_type": "code",
319+
"execution_count": null,
320+
"metadata": {},
321+
"outputs": [],
322+
"source": [
323+
"imshow_grid(mg2, 'topographic__elevation')"
324+
]
325+
},
326+
{
327+
"cell_type": "markdown",
328+
"metadata": {},
329+
"source": []
330+
},
331+
{
332+
"cell_type": "code",
333+
"execution_count": null,
334+
"metadata": {},
335+
"outputs": [],
336+
"source": [
337+
"%opts Image (cmap='viridis')\n",
338+
"\n",
339+
"hvds2 = hv.Dataset(ds2)\n",
340+
"topo2 = hvds2.to(hv.Image, ['x', 'y'], ['topographic__elevation'])\n",
341+
"rock2 = hvds2.to(hv.Image, ['x', 'y'], ['rock_type__id'])\n",
342+
"\n",
343+
"topo2 + rock2"
344+
]
345+
},
346+
{
347+
"cell_type": "code",
348+
"execution_count": null,
349+
"metadata": {
350+
"collapsed": true
351+
},
352+
"outputs": [],
353+
"source": []
354+
},
355+
{
356+
"cell_type": "code",
357+
"execution_count": null,
358+
"metadata": {
359+
"collapsed": true
360+
},
361+
"outputs": [],
362+
"source": []
363+
}
364+
],
365+
"metadata": {
366+
"kernelspec": {
367+
"display_name": "Python 3",
368+
"language": "python",
369+
"name": "python3"
370+
},
371+
"language_info": {
372+
"codemirror_mode": {
373+
"name": "ipython",
374+
"version": 3
375+
},
376+
"file_extension": ".py",
377+
"mimetype": "text/x-python",
378+
"name": "python",
379+
"nbconvert_exporter": "python",
380+
"pygments_lexer": "ipython3",
381+
"version": "3.6.2"
382+
}
383+
},
384+
"nbformat": 4,
385+
"nbformat_minor": 2
386+
}

0 commit comments

Comments
 (0)