-
Notifications
You must be signed in to change notification settings - Fork 0
/
admtools.Rmd
330 lines (235 loc) · 10.9 KB
/
admtools.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
---
title: "admtools: Getting Started"
output: rmarkdown::html_vignette
author: "Niklas Hohmann"
vignette: >
%\VignetteIndexEntry{admtools: Getting Started}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
```{r setup, echo=FALSE}
library(admtools)
```
# Introduction
This vignette is an introduction to the *admtools* package.
## Installation
### From GitHub
To install the package from *GitHub*, first install the *remotes* package:
```{r, eval=FALSE}
install.packages("remotes")
```
Then run
```{r, eval=FALSE}
remotes::install_github(repo = "MindTheGap-ERC/admtools",
build_vignettes = TRUE,
ref = "HEAD",
dependencies = TRUE)
```
to install the latest stable release. To install the version under development, use
```{r, eval=FALSE}
remotes::install_github(repo = "MindTheGap-ERC/admtools",
build_vignettes = TRUE,
ref = "dev",
dependencies = TRUE)
```
### From CRAN
To install the package from CRAN, run
```{r, eval=FALSE}
install.packages("admtools")
```
## Loading the package
Load the package using
```{r, eval=FALSE}
library(admtools)
```
This makes all functions in the package available.
## Getting help
Use
```{r, eval=FALSE}
help(package = "admtools")
```
to get an overview of the available help pages of the package, and
```{r, eval=FALSE}
?admtools
```
to view a simple help page for the package.
Vignettes are a long form of package documentation that provide more detailed examples. To list the available vignettes, use
```{r, eval=FALSE}
browseVignettes(package = "admtools") # opens in Browser
#or
vignette(package = "admtools")
```
## The adm and multiadm classes
The *admtools* package defines three main classes: `adm,` `sac` and `multiadm`. The class `adm` represent a single age-depth model, from which information can be extracted (e.g. completeness, number of hiatuses, etc.) and that can be used to transform data between the stratigraphic domain and time domain. The `multiamd` class is a list of `adm` objects. `multiadm` objects are used to represent uncertainties of age-depth models. Objects of type `sac` are sediment accumulation curves that can contain information on erosion, and can be turned into age-depth models.
## Conventions
In contrast to its name, the *admtools* package currently deals with time and height instead of age and depth. In this sense, the age-depth models are time-height models. Both time and height can be negative values. To handle ages, use time before the present. To handle depths, use height below a point of reference (e.g., the sediment surface).
## Example
This example explains the construction and application of `adm` objects. As example data we use outputs from CarboCAT Lite, a model of carbonate platform growth (Burgess 2013, 2023). This data is automatically loaded in the background by the package. To get some info of the data use
```{r, eval=FALSE}
?CarboCATLite_data
```
### Defining age-depth models
The standard constructor for age-depth models is `tp_to_adm` ("tiepoint to age-depth model"). It returns an objecto of class `adm`. This object combines information of stratigraphic heights and times and erosive interval. It allows to transform data between the stratigraphic and the time domain, and identify which data is destroyed due to hiatuses.
As example, I use the timing and stratigraphic positions of tie points taken from CarboCAT Lite to construct an age-depth model, and use the option to directly associate length and time units with it.
```{R}
# see ?tp_to_adm for detailed documentation
my_adm = tp_to_adm(t = CarboCATLite_data$time_myr,
h = CarboCATLite_data$height_2_km_offshore_m,
L_unit = "m",
T_unit = "Myr")
```
This age-depth model represents the relationship between elapsed model time and accumulated sediment 2 km offshore in a synthetic carbonate platform.
## Representation
Typing the name `my_adm` in the console will only tell that the generated variable is an age-depth model
```{r}
my_adm
```
To get a quick overview of the properties of `my_adm`, use `summary`:
```{R}
summary(my_adm)
```
If you want to inspect the insides of the object, use `str`:
```{r}
str(my_adm)
```
You can manually manipulate the fields of the `adm` object by treating it like a list. I do not recommend doing so, as it might result in unexpected downstream behavior.
You can plot `adm` objects via the standard `plot` function. Here, I use the option to highlight hiatuses in red, and increase the linw width of the conservative ( = non-destructive) intervals.
```{R}
# see ?plot.adm for plotting options for adm objects
plot(my_adm,
col_destr = "red",
lwd_acc = 2)
T_axis_lab() # plot time axis label, see ?T_axis_lab for details
L_axis_lab() # plot height axis label, see ?L_axis_lab for details
```
You can also plot sedimentation rates in the time domain using `plot_sed_rate_t`":
```{r}
plot_sed_rate_t(my_adm)
```
For more information on the extraction of sedimentation rates xin the time domain see the functions `sed_rate_t` and `sed_rate_t_fun`. Sedimentation rates in the length domain can be extracted using `sed_rate_l` and `sed_rate_l_fun` and plotted via `plot_sed_rate_l`. In addition, condensation (time preserved per stratigraphic increment) can be examined using the functions `condensation`, `condensation_fun` and `plot_condensation`.
### Extracting data from age-depth models
Use the functions `get_total_duration`, `get_total_thickness`, `get_completeness`, and `get_hiat_no` to extract information:
```{r}
get_total_duration(my_adm) #total time covered by the age-depth model
get_total_thickness(my_adm) # total thickness of section represented by the adm
get_completeness(my_adm) # stratigraphic completeness as proportion
get_incompleteness(my_adm) # stratigraphic incompleteness (= 1- strat. incompleteness)
get_hiat_no(my_adm) # number of hiatuses
```
For more detailed information, you can use
- `get_hiat_duration` to get a vector of hiatus durations
- `get_hiat_list` to get a list of hiatus positions and duration.
For example, to plot a histogram of hiatus durations, use
```{r}
hist(x = get_hiat_duration(my_adm),
freq = TRUE,
xlab = "Hiatus duration [Myr]",
main = "Hiatus duration 2 km offshore")
```
The function `is_destructive` can be used to examine whether points in time coincide with hiatuses:
```{R}
is_destructive(my_adm,
t = c(0.1,0.5))
```
### Transforming data between time and stratigraphic domain
#### Heights and times
The functions `get_height` and `get_time` are the workhorses to transform data using age-depth models.
- `get_time` takes and `adm` object and vector of heights `h` (stratigraphic positions), and returns a vector of times
- `get_height` takes an `adm` object and vector of times `t` and returns a vector of associated stratigraphic positions
As example, say we want to know the time of deposition of the following stratigraphic positions:
```{r}
h = c(30,120) # stratigraphic positions
get_time(my_adm,
h = h)
```
Conversely, to determine what parts of the section are deposited as a specific time, use
```{r}
t = c(0.2,1.4)
get_height(my_adm,
t = t)
```
Here, the `NA` indicates that the time 1.4 coincides with erosion. If you want to know the stratigraphic position of the hiatus that coincides with that time, use the option `destructive = FALSE`:
```{r}
t = c(0.2,1.4)
get_height(my_adm,
t = t,
destructive = FALSE)
```
#### Phylogenetic trees
The `admtools` package can transform complex objects between the time and stratigraphic domain. This is done using the functions `strat_to_time` and `time_to_strat`.
As an example, we transform a chronogram (a phylogenetic tree where branch length represents time). For this, we first generate a tree following the birth-death model using the `ape` package.
```{r}
#install.packages("ape") Package for analyses of phylogenetics and evolution
# see ?ape::rlineage for help
set.seed(1)
tree_in_time = ape::rlineage(birth = 1.8,
death = 0.2,
Tmax = 2)
plot(tree_in_time) # see also ?ape::plot.phylo
axis(1)
mtext("Time [Myr]", side = 1, line = 2.5)
```
You can transform the tree using `time_to_strat`:
```{r}
tree_in_strat_domain = time_to_strat(obj = tree_in_time,
x = my_adm)
```
Plotting the resulting tree along the tratigraphic column shows how the evolutinoary relationships would be preserved 2 km offshore in the simulated carbonate platform:
```{r}
plot(tree_in_strat_domain, direction = "upwards")
axis(side = 2)
mtext("Stratigraphic Height [m]",
side = 2,
line = 2)
```
#### Lists and time/stratigraphic series
*admtools* can transfrom lists from the time to the height domain and vice versa given they have elements with names `h` or `t`. These lists can be interpreted as time/stratigraphic series, where times and stratigraphic positions are associated with measured values. Note that these are not `t` objects as used by the `stats` package, as they will be generally heterodistant due to hte irregular nature of the age-depth transformation.
As example, we simulate trait evolution over 2 Myr using a Brownian motion, and transform the simulation into the stratigraphic domain.
```{r}
t = seq(0, 2, by = 0.001) # times
BM = function(t){
#" Simulate Brownian motion at times t
li = list("t" = t,
trait_val = cumsum(c(0, rnorm(n = length(t) - 1, mean = 0, sd = sqrt(diff(t))))))
return(li)
}
evo_list = BM(t)
plot(x = evo_list$t,
y = evo_list$trait_val,
xlab = "Time [Myr]",
ylab = "Trait Value",
type = "l")
```
```{r}
strat_list = time_to_strat(obj = evo_list,
x = my_adm)
plot(x = strat_list$h,
y = strat_list$trait_val,
type = "l",
xlab = "Stratigraphic Height [m]",
ylab = "Trait Value",
main = "Trait Evolution 2 km Offshore")
```
Note the jump in traits generated by the erosional interval in `my_adm`.
## Further information
For an overview of the structure of the `admtools` package and the classes used therein see
```{r, eval=FALSE}
vignette("admtools_doc")
```
For information on estimating age-depth models from sedimentation rates, see
```{r, eval=FALSE}
vignette("adm_from_sedrate")
```
For information on estimating age-depth models from tracer contents of rocks and sediments, see
```{r, eval=FALSE}
vignette("adm_from_trace_cont")
```
## References
- Burgess, Peter. "CarboCAT: A cellular automata model of heterogeneous carbonate strata." Computers & geosciences 53 (2013): 129-140. [DOI: 10.1016/j.cageo.2011.08.026](https://doi.org/10.1016/j.cageo.2011.08.026)
- Burgess, Peter. "CarboCAT Lite v1.0.1". Zenodo 2023. [DOI: 10.5281/zenodo.8402578](https://doi.org/10.5281/zenodo.8402578)