---
title: "FLightR with black-tailed godwit vignette from MEE 2017"
author: "Eldar Rakhimberdiev & Julia Karagicheva"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{FLightR with black-tailed godwit vignette from MEE 2017}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


# FLightR analysis example of Black tailed godwit geolocator
_this vignette is supplementary material to_
**Rakhimberdiev E., Saveliev A., Piersma, T., Karagicheva J. 2017 FLightR: An R package for reconstructing animal paths from solar geolocation loggers. Methods in Ecology and Evolution. [DOI](https://doi.org/10.1111/2041-210X.12765)**


## Install package
Note that we used FLightR 0.4.5 version, and running this example in later versions may provide somehow different results.
To install the latest CRAN version try
```{r, eval = F}
install.packages('FLightR')
```
To run version 0.4.5 which we use in this example try 
```{r, eval = F}
library(devtools)
install_github('eldarrak/FLightR@0.4.5')
```
The latest version is available here:
```{r, eval = F}
install_github('eldarrak/FLightR')
```
Load package
```{r}
library(FLightR)
```
## 1. Data preparation and import
The procedure of the definition of the twilight events, i. e. sunrises and sunsets, is not implemented in FLightR. Please, use appropriate functions in other R packages, e.g. BAStag, twGeos or GeoLight [Lisovski et al. 2012a](#lisovski_2012a). In this software, user guides semi-automatically search for the twilight times and then visually explores the light data, manually removing twilights defined apparently wrong or around which the light pattern is atypical or systematically biased, e.g. when a bird enters or exits a cavity during twilight ([Rakhimberdiev et al. 2016](#rakhimberidev_jab_2016)). An example of the routine for the twilight definition in BAStag is available as [appendix A4](https://github.com/eldarrak/FLightR/blob/master/examples/Black-Tailed_Godwit_JAB_example/A4_BAStag_routine.Rmd) to [Rakhimberdiev et al. 2016](#rakhimberdiev_jab_2016).

FLightR works with the ‘TAGS’ format, which is an agreed general annotation of twilight data by National Centre for Ecological Analysis and Synthesis working group ‘Establishing an open-source animal-tracking analysis platform for archival geolocators https://www.nceas.ucsb.edu/workinggroups/establishing-open-source-animal-tracking-analysis-platform-archival-geolocators. A TAGS file is a CSV file containing the following fields: 

* `datetime` – date and time in ISO 8601 format e.g. 2013-06-16T00:00:11.000Z;
* `light` – light value measured by tag;
* `twilight` – assigned by the software numeric indication of whether the record belongs to sunrise (1), sunset (2) or none of those (0);
* `excluded` – indication of whether a twilight was excluded during manual inspection (logical, `TRUE | FALSE`);
* `interp` - indication of whether the light value at twilight was interpolated (logical, `TRUE | FALSE`). 
The fields `excluded` and `interp` may have values of `TRUE` only for `twilight > 0`. 
The online TAGS service saves data in the TAGS format. In the R packages GeoLight and BAStag or twGeos, the annotated twilight data need to be exported to TAGS, for which the functions in the FLightR (`GeoLight2TAGS`,  `BAStag2TAGS` or `twGeos2TAGS`) can be used.

We will use the TAGS formatted data provided as appendix A3 to [Rakhimberdiev et al. 2016](#rakhimberdiev_jab_2016) also available on the authors’ GitHub page:
```{r, eval=FALSE, tidy=FALSE}
download.file('https://git.io/vrJgv', 'example_TAGS_format.csv')
```
The function `get.tags.data` reads comma separated file in the TAGS format, detects the tag type, checks whether the light data are log-transformed, transforms them back from the log scale if needed and creates an object, containing 

1. the recorded light data,
2. the detected twilight events,
3. light level data at the moment of each determined sunrise and sunset and around them (24 fixes before and 24 after it events into an object of two lists
4. technical parameters of the tag, added automatically, unless preset by the user.

The function works with all the common tag types: mk tags (produced by British Antarctic Survey, Lotek, and Migrate Technology Ltd.), Intigeo tags (Migrate technology Ltd.) and GDL tags (Swiss Ornithological Institute).
```{r, eval=FALSE, tidy=FALSE}
Proc.data<-get.tags.data("example_TAGS_format.csv")
```
 
## 2. Calibration

Geolocators measure light levels with different precision, and calibration is needed to establish the relationship between the observed and the expected light levels for each device. This relationship is depicted by the calibration parameters (slopes), calculated using the data recorded in known (calibration) geographic positions, e. g. where the animal was tagged, recaptured or observed. FLightR uses a 'template fit' for calibration [Ekstrom 2004, 2007](#ekstrom_2004). For each tag it finds the linear (on a log-log scale) relationship between the light levels measured in known locations and the theoretical light levels, estimated from current sun angle in these locations with the deterministic equation developed by Ekstrom [Rakhimberdiev et al. 2015](#rakhimberdiev_movecol_2015).

To calculate the calibration parameters user needs to create a data frame where the geographic coordinates of the calibration location, and the start and end dates of the calibration period, i. e. the period of residence in the known location, are specified:
* `calibration.start` (POSIXct format)
* `calibration.stop` (POSIXct format)
* `lon` (numeric)
* `lat` (numeric)
The data frame contains as many rows as many distinct calibration periods the track contains. 

```{r, eval = F}
Calibration.periods<-data.frame(
        calibration.start=as.POSIXct(c(NA, "2014-05-05"), tz='GMT'),
        calibration.stop=as.POSIXct(c("2013-08-20", NA), tz='GMT'),
        lon=5.43, lat=52.93) 
		#use c() also for the geographic coordinates, 
		#if you have more than one calibration location
		# (e. g.,  lon=c(5.43, 6.00), lat=c(52.93,52.94))
print(Calibration.periods)
```
In this example, we have two calibration periods in the same location, at the beginning and at the end of the track. This is a common case, as the birds are often recaptured at the same location, where they were tagged. 

When multiple calibration locations are available, each of them has to be processed with the function `plot_slopes_by_location`.  In this case, in the `Calibration periods` data frame, each row should refer to one calibration period. Compiling the data frame with multiple calibration locations, use `c()` also for the geographic coordinates (e. g., `lon=c(5.43, 6.00), lat=c(52.93,52.94)`).

A 'Calibration' object is compiled with the function`make.calibration` from the created `Calibration periods` data frame and the `Proc.data` object. 
```{r, eval = F}
Calibration<-make.calibration(Proc.data, Calibration.periods)
```
This object contains all the calibration parameters for the tag, and it will be further used for calculation of geographic positions across the track. When there are more than one calibration periods, the parameter `model.ageing` can be set `TRUE` to account for the tag ageing. In this case, the calibration parameters are calculated, based on the assumption that the calibration slope changes linearly with time. 

### Find calibration periods for a known calibration location
The exact period of a tagged animal's stay in a known location is usually unknown, but it can be derived from the data. For this, calibration slopes for the sunset and sunrise of each day of the tracking period are calculated, based on the assumption that the tag remained in the same known position all the time. The slopes are calculated and plotted with the function `plot_slopes_by_location`.

```{r, eval = F}
plot_slopes_by_location(Proc.data=Proc.data, location=c(5.43, 52.93))
```

Looking at the plot, we can define the time periods, during which the tag resided in the calibration location (recall, that we assume that the tag remained in this location all the time). Because calibration slopes reflect the adequacy of the light level measured by the device, they vary little, in time and between sunsets and sunrises, as long as the tagged animal stays in the calibration location, but become apparently diverse, when it moves away from it. Both patterns are clearly distinguishable at the plot. 

Play with `abline()` to find the proper boundaries for the calibration.
```{r, eval = F}
abline(v=as.POSIXct("2013-08-20", tz='GMT')) # end of first calibration period
abline(v=as.POSIXct("2014-05-05", tz='GMT')) # start of the second calibration period
```
### Find a calibration location for a known calibration period
It may happen that an animal was tagged in the High Arctic under polar day conditions or that it moved far away from the capture site immediately after tagging and the roof-top calibration data are not available. Even in such cases it is still possible to obtain calibration parameters for a resident period at unknown location. FLightR approach to this problem is similar to Hill-Ekstrom calibration [
[Lisovski et al. 2012b](#lisovski_2012b) implemented in GeoLight [Lisovksi et al. 2012a](#lisovski_2012a). If bird is assumed to be resident at some period one can try:
```{r, eval=FALSE}
# ~ 15 min run time
Location<-find.stationary.location(Proc.data, '2013-07-20', '2013-08-20',
                                   initial.coords=c(10, 50))
```
The function will return geographic coordinates of the location for which the range of errors in slopes is minimal. User has to provide the initial coordinates, which should be within a few thousand kilometers from the hypothetical real location.

## 3. Assign spatial extent

The function `make.grid` sets up a spatial grid (50 X 50 km on default) delimited by user-defined boundaries: `left`, `right`. `bottom` and `top`. When the tagged animal cannot occur or stay between two subsequent twilights over particular areas, e.g. open water for a landbird or deep inland for a marine animal, additional parameters `distance.from.land.allowed.to.use` and `distance.from.land.allowed.to.stay` can be specified. Each of the parameters require a vector of two numbers: the minimal and the maximal distances (in km) from shoreline, at which the animal is allowed to occur/stay.
```{r, eval = F}
Grid<-make.grid(left=-14, bottom=30, right=13, top=57,
   distance.from.land.allowed.to.use=c(-Inf, Inf),
   distance.from.land.allowed.to.stay=c(-Inf, Inf))
```
The resulting `Grid` is a matrix with the columns: `lon` (longitude), `lat` (latitude) and `Stay` (probability of stay). The grid cells, which the animal presumably cannot use, are excluded from the data, while the locations at which an animal cannot be stationary are given a low probability of stay. Using masks can side track model estimation to the local minima, and we recommend to initially run model without a mask, enable them for the second run and visually compare the results, to see if the model converges to a similar track.

## 4. Prepare the model for run

The function `make.prerun.object` creates a complex object, which will be used in the main run. It incorporates all the objects, created at earlier steps: the light data with the detected twilight events (`Proc.data`), the spatial parameters (`Grid`), geographic coordinates of the initial location, where the tracking has started (`start`), and the calibration parameters (`Calibration`).

```{r, eval = F}
# ~ 15 min run time
all.in<-make.prerun.object(Proc.data, Grid, start=c(5.43, 52.93), Calibration=Calibration)
```
## 5. Particle filter run

At this stage, the model output is generated. It contains: a table of positions at each twilight (`$Results$Quantiles`) and their statistics (mean, median values and credible intervals), a table of parameters of the movement model (`$Results$Movement.results`) and posterior distribution at every twilight (`$Results$Points.rle`) and at every transition between twilights (`$Results$Transitions.rle`). 
Within `run.particle.filter`, the following parameters can be defined: 
* `nParticles` - number of particles (1e4 is recommended for test and 1e6 for the analysis);
* `threads` - amount of parallel threads to use for the run default is -1 that means all available except one;
* `known.last` - TRUE if you know that in the end of the logging period tag occurred in a known place (FALSE is the default option);
* `check.outliers` – FALSE by default. Set it TRUE if you wish on a fly outliers detection, we recommend to use it if the results have strong outliers.

```{r, eval = F}
nParticles=1e6
# ~ 45 min run time
Result<-run.particle.filter(all.in, threads=-1,
            nParticles=nParticles, known.last=TRUE,
		    precision.sd=25, check.outliers=F)
save(Result, file="Result.bltg.ageing.model.noOD.RData")
```
Now the results are saved in your work directory as an RData object.

## 6. Estimation of arrival and departure dates

The function `find.times.distribution` derives the time at which an animal arrived or departed from the area and provides the measure of its uncertainty. First, select grid points of interest. For example in the current data we are interested in the date when our bird left the Netherlands. We will make  a boundary at 2&deg; longitude:
```{r, eval = F}
Index<-which(Result$Spatial$Grid[,1]>(2))
```
Estimate probabilities of occurrence within the area at each twilight:
```{r, eval = F}
Arrivals.NL<-find.times.distribution(Result,Index)
Arrivals.NL
```
Print a simple summary of the periods of residence and movement:
```{r, eval = F}
stationary.migration.summary(Result, prob.cutoff = 0.1, min.stay = 3)
```
## 7. Visualisation of the results

### Plot a simple map

Plot a map with the most probable positions, i. e. combinations of the most probable latitude and longitude for each twilight:
```{r, eval = F}
map.FLightR.ggmap(Result)
```
### Plot lon lat graph

Plot the estimated for each twilight probabilities of an animal to occur at particular latitude and longitude:
```{r, eval = F}
plot_lon_lat(Result)
```
### Plot utilization distribution

Plot space utilisation distribution for the wintering range:
```{r, eval = F}
plot_util_distr(Result, 
    dates=data.frame(as.POSIXct('2013-12-01', tz='GMT'), as.POSIXct('2014-01-31', tz='GMT')),
    add.scale.bar=TRUE, percentiles=0.5)
```

## References

1. Ekstrom, P. (2004). An advance in geolocation by light. Memoirs of the National Institute of Polar Research, Special Issue, 58, 210–226.<a id="ekstrom_2004"></a>
1. Ekstrom, P. (2007). Error measures for template-fit geolocation based on light. Deep Sea Research Part II: Topical Studies in Oceanography, 54, 392–403.<a id="ekstrom_2007"></a>
1. Lisovski, S., Hahn, S. & Hodgson, D. (2012a). GeoLight - processing and analysing light-based geolocator data in R. Methods in Ecology and Evolution, 3, 1055–1059.<a id="lisovski_2012b"></a>
1. Lisovski, S., Hewson, C.M., Klaassen, R.H.G., Korner-Nievergelt, F., Kristensen, M.W. & Hahn, S. (2012b). Geolocation by light: accuracy and precision affected by environmental factors. Methods in Ecology and Evolution, 3, 603–612.<a id="lisovski_2012a"></a>
1. Rakhimberdiev, E., Winkler, D.W., Bridge, E., Seavy, N.E., Sheldon, D., Piersma, T. & Saveliev, A. (2015). A hidden Markov model for reconstructing animal paths from solar geolocation loggers using templates for light intensity. Movement Ecology, 3, 25.<a id="rakhimberdiev_movecol_2015"></a>
1. Rakhimberdiev, E., Senner, N. R., Verhoeven, M. A., Winkler, D. W., Bouten, W. and Piersma T. (2016) Comparing inferences of solar geolocation data against high-precision GPS data: annual movements of a double-tagged Black-Tailed Godwit. - Journal of Avian Biology 47: 589-596.<a id="rakhimberdiev_jab_2016"></a>
