### Much of the code below was modified from code originally provided by S. Lisovski

library(GeoLight) 
library(SGAT) 
library(BAStag) 
library(MASS) 
library(maptools)
library(geosphere)
library(raster)
library(geosphere)

###################################
########## LOADING IN DATA ########
###################################

### Function written by Simon Wotherspoon for reading .lig files from British Antarctic Survey tags. The file 
### "read_lig.R" needs to be in your folder linked to R, and can be found at ###"https://github.com/bsevansunc/MigratoryBirds/blob/master/R_scripts/read_lig.R"
### This step will vary depending on geolocator manufacturer
source("read_lig.R")

### Import .lig file using read.lig fucntion to convert dates that GeoLight recognizes. 
d.lux<-read.lig("GL14_000.lig")

### Subset the Date and Light columns needed for analysis
d.lux<-subset(d.lux, select=c("Date","Light"))

head(d.lux)           
#                 Date Light
#1 2014-06-11 02:17:33    44
#2 2014-06-11 02:19:33     0
#3 2014-06-11 02:21:33     1
#4 2014-06-11 02:23:33     1
#5 2014-06-11 02:25:33     0
#6 2014-06-11 02:27:33     0

##########################################################################################
########## SETTING UP CALIBRATION PERIOD AND OUTER BOUNDARIES OF ANALYZABLE DATA #########
##########################################################################################

### Enter the lat/lon of each bird's location of geolocator deployment
lat.calib<- 44.655523 
lon.calib<- -84.647636

### Setting up data-frame for use below
tm <- seq(d.lux[1,1], d.lux[nrow(d.lux),1], by = "day")
rise <- rep(c(TRUE, FALSE), length(tm))

### Making predicted twilight times given location and zenith for graphs below
c.dat <- data.frame(Twilight = twilight(rep(tm, each = 2), lon = lon.calib, lat = lat.calib, rise = rise, zenith = 93.17), Rise = rise) 
offset = 19

### Plot your light data with predicted twilight times assuming bird was stationary at deployment site for entire year
lightImage(d.lux, offset = offset, zlim = c(0,12))
tsimagePoints(c.dat$Twilight, offset = offset, pch = 16, cex = 0.5, col = ifelse(c.dat$Rise, "dodgerblue", "firebrick"))

### Adds line at equinoxes
eqnx<-as.POSIXct(c("2014-09-23", "2015-03-20"), tz = "GMT") 
abline(v = eqnx, lwd=3, lty=3, col="purple")

### Leave plot from above open. The locator function can be used click on the outer boundaries of the calibration 
### period. Once you run code below, you get a target cursor. First click defines the start 
### of the calibration period, second click defines the end. User can define calibration period however they choose. 
### Choice will affect resulting twilight error distribution.
tm.calib <- as.POSIXct(locator(n=2)$x, origin = "1970-01-01")

### See results of your selection
tm.calib

### Take the dates from above and enter into below
tm.calib <- as.POSIXct(c("2014-06-13", "2014-07-29"), tz = "GMT") 

### Visualize your calibration period
rect(tm.calib[1], 26, tm.calib[2], 39, col = rgb(1, .5, 0, alpha = 0.4))

### In the same way, we can define the start/end dates of usable data (must still have figure from above open) 
### Alternatively you can just enter the day after deployment at the start and the day before capture as the end
tm.track <- as.POSIXct(locator(n=2)$x, origin = "1970-01-01") 
tm.track

### Enter known start/end dates or dates from tm.track above
tm.track <- as.POSIXct(c("2014-06-13", "2015-05-11"), tz = "GMT")

### Visualize start/end dates of usable data
abline(v = tm.track, lwd = 3, lty = 3, col = "firebrick")

###########################################
########## DEFINING TWILIGHT TIMES ########
###########################################

### As of present, this next step unfortunately does not work on Max OSX machines due to the need to work with 
### interactive plots

### The function preprocessLight opens two interactive plots and you have to go through 4 different steps: 

###	-1. Define the desired period with left click on the left boundary and right clock on the right boundary. Press "a" 	
###		to go to the next step. 
###	-2. Define seeds (one or multible) - a position during the night that shows clear sunrise/sunset boundaries. You can 		
###		click on plot 1 to zoom into a region that will be shown in plot 2 (seeds need to be defined in plot 2). Press ###     "a" to go to the next step. 
###	-3. Add non-defined twilights if necessary. Press "a" to go to the next step. 
###	-4. You can go through each twilight using the forward and backward arrows or click on specific twilight times. The 
###		selected twilight times will be shown in the second plot. You can change the twilight time by clickinga at the 
###		desired position and press "a". Press "q" to close plots and return to R.
### See ?preprocessLight for other functions

### Sets user-defined light-level threshold
threshold <- 1

### Brings up interactive process. Can define minimum dark period. Here it is set at 4 hours (240 minutes)
twl <- preprocessLight(d.lux, threshold, offset = offset, zlim = c(0, 12), dark.min = 240)

### Saves results of above process
write.csv(twl, paste("GL14_twl.csv", sep = "?/"), row.names = F)

### Reading in data transfered from PC. All code below runs on Mac OSX or PC
twl <- read.csv(paste("GL14_twl.csv", sep = "/"))

### Converts twilights into appropriate format
twl$Twilight <- as.POSIXct(strptime(twl$Twilight, format="%Y-%m-%d %H:%M"), tz = "GMT")

### Account for storing maximum values for each 2 minute period. This step will depend on tag type. 
twl <- twilightAdjust(twl, 60)

### Show results of defining twilights against raw light data
lightImage(d.lux, offset = offset, zlim = c(0,12))
tsimagePoints(twl$Twilight, offset = offset, pch = 16, cex = 0.5,
               col = ifelse(twl$Rise, "dodgerblue", "firebrick"))

################################################
########## DETERMINING THE ZENITH ANGLE ########
################################################

### Subsetting data based on calibration period determined above and takes into account any twilights deleted.
d.calib <- subset(twl, Twilight>=tm.calib[1] & Twilight<=tm.calib[2] & !Deleted)

### Calculates solar information needed for determining zenith angle.
sun  <- solar(d.calib[,1]) 

### Estimates zenith angle for each date of calibration period.
z    <- refracted(zenith(sun, lon.calib, lat.calib))

### Zenith angle for simple threshold estimates.
zenith0 <- median(z) 
zenith0

### Zenith angle for estelle model.
zenith  <- quantile(z, prob = .95) 
zenith

###########################################################
########## DETERMINING TWILIGHT ERROR DISTRIBUTION ########
###########################################################

### Calculates known twilight times at deployment site for all dates during calibration period.
twl_t   <- twilight(d.calib[,1], lon.calib, lat.calib, rise = d.calib[,2], zenith = max(z)+.1) 

### Calculates difference between twilight times for deployment site and twilight times estimated from light-level data.
twl_dev <- ifelse(d.calib$Rise, as.numeric(difftime(d.calib[,1], twl_t, units = "mins")), as.numeric(difftime(twl_t, d.calib[,1], units = "mins")))

### Plots histogram of differences between estimated and known twilight times and then attempts to fit the distribution
### with a log-Normal distribution. Here you are looking for a good fit between the dashed red line (alpha) and the 
### historgram. If you do not find a find a good fit, you can change "zenith = max(z)+.1" above, change calibration 
### period, distribution type, or light-level threshold.
hist(twl_dev, freq = F, ylim = c(0, 0.12), xlim = c(0, 35), main = "", xlab = "Twilight error (mins)") 
seq <- seq(0,35, length = 100)
fitml <- fitdistr(twl_dev, "log-Normal")
alpha <- c(fitml$estimate[1], fitml$estimate[2])
lines(seq, dlnorm(seq, alpha[1], alpha[2]), col = "firebrick", lwd = 3, lty = 2)

#################################################################
########## CREATING LOGICAL PRIOR TO BEGIN MCMC SAMPLING ########
#################################################################

### This process estimates positions based on light-level data using threshold method. Estelle model below needs a 
### a logical starting point to begin the sampling process

### Subsets data by previously chosen boundaries, and takes into account any twilight events you may have deleted during
### interactive twilight estimation process
d.twl <- subset(twl, Twilight>=tm.track[1] & Twilight<=tm.track[2] & !Deleted)

### Creates path from simple threshold estimates. Varying the "tol" value will change the duration of period surrounding
### the equinoxes that latitudes will not be estimated from light-level data. Instead latitudes will be estimated
### in a linear progression from the first twilight before the equinox period until the first twilight after the equinox ### period
path <- thresholdPath(d.twl$Twilight, d.twl$Rise, zenith = zenith0, tol = 0.1)

### Loads in basic world map
data(wrld_simpl)

### Plots estimated longitudes from light-level data against time, adding a line indicating the longitude from 
### deployment site for reference
opar <- par(mfrow = c(2, 1), mar = c(2,4,1,1)+0.1) 
plot(path$time, path$x[, 1], type = "b", pch = 16, cex = 0.5, ylab = "Lon", xlab = '') 
abline(h = lon.calib)

### Adds line at two equinoxes for 2014/2015
eqnx<-as.POSIXct(c("2014-09-23", "2015-03-20"), tz = "GMT") 
abline(v = eqnx, lwd=3, lty=3, col="purple")

### Adds line indicating end of calibration period
abline(v = tm.calib[2])

### Plots estimated latitudes from light-level data against time, adding a line indicating the latitude from 
### deployment site for reference
plot(path$time, path$x[, 2], type = "b", pch = 16, cex = 0.5, ylab = "Lat", xlab = '') 
abline(h = lat.calib)

### Adds line at two equinoxes for 2014/2015 (dates can vary by year) 
eqnx<-as.POSIXct(c("2014-09-23", "2015-03-20"), tz = "GMT") 
abline(v = eqnx, lwd=3, lty=3, col="purple")
abline(v = tm.calib[2])

par(opar)

### Plots estimated movement path of bird and lat/lon of deployment site
plot(path$x, type = "n", ylim = c(15, 50), xlim = c(-100,-60)) 
plot(wrld_simpl, add = T, col = "grey95")
box() 
lines(path$x, col = "blue")
points(path$x, pch = 16, cex = 0.5, col = "blue")
abline(h = lat.calib)
abline(v = lon.calib)

### Above process may vary quite a bit depending on your data. If possible, you need to confirm that the zenith angle is 
### appropriate throughout the year. For our males, the latitudes just after the spring equinox were farther to the 
### north than the rest of the wintering period, and often located in open ocean, indicating that the zenith angle 
### estimated from the calibration period was not appropriate for all periods of the annual cycle. To estimate a zenith ### angle for the wintering period, we first needed to estimate the dates for the breeding period, wintering period, and ### fall and spring migration. In order to do this we calculated the distance from the estimated location on each day to ### both the deployment site and the centroid of the locations between November 15th and February 27th (i.e., the 
### wintering period). We calculated these distances and estimated the dates of each migration began by determining the ### first day each male moved greater than 150 km away from the breeding or wintering period. The ends of migration were ### estimated as the first day each male was within 150 km of the breeding or wintering period. These dates are not the ### final migration dates and are just used to determine the wintering zenith angle. This process may or may not be 
### necessary depending upon your data. For simplicity code above and below assumes single zenith angle for entire year. ### However, multiple zenith angles can be used with the following code:

### sets up empty vector where length equals number of positions you have for each bird
zenithlist<-vector("numeric" , 664)

### Here zenith0 is used when bird on breeding grounds, 94.3 is used during wintering period, and average of two used
### during migration periods

zenithlist[1:155]<-zenith0
zenithlist[156:259]<- (zenith0+94.3)/2
zenithlist[257:632]<- 94.3
zenithlist[633:658]<- (zenith0+94.3)/2
zenithlist[659:664]<- zenith0

### Creating initial path using multiple zenith angles. Ignores warnings
suppressWarnings(path <- thresholdPath(d.twl$Twilight, d.twl$Rise, zenith = zenithlist, tol = 0.20))

######################################################################
########## CREATING FLIGHT SPEED DISTRIBUTION ########################
######################################################################

### Define and plot flight speed distribution. May need to customize depending upon species. We assumed most likely 
### speed was close to zero given that birds are mostly stationary throughout the year. Speeds up to ~50 km were 
### allowed
beta = c(0.7, 0.08)
matplot(0:100, dgamma(0:100, beta[1], beta[2]), type = "l", col="red", xlab = "Speed (km/h)", ylab = "Probability", xlim = c(0,100))
abline(h = 0)

######################################################################
########## CREATING THE SPATIAL PROBABILITY MASK #####################
######################################################################

### We created a three-part spatial probability mask. The first part covers from Florida to just below Hudson's Bay
### The second part covers the Gulf of Mexico region and the third part covers the Caribbean. 

### Load shapefile for the Americas
Americas<-shapefile("Americas.shp")

### Defining the boundaries of the first part of the mask
xlim1 = c(-100,-60)
ylim1 = c(31, 51)

### Creates empty raster
r1 <- raster(extent(xlim1[1], xlim1[2], ylim1[1], ylim1[2]), resolution = 0.15)

### Rasterizes Americas shapefile
wrld1 <- rasterize(Americas, r1)
  wrld1[] <- ifelse(wrld1[]>1, 1, NA)

d.mask1 <- distance(wrld1)

### Change according to your preferences. In this case, the locations within 100 km of land are 5 times more likely
### than those over water
sX1 <- seq(0, max(d.mask1[], na.rm = T), by = 200) 
sY1 <- 1 + 5*exp(-((sX1)/100000)^1.5) 
 
prob1 <- approxfun(x=sX1, y=sY1)
 
prob.mask1 <- d.mask1
prob.mask1[] <- prob1(d.mask1[])

### Defining the boundaries of the second part of the mask
xlim2 = c(-100,-82)
ylim2 = c(15, 30.99999)

### Creates empty raster
r2 <- raster(extent(xlim2[1], xlim2[2], ylim2[1], ylim2[2]), resolution = 0.15)

wrld2 <- rasterize(Americas, r2)
  wrld2[] <- ifelse(wrld2[]>1, 1, NA)

d.mask2 <- distance(wrld2)

### Change according to your preferences. In this case, the locations within 100 km of land are 5 times more likely
### than those over water.
sX2 <- seq(0, max(d.mask2[], na.rm = T), by = 200) 
sY2 <- 1 + 5*exp(-((sX2)/100000)^1.5)
 
prob2 <- approxfun(x=sX2, y=sY2)
 
prob.mask2 <- d.mask2
prob.mask2[] <- prob2(d.mask2[])

### Defining the boundaries of the third part of the mask
xlim3 = c(-81.99999999,-60)
ylim3 = c(15, 30.99999)

### Creates empty raster
r3 <- raster(extent(xlim3[1], xlim3[2], ylim3[1], ylim3[2]), resolution = 0.15)

wrld3 <- rasterize(Americas, r3)
  wrld3[] <- ifelse(wrld3[]>1, 1, NA)

d.mask3 <- distance(wrld3)

### Change according to your preferences. In this case, the locations within 300 km of land are 5 times more likely
### than those over water. We found that when using only 100 km for the Caribbean, the large landmass of Cuba acted as
### a magnet and nearly all birds ended up wintering there even though the raw positions of only one male were on Cuba
### This mask, basically treats the entire Caribbean as a land mass to avoid this problem.
sX3 <- seq(0, max(d.mask3[], na.rm = T), by = 200) 
sY3 <- 1 + 5*exp(-((sX3)/300000)^15)
 
prob3 <- approxfun(x=sX3, y=sY3)
 
prob.mask3 <- d.mask3
prob.mask3[] <- prob3(d.mask3[])

prob.mask<-merge(prob.mask1, prob.mask2, prob.mask3, tolerance = 0.4)

### Plots final mask
plot(prob.mask)
map('state', fill = T, col = "transparent",border = "grey80", add = T)
plot(Americas,xlim=xlim,ylim=ylim, col = "transparent",border = "grey10", add = T)

### Reset xlim and ylim to include both rasters
xlim = c(-100,-60)
ylim = c(15, 51)

### Below code sets up the spatial probability as a prior for use during MCMC sampling later
lookup <- function(r, xlim, ylim) {
  r <- as.matrix(r)[nrow(r):1,]
  
  xbin <- seq(xlim[1], xlim[2], length = ncol(r) + 1)
  ybin <- seq(ylim[1], ylim[2], length = nrow(r) + 1)
  
  function(p) {
    r[cbind(.bincode(p[, 2], ybin), .bincode(p[, 1], xbin))]
  }
}

prob <- lookup(prob.mask, xlim = xlim, ylim = ylim)

log.prior <- function(p) {
  f <- prob(p)
  ifelse(is.na(f), -1000, log(f))
}


######################################################################
########## RUNNING THE ESTELLE MODEL #################################
######################################################################

### Note that below code can take hours to run for each bird depending upon the number of iterations

### Sets your initial movement path from above as a starting point for MCMC sampling
x0 <- path$x

### Sets the first two and last two locations as fixed at the deployment and re-capture site. This assumes deployment
### and re-capture occurred in same place but you can alter code if that is not the case.
fixedx <- rep(F, nrow(x0))
fixedx[1:2] <- T
fixedx[nrow(x0) - 1] <- T
fixedx[nrow(x0)] <- T

x0[fixedx, 1] <- lon.calib
x0[fixedx, 2] <- lat.calib

### Sets intermediate points along path
z0 <- trackMidpts(x0)

### Defines estelle model. Note that estelle model requires different zenith angle. Calculated above as zenith instead
### of zenith0. Can use list of multiple zenith angles for different periods of annual cycle if desired.
model <- thresholdModel(d.twl$Twilight,d.twl$Rise, 
                        twilight.model="ModifiedLogNormal",
                        alpha=alpha,beta=beta,
                        logp.x = log.prior, logp.z = log.prior,
                        x0=x0,z0=z0,zenith= zenith,fixedx=fixedx)

proposal.x <- mvnorm(S=diag(c(0.005,0.005)),n=nlocation(x0))
proposal.z <- mvnorm(S=diag(c(0.005,0.005)),n=nlocation(z0))

### Burn-in sampling
CfitGL14 <- estelleMetropolis(model,proposal.x,proposal.z,iters=10000,thin=15,chains=3)

### Tuning the model
x0 <- chainLast(CfitGL14$x)
z0 <- chainLast(CfitGL14$z)

model <- thresholdModel(d.twl$Twilight, d.twl$Rise, twilight.model = "LogNormal", 
    alpha = alpha, beta = beta, x0 = x0, z0 = z0, zenith = zenith, 
    logp.x = log.prior, logp.z = log.prior,
    fixedx = fixedx)

proposal.x <- mvnorm(S=diag(c(0.005,0.005)),n=nlocation(x0))
proposal.z <- mvnorm(S=diag(c(0.005,0.005)),n=nlocation(z0))


CfitGL14 <- estelleMetropolis(model,proposal.x,proposal.z,
                          iters=10000,thin=15,chains=3)

for(k in 1:3) {
  proposal.x <- mvnorm(chainCov(CfitGL14$x),s=0.2)
  proposal.z <- mvnorm(chainCov(CfitGL14$z),s=0.2)
  CfitGL14 <- estelleMetropolis(model,proposal.x,proposal.z,
                            x0=chainLast(CfitGL14$x),
                            z0=chainLast(CfitGL14$z),
                            iters=10000,thin=15,chains=3)
}
                            
### Final run of the model
proposal.x <- mvnorm(chainCov(CfitGL14$x),s=0.25)
proposal.z <- mvnorm(chainCov(CfitGL14$z),s=0.25)

CfitGL14 <- estelleMetropolis(model,proposal.x,proposal.z,
                          x0=chainLast(CfitGL14$x),
                          z0=chainLast(CfitGL14$z),
                          iters=5000,thin=20,chains=3)
                          
###Save model for later analysis/exploration
save(CfitGL14, file = "GL14Cfit.rda")

######################################################################
########## MODEL DIAGNOSTICS #########################################
######################################################################

### Load model
load("GL14Cfit.rda")
Cfit<-CfitGL14

### Adds some colors for graphs below
sr.col <- "#FF7F00"
ss.col <- "#377EB8"
mn.col <- "firebrick1"
ci.col <- "dodgerblue1"
asm.col <- "grey70"
trk.col <- "steelblue"
map1.col <- "honeydew3"
map2.col <- "honeydew4"
grp.pal <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3")

### Summarize posterior distribution
s <- locationSummary(Cfit$x,time=model$time,collapse=F)

### Plots chain convergence graphs for longitude. Visually inspect for convergence. Convergence may not occur near
### equinoxes regardless of number of iterations.
opar <- par(mfrow=c(1,3),mar=c(3,5,2,1)+0.1)
for(k in 1:length(s)) {
  ci <- as.matrix(s[[k]][,c("Lon.2.5%","Lon.97.5%")]-s[[k]][,"Lon.mean"])
  ci[fixedx,] <- NA
  matplot(s[[k]]$Time,ci,type="l",lty=1,ylab="Lon",xaxt="n",col=ci.col,cex.axis=0.7,yaxp = c(-6,6,24))
  axis.POSIXct(1,at=seq(round(min(s[[k]]$Time),"day"),
                   round(max(s[[k]]$Time),"day"),"weeks"),
               format="%d %b",cex.axis=0.7,las=2)
  for(j in 1:length(s))
    lines(s[[j]]$Time,s[[j]][,"Lon.mean"]-s[[k]][,"Lon.mean"],col=grp.pal[j])
}

par(opar)

### Plots chain convergence graphs for latitude. Visually inspect for convergence. Convergence may not occur near
### equinoxes regardless of number of iterations.
opar <- par(mfrow=c(1,3),mar=c(3,5,2,1)+0.1)
for(k in 1:length(s)) {
  ci <- as.matrix(s[[k]][,c("Lat.2.5%","Lat.97.5%")]-s[[k]][,"Lat.mean"])
  ci[fixedx,] <- NA
  matplot(s[[k]]$Time,ci,type="l",lty=1,ylab="Lat",xaxt="n",col=ci.col,cex.axis=0.7, yaxp = c(-6,6,24))
  axis.POSIXct(1,at=seq(round(min(s[[k]]$Time),"day"),
                   round(max(s[[k]]$Time),"day"),"weeks"),
               format="%d %b",cex.axis=0.7,las=2)
  for(j in 1:length(s))
    lines(s[[j]]$Time,s[[j]][,"Lat.mean"]-s[[k]][,"Lat.mean"],col=grp.pal[j])
}

### Plot movement paths estimated by each chain. Compare for convergence.
xlim = c(-90, -65)
ylim = c(15, 50)

plot(wrld_simpl,xlim=xlim,ylim=ylim,col=map1.col,border=map2.col)
plot(elide(wrld_simpl,shift=c(360,0)),xlim=xlim,ylim=ylim,add=T,col=map1.col,border=map2.col)
for(k in 1:length(s)) {
  lines(s[[k]][,"Lon.mean"],s[[k]][,"Lat.mean"],
        col=rgb(t(col2rgb(grp.pal[k]))/255,alpha=0.4))
}
box()

######################################################################
########## BASIC PLOTTING ############################################
######################################################################

load("GL14Cfit.rda")
Americas<-shapefile("Americas.shp")

### Set boundaries of raster for mapping
xlim = c(-95, -65)
ylim = c(10, 52)

### Sets up raster
r <- raster(ncol = sum(diff(xlim))*4, nrow = sum(diff(ylim))*4, xmn = xlim[1], xmx = xlim[2],

            ymn = ylim[1], ymx = ylim[2])

### Slices raster into sections so that you can plot any time period. Default slices by each location
s  <- slices(type = "intermediate", mcmc = CfitGL14, grid = r) ## intermediate == z

### Select which slices you want. Depends on desired output. If you want to display full year at once select all 
### locations, otherwise can display just migration period or just wintering period, etc. by changing k = 
### Important to remember that probability of residency is relative to the time period selected. If full year is
### selected highest probability will occur where most time spent (usually wintering grounds). Depending on question
### you may want to inspect figures of just spring or just fall migration or other periods of time. You can also exclude
### the equinox periods if desired by using multiple slices (e.g., k = c(1:200, 230:500)
sk <- slice(s, k = 1:663)

### Plot base map
plot(Americas,xlim=xlim,ylim=ylim,col = "grey90",border = "grey10")

### Plot probability of residency from posterior distribution. Breaks are user defined and depend on what you want to 
### do. Here the 95th quantile is used.
plot(sk, useRaster = F, legend = F, add = T, breaks = c(seq(0, quantile(sk, prob = 0.95), length = 100), maxValue(sk)),

     col = c(rep("transparent", 1), rev(topo.colors(99))))
 
### Add state boundaries     
map('state', fill = T, col = "transparent",border = "grey80", add = T)
 
### Re-plot Americas boundaries on top    
plot(Americas,xlim=xlim,ylim=ylim, col = "transparent",border = "grey10", add = T)

### Take mean from posterior distribution to estimate movement path
Cxm <- locationMean(CfitGL14$x)

### Define spring migration
springmig<-Cxm[632:659,]

### Define Fall migration
fall2<-Cxm[209:259,]

### Plot spring and fall migration lines
lines(springmig, lwd = 1.2, col=rgb(t(col2rgb("red"))/255,alpha=0.8))
lines(fall2, lwd = 1.2, col=rgb(t(col2rgb("black"))/255,alpha=0.4))


######################################################################
########## DETERMINING SPRING STOPOVER PERIODS #######################
######################################################################

### Load in posterior distribution
load("GL14Cfit.rda")
fit<-CfitGL14

### Summarize posterior distribution
zm <- locationSummary(fit$z,time=fit$model$time,collapse=T)

### Convert mean of posterior distribution (i.e., mean migration path) into twilight times
twl.back <- data.frame(Twilight=twilight(fit$model$time[-length(fit$model$time)], zm$Lon.mean, zm$Lat.mean, fit$model$rise[-length(fit$model$time)], zenith=zenith0, iters=5), Rise=fit$model$rise[-length(fit$model$time)])

### Visualize light data
lightImage(d.lux,offset=offset,zlim=c(0,12))
tsimagePoints(twl.back$Twilight,offset=offset,pch=16,cex=0.8,col=ifelse(twl.back$Rise,"dodgerblue","firebrick"))

### Re-format data for GeoLight
twl.gl <- data.frame(tFirst=twl.back[-nrow(twl.back),1], tSecond=twl.back[-1,1], type=ifelse(twl.back[,2],1,2)[-nrow(twl.back)])

### Fix first and last two locations at site of deployment/re-capture
fixed <- matrix(FALSE, nrow = nrow(twl.gl), ncol = 2)
fixed[c(1:2, (nrow(twl.gl)-1):nrow(twl.gl)),1] <- TRUE
fixed[c(1:2, (nrow(twl.gl)-1):nrow(twl.gl)),2] <- TRUE

### changeLight analysis. User can define quantile and minimum stopover period
cL <- changeLight(twl = twl.gl, fixed = fixed, quantile = .95, days = 1, summary = T) 

### Set xlim/ylim for maps below
xlim = c(-95, -65)
ylim = c(10, 52)

### Create map of possible stopover locations. You may need to vary quantile and minimum stopover period above to 
### achieve good results. For our analysis, this process did not work for fall migration, presumably due to overlap
### with the autumnal equinox
siteMap(cbind(zm$Lon.mean,zm$Lat.mean),site=cL$site,xlim=xlim,ylim=ylim,type='cross',hull=T)
map('state', fill = T, col = "transparent",border = "grey80", add = T)

### Merge sites based on distance between sites. Sometimes the above process will identify multiple sites in close
### proximity. Changing distThreshold will allow you to set this distance. Note - can take a few minutes of processing
mS200 <- mergeSites(twl = twl.gl, fixed = fixed, site = cL$site, degElevation = 90-zenith0, distThreshold = 200)

### Create new map with merged sties.
siteMap(cbind(zm$Lon.mean, zm$Lat.mean),site=mS200$site,xlim=xlim,ylim=ylim,type='cross',hull=T)
map('state', fill = T, col = "transparent",border = "grey80", add = T)

### View departure/arrival schedule. One could use this to identify departure and arrival dates. We did not, but 
### estimated dates were similar in general. 
schedule(twl.gl$tFirst, twl.gl$tSecond, site = mS200$site)

######################################################################
########## DETERMINING SPRING MIGRATION DISTANCE #####################
######################################################################

### The basic idea here is that if you calculate the distance the along mean path you are incorporating a lot of
### error into the estimate because of the inherent error in light-level geolocation. To compensate for this we assumed
### each bird was stationary during its stopover periods identified above. We then calculated distances from the 
### wintering entroid along the mean path to each stopover centroid, and finally to the location of geolocator recovery.

### Calculate mean migration path
Cxm <- locationMean(fit$z)
springmig<-Cxm[633:658,]

### Calculate winter centroid
WinterCentroid <- colMeans(Cxm[260:518,])

### Calculate Centroid of first stopover
Stop1Centroid<- colMeans(Cxm[632:639,])

### Calculate Centroid of second stopover
Stop2Centroid<- colMeans(Cxm[641:653,])

### Enter location of geolocator recovery
BreedingLoc<-c(-84.647636,44.655523)

### Here you add in the segments of the mean migration path that are in between the stationary periods
DistMigPath<-WinterCentroid
DistMigPath<-rbind(DistMigPath,Stop1Centroid)
DistMigPath<-rbind(DistMigPath,Cxm[640,])
DistMigPath<-rbind(DistMigPath,Stop2Centroid)
DistMigPath<-rbind(DistMigPath,Cxm[654:658,])
DistMigPath<-rbind(DistMigPath,BreedingLoc)

### Visualize difference between path for migration distance and mean migration path. They should be very similar, but
### path for migration distance should smooth out some of the error in the mean migration path
plot(Americas,xlim=xlim,ylim=ylim,col = "grey90",border = "grey10")
lines(DistMigPath, lwd = 1.5, col=rgb(t(col2rgb("black"))/255,alpha=0.8))
lines(springmig, lwd = 1.5, col=rgb(t(col2rgb("red"))/255,alpha=0.8))

### Calculate great circle distance between points along the path
GL14Dists <- (sapply(2:nrow(DistMigPath),function(i){distm(DistMigPath[i-1,],DistMigPath[i,], fun = distVincentySphere)}))/1000

### Sum distance to calculate final migration distance
sum(GL14Dists)