Skip to content

Commit b0f05eb

Browse files
kbarnhartkbarnhart
committed
adding some additional information regarding xarray and annotated data
1 parent e1a7c6d commit b0f05eb

1 file changed

Lines changed: 96 additions & 24 deletions

File tree

‎lithology/introduction_to_lithology.ipynb‎

Lines changed: 96 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,8 @@
1717
"\n",
1818
"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",
1919
"\n",
20+
"We will also use [xarray](https://xarray.pydata.org/en/stable/) to store and annotate our model output. \n",
21+
"\n",
2022
"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."
2123
]
2224
},
@@ -115,9 +117,7 @@
115117
"\n",
116118
"Next, lets instantiate a FlowAccumulator and a FastscapeEroder to create a simple landscape evolution model. \n",
117119
"\n",
118-
"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",
119-
"\n",
120-
"We will also instatiate an xarray dataset used to store the output of our model through time for visualization. "
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. "
121121
]
122122
},
123123
{
@@ -133,39 +133,95 @@
133133
"dt = 1000\n",
134134
"\n",
135135
"fa = FlowAccumulator(mg)\n",
136-
"sp = FastscapeEroder(mg, K_sp='K_sp')\n",
136+
"sp = FastscapeEroder(mg, K_sp='K_sp')"
137+
]
138+
},
139+
{
140+
"cell_type": "markdown",
141+
"metadata": {},
142+
"source": [
143+
"Before we run the model we will also instatiate an xarray dataset used to store the output of our model through time for visualization. \n",
137144
"\n",
138-
"out_fields = ['topographic__elevation',\n",
139-
" 'rock_type__id']\n",
145+
"The next block may look intimidating, but I'll try and walk you through what it does. \n",
146+
"\n",
147+
"[xarray](https://xarray.pydata.org/en/stable/) allows us to create a container for our data and label it with information like units, dimensions, short and long names, etc. xarray gives all the tools for dealing with N-dimentional data provided by python packages such as [numpy](http://www.numpy.org), the labeling and named indexing power of the [pandas](https://pandas.pydata.org) package, and the data-model of the [NetCDF file](https://www.unidata.ucar.edu/software/netcdf/).\n",
140148
"\n",
141-
"ds = xr.Dataset(data_vars={'topographic__elevation' : (('time', 'y', 'x'), \n",
142-
" np.empty((nts, mg.shape[0], mg.shape[1])),\n",
143-
" {'units' : 'meters',\n",
149+
"This means that we can use xarray to make a \"self-referetial\" dataset that contains all of the variables and attributes that describe what each part is and how it was made. In this application, we won't make a fully self-referential dataset, but if you are interested in this, check out the [NetCDF best practices](https://www.unidata.ucar.edu/software/netcdf/docs/BestPractices.html). \n",
150+
"\n",
151+
"Important for our application is that later on we will use the [HoloViews package](http://holoviews.org) for visualization. This package is a great tool for dealing with multidimentional annotated data and will do things like automatically create nice axis labels with units. However, in order for it to work, we must first annotate our data to include this information.\n",
152+
"\n",
153+
"Here we create an xarray Dataset with two variables `'topographic__elevation'` and `'rock_type__id'` and three dimensions `'x'`, `'y'`, and `'time'`. \n",
154+
"\n",
155+
"We pass xarray two dictionaries, one with information about the data variabiables (`data_vars`) and one with information about the coordinate system (`coords`). For each data variable or coordinate, we pass a tuple of three items: `(dims, data, atts)`. The first element is a tuple of the name of the dimensions, the second element is the data, an the third is a dictionary of attributes. "
156+
]
157+
},
158+
{
159+
"cell_type": "code",
160+
"execution_count": null,
161+
"metadata": {
162+
"collapsed": true
163+
},
164+
"outputs": [],
165+
"source": [
166+
"ds = xr.Dataset(data_vars={'topographic__elevation' : (('time', 'y', 'x'), # tuple of dimensions\n",
167+
" np.empty((nts, mg.shape[0], mg.shape[1])), # n-d array of data\n",
168+
" {'units' : 'meters', # dictionary with data attributes\n",
144169
" 'long_name': 'Topographic Elevation'}),\n",
145170
" 'rock_type__id': (('time', 'y', 'x'), \n",
146171
" np.empty((nts, mg.shape[0], mg.shape[1])),\n",
147172
" {'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",
173+
" 'long_name' : 'Rock Type ID Code'})},\n",
174+
" coords={'x': (('x'), # tuple of dimensions\n",
175+
" mg.x_of_node.reshape(mg.shape)[0,:], # 1-d array of coordinate data\n",
176+
" {'units' : 'meters'}), # dictionary with data attributes\n",
153177
" 'y': (('y'), \n",
154178
" mg.y_of_node.reshape(mg.shape)[:, 1],\n",
155179
" {'units' : 'meters'}),\n",
156180
" 'time': (('time'), \n",
157181
" dt*np.arange(nts)/1e6,\n",
158182
" {'units': 'millions of years since model start',\n",
159-
" 'standard_name' : 'time'})\n",
160-
" }\n",
161-
" )\n"
183+
" 'standard_name' : 'time'})})"
184+
]
185+
},
186+
{
187+
"cell_type": "markdown",
188+
"metadata": {},
189+
"source": [
190+
"We can print the data set to get some basic information about it."
191+
]
192+
},
193+
{
194+
"cell_type": "code",
195+
"execution_count": null,
196+
"metadata": {},
197+
"outputs": [],
198+
"source": [
199+
"print(ds)"
162200
]
163201
},
164202
{
165203
"cell_type": "markdown",
166204
"metadata": {},
167205
"source": [
168-
"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",
206+
"We can also print a single variable to get more detailed information about it. \n",
207+
"\n",
208+
"Since we initialized the datset with empty arrays for the two data variables, we just see zeros for the data values. "
209+
]
210+
},
211+
{
212+
"cell_type": "code",
213+
"execution_count": null,
214+
"metadata": {},
215+
"outputs": [],
216+
"source": [
217+
"ds.topographic__elevation"
218+
]
219+
},
220+
{
221+
"cell_type": "markdown",
222+
"metadata": {},
223+
"source": [
224+
"Next, 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",
169225
"\n",
170226
"`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",
171227
"\n",
@@ -180,6 +236,9 @@
180236
},
181237
"outputs": [],
182238
"source": [
239+
"out_fields = ['topographic__elevation',\n",
240+
" 'rock_type__id']\n",
241+
"\n",
183242
"for i in range(nts):\n",
184243
" fa.run_one_step()\n",
185244
" sp.run_one_step(dt = dt)\n",
@@ -202,7 +261,9 @@
202261
{
203262
"cell_type": "code",
204263
"execution_count": null,
205-
"metadata": {},
264+
"metadata": {
265+
"collapsed": true
266+
},
206267
"outputs": [],
207268
"source": [
208269
"imshow_grid(mg, 'topographic__elevation', cmap='viridis')"
@@ -222,7 +283,9 @@
222283
{
223284
"cell_type": "code",
224285
"execution_count": null,
225-
"metadata": {},
286+
"metadata": {
287+
"collapsed": true
288+
},
226289
"outputs": [],
227290
"source": [
228291
"hvds_topo = hv.Dataset(ds.topographic__elevation)\n",
@@ -242,7 +305,9 @@
242305
{
243306
"cell_type": "code",
244307
"execution_count": null,
245-
"metadata": {},
308+
"metadata": {
309+
"collapsed": true
310+
},
246311
"outputs": [],
247312
"source": [
248313
"%opts Image style(interpolation='bilinear', cmap='viridis') plot[colorbar=True]\n",
@@ -362,7 +427,9 @@
362427
{
363428
"cell_type": "code",
364429
"execution_count": null,
365-
"metadata": {},
430+
"metadata": {
431+
"collapsed": true
432+
},
366433
"outputs": [],
367434
"source": [
368435
"imshow_grid(mg2, 'topographic__elevation', cmap='viridis')"
@@ -380,7 +447,9 @@
380447
{
381448
"cell_type": "code",
382449
"execution_count": null,
383-
"metadata": {},
450+
"metadata": {
451+
"collapsed": true
452+
},
384453
"outputs": [],
385454
"source": [
386455
"volcanic_deposits = np.zeros(mg2.size('node'))\n",
@@ -436,6 +505,7 @@
436505
"cell_type": "code",
437506
"execution_count": null,
438507
"metadata": {
508+
"collapsed": true,
439509
"scrolled": true
440510
},
441511
"outputs": [],
@@ -453,7 +523,9 @@
453523
{
454524
"cell_type": "code",
455525
"execution_count": null,
456-
"metadata": {},
526+
"metadata": {
527+
"collapsed": true
528+
},
457529
"outputs": [],
458530
"source": [
459531
"hvds_topo2 = hv.Dataset(ds2.topographic__elevation)\n",

0 commit comments

Comments
 (0)