###     Determine first and last lat and long on migration        
##############################################################
library(rgdal)
data3 <- subset(data2,abs(data2$long.m.spd) < 50)

## First and last fix per bird
first.dts <- aggregate(data3$dt,by=list(data3$name),FUN='min')
colnames(first.dts) <- c("name","start")
attributes(first.dts$start)$tzone <- "UTC"  

last.dts <- aggregate(data3$dt,by=list(data3$name),FUN='max')
colnames(last.dts) <- c("name","end")
attributes(last.dts$end)$tzone <- "UTC"  

# Merge starts and ends to df
startends <- merge(first.dts,last.dts)
data3 <- merge(data3,startends,all.x=T)

## Determine end lat in degrees
end.lats <- data3[which(data3$dt == data3$end),c("name","lat")]
colnames(end.lats)[2] <- c("end.lat")
data3 <- merge(data3,end.lats,all.x=T)

## Determine end long in degrees
end.longs <- data3[which(data3$dt == data3$end),c("name","long")]
colnames(end.longs)[2] <- c("end.long")
data3 <- merge(data3,end.longs,all.x=T)

## Determine start long in degrees
st.longs <- data3[which(data3$dt == data3$start),c("name","long")]
colnames(st.longs)[2] <- c("st.long")
data3 <- merge(data3,st.longs,all.x=T)

#subset data to include only flight: speed threshold 5km per hour over sea and during daylight hours
data3 <- data3[which(data3$spd >=5 &((data3$hr >= 7 & data3$hr < 21)| data3$landmask == 'sea')),]

# set spd in scale m/s instead of km/h
data3$spd.m <- data3$spd / 3.6

# remove outliers in terms of spd
data3 <- subset(data3,data3$spd.m <=25)

data3 <- data3[order(data3$migr,data3$dt),]

#### Run Multiple Linear Regression Models        
####################################################

m1 <- lm(long.m.spd~ u ,data=data3)
m2 <- lm(long.m.spd~ u + v,data=data3)
m3 <- lm(long.m.spd~ u * v,data=data3)

(m<-anova(m1,m2))
(m<-anova(m1,m3))
(m<-anova(m2,m3))

data3$res <- data3$long.m.spd - predict(m3)
outliers <- which(abs(data3$res)>14)
data3 <- data3[-outliers,]

m1 <- lm(long.m.spd~ u ,data=data3)
m2 <- lm(long.m.spd~ u + v,data=data3)
m3 <- lm(long.m.spd~ u * v,data=data3)

(m<-anova(m1,m2))
(m<-anova(m1,m3))
(m<-anova(m2,m3))

#### PLOT 3D RESPONSE SURFACE AND RESIDUAL VALUES
#######################################################
library(rgl)

myColorRamp <- function(colors, values) {
    v <- (values - min(values))/diff(range(values))
    x <- colorRamp(colors)(v)
    rgb(x[,1], x[,2], x[,3], maxColorValue = 255)
}

cols <- myColorRamp(c("blue", "yellow","red"), data3$res) 

newdat <- expand.grid(u=seq(-20,20,by=.5),v=seq(-15,15,by=.5))
newdat$long.m.spd <- predict(m3,newdat)
dev.new()
with(data3,plot3d(u,v,long.m.spd, col=cols, size=.5,type="s", main=""))
with(newdat,surface3d(unique(u),unique(v),long.m.spd,
                      alpha=0.3,front="line", back="line"))

rgl.snapshot("Fig2A.png", fmt = "png", top = TRUE ) 

### PLOT FIGURE 2B
####################################################
## To do this replace U-wind component in script for figure 1 by model residuals

### MODEL DELTA-LONG AS A FUNCTION OF WIND EN ROUTE
#######################################################
meanu <- aggregate(data3$u,by=list(data3$name),FUN='mean')
colnames(meanu)[1:2] <- c("name","meanu")

meanv <- aggregate(data3$v,by=list(data3$name),FUN='mean')
colnames(meanv)[1:2] <- c("name","meanv")

finlong <- unique(data3[,c("name","st.long","end.long","end.lat","cat")])

birds <- merge(meanu,finlong)
birds <- merge(birds,meanv)
# Subset individuals that survived until end of migration based on final lat
survivors <- subset(birds,birds$end.lat.d < 20)

m1 <- lm((end.long-st.long)~meanu,data=survivors)
m2 <- lm((end.long-st.long)~meanu+meanv,data=survivors)
m3 <- lm((end.long-st.long)~meanu*meanv,data=survivors)

m4 <- lm((end.long-st.long)~meanu+meanv+cat,data=survivors)
m5 <- lm((end.long-st.long)~meanu*meanv+cat,data=survivors)

(m<-anova(m1,m2))
(m<-anova(m2,m3))

(m<-anova(m2,m4))
(m<-anova(m3,m5))

survivors$pred <- predict(m2)
survivors$obs <- survivors$end.long - survivors$st.long

#### PLOT FIGURE 3
###############################################
# determine which outleirs need to be labelled
outliers <- c("Matti","Hans","Venus")

# Generate plot
ggplot() + geom_point(data=survivors,aes(x=obs,y=pred))+
	theme_bw() + ylab('Predicted ?long') + xlab('Observed ?long') +
	geom_abline(intercept=0,slope=1)+
	geom_text(data=survivors[which(survivors$name %in% outliers),],aes(obs,pred,label=name),hjust = 0, nudge_x = 0.25)+
  	theme(legend.position = 'bottom',
	  	legend.direction= 'horizontal',
		plot.background = element_rect(fill='white'),
		panel.grid.major	= element_line(colour='grey50',linetype='dashed'),
		panel.grid.minor	= element_line(colour='grey70',linetype='dashed'),
		legend.background = element_rect(fill='white'),
		panel.background	= element_rect(fill='white'),
		strip.background	= element_rect(fill='white'),
		legend.title	= element_text(size=10,colour='black'),
		legend.text		= element_text(size=10,colour='black'),
		strip.text		= element_text(size=14,colour='black',face='bold'),
		axis.text		= element_text(size=12,colour='black',face='italic'),
		axis.title		= element_text(size=12,colour='black',face='bold'))
ggsave('Fig3.tiff',width=4,height=4,dpi=300)



