---
title: "Learning to `crawl`"
author: "Kenady Wilson, Josh London, and Devin Johnson"
date: '`r Sys.Date()`'
output:
  html_document:
  number_sections: yes
pdf_document: default
subtitle: |
  Modeling Animal Movement by Example
vignette: |
  %\VignetteIndexEntry{Intro to crawl} 
  %\VignetteEncoding{UTF-8} 
  %\VignetteEngine{knitr::rmarkdown}
---

<style type="text/css">

.main-container {
max-width: 75rem;
font-size: 17px;
}
</style>

# Introduction

Electronic data collection via telemetry instruments is common for many species. Location estimates for the animals are not evenly-spaced in time due to the movement of satellites and the behavior of the animals themselves. In order to view a more continuous representation of the animals' movement, one must re-create the movement path. Several methods are available to interpolate the movement path, in the Correlated RAndom Walk Library (crawl) we use continuous-time correlated random walk models with time-indexed covariates. The model is fit using a Kalman filter on a state-space version of the continuous-time stochastic movement process. 

Telemetry data transmitted via Argos typically have error associated with each location estimate. `crawl` is now able to incorporate the new error ellipses as well as the older location class designations. 

# Preliminary Procedures
To use the package, you must first load the library:

```{r,message=FALSE}
library(crawl)
```

## Data Preparation
Telemetry data are often downloaded in the proprietory format of the tag manufacturer. Consequently, to run `crawl`, the data might need to be adjusted. Dataframes must be in a specific format in order for the functions to read the data properly. 

Run the following checks on your data prior to building a `crawl` model:

* **Time**

The date/time stamp must not contain any missing values, it should be ordered, and the class should be either numeric or POSIXct. To convert a timestamp of the form ("2011-03-27 01:30:00") to POSIX the following code may be helpful:

`as.POSIXct(strptime("2011-03-27", "%Y-%m-%d %H:%M:%S"))`

Users may also find the functions within the [lubridate](https://cran.r-project.org/package=lubridate) package also useful. Additionally, the date-time values should all be unique for a given animal. Duplicate date-time values are a known aspect of Argos data and duplicates should either be filtered from the data set or adjusted. Users can consult the `make.time.unique()` function within the [xts](https://cran.r-project.org/package=xts) package or the `adjust.duplicateTimes()` function wtihin the [trip](https://cran.r-project.org/package=trip) package in order to adjust any duplicate date-time values by 1 second.

* **Latitude and Longitude (X and Y values)**

The lat/long coordinates of your locaiton estimates should be numeric or defined as a spatial points data frame. If they are already spatial points, then check the projection to verify it's accuracy. If they are numeric, then `crawl` will **assume** the data are projected. 

* **Covariates**

If you are using covariates in your model, check to see that there are no missing values. You **can** have missing locations, but not missing covariates. Covariates that will be included in the model itself must be time-indexed.

You can also use the included Shiny App to check your data for any issues. The app will run basic checks on your Time and lat/long columns to make sure they comply with `crawl` specifications.

------------------------------------------------

## Northern Fur Seal Demo
This is an example using crawl to model the movements of one animal using Argos location classes for location error estimates, and incorporating drift into the movement process. These data represent Northern fur seal pup relocation data that were used in Johnson et al. (2008). The data are for one seal, including 795 observations with 4 variables. Northern fur seal pups travel long distances and may exhibit both directed travel and movement within large-scale ocean currents. Therefore, a varying drift model is included in this example for the mean velocity. 

### Load and clean the data
```{r}
data("northernFurSeal")
head(northernFurSeal)
```
Define the location classes as factors
```{r}
northernFurSeal$Argos_loc_class <- factor(northernFurSeal$Argos_loc_class,
                                          levels=c("3", "2", "1","0","A"))
```

### Make sure your data is projected
First, tell R which columns represent the coordinates for the data
```{r, message=FALSE}
library(sp)
library(rgdal)
coordinates(northernFurSeal) = ~longitude+latitude
```
Get the projection information from the data
```{r}
proj4string(northernFurSeal) <- CRS("+proj=longlat")
```
Run a spatial transform for map projection. In this example we define a custom projection for the data based on the location of the animals.
```{r,message=FALSE}
northernFurSeal <- spTransform(northernFurSeal, 
                               CRS(paste("+proj=aea +lat_1=30 +lat_2=70",
                                         "+lat_0=52 +lon_0=-170 +x_0=0 +y_0=0",
                                         "+ellps=GRS80 +datum=NAD83",
                                         "+units=m +no_defs"))
)
```

### Set initial parameters and priors for the model

`initial` is a list of starting values for the mean and variance-covariance for the initial state of the model. When choosing the initial parameters, it is typical to have the mean centered on the first observation with zero velocity. `a` is the starting location for the model -- the first known coordinate; and `P` is the initial velocity -- a 4x4 var-cov matrix. For these data, `a` should correspond to the location where the fur seal was instrumented.

```{r,message=FALSE}
initial = list(a=c(coordinates(northernFurSeal)[1,1],0,
                   coordinates(northernFurSeal)[1,2],0),
               P=diag(c(10000^2,54000^2,10000^2,5400^2)))
```

`fixPar` is used if you want to fix values in the model (usually to 0 or 1). Here are some examples of why you may want to fix parameters:

1. **You have Argos data with errors 1, 2, 3, A, B.** 

You want to fix the errors for location classes 1, 2, and 3 to what Argos 
suggests, but want to estimate the error associated with location classes 
A and B. In this example, we are fixing the known Argos errors. 

2. **You have an animal that exhibits limited movement.**

This scenario makes it difficult to estimate the autocorrelation parameter 
for the model. In this case, you can fix the autocorrelation parameter to 
improve optimization. 

3. **Your animal hauls out on land.** 

The activity model specifies the autocorrelation parameter $\beta_i = \beta/A^\phi$ and $\sigma_i^2 = \sigma_i^2*A^\phi$ for activity index of an animal, $A \in [0,1]$. E.g., See the harbor seal example in the next section. If the researcher only has binary data $A = 0$ or 1, then $\phi$ is unidentafiable and can be set to $\phi = 1$ 

```{r,message=FALSE}
fixPar = c(log(250), log(500), log(1500), rep(NA,3), NA)
```
To make sure everything looks as it should, view the parameters that will be used in the crawl model
```{r,message=FALSE}
displayPar(mov.model=~1,
           err.model=list(x=~Argos_loc_class-1),
           data=northernFurSeal,
           fixPar=fixPar)
```

The constraint parameter is a list of parameters with vectors for the upper and lower limits. This parameter can be used to estimate Argos error. You don't want to generate estimates lower than the values for error classes 3, 2, or 1 so you can set a lower bound. You can also constrain the autocorrelation paramater to between -4 and 4. (It is highly unlikely that this value will be outside this range.)

```{r, message=FALSE}
constr=list(lower=c(rep(log(1500),2), rep(-Inf,2)),
            upper=rep(Inf,4))
```

#### Set a prior. 

The prior used in this example is a Laplace prior (double exponential). To set your own Laplace prior or to get more information on the distribution itself, look into the `ddoublex()` function in the R package [smoothmest](https://cran.r-project.org/package=smoothmest). This function has the foloowing density: $exp(-abs(x-mu)/lambda)/(2*lambda)$. The prior below is on the log scale. If you need to adjust the prior, we recommend first adjusting the lambda value (0.5).  For more information on using a Laplace prior in general see: [Hooten and Hobbs (2015) Ecological Monographs 85:3-28] (http://www.esajournals.org/doi/full/10.1890/14-0661.1).

```{r, message=FALSE}
ln.prior = function(theta){-abs(theta[4]-3)/0.5}
```

### Fit the model

In this example we used an intercept model. You can add more covariates to the model if you have them, but they must be time-indexed. (ie. the number of dives per hour or whether it was day or night)
```{r,message=FALSE}
set.seed(123)
fit1 <- crwMLE(mov.model=~1, 
               err.model=list(x=~Argos_loc_class-1),
               data=northernFurSeal, 
               Time.name="Time",
               initial.state=initial,
               fixPar=fixPar, 
               constr=constr, 
               prior=ln.prior,
               control=list(maxit=30, trace=0,REPORT=1),
               initialSANN=list(maxit=200, trace=0, REPORT=1))
```

View the model output:
```{r}
fit1
```

### Predict regularly-spaced locations.

In order to standardize the data and make hourly locaiton predictions, use the `crwPredict` function. This function predicts the regular-timed locations along the movement path using the posterior mean and variance of the track.  

The speed estimates from this function are measures of instantaneous speed. The data must be projected and the output is in meters per whatever time unit you have specified. If the time unit is projected and POSIXct, then divide the estimate by 3600 to get a speed estimate in m/s. 

First, define the min and max (floor and cieling) times in your data
```{r,message=FALSE}
predTime <- seq(ceiling(min(northernFurSeal$Time)), 
                floor(max(northernFurSeal$Time)), 1)
```
Next, using the MLE model you fit above, predict locations using the time range specified by predTime
```{r,message=FALSE}
predObj <- crwPredict(object.crwFit=fit1, 
                      predTime, 
                      speedEst=TRUE, 
                      flat=TRUE)
```
Now, view the predicted movement path
```{r,message=FALSE}
crwPredictPlot(predObj, "map")
```

### Simulation 

Create a simulation object with 100 parameter draws. The simulator function is different from `crwPredict` in that you get a distribution of distances traveled. 
```{r,message=FALSE}
set.seed(123)
simObj <- crwSimulator(fit1, 
                       predTime, 
                       method="IS", 
                       parIS=100, 
                       df=5, 
                       scale=18/20)
```

### Examine the simulation

First, look at the importance sampling weight distribution. You want to have more weights near 1. If weights are not near one, you may want to adjust your prior.
```{r,message=FALSE}
w <- simObj$thetaSampList[[1]][,1]
hist(w*100, main='Importance Sampling Weights', sub='More weights near 1 is desirable')
```

Next, look at the approximate number of independent samples
```{r,message=FALSE}
round(100/(1+(sd(w)/mean(w))^2))
```

### Sample tracks and make maps

If the simulation looks good, sample 20 tracks from the posterior predictive distribution.

First, define your color ramp:
```{r}
my.colors <-colorRampPalette(c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a'))
```

Set the number of tracks you want to sample and define your colors:
```{r}
iter <- 20
cols <- my.colors(iter)
```

Next, sample from the posterior using the simulation object you created above.
```{r,message=FALSE}
crwPredictPlot(predObj, 'map')
for(i in 1:iter){
  samp <- crwPostIS(simObj)
  lines(samp$alpha.sim[,'mu.x'], samp$alpha.sim[,'mu.y'],col=cols[i]) 
}
```

Compare the sampled data to the predicted movement path from step 5

--------------------------------------------------

## Harbor Seal Demo

This example uses crawl to model the movements of one animal using location classes for locaiton error estimates, and incorporating a haul-out model into the movement process. The basic CTCRW model assumes the animal is in continuous motion for the length of the track, which is not the case for harbor seals and other animals that haul-out between foraging trips. This demo includes a continuously valued covariate which produces a model where a smooth range of haul-out behavior is allowed to act on the movement of the seal. This example also provides details on how to handle missing values that occur when locations and dry-time values are measured at different rates. 

### Load and clean the data
```{r}
library(crawl)
library(sp)
library(rgdal)
library(ggplot2)
data("harborSeal")
head(harborSeal)
```
Define the location classes as factors
```{r}
harborSeal$Argos_loc_class = factor(harborSeal$Argos_loc_class, levels=c("3","2","1","0","A","B"))
```

### Make sure your data is projected
Locations and dry-times are measured at different times; therefore, these data include times that do not have an associated location estimate. The lat/long for these times is set as a missing value. The rgdal package will not work with NAs in the data so projection will involve a few more steps than what was done with the Northern fur seal example above. 

First, get rid of the rows with NAs in the lat/long. (We are assuming that if a latitude value is missing, the longitude value is also missing.) 
```{r}
toProj = harborSeal[!is.na(harborSeal$latitude),c("Time","latitude","longitude")]
```
Tell R which columns represent the coordinates for the data
```{r, message=FALSE}
coordinates(toProj) = ~longitude+latitude
```
Convert to a Spatial Points Data Frame and project. This example uses a library of defined projections (epsg:3338). To see if your region has already been defined serach for it on [spatialreference.org] (spatialreference.org). 
```{r}
proj4string(toProj) <- CRS("+proj=longlat")
toProj <- spTransform(toProj, CRS("+init=epsg:3338"))
```
Now that the data are projected, convert the spdf back into a data.frame and merge the projected coordinates with the missing values we removed in step 1. 
```{r}
toProj = as.data.frame(toProj)
colnames(toProj)[2:3] = c("x","y")
harborSeal = merge(toProj, harborSeal, by="Time", all=TRUE)
harborSeal = harborSeal[order(harborSeal$Time),]
```

### Set initial parameters and priors for the model

Our initial locations were unknown because we allowed the animals to range freely, thus reducing capture effects on the data. Therefore, we used `a` and `P` values that were large enough to represent an unknown initial state when recording lat/long. 
```{r}
initial = list(a=c(harborSeal$x[1],0,harborSeal$y[1],0),
               P=diag(c(10000^2,5400^2,10000^2,5400^2)))
```
Here, we are fixing the known Argos location errors and the haul-out parameter.
```{r}
fixPar = c(log(250), log(500), log(1500), rep(NA,5), 0)
displayPar( mov.model=~1, err.model=list(x=~Argos_loc_class-1),data=harborSeal,activity=~I(1-DryTime),fixPar=fixPar)
```
Set the constraint parameter
```{r}
constr=list(lower=c(rep(log(1500),3), rep(-Inf,2)),
            upper=rep(Inf,5))
```

### Fit the model
This example uses an intercept model, with the haul-out model represented in the 'activity' parameter, and the error model defined using the Argos location classes. We did not define a prior for this model. A prior is generally not necessary for Maximum Likelihood Estimation; however, you may want to include one if you have optimization problems or you want to 'select' for fixed parameters via regularization.  
```{r,message=FALSE}
set.seed(123)
fit1 <- crwMLE(
  mov.model=~1, err.model=list(x=~Argos_loc_class-1), activity=~I(1-DryTime),
  data=harborSeal, coord=c("x","y"), Time.name="Time", 
  initial.state=initial, fixPar=fixPar, theta=c(rep(log(5000),3),log(3*3600), 0),
  constr=constr,
  control=list(maxit=2000, trace=1, REPORT=1)
)
```
View the model output:
```{r}
print(fit1)
```

###Predict regularly-spaced locations.
In this example `predTime=NULL` because there are already hourly times in the data due to the activity parameter 
```{r}
pred1 = crwPredict(fit1, predTime=NULL, flat=TRUE)
```
View the predicted movement path using ggplot2
```{r}
p1=ggplot(aes(x=mu.x, y=mu.y), data=pred1) + geom_path(col="red") + geom_point(aes(x=x, y=y), col="blue") + coord_fixed()

p2=ggplot(aes(x=Time, y=mu.x), data=pred1) + geom_ribbon(aes(ymin=mu.x-2*se.mu.x,ymax=mu.x+2*se.mu.x),fill="green", alpha=0.5)  + geom_path(col="red") + geom_point(aes(x=Time, y=x), col="blue", size=1)

p3=ggplot(aes(x=Time, y=mu.y), data=pred1) + 
  geom_ribbon(aes(ymin=mu.y-2*se.mu.y,ymax=mu.y+2*se.mu.y), fill="green", alpha=0.5)  + 
  geom_path(col="red") + geom_point(aes(x=Time, y=y), col="blue", size=1)

suppressWarnings(print(p1))
suppressWarnings(print(p2))
suppressWarnings(print(p3))
```

--------------------------------------------------

## Bearded Seal Example

This example builds off of the details listed in the Northern fur seal example and condenses the code to model the movement for multiple animals. We also run the data through an initial, course speed filter to remove extreme outliers. This example differs from the Northern fur seal example in that it:

* demonstrates some data preparation steps
* employs a course speed filter using the [argosfilter](https://cran.r-project.org/package=argosfilter) package
* demonstrates modeling of multiple animals, 
* Argos location errors are specified as ellipses instead of location classes,

The bearded seal movement data from three bearded seals are included with `crawl'. These data are generally representative of data delivered directly from Argos.

### Load the necessary packages and examine the data

```{r message=FALSE}
library(sp)
library(rgdal)
library(tidyr)
library(dplyr)
library(lubridate)
library(ggplot2)
```

Some initial tidying of the data has been done prior to inclusion with the package [dplyr](https://cran.r-project.org/package=dplyr) was used for this process and this resulted in a `tbl_df` object. 

```{r message=FALSE}
data("beardedSeals")
beardedSeals
```

The key columns to note are

* __deployid__ -- unique identifier for each seal
* __date_time__ -- POSIXct date-time with time zone set to "UTC"
* __quality__ -- categorical Argos location quality class
* __latitude__ -- y coordinate in decimal degrees; unprojected
* __longitude__ -- x coordinate in decimal degrees; unprojected
* __error_semimajor_axis__ -- length of the semi-major error axis (meters)
* __error_semiminor_axis__ -- length of the semi-minor error axis (meters)
* __error_ellipse_orientation__ -- error ellipse orientation (degrees)

### Adjust Duplicate Times

We have already confirmed that our `date_time` column is of class `POSIXct` and that the time zone is set to "UTC". As was stated earlier, Argos data are known to occasionally produce records with duplicate date-time values -- this dataset is no exception.

```{r message=FALSE}
beardedSeals %>% 
  group_by(deployid,date_time) %>% 
  filter(n()>1)
```

Note that while the `date_time` columns are duplicated, values for `quality`, `latitude`, and `longitude` differ. Since we have no objective method for selecting one record over the other, we will opt to adjust one of the duplicate times by increasing the `date_time` by 1 second. In this example, we use the `make.time.unique()` function from the [xts](https://cran.r-project.org/package=xts) package (along with some dplyr and [tidyr](https://cran.r-project.org/package=tidyr) magic) to create a new column, `unique_posix`.

```{r message=FALSE}
library(xts)
date_unique <-beardedSeals %>% 
  group_by(deployid) %>%
  do(unique_date = xts::make.time.unique(.$date_time,eps=1)) %>%
  tidyr::unnest(unique_date) %>%
  mutate(unique_posix = as.POSIXct(.$unique_date,origin='1970-01-01 00:00:00',tz='UTC')) %>%
  dplyr::arrange(deployid,unique_posix) %>% 
  dplyr::select(unique_posix)

beardedSeals <- beardedSeals %>% arrange(deployid,date_time) %>%
  bind_cols(date_unique)
```

Now, if we return to our previous code to check for duplicates but examine the `unique_posix` column, we see that no duplicate date-time values are found.

```{r message=FALSE}
beardedSeals %>% 
  group_by(deployid,unique_posix) %>% 
  filter(n()>1)
```

### Remove Obviously Erroneous Locations

Prior to modeling the movement, we will use the `argosfilter` package to pass the Argos locations through a speed filter algorithm. The nature of the Argos satellite system and deployments on wild, marine animals inevitably leads to obviously erroneous locations in the dataset. While the movement models can handle inclusion of these locations, they can sometimes lead to non-convegence or biologically unreasonable movement parameters. Users should use this approach with some caution, however, as we do not want to remove informative observations. To mitigate these issues, here, we pass the data through a speed filter with a speed greater than twice the expected max speed for bearded seals (vmax=5.0 m/s). The `sdafilter` function also has the capability to remove spikes along the observed track. We set `ang = -1` to disable this functionality.

Also notice that we emply the `doParallel` and `foreach` function in order to perform this procedure in parallel and save some time.

```{r message=FALSE}
beardedSeals <- beardedSeals %>%  
  dplyr::arrange(deployid,unique_posix)

library(doParallel)
library(argosfilter)

split_data <- split(beardedSeals,beardedSeals$deployid)

registerDoParallel(cores=2)
beardedSeals$filtered <- foreach(i = 1:length(split_data), .combine = c) %dopar% {
  argosfilter::sdafilter(
    lat=split_data[[i]]$latitude, 
    lon=split_data[[i]]$longitude, 
    dtime=split_data[[i]]$unique_posix,
    lc=split_data[[i]]$quality, 
    ang=-1,
    vmax=5)
}
stopImplicitCluster()

beardedSeals <- beardedSeals %>% 
  dplyr::filter(., filtered=="not" & !is.na(error_semimajor_axis)) %>%
  arrange(.,deployid,unique_posix)
```

### Convert to SpatialPointsDataFrame and Project

Passing a properly projected _SpatialPointsDataFrame_ to `crawl::crwMLE()` is the preferred practice. The _SpatialPointsDataFrame_ insures the spatial structure of the data is properly specified. By projecting the coordinates from geographic, the user can avoid common pitfalls associated with geographic coordinate systems (e.g., crossing 180, converting decimal degrees to meters away from the equator). Here, since our data range from the Bering and Chukchi seas, we specify the [North Pole Lambers Albert Equal Area Bering Sea](http://spatialreference.org/ref/epsg/wgs-84-north-pole-laea-bering-sea/) projection by refering the shorthand EPSG code of 3571.


```{r message=FALSE}
beardedSeals <- as.data.frame(beardedSeals)
coordinates(beardedSeals) = ~longitude+latitude
proj4string(beardedSeals) = CRS("+proj=longlat +datum=WGS84")

beardedSeals <- spTransform(beardedSeals, CRS("+init=epsg:3571"))
```

### Fit the model

It is not uncommon for a user to have multiple animals from the same study that need to be modeled as part of the same analysis. This approach is well suited for a parallelized programming approach that takes advantage of modern computer hardware (often with 8 or more CPU cores) and the functionality provided by the `foreach` and `doParallel` packages. This is not intended to be a detailed tutorial on running parallel analysis within R, but should provide most users a basic example. The code below will loop over the `deployid` identifiers and fit the CTCRW models for each seal and returned object, `model_fits` is a list of model fits. Refer to Step 3 in the Northern fur seal model for details on the initial parameters. 

Unlike with the Northern fur seal model, the data in this example contain parameters for the ellipse errors. 

```{r message=FALSE}
ids = unique(beardedSeals@data$deployid)      #define seal IDs

registerDoParallel(cores=2)
model_fits <-
  foreach(i = 1:length(ids)) %dopar% {
    id_data = subset(beardedSeals,deployid == ids[i])
    diag_data = model.matrix(
      ~ error_semimajor_axis + error_semiminor_axis + error_ellipse_orientation,
      model.frame( ~ ., id_data@data, na.action = na.pass)
    )[,-1]
    
    id_data@data = cbind(id_data@data, 
                         crawl::argosDiag2Cov(
                           diag_data[,1], 
                           diag_data[,2], 
                           diag_data[,3]))
    
    init = list(a = c(sp::coordinates(id_data)[1,1],0,
                      sp::coordinates(id_data)[1,2],0),
                P = diag(c(5000 ^ 2,10 * 3600 ^ 2, 
                           5000 ^ 2, 10 * 3600 ^ 2)))
    
    fit <- crawl::crwMLE(
      mov.model =  ~ 1,
      err.model = list(
        x =  ~ ln.sd.x - 1, 
        y =  ~ ln.sd.y - 1, 
        rho =  ~ error.corr
      ),
      data = id_data,
      Time.name = "unique_posix",
      initial.state = init,
      fixPar = c(1,1,NA,NA),
      theta = c(log(10), 3),
      initialSANN = list(maxit = 2500),
      control = list(REPORT = 10, trace = 1)
    )
    fit
  }
stopImplicitCluster()

names(model_fits) <- ids

print(model_fits)
```

### Predict regularly-spaced locations
This function predicts the regular-timed -- in this case, hourly -- locations along the movement path using the posterior mean and variance of the track. 

```{r message=FALSE}
registerDoParallel(cores=2)
predData <- foreach(i = 1:length(model_fits), .combine = rbind) %dopar% {
  
  model_fits[[i]]$data$unique_posix <- lubridate::with_tz(
    model_fits[[i]]$data$unique_posix,"GMT")
  predTimes <- seq(
    lubridate::ceiling_date(min(model_fits[[i]]$data$unique_posix),"hour"),
    lubridate::floor_date(max(model_fits[[i]]$data$unique_posix),"hour"),
    "1 hour")
  tmp = crawl::crwPredict(model_fits[[i]], predTime=predTimes)
}
stopImplicitCluster()

predData$predTimes <- intToPOSIX(predData$TimeNum)
```

While a projected _SpatialPointsDataFrame_ was passed to the `crwMLE()` funtion, the prediction object returned from `crwPredict()` is not a _SpatialPointsDataFrame_. The columns `mu.x` and `mu.y` represent the predicted coordinates which we can coerce into a _SpatialPointsDataFrame_ with the `coordinates()` function and then specify our projection with the `proj4string` function.

```{ message=FALSE}
predData_sp <- predData
coordinates(predData_sp) <- ~mu.x+mu.y
proj4string(predData_sp) <- CRS("+init=epsg:3571")
```

### Plot the output
Using ggplot, we can examin the predicted movement path for each bearded seal. We will create a custom ggplot theme to provide a mapping look and feel.

```{r plot-1}
theme_map = function(base_size=9, base_family="")
{
  require(grid)
  theme_bw(base_size=base_size, base_family=base_family) %+replace%
    theme(axis.title.x=element_text(vjust=0),
          axis.title.y=element_text(angle=90, vjust=1.25),
          axis.text.y=element_text(angle=90),
          axis.ticks=element_line(colour="black", size=0.25),
          legend.background=element_rect(fill=NA, colour=NA),
          legend.direction="vertical",
          legend.key=element_rect(fill=NA, colour="white"),
          legend.text=element_text(),
          legend.title=element_text(face="bold", hjust=0),
          panel.border=element_rect(fill=NA, colour="black"),
          panel.grid.major=element_line(colour="grey92", size=0.3, linetype=1),
          panel.grid.minor=element_blank(),
          plot.title=element_text(vjust=1),
          strip.background=element_rect(fill="grey90", colour="black", size=0.3),
          strip.text=element_text()
    )
}

p1 <- ggplot(data=predData,aes(x=mu.x,y=mu.y)) + 
  geom_path(aes(colour=deployid)) + xlab("easting (meters)") +
  ylab("northing (meters)") + theme_map()
p1
```
