-
Notifications
You must be signed in to change notification settings - Fork 0
/
README.Rmd
327 lines (258 loc) · 14.5 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
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
---
output:
md_document:
variant: markdown_github
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r setup, echo=FALSE}
library(knitr)
opts_chunk$set(dev="svg", fig.path='README_files/')
```
# mptools
[![Travis-CI Build Status](https://travis-ci.org/johnbaums/mptools.svg?branch=master)](https://travis-ci.org/johnbaums/mptools)
[![CRAN](http://www.r-pkg.org/badges/version/mptools)](https://cran.r-project.org/package=mptools)
[![Downloads](http://cranlogs.r-pkg.org/badges/mptools?color=brightgreen)](http://www.r-pkg.org/pkg/mptools)
In recent years, spatially-explicit, coupled habitat-demographic models have
been developed to acknowledge the importance of population processes when
evaluating the impacts of habitat change (Keith *et al.* 2008, 2014; Southwell
*et al.* 2008; Fordham *et al.* 2012; Cadenhead *et al.* 2015). Frequently,
these coupled models use the Metapop module of RAMAS GIS (Akçakaya & Root 2005)
to relate spatially-dynamic habitat suitability to carrying capacity, and then
simulate metapopulation dynamics. Metapop saves settings, parameters values, and
simulation results to a text file (with .mp extension). Unfortunately, there is
limited flexibility within Metapop regarding presentation of results, and
extraction of results via the GUI is awkward.
**The `mptools` package provides several useful functions for extracting and
processing data from RAMAS Metapop .mp files.**
## Function overview:
Function | Description
------------- | -------------
`actions()` | Extract management action details from .mp files
`kch()` | Extract carrying capacity times series from .kch files
`knt()` | Plot the carrying capacity and abundance dynamics of one or more populations
`meta()` | Extract population details from .mp files
`mp2sp()` | Create a `SpatialPointsDataFrame` representing the centroid of each population, with attributes `pop` (population name), `time` (the time step), and `N` (the mean population size)
`mp2xy()` | Extract coordinates for populations from .mp files, and transform them to the original coordinate reference system
`mp_animate()`| Animate temporal habitat and abundance dynamics on a gridded landscape
`results()` | Extract simulation results (mean, sd, min and max of population size), as well as expected minimum abundance (EMA) and its standard deviation, from .mp files
`ths()` | Extract total habitat suitability from .ptc files.
## Installation
First of all, we install and load the `mptools` package. This requires the
`devtools` package to be present, which can be installed with
`install.packages(devtools)`.
```{r load_pkg_echo, eval=FALSE}
library(devtools)
install_github('johnbaums/mptools')
library(mptools)
```
```{r load_pkg, eval=3, echo=FALSE}
library(devtools)
install_github('johnbaums/mptools')
library(mptools)
```
## Examples
`mptools` comes packaged with some data relating to a RAMAS Metapop model. The
model is of a metapopulation of a hypothetical species, comprising 263
populations, the population dynamics of which have been simulated over 100
years. In reality, the number of populations, iterations, and time steps are
only limited by RAMAS Metapop itself, and the amount of memory available to R.
Most functions in the mptools package operate on RAMAS Metapop .mp files that
have been run and contain simulation results. An example .mp file representing
our hypothetical metapopulation is included in the package. The path to this
file is given by `system.file('example.mp', package='mptools')`. Here we assign
that path to the object `mp`.
```{r mpfile}
mp <- system.file('example.mp', package='mptools')
```
### `results()`: extract simulation results
To extract the simulation results, we use the `results` function, pointing it to
the .mp file.
```{r res, message=FALSE}
res <- results(mp=mp)
res
```
Summary information is printed by default, but the resulting list comprises
elements containing: the mean, minimum, maximum and standard deviation of
population size for each time step, given for each patch and for the
metapopulation as a whole (in element `results`); the minimum (`iter_min`),
maximum (`iter_max`) and final (`iter_terminal`) abundance across the time
horizon, for each iteration; the probability of quasi-extinction (i.e., the
proportion of iterations during which total abundance fell beneath a
predetermined threshold, `qe_prob`); the mean of minimum abundance (`EMA`;
McCarthy & Thompson 2001) and standard deviation (`SDMA`) of minimum abundance.
Additionally, some metadata about the simulation is provided: the date and time
that the model was run (`timestamp`); and the number of populations (`n_pops`),
time steps (`duration`), and iterations (`n_iters`).
For example, the mean, minimum, maximum and standard deviation of total
population size at each time step (first six time steps shown here) is given by:
```{r res_head}
head(res$results[,, 'ALL'])
```
The same could be extracted for population 40, referring to it by name:
```{r res_head40}
head(res$results[,, 'Pop 40'])
```
### `meta()`: extract Metapop settings
The `meta` function returns a `data.frame` containing information about RAMAS
Metapop settings for each population, including their names, initial abundances,
density dependence types, carrying capacities, and maximum growth rates. Here
the first few rows are shown for brevity.
```{r meta, message=FALSE}
met <- meta(mp=mp)
head(met)
```
### `mp2xy()`: extract patch coordinates
RAMAS Metapop simulations are often based on spatial grids describing the
distribution of habitat (i.e., when using the Spatial Data module to identify
patch structure). In these cases, spatial coordinates are converted by RAMAS
such that they describe the position relative to the top left corner of the
grid. Such coordinates are returned by the `meta` function in the columns
`xMetapop` and `yMetapop`. In order to relate simulation results to the true
landscape, the original (untransformed) coordinates can be recovered with
`mp2xy`. This requires one of the original grids used by Spatial Data (or its
`Raster*` object representation), and knowledge of the cell length setting
passed to that module. By default, `mp2xy` creates a plot of the points,
overlaid upon the provided raster data.
The raster grids that were originally used to define the patch structure are
included with `mptools`. We identify the file path of one of these, and pass it
to `mp2xy`, below.
```{r xy, echo=c(1, 2, 4), results='hide'}
library(raster)
r <- system.file('example_001.tif', package='mptools')
svg('README_files/xy-1.svg')
xy <- mp2xy(mp=mp, r=r, cell.length=9.975)
dev.off()
```
![](https://rawgit.com/johnbaums/mptools/master/README_files/xy-1.svg)
```{r xy_head}
head(xy)
```
Above we see that the `data.frame` returned by `mp2xy` includes the populations'
names, their RAMAS coordinates, and their original coordinates. In this case,
the original coordinates were defined by the Australian Albers equal area
coordinate system.
### `mp2sp()`: create a `SpatialPointsDataFrame` describing population dynamics
The coordinates returned by `mp2xy`, and the abundance at the corresponding
patches at each time step of the simulation, can be used to create a spatial
object that enables further analysis and visualisation. This is done with
`mp2sp`, which creates a `SpatialPointsDataFrame` with attributes containing the
average (across iterations) abundance at each patch, for each time step of the
simulation.
In the following example, we set `start` to 2000, indicating that the first time
step of the simulation corresponds to the year 2000. This controls labelling of
the shapefile's attributes---subsequent time steps are labelled with sequential
integers. Below we also pass the source and target [PROJ.4
strings](https://en.wikipedia.org/wiki/PROJ.4) describing the coordinate
reference system (CRS) for the input coordinates, and the output
`SpatialPointsDataFrame`, respectively.
```{r mp2sp, message=FALSE}
spdf <- mp2sp(mp=mp, coords=xy, start=2000,
s_p4s='+init=epsg:3577', t_p4s='+init=epsg:4326')
```
Above, we indicated that the source CRS is Australian Albers (EPSG 3577) and the
target CRS is WGS 84 (EPSG 4326). (These could also be passed as [`CRS`
objects](http://www.inside-r.org/packages/cran/sp/docs/CRS) rather than as
PROJ.4 strings.)
We can plot the output with `sp::spplot`, which can colour points according to
the value of an attribute:
```{r plot_spplot, echo=c(1, 2, 4), results='hide'}
library(sp)
library(viridis)
svg('README_files/plot_spplot-1.svg', height=3.75)
spplot(subset(spdf, time==2000), zcol='N',
cuts=c(-Inf, 0, 10^(0:6)), key.space='right',
col.regions=c('gray80', viridis(100)))
dev.off()
```
![](https://rawgit.com/johnbaums/mptools/master/README_files/plot_spplot-1.svg)
We can then write these data out to an ESRI shapefile or KML for use in, e.g., a
GIS or Google Earth. For example, to write out as a shapefile, we can use:
```{r write_shp, message=FALSE}
library(rgdal)
writeOGR(spdf, dsn=tempdir(), layer='mp', driver='ESRI Shapefile')
```
The shapefile was written out to a file called "mp.shp" (along with its
accessory files) within the current temporary directory, which can be viewed
with `browseURL(tempdir())`.
To write out to KML, we include the output file name in the string passed to the
`dsn` argument, and specify the KML driver. Note that below, we subset the data
to the first time step (2000), to reduce the size of the output.
```{r write_kml}
library(rgdal)
writeOGR(subset(spdf, time==2000),
dsn=file.path(tempdir(), 'mp.kml'),
layer='', driver='KML')
```
If Google Earth is installed on your system, and if .kml files are associated
with it, you can open the new .kml file with:
```{r show_kml, eval=FALSE}
file.show(file.path(tempdir(), 'mp.kml'))
```
Note that to display correctly in Google Earth, the data should be in the WGS 84
datum. When using `mp2xy` above, we ensured that `spdf` was transformed to WGS
84.
### `kch()`: extract carrying capacity dynamics from .kch files
If a model involves carrying capacity dynamics, RAMAS stores information about
these in .kch files. Changes in carrying capacity can be extracted from these
files with the `kch` function. Rather than show this as text output, we see
below how these carrying capacity sequences are used by `knt`, along with the
mean abundances returned by `results`.
```{r kch}
k <- kch(meta=met, path=dirname(mp))
str(k)
```
### `knt()`: visualise carrying capacity and abundance dynamics
The `knt` function plots the change in carrying capacity and mean abundance
through time, for each population, or for a given subset of populations. Below,
we show these dynamics for four chosen populations. Black lines show mean total
abundance, and grey lines show carrying capacity. Populations whose names are
shown in bold on a blue background had positive abundance at the first time
step.
```{r knt, message=FALSE, echo=2, results='hide'}
svg('README_files/knt-1.svg')
knt(meta=met, kch=k, pops=c('Pop 169', 'Pop 170', 'Pop 174', 'Pop 175'),
show_N=TRUE, results=res, samelims=TRUE, layout=c(2, 2))
dev.off()
```
![](https://rawgit.com/johnbaums/mptools/master/README_files/knt-1.svg)
### `mp_animate()`: animate habitat and abundance dynamics
With `mp_animate`, we can create a gif animation showing temporal dynamics in
habitat suitability and mean abundance. This can reveal lags in response to
changing habitat suitability.
The function requires a `RasterStack` or `RasterBrick` with layers, describing
habitat change, in temporal order. The interval between successive rasters (i.e.,
the interval between time steps) is assumed to be constant. Below, we create a
`RasterStack` rasters included with `mptools`, and pass it to `mp_animate`.
```{r animate, message=FALSE, echo=FALSE}
library(raster)
tifs <- list.files(system.file(package='mptools'), '\\.tif$', full.names=TRUE)
spdf <- mp2sp(mp=mp, coords=xy, start=2000)
mp_animate(spdf, habitat=stack(tifs), outfile='README_files/dynamics.gif',
zlim=c(0, 800), width=630, height=615, overwrite=TRUE)
```
```{r animate_echo, message=FALSE, eval=FALSE}
library(raster)
tifs <- list.files(system.file(package='mptools'), '\\.tif$', full.names=TRUE)
spdf <- mp2sp(mp=mp, coords=xy, start=2000)
mp_animate(spdf, habitat=stack(tifs), outfile='README_files/dynamics.gif',
zlim=c(0, 800), width=630, height=615)
```
![](README_files/dynamics.gif)
Points indicate populated patches, with colours scaled linearly from white
(mean abundance = 1) to black (mean abundance is equal to the maximum mean
abundance across populations and years). Points are not shown when mean
abundance is zero. The raster is coloured according to the carrying capacity of
the cell.
### References
Akçakaya, H. & Root, W. (2005). *RAMAS GIS: Linking spatial data with population viability analysis (version 5.0)*. Applied Biomathematics, Setauket, New York.
Cadenhead, N.C., Kearney, M.R., Moore, D., McAlpin, S. & Wintle, B.A. (2015). Climate and fire scenario uncertainty dominate the evaluation of options for conserving the great desert skink. *Conservation Letters*.
Fordham, D.A., Resit Akçakaya, H., Araújo, M.B., Elith, J., Keith, D.A., Pearson, R., Auld, T.D., Mellin, C., Morgan, J.W. & Regan, T.J. (2012). Plant extinction risk under climate change: Are forecast range shifts alone a good indicator of species vulnerability to global warming? *Global Change Biology*, **18**, 1357–1371.
Keith, D.A., Akçakaya, H.R., Thuiller, W., Midgley, G.F., Pearson, R.G., Phillips, S.J., Regan, H.M., Araújo, M.B. & Rebelo, T.G. (2008). Predicting extinction risks under climate change: Coupling stochastic population models with dynamic bioclimatic habitat models. *Biology Letters*, **4**, 560.
Keith, D.A., Mahony, M., Hines, H., Elith, J., Regan, T.J., Baumgartner, J.B., Hunter, D., Heard, G.W., Mitchell, N.J. & Parris, K.M. (2014). Detecting extinction risk from climate change by IUCN red list criteria. *Conservation Biology*, **28**, 810–819.
McCarthy, M.A. & Thompson, C. (2001). Expected minimum population size as a measure of threat. *Animal Conservation*, **4**, 351–355.
Southwell, D., Lechner, A., Coates, T. & Wintle, B. (2008). The sensitivity of population viability analysis to uncertainty about habitat requirements: Implications for the management of the endangered southern brown bandicoot. *Conservation Biology*, **22**, 1045–1054.
[//]: # (Citations styled using the Methods in Ecology and Evolution style
written by Xiaodong Dang and provided by the Citation Style Language project
under the Creative Commons Attribution-ShareAlike 3.0 Unported license
[http://creativecommons.org/licenses/by-sa/3.0/ license] -
http://citationstyles.org/)