##########################################################################
####    CODE ACCOMPANYING VANSTEELANT ET AL 2017 - PROC R SOC B       ####
###                     by WMG VANSTEELANT                              ##
###                (vansteelant.wouter@gmail.com)                       ##
##                                                                      ##                                                                      ##                                                                      ##
##                                                                      ##                                                                      ##                                                                      ##
##   LAST UPDATE: 23.02.2017                                            ##
##########################################################################
setwd('C:/Users/wvanste1/Documents/2015 HB EOU')

### READ TRACKING DATA ANNOTATED WITH WIND AND ENVIRONMENTAL DATA FOR MAPPING
#############################################################################
# Set working directory and read data downloaded as online supplementary materials with this article
setwd('...')
data2 <- read.csv('Vansteelant-etal-ProcRSocB-fin.csv')
data2$dt <- as.POSIXct(strptime(data2$dt, format="%Y-%m-%d %H:%M:%S"), tz='UTC')

### PROCESS ENVIRONMENTAL DATA FOR MAPPING
##########################################
### SPECIFY STUDY RANGE FOR MAPPING
latlimits <- c(-10, 65) 
longlimits <- c(-15, 45) 

### OBTAIN MAP 
library(raster)
library(tiff)
library(maptools)
# Read any DEM stored in your working directory
dem <- raster('...')

# Crop DEM to study range
dem2 <- crop(dem, extent(longlimits[1],longlimits[2],latlimits[1],latlimits[2]))
dem3 <- aggregate(dem2,5,FUN='mean')
map2 <- rasterToPoints(dem3)
map2 <- data.frame(map2)
colnames(map2) <- c("long","lat","alt")

### ADD COASTLINES
# read high resolution coastlines downloaded from
# gshhs.h.b <- "./gshhg-bin-2.3.0/gshhs_h.b"
# create spatial object out of coastlines
shore.hr <- getRgshhsMap(gshhs.h.b, xlim = longlimits, ylim = latlimits)
shore.high <- fortify(shore.hr)

### EXTRACT DEM within coastlines only (to plot borders Caspian and Black Sea)
map.hr.sub <- mask(dem3, shore.hr)
map.hr.suba <- rasterToPoints(map.hr.sub)
map.hr.suba <- data.frame(map.hr.suba)
colnames(map.hr.suba) <- c("long","lat","alt")

### PRODUCE FIGURE 1 
#########################
# Create colour palette for colouring DEM
cols <- colorRampPalette(c("grey40","grey50","grey60","grey70","grey80","grey90"))(266)
breaks = quantile(map2$alt, seq(0, 1, length.out=256))

# Set colour scale for zonal wind
colsu <- colorRampPalette(c("darkblue","blue","cornflowerblue","yellow","orange","red","darkred"))(32)
breaksu = seq(-22, 22, length.out=31)

# Plot using ggplot2
library(ggplot2)
data2 <- data2[order(data2$migr,data2$dt),]

# Produce map
dev.new()
ggplot() + 
	geom_polygon(data = shore.high, aes(x=long, y = lat, group = group), fill= 'grey10', color = "black") +
	geom_raster(map.hr.suba,mapping=aes(long,lat,fill=alt,alpha=alt))+
	scale_fill_gradientn(colours=cols,breaks=breaks)+
	scale_alpha_continuous(range = c(0.6, 1.0))+
	coord_cartesian(xlim=c(-18,36),ylim=c(0,65)) +
	scale_x_continuous(limits = c(-18,36), expand = c(0, 0)) +
  	scale_y_continuous(limits = c(0,65), expand = c(0, 0)) +
	geom_path(data=data2,aes(long,lat,group=migr,colour=u),size=.9,alpha=.7)+
  	theme(legend.position = 'none')+
	scale_colour_gradientn(colours=colsu,breaks=breaksu,limits=c(-16,16))+
	facet_grid(~season2)+
  	theme(legend.position = 'none',
	  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'))

# Save outputted graph
#ggsave('Figure1.tiff',dpi=300,width=15,height=7)

###     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)

#### Run Multiple Linear Regression Models        
####################################################
#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),]

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'))

## Save output
#ggsave('Fig3.tiff',width=4,height=4,dpi=300)
