-
Notifications
You must be signed in to change notification settings - Fork 7
/
README.Rmd
262 lines (204 loc) · 10.3 KB
/
README.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
---
output: github_document
editor_options:
chunk_output_type: console
---
paleocar
========
`paleocar` is an *R* package implementing functions to perform spatio-temporal paleoclimate reconstruction from tree-rings using the CAR (Correlation Adjusted corRelation) approach of Zuber and Strimmer as implemented in the [`care` package](https://CRAN.R-project.org/package=care) for *R*. It is optimized for speed and memory use.
This is based on the approach used in Bocinsky and Kohler (2014):
Bocinsky, R. K. and Kohler, T. A. (2014). A 2,000-year reconstruction of the rain-fed maize agricultural niche in the US Southwest. *Nature Communications*, 5:5618. doi: [10.1038/ncomms6618](http://www.nature.com/ncomms/2014/141204/ncomms6618/full/ncomms6618.html).
The primary difference between the latest version of `paleocar` and that presented in Bocinsky and Kohler (2014) is, here, model selection is performed by minimizing the corrected Akaike's Information Criterion.
A more recent reference would be Bocinsky et al. (2016):
Bocinsky, R. K., Rush, J., Kintigh, K. W., and Kohler, T. A. (2016). Exploration and exploitation in the macrohistory of the pre-Hispanic Pueblo Southwest. *Science Advances*, 2:[e1501532](http://advances.sciencemag.org/content/2/4/e1501532).
This package has been built and tested on a source (Homebrew) install of *R* on macOS 10.12 (Sierra), and has been successfully run on Ubuntu 14.04.5 LTS (Trusty), Ubuntu 16.04.1 LTS (Xenial) and binary installs of *R* on Mac OS 10.12 and Windows 10.
### Development
+ [Kyle Bocinsky](http://bocinsky.io) - Crow Canyon Archaeological Center, Cortez, CO
### Install `paleocar`
+ Development version from GitHub:
```r
install.packages("devtools")
devtools::install_github("bocinsky/paleocar")
library(paleocar)
```
+ Linux (Ubuntu 14.04.5 or 16.04.1):
First, in terminal:
```bash
sudo add-apt-repository ppa:ubuntugis/ppa -y
sudo apt-get update -q
sudo apt-get install libssl-dev libcurl4-openssl-dev netcdf-bin libnetcdf-dev gdal-bin libgdal-dev
```
Then, in R:
```r
update.packages("survival")
install.packages("devtools")
devtools::install_github("bocinsky/paleocar")
library(paleocar)
```
### Demonstration
This demo script is available in the `/inst` folder at the location of the installed package.
#### Load `paleocar` and set a working directory
```{r, echo = TRUE}
library(paleocar)
library(magrittr) # The magrittr package enables piping in R.
library(ggplot2)
# Set a directory for testing
testDir <- "./paleocar_test/"
# and create it if necessary
dir.create(testDir, showWarnings=F, recursive=T)
```
#### Load test datasets
`paleocar` ships with test files defining a study area (Mesa Verde National Park), and pre-extracted data from the International Tree Ring Databank using the [`FedData` package](https://github.com/bocinsky/FedData). See the `data-raw/data.R` script (or the documentation for `FedData`) to learn how to download these data.
```{r, echo = TRUE}
# Load spatial polygon for the boundary of Mesa Verde National Park (MVNP) in southwestern Colorado:
data(mvnp)
# Get Tree-ring data from the ITRDB for 10-degree buffer around MVNP
data(itrdb)
# Get 1/3 arc-second PRISM gridded data for the MVNP north study area (water-year [October--September] precipitation, in millimeters)
data(mvnp_prism)
```
#### Run `paleocar`
`paleocar` can be run for either single location given by a vector of annualized climate data, a matrix of locations, or over gridded climate data such as PRISM in raster format. There are three primary functions:
- `paleocar_models()` calculates the CAR-ranked linear models for all reconstructions
- `predict_paleocar_models()` generates climate predictions over a specified prediction period, and
- `uncertainty_paleocar_models()` generates an estimate of model uncertainty over a specified prediction period.
Finally, the `paleocar()` method is a convenience wrapper that runs all three of these functions and returns a list with their output. See the documentation for each function for details.
##### `paleocar` reconstruction for a single location
`paleocar` may be run for a single location by providing a vector of annualized values to be reconstructed. Simply provide a numeric vector the same length as your calibration years as the `predictands` parameter.
```{r, echo = TRUE}
# Extract a vector of annualized climate data (the first cell in the raster)
mvnp_prism.vector <- mvnp_prism[1][1,]
test.vector <- paleocar_models(predictands = mvnp_prism.vector,
chronologies = itrdb,
calibration.years = 1924:1983,
prediction.years = 1:2000,
verbose = T)
# Generate predictions and uncertainty (and plot timeseries of each)
test.prediction <- predict_paleocar_models(models = test.vector,
prediction.years = 600:1299)
test.prediction %>%
ggplot(aes(x = year,
y = Prediction)) +
geom_ribbon(aes(ymin = Prediction - `PI Deviation`,
ymax = Prediction + `PI Deviation`),
color = NA,
fill = "dodgerblue") +
geom_line(size = 0.2)
```
##### `paleocar` reconstruction for multiple locations using the same set of predictors (in this case, tree-ring chronologies)
Running `paleocar` on a matrix of locations (`predictands`) will generate reconstructions that select from
the same set of predictors (`chronologies`). The matrix must be formatted such that each location is in a column, and each row is a year of data. Note that the number of rows of the matrix must be the same as the
number of years provided to `calibration.years`.
```{r, echo = TRUE}
# Extract a matrix of annualized climate data (all cells in the raster)
mvnp_prism.matrix <- mvnp_prism %>%
raster::as.matrix() %>%
t()
test.matrix <- paleocar_models(predictands = mvnp_prism.matrix,
chronologies = itrdb,
calibration.years = 1924:1983,
prediction.years = 1:1985,
verbose = T)
# Generate predictions and uncertainty (and plot location means in uncertainty)
test.prediction <- predict_paleocar_models(models = test.matrix,
prediction.years = 600:1299)
test.prediction %>%
dplyr::mutate(cell = as.factor(cell)) %>%
dplyr::filter(cell %in% c(1,200,400,600)) %>%
ggplot(aes(x = year,
y = `Prediction (scaled)`)) +
geom_ribbon(aes(ymin = `Prediction (scaled)` - `PI Deviation (scaled)`,
ymax = `Prediction (scaled)` + `PI Deviation (scaled)`,
fill = cell),
color = NA) +
geom_line(size = 0.2) +
facet_wrap(~cell, nrow = 2) +
xlab("Year CE")
```
##### `paleocar` reconstruction over a grid
Paleocar can also be performed over a gridded climate dataset such as PRISM, so long as it is a `RasterStack` or `RasterBrick` as defined in the [`raster` package for *R*](https://CRAN.R-project.org/package=raster). Results will be returned in `RasterBrick` format.
```{r, echo = TRUE}
# Print to show format
mvnp_prism
test.raster <- paleocar_models(predictands = mvnp_prism,
chronologies = itrdb,
calibration.years = 1924:1983,
prediction.years = 600:1299,
verbose = T)
# Generate predictions and errors
test.raster.predictions <- predict_paleocar_models(models = test.raster,
prediction.years = 600:1299)
test.raster.predictions$`Prediction (scaled)` %>%
raster::mean() %>%
raster::plot()
# test.raster.predictions$`PI Deviation (scaled)` %>%
# raster::mean() %>%
# raster::plot()
```
##### `paleocar()` convenience wrapper
The `paleocar()` convenience wrapper returns a list containing the `models`, `reconstructions`, and `uncertainty`. The `paleocar()` method also automatically saves the output of `predict_paleocar_models()` and `errors_paleocar_models()`. Pass variables through this function to other ones (e.g., `meanVar = "chained"`).
```{r, echo = TRUE}
# Generate models and perform the reconstruction and error predictions.
mvnp_models <- paleocar_models(predictands = mvnp_prism,
label = "mvnp_prism",
chronologies = itrdb,
calibration.years = 1924:1983,
prediction.years = 600:1299,
out.dir = testDir,
force.redo = T,
verbose = T)
mvnp_recon <- paleocar(models = mvnp_models,
predictands = mvnp_prism,
label = "mvnp_prism",
chronologies = itrdb,
calibration.years = 1924:1983,
prediction.years = 600:1299,
out.dir = testDir,
force.redo = T,
verbose = T)
mvnp_recon <- paleocar(predictands = mvnp_prism,
label = "mvnp_prism",
chronologies = itrdb,
calibration.years = 1924:1983,
prediction.years = 600:1299,
out.dir = testDir,
force.redo = T,
verbose = T)
# Examine the structure of the output
str(mvnp_recon,
max.level = 2)
```
You can quickly load a prior reconstruction by setting `force.redo = FALSE`:
```{r, echo = TRUE}
# Generate models and perform the reconstruction and error predictions.
mvnp_recon <- paleocar(predictands = mvnp_prism,
label = "mvnp_prism",
chronologies = itrdb,
calibration.years = 1924:1983,
prediction.years = 600:1299,
out.dir = testDir,
force.redo = F,
verbose = T)
```
#### Plot results
```{r, echo = TRUE}
mvnp_recon$predictions$Prediction %>%
raster::mean() %>%
raster::plot()
mvnp_recon$predictions$`PI Deviation` %>%
raster::mean() %>%
raster::plot()
mvnp_recon$predictions$`Prediction (scaled)` %>%
raster::mean() %>%
raster::plot()
mvnp_recon$predictions$`PI Deviation (scaled)` %>%
raster::mean() %>%
raster::plot()
```
```{r, echo = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "README-"
)
```