#Rcodes to test the influence of the moon on the diet of Bulwer´s petrel #Full_Occurrences rm(list=ls()) require(vegan) diet<-read.table("/home/granadeiro/Dados/RotinasR/SilkeMoon/Scientific_Reports_Data.csv", sep=";", header=T) str(diet) diet$dateok<-as.POSIXct(diet$DATA.CERTA, format="%d-%m-%y", tz="GMT") ###write.table(diet, file="/home/granadeiro/Dados/RotinasR/SilkeMoon/DadosDieta.csv", sep=";", row.names=F) cloud<-read.table("/home/granadeiro/Dados/RotinasR/MODISClouds/dadosLua.csv", sep=";", header=T) cloud$dateok<-as.POSIXct(cloud$date, fomat="%Y-%m-%d", tz="GMT") #l<-merge(diet[, 1:8], cloud, all.x=F, all.y=T, by.x="dateok", by.y="dateok") l<-merge(diet, cloud, all.x=T, all.y=F, by.x="dateok", by.y="dateok") l$Yearf<-factor(l$Year) l<-na.omit(l) ##Data of Deserta 2013 des<-which(l$Local=="DS"&l$Year==2013) pres<-l[des, 9:83] var1<-data.frame(cloud=l$cloudpD[des], yearf=l$Yearf[des], Moon=l$Moon_Phase[des]) ##Checks for infrequent prey selc<-which((apply(pres, 2, sum)/dim(l)[1])<.05) pres1<-pres[,-selc] ##checks for empty samples sell<-which((apply(pres1, 1, sum)<1)) pres2<-pres1 if(length(sell)>0) { pres2<-pres1[-sell,] var1<-var1[-sell,] } presenca<-pres2 variaveis<-var1 row.names(presenca)<-1:(dim(presenca)[1]) ad<-adonis(presenca~cloud*Moon, data=variaveis) ##Data from Deserta and Selvagem (both sites combined) pres<-l[, 9:83] var1<-data.frame(cloud=l$cloudpD, yearf=l$Yearf, Moon=l$Moon_Phase, site=l$Local) ##Checks for infrequent prey selc<-which((apply(pres, 2, sum)/dim(l)[1])<.05) pres1<-pres[,-selc] ##checks for empty samples sell<-which((apply(pres1, 1, sum)<1)) pres2<-pres1 if(length(sell)>0) { pres2<-pres1[-sell,] var1<-var1[-sell,] } presenca<-pres2 variaveis<-var1 row.names(presenca)<-1:(dim(presenca)[1]) ad<-adonis(presenca~site+yearf+cloud*Moon, data=variaveis) ############ Figure for Supplementary material dados<-read.table("/home/granadeiro/Dados/RotinasR/MODISClouds/dadosLua.csv", sep=";", header=T) dados$dateok<-as.POSIXct(paste(dados$date, "00:00:01"), tz="UTC") grafpdf<-TRUE if(grafpdf) pdf("/home/granadeiro/Dados/RotinasR/SilkeMoon/FigureSupMat.pdf", width=10, height=10) par(las=2, mfrow=c(2,2)) Sys.setlocale("LC_TIME", "en_GB.UTF-8") ### 2012 ### DESERTA with(dados[year(dados$dateok)==2012,], plot(dateok, cloudpD, type="n", ylim=c(0,1), ylab="Cloud cover", xlab=NA)) dias<-unique(l$dateok[l$Local=="DS"&l$Year==2012]) for (x in 1:length(dias)) { rect(dias[x]-60*60*12, -.5, dias[x]+60*60*12, 1.5, col="lightgrey", border="lightgrey") } box() ptcol<-(dados[year(dados$dateok)==2012,"dateok"]-1)%in%dias with(dados[year(dados$dateok)==2012,], lines(dateok, cloudpD)) with(dados[year(dados$dateok)==2012,], points(dateok, cloudpD, pch=21, bg=c("white", "black")[ptcol+1])) title("Deserta island 2012") with(dados[year(dados$dateok)==2012,], lines(dateok, moon, col="red")) lua<-dados[year(dados$dateok)==2012, "moon"]*2*pi cover<-dados[year(dados$dateok)==2012, "cloudpD"] #cover<-rnorm(length(cover)) p<-cor.test(sin(lua), cover) print(paste("Correlation between cover and sin(moon) =", round(p$estimate,2), ", p =", round(p$p.value,3), ", n =", p$parameter)) ### 2012 ### SELVAGEM with(dados[year(dados$dateok)==2012,], plot(dateok, cloudpS, type="n", ylim=c(0,1), ylab="Cloud cover", xlab=NA)) dias<-unique(l$dateok[l$Local=="SG"&l$Year==2012]) for (x in 1:length(dias)) { rect(dias[x]-60*60*12, -.5, dias[x]+60*60*12, 1.5, col="lightgrey", border="lightgrey") } box() ptcol<-(dados[year(dados$dateok)==2012,"dateok"]-1)%in%dias with(dados[year(dados$dateok)==2012,], lines(dateok, cloudpS)) with(dados[year(dados$dateok)==2012,], points(dateok, cloudpS, pch=21, bg=c("white", "black")[ptcol+1])) title("Selvagem Grande 2012") with(dados[year(dados$dateok)==2012,], lines(dateok, moon, col="red")) lua<-dados[year(dados$dateok)==2012, "moon"]*2*pi cover<-dados[year(dados$dateok)==2012, "cloudpS"] #cover<-rnorm(length(cover)) p<-cor.test(cos(lua), cover) print(paste("Correlation between cover and sin(moon) =", round(p$estimate,2), ", p =", round(p$p.value,3), ", n =", p$parameter)) ### 2013 ### DESERTA with(dados[year(dados$dateok)==2013,], plot(dateok, cloudpD, type="n", ylim=c(0,1), ylab="Cloud cover", xlab=NA)) dias<-unique(l$dateok[l$Local=="DS"&l$Year==2013]) for (x in 1:length(dias)) { rect(dias[x]-60*60*12, -.5, dias[x]+60*60*12, 1.5, col="lightgrey", border="lightgrey") } box() ptcol<-(dados[year(dados$dateok)==2013,"dateok"]-1)%in%dias with(dados[year(dados$dateok)==2013,], lines(dateok, cloudpD)) with(dados[year(dados$dateok)==2013,], points(dateok, cloudpD, pch=21, bg=c("white", "black")[ptcol+1])) title("Deserta island 2013") with(dados[year(dados$dateok)==2013,], lines(dateok, moon, col="red")) #mod<-activity::fitlincirc(dados$moon*2*pi, dados$cloudpD, pCI = 0.95, reps = 1000, res = 512) #plot(mod, ) if(grafpdf) dev.off()