Skip to content

Commit ac1e8b1

Browse files
kbarnhartkbarnhart
committed
better x array documentation
1 parent f1236aa commit ac1e8b1

1 file changed

Lines changed: 77 additions & 47 deletions

File tree

‎lithology/introduction_to_lithology.ipynb‎

Lines changed: 77 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -70,7 +70,7 @@
7070
"attrs = {'K_sp': {0: 0.0003,\n",
7171
" 1: 0.0001}}\n",
7272
"\n",
73-
"mg = RasterModelGrid((70, 70), 200)\n",
73+
"mg = RasterModelGrid((100, 60), 200)\n",
7474
"z = mg.add_zeros('node', 'topographic__elevation') \n",
7575
"random_field = 0.01*np.random.randn(mg.size('node'))\n",
7676
"z += random_field - random_field.min()\n",
@@ -80,13 +80,13 @@
8080
"ids = np.tile([0,1], 20) \n",
8181
"\n",
8282
"# Anticline\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",
83+
"lith = LithoLayers(mg, z0s, ids, x0=6000, y0=6000, function = lambda x, y : ((0.002*x)**2+(0.001*y)**2), attrs=attrs)\n",
8484
"\n",
8585
"# Shallow dips\n",
86-
"#lith = LithoLayers(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, function = lambda x, y : ((0.001*x)+(0.003*y)), attrs=attrs)\n",
8787
"\n",
8888
"# Steeper dips\n",
89-
"#lith = LithoLayers(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, function = lambda x, y : ((0.01*x)+(0.01*y)), attrs=attrs)"
9090
]
9191
},
9292
{
@@ -101,9 +101,7 @@
101101
{
102102
"cell_type": "code",
103103
"execution_count": null,
104-
"metadata": {
105-
"collapsed": true
106-
},
104+
"metadata": {},
107105
"outputs": [],
108106
"source": [
109107
"imshow_grid(mg, 'K_sp')"
@@ -140,11 +138,26 @@
140138
"out_fields = ['topographic__elevation',\n",
141139
" 'rock_type__id']\n",
142140
"\n",
143-
"ds = xr.Dataset(data_vars={'topographic__elevation' : (('t', 'y', 'x'), np.empty((nts, mg.shape[0], mg.shape[1]))),\n",
144-
" 'rock_type__id': (('t', 'y', 'x'), np.empty((nts, mg.shape[0], mg.shape[1])))},\n",
145-
" coords={'x': mg.x_of_node.reshape(mg.shape)[0,:],\n",
146-
" 'y': mg.y_of_node.reshape(mg.shape)[:, 1],\n",
147-
" 't': dt*np.arange(nts)/1e6})"
141+
"ds = xr.Dataset(data_vars={'topographic__elevation' : (('t', 'y', 'x'), \n",
142+
" np.empty((nts, mg.shape[0], mg.shape[1])),\n",
143+
" {'units' : 'meters',\n",
144+
" 'long_name': 'Topographic Elevation'}),\n",
145+
" 'rock_type__id': (('t', 'y', 'x'), \n",
146+
" np.empty((nts, mg.shape[0], mg.shape[1])),\n",
147+
" {'units' : '-',\n",
148+
" 'long_name' : 'Rock Type ID Code'})\n",
149+
" },\n",
150+
" coords={'x': (('x'), \n",
151+
" mg.x_of_node.reshape(mg.shape)[0,:],\n",
152+
" {'units' : 'meters'}),\n",
153+
" 'y': (('y'), \n",
154+
" mg.y_of_node.reshape(mg.shape)[:, 1],\n",
155+
" {'units' : 'meters'}),\n",
156+
" 'time': (('t'), \n",
157+
" dt*np.arange(nts)/1e6,\n",
158+
" {'units': 'millions of years since model start'})\n",
159+
" }\n",
160+
" )\n"
148161
]
149162
},
150163
{
@@ -161,9 +174,7 @@
161174
{
162175
"cell_type": "code",
163176
"execution_count": null,
164-
"metadata": {
165-
"collapsed": true
166-
},
177+
"metadata": {},
167178
"outputs": [],
168179
"source": [
169180
"for i in range(nts):\n",
@@ -188,9 +199,7 @@
188199
{
189200
"cell_type": "code",
190201
"execution_count": null,
191-
"metadata": {
192-
"collapsed": true
193-
},
202+
"metadata": {},
194203
"outputs": [],
195204
"source": [
196205
"imshow_grid(mg, 'topographic__elevation', cmap='viridis')"
@@ -210,9 +219,7 @@
210219
{
211220
"cell_type": "code",
212221
"execution_count": null,
213-
"metadata": {
214-
"collapsed": true
215-
},
222+
"metadata": {},
216223
"outputs": [],
217224
"source": [
218225
"hvds_topo = hv.Dataset(ds.topographic__elevation)\n",
@@ -232,9 +239,7 @@
232239
{
233240
"cell_type": "code",
234241
"execution_count": null,
235-
"metadata": {
236-
"collapsed": true
237-
},
242+
"metadata": {},
238243
"outputs": [],
239244
"source": [
240245
"%opts Image style(interpolation='bilinear', cmap='viridis') plot[colorbar=True]\n",
@@ -261,19 +266,17 @@
261266
"\n",
262267
"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",
263268
"\n",
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. "
269+
"We also are handling the model grid boundary conditions differently than in the last example, setting the boundaries on the top and bottom to closed. "
265270
]
266271
},
267272
{
268273
"cell_type": "code",
269274
"execution_count": null,
270-
"metadata": {
271-
"collapsed": true
272-
},
275+
"metadata": {},
273276
"outputs": [],
274277
"source": [
275-
"mg2 = RasterModelGrid((70, 70), 200)\n",
276-
"mg2.set_closed_boundaries_at_grid_edges(True, False, True, False)\n",
278+
"mg2 = RasterModelGrid((100, 60), 200)\n",
279+
"mg2.set_closed_boundaries_at_grid_edges(False, True, False, True)\n",
277280
"z2 = mg2.add_zeros('node', 'topographic__elevation') \n",
278281
"random_field = 0.01*np.random.randn(mg2.size('node'))\n",
279282
"z2 += random_field - random_field.min()\n",
@@ -302,11 +305,33 @@
302305
"out_fields = ['topographic__elevation',\n",
303306
" 'rock_type__id']\n",
304307
"\n",
305-
"ds2 = xr.Dataset(data_vars={'topographic__elevation' : (('t', 'y', 'x'), np.empty((nts, mg2.shape[0], mg2.shape[1]))),\n",
306-
" 'rock_type__id': (('t', 'y', 'x'), np.empty((nts, mg2.shape[0], mg2.shape[1])))},\n",
307-
" coords={'x': mg2.x_of_node.reshape(mg2.shape)[0,:],\n",
308-
" 'y': mg2.y_of_node.reshape(mg2.shape)[:, 1],\n",
309-
" 't': dt*np.arange(nts)/1e6})\n",
308+
"nts = 500\n",
309+
"U = 0.001\n",
310+
"dt = 1000\n",
311+
"\n",
312+
"\n",
313+
"ds2 = xr.Dataset(data_vars={'topographic__elevation' : (('t', 'y', 'x'), \n",
314+
" np.empty((nts, mg2.shape[0], mg2.shape[1])),\n",
315+
" {'units' : 'meters',\n",
316+
" 'long_name': 'Topographic Elevation'}),\n",
317+
" 'rock_type__id': (('t', 'y', 'x'), \n",
318+
" np.empty((nts, mg2.shape[0], mg2.shape[1])),\n",
319+
" {'units' : '-',\n",
320+
" 'long_name' : 'Rock Type ID Code'})\n",
321+
" },\n",
322+
" coords={'x': (('x'), \n",
323+
" mg2.x_of_node.reshape(mg2.shape)[0,:],\n",
324+
" {'units' : 'meters'}),\n",
325+
" 'y': (('y'), \n",
326+
" mg2.y_of_node.reshape(mg2.shape)[:, 1],\n",
327+
" {'units' : 'meters'}),\n",
328+
" 'time': (('t'), \n",
329+
" dt*np.arange(nts)/1e6,\n",
330+
" {'units': 'millions of years since model start'})\n",
331+
" }\n",
332+
" )\n",
333+
"\n",
334+
"\n",
310335
"half_nts = int(nts/2)\n",
311336
"for i in range(half_nts):\n",
312337
" fa2.run_one_step()\n",
@@ -331,9 +356,7 @@
331356
{
332357
"cell_type": "code",
333358
"execution_count": null,
334-
"metadata": {
335-
"collapsed": true
336-
},
359+
"metadata": {},
337360
"outputs": [],
338361
"source": [
339362
"imshow_grid(mg2, 'topographic__elevation', cmap='viridis')"
@@ -351,9 +374,7 @@
351374
{
352375
"cell_type": "code",
353376
"execution_count": null,
354-
"metadata": {
355-
"collapsed": true
356-
},
377+
"metadata": {},
357378
"outputs": [],
358379
"source": [
359380
"volcanic_deposits = np.zeros(mg2.size('node'))\n",
@@ -408,9 +429,7 @@
408429
{
409430
"cell_type": "code",
410431
"execution_count": null,
411-
"metadata": {
412-
"collapsed": true
413-
},
432+
"metadata": {},
414433
"outputs": [],
415434
"source": [
416435
"imshow_grid(mg2, 'topographic__elevation', cmap='viridis')"
@@ -426,9 +445,7 @@
426445
{
427446
"cell_type": "code",
428447
"execution_count": null,
429-
"metadata": {
430-
"collapsed": true
431-
},
448+
"metadata": {},
432449
"outputs": [],
433450
"source": [
434451
"hvds_topo2 = hv.Dataset(ds2.topographic__elevation)\n",
@@ -441,6 +458,19 @@
441458
"topo2 + rock2"
442459
]
443460
},
461+
{
462+
"cell_type": "code",
463+
"execution_count": null,
464+
"metadata": {
465+
"collapsed": true
466+
},
467+
"outputs": [],
468+
"source": [
469+
"# if you wanted to output to visualize in something like ParaView, the following commands can be used\n",
470+
"#ds.to_netcdf('anticline.nc')\n",
471+
"#ds2.to_netcdf('inversion.nc')"
472+
]
473+
},
444474
{
445475
"cell_type": "markdown",
446476
"metadata": {

0 commit comments

Comments
 (0)