---
title: "ausplotsR: quickstart guide to basic analysis of TERN AusPlots vegetation data"
author: "Greg Guerin & Bernardo Blanco-Martin"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{ausplotsR: quickstart guide to basic analysis of TERN AusPlots vegetation data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  tidy.opts=list(width.cutoff=80), 
  tidy = 'formatR', 
  comment = "#>"
)
```

## Introduction

TERN AusPlots is a national plot-based terrestrial ecosystem surveillance monitoring method and dataset for Australia (Sparrow et al. 2020). Through **ausplotsR**, users can directly access AusPlots data collected by on-ground observers on vegetation and soils, including physical sample/voucher details and barcode numbers. The dataset can be downloaded in its entirety or as individual modules, or as a subset by geographic bounding box or species name or site search. The package also includes a series of bespoke functions for working with AusPlots data, including visualisation, creating tables of species composition, and calculation of tree basal area, fractional cover or vegetation cover by growth form/structure/strata and so on.

This is a short guide for getting started with analysis of AusPlots data through the **ausplotsR** R package. More information on making use of AusPlots data in **ausplotsR** is available through the package help files and manual. Below, we demonstrate installing the package, accessing some AusPlots data, generating matrices and running simple example analyses.

More comprehensive (but not necessarily updated)  tutorials on accessing and analysing AusPlots data (Blanco-Martin 2019) are available at:
[https://github.com/ternaustralia/TERN-Data-Skills/tree/master/EcosystemSurveillance_PlotData](https://github.com/ternaustralia/TERN-Data-Skills/tree/master/EcosystemSurveillance_PlotData)


## Installing the package and accessing raw data

The latest version of **ausplotsR** can be installed from CRAN or directly from github using the `devtools` package, which must be installed first.

```
library(devtools)
install_github("ternaustralia/ausplotsR", build_vignettes = TRUE, dependencies = TRUE)
```

Once installed, load the package as follows. Note, packages ggplot2 and mapdata are required for **ausplotsR** to load, and functions are also imported from packages: vegan, plyr, R.utils, httr, jsonlite, gtools, jose, curl, r2r, stringr, and betapart, while knitr and rmarkdown are required to build this package vignette (i.e., if 'build_vignettes' is set to TRUE above).

```{R, warning=FALSE, message=FALSE, error=FALSE}
library(ausplotsR)
```

We can now access live data, starting here with basic site information and vegetation point-intercept modules and using a bounding box to spatially filter the dataset to central Australia. All data modules are extracted via a single function, `get_ausplots`:

```{r, echo=FALSE}
oldpar <- par(no.readonly = TRUE)
```

```{r}
#See ?get_ausplots to explore all data modules available
```


```{r}
my.ausplots.data <- try(get_ausplots(veg.PI=TRUE, veg.vouchers=TRUE, bounding_box = c(125,140,-40,-10)))
```

```{r, echo=FALSE}
if(!inherits(my.ausplots.data, "list")) {
  message("Vignette aborted due to database connection issue.")
  knitr::knit_exit()
}

```

The output of the above call is a list with the following `$elements`:

```{r}
names(my.ausplots.data)
```

The 'site.info' table contains basic site and visit details. Here are a selected few of the many fields:
```{r}
head(my.ausplots.data$site.info[,c("site_location_name", "site_unique", "longitude", "latitude", "bioregion_name")])
```

**Each survey is identified by the *'site_unique'* field, which is unique combination of site ID ('site_location_name') and visit ID ('site_location_visit_id'). The *'site_unique'* field therefore links all tables returned from the `get_ausplots` function**.

The 'site.info' table and can be used to identify, subset or group surveys in space and time, for example:

```{r}
#count plot visits per Australian States:
summary(as.factor(my.ausplots.data$site.info$state))
```


## Map AusPlots sites and visualise data

The package has an in-built function - see `?ausplots_visual` - to rapidly map AusPlots over Australia and to visualise the relative cover/abundance of green vegetation, plant growth forms and species. Maps can also be generated manually using the longitude and latitude fields in the `$site.info` table.

```{r fig1, fig.height = 4, fig.width = 6}
#Sites are coded by IBRA bioregion by default. 
map_ausplots(my.ausplots.data)
```

Alternatively, the following call generates a pdf with a map of all sites and attribute graphics for selected AusPlots:
`ausplotsR::ausplots_visual()`


Here is a snippet of the raw point-intercept data that will be used in the following examples to derive vegetation attributes:

```{r}
head(subset(my.ausplots.data$veg.PI, !is.na(herbarium_determination)))
```

Note that 'veg_barcode' links species hits to the vegetation vouchers module, while the 'hits_unique' field identifies the individual point-intercept by transect and point number (see `help(ausplotsR)` and references for more details on the plot layout and survey method). At each point, plant species (if any), growth form and height are recorded along with substrate type.

## Example 1: latitudinal pattern in proportional vegetation cover

Let's visualise basic vegetation cover as a function of latitude. First, we call the `fractional_cover` function on the extracted point-intercept data (`$veg.PI`). The function converts the raw data to proportional cover of green/brown vegetation and bare substrate. Note the calculation may take a few minutes for many AusPlots, so for this example we will pull out a subset of 100 randomly drawn sites to work with.

```{r, warning=FALSE}
sites100 <- my.ausplots.data$veg.PI[which(my.ausplots.data$veg.PI$site_unique  %in% sample(my.ausplots.data$site.info$site_unique, 100)), ]
my.fractional <- fractional_cover(sites100)

head(my.fractional)
```

Next, we need to merge the fractional cover scores with longlat coordinates from the site information table. We use the 'site_unique' field (unique combination of site and visit IDs) to link tables returned from the `get_ausplots` function:

```{r}
my.fractional <- merge(my.fractional, my.ausplots.data$site.info, by="site_unique")[,c("site_unique", "bare", "brown", "green", "other", "longitude", "latitude")]

my.fractional <- na.omit(my.fractional)

head(my.fractional)
```

Now we can plot out the continental relationship, e.g., between the proportion of bare ground with no kind of vegetation cover above and latitude.

```{r, fig.height = 4, fig.width = 6}
plot(bare ~ latitude, data=my.fractional, pch=20, bty="l")
```

There appears to be a hump-backed relationship, with a higher proportion of bare ground in the arid inland at mid-latitudes. We can add a simple quadratic model to test/approximate this:

```{r, fig.height = 4, fig.width = 6}
my.fractional$quadratic <- my.fractional$latitude^2

LM <- lm(bare ~ latitude + quadratic, data=my.fractional)
summary(LM)

#generate predicted values for plotting:
MinMax <- c(min(my.fractional$latitude), max(my.fractional$latitude))
ND <- data.frame(latitude=seq(from=MinMax[1], to=MinMax[2], length.out=50), quadratic=seq(from=MinMax[1], to=MinMax[2], length.out=50)^2)
ND$predict <- predict(LM, newdata=ND)
#
plot(bare ~ latitude, data=my.fractional, pch=20, bty="n")
points(ND$latitude, ND$predict , type="l", lwd=2, col="darkblue")
```

## Example 2: Species by sites table

Aside from 'gross' values from plots such as fractional cover, many analyses in community ecology begin with species abundance information. With **ausplotsR** you can generate this easily from the more complex vegetation point-intercept data. The first step to work with species-level AusPlots data is to create a species occurrence matrix. The `species_table` function in the **ausplotsR** package can be used to create this type of matrix. This function takes a data frame of individual raw point-intercept hits (i.e. a `$veg.PI` data frame) generated using the `get_ausplots` function and returns a ‘species against sites’ matrix:

```{r}
#The species_table function below can also take the `$veg.voucher` module as input, but `m_kind="PA"` must be specified to get a sensible presence/absence output.
```


```{r}
#The 'species_name' argument below specifies use of the "standardised_name" field to identify species, which is based on herbarium_determination names (i.e., "HD" option in species_name) matched to accepted scientific name according to a standard (APC: https://www.anbg.gov.au/cpbr/program/hc/hc-APC.html).
```


```{r}
my.sppBYsites <- species_table(my.ausplots.data$veg.PI, m_kind="percent_cover", cover_type="PFC", species_name="SN")

#check the number of rows (plots) and columns (species) in the matrix
dim(my.sppBYsites)

#look at the top left corner (as the matrix is large)
my.sppBYsites[1:5, 1:5] 
```

We can crudely pull out the 10 highest ranking species in terms of their percent cover cumulative across all plots they occur in:

```{r}
rev(sort(colSums(my.sppBYsites)))[1:10]
```

A simple example of downstream visualisation and analysis of species-level AusPlots data is Rank-Abundance Curves (also known as Whittaker Plots). Rank-Abundance Curves provide further information on species diversity. They provide a more complete picture than a single diversity index. Their x-axis represents the abundance rank (from most to least abundant) and in the y-axis the species relative abundance. Thus, they depict both Species Richness and Species Evenness (slope of the line that fits the rank; steep gradient indicates low evenness and a shallow gradient high evenness).


```{r, fig.height = 4, fig.width = 6}
#Whittaker plots for some selected AusPlots with alternative relative abundance models fitted to the plant community data:
par(mfrow=c(2,2), mar=c(4,4,1,1))
for(i in c(1:4)) {
  plot(vegan::radfit(round(my.sppBYsites[9+i,], digits=0), log="xy"), pch=20, legend=FALSE, bty="l")
  legend("topright", legend=c("Null", "Preemption", "Lognormal", "Zipf", "Mandelbrot"), lwd=rep(1, 5), col=c("black", "red", "green", "blue", "cyan"), cex=0.7, bty="n")
}
```


## Example 3: Quick species lists

Perhaps you simply want to browse which plant species have been recorded in AusPlots, without all the associated raw data? Here, the `species_list` function is your friend:

```{r}
#The species_list function is designed to take $veg.voucher as input but can also take $veg.PI
```


```{r}
#print a list of genus_species-only records from selected plots (for demonstration we print only part):
species_list(my.ausplots.data$veg.vouch, grouping="by_site", species_name="GS")[1:2]

#overall species list ordered by family (for demonstration we print only part):
species_list(my.ausplots.data$veg.vouch, grouping="collapse", species_name="SN", append_family=TRUE)[1:50]
```

```{r, echo=FALSE}
par(oldpar)
```

## Explore TERN AusPlots

In addition to the key site info and vegetation point-intercept modules introduced above, `get_ausplots` is your gateway to raw data modules for vegetation structural summaries, vegetation vouchers (covers the full species diversity observed at the plot and includes tissue sample details), basal wedge, and soils subsites, bulk density and pit/characterisation (including bulk and metagenomics soil samples).


## References

Blanco-Martin, B. (2019) Tutorial: Understanding and using the 'ausplotsR' package and AusPlots data. Terrestrial Ecology Research Network. Version 2019.04.0, April 2019. https://github.com/ternaustralia/TERN-Data-Skills/

Sparrow, B., Foulkes, J., Wardle, G., Leitch, E., Caddy-Retalic, S., van Leeuwen, S., Tokmakoff, A., Thurgate, N., Guerin, G.R. and Lowe, A.J. (2020) A vegetation and soil survey method for surveillance monitoring of rangeland environments. *Frontiers in Ecology and Evolution*, 8:157.
