##Contents ##### 1. Data upload and preparation - line 26 ##### 2. Seasonality - line 384 ##### 3. Fourier spectra - line 631 ##### 4. Long-term trends - line 706 ##### 5. Periodicity over time - line 927 ##### 6. Influence of oceans - line 1003 ## - Load libraries library(plyr) library(ggplot2) library(zoo) library(lubridate) library(timeDate) library(scales) library(biwavelet) library(ggpubr) library(corrplot) library(lme4) library("itsadug") library("tidyr") library("broom") library("gridExtra") ########################################## ##### 1. Data upload and preparation ##### ########################################## ### Rainfall ### ## - Data available for download at the University of Stirling's DataSTORRE (http://hdl.handle.net/11667/133) ## - Upload daily rainfall data (Dataset C - a combination of datasets A and B with adjusted data) Rainfall_daily_C <- read.csv("Rainfall_daily_v.2017-12-31.csv") Rainfall_daily_C$Date<-as.Date(Rainfall_daily_C$Date) Rainfall_daily_C<-Rainfall_daily_C[order(Rainfall_daily_C$Date),] ## - Calculate mean value for each day in the time series and interpolate missing values using the 10 day rolling mean (to create Dataset D) Rainfall_daily_D<-ddply(Rainfall_daily_C,.(Date),summarise,Rainfall_daily_mm=mean(Rainfall_adjusted_mm,na.rm=T)) # - Create full calendar of dates Rainfall_daily_D<-merge(Rainfall_daily_D, data.frame(Date=seq.Date(from=min(Rainfall_daily_D$Date,na.rm=T),to=max(Rainfall_daily_D$Date,na.rm=T),by="day")) ,"Date",all.y=T) Rainfall_daily_D$Data<-as.factor(ifelse(is.na(Rainfall_daily_D$Rainfall_daily_mm),"Missing","Data")) # - Interpolate missing data Rainfall_daily_D$MeanRain_10day<-rollapply(Rainfall_daily_D$Rainfall_daily_mm,10,mean,na.rm=T,fill=NA) Rainfall_daily_D$Rainfall_interp_mm<-Rainfall_daily_D$Rainfall_daily_mm Rainfall_daily_D$Rainfall_interp_mm[is.na(Rainfall_daily_D$Rainfall_interp_mm)]<-Rainfall_daily_D$MeanRain_10day[is.na(Rainfall_daily_D$Rainfall_interp_mm)] Rainfall_daily_D$Comments<-ifelse(is.na(Rainfall_daily_D$Rainfall_daily_mm),"Interpolated from 10 day running mean", "NA") Rainfall_daily_D<-data.frame(Date=Rainfall_daily_D$Date, Rainfall_mm=Rainfall_daily_D$Rainfall_daily_mm,Rainfall_interp_mm=Rainfall_daily_D$Rainfall_interp_mm, Comments=Rainfall_daily_D$Comments, Site="SEGC") Rainfall_daily_D$DOY<-as.integer(format(Rainfall_daily_D$Date,"%j")) Rainfall_daily_D$Month<-month(Rainfall_daily_D$Date) Rainfall_daily_D$Season<-mapvalues(Rainfall_daily_D$Month,from=c(1:12),to=c("DJF","DJF","MAM","MAM","MAM","JJAS","JJAS","JJAS","JJAS","ON","ON","DJF")) Rainfall_daily_D$Year<-year(Rainfall_daily_D$Date) Rainfall_daily_D$Year_r<-scale(Rainfall_daily_D$Year) ## - Sum rainfall for each month in the time series from interpolated daily data and interpolate missing months (to create Dataset E) Rainfall_monthly_E<-ddply(Rainfall_daily_D,.(Year,Season,Month,YearMonth=as.yearmon(Date)),summarise, Rainfall_daily_mm=mean(Rainfall_mm,na.rm=T), Rainfall_mm=sum(Rainfall_interp_mm,na.rm=T), sample=length(Date[complete.cases(Rainfall_interp_mm)]), sample_max=length(Date)) # - Remove monthly values with more than 10% missing days Rainfall_monthly_E$sample_prop<-Rainfall_monthly_E$sample/Rainfall_monthly_E$sample_max Rainfall_monthly_E$Rainfall_mm[Rainfall_monthly_E$sample_prop<0.9]<-NA Rainfall_monthly_E$Rainfall_daily_mm[Rainfall_monthly_E$sample_prop<0.9]<-NA # - Complete missing values with mean value for corresponding calendar month Rainfall_monthly_E<-merge(Rainfall_monthly_E, ddply(Rainfall_monthly_E,.(Month=month(YearMonth)),summarise, Rainfall_month_mm=mean(Rainfall_mm,na.rm=T)), "Month") Rainfall_monthly_E$Rainfall_interp_mm<-Rainfall_monthly_E$Rainfall_mm Rainfall_monthly_E$Rainfall_interp_mm[is.na(Rainfall_monthly_E$Rainfall_interp_mm)]<-Rainfall_monthly_E$Rainfall_month_mm[is.na(Rainfall_monthly_E$Rainfall_interp_mm)] Rainfall_monthly_E<-Rainfall_monthly_E[order(Rainfall_monthly_E$YearMonth),] Rainfall_monthly_E$Rainfall_standard_mm<-scale(Rainfall_monthly_E$Rainfall_interp_mm) ## - Sum rainfall for each year in the timeseries from interpolated daily data (to create Dataset F) Rainfall_yearly_F<-ddply(Rainfall_daily_D,.(Year=year(Rainfall_daily_D$Date)),summarise, Rainfall_daily_mm=mean(Rainfall_interp_mm,na.rm=T), Rainfall_yearly_mm=sum(Rainfall_interp_mm)) Rainfall_yearly_F$Year_r<-scale(Rainfall_yearly_F$Year) ## - Summary of rainfall dataframes (match flowchart in supporting information) summary(Rainfall_daily_C) #Lopé daily rainfall dataset including original and adjusted data (Dataset C) summary(Rainfall_daily_D) #Lopé mean daily rainfall time series with non-interpolated and interpolated data (Dataset D) summary(Rainfall_monthly_E) #Lopé monthly rainfall time series with non-interpolated and interpolated data (Dataset E) summary(Rainfall_yearly_F) #Lopé yearly rainfall time series (Dataset F) ### Temperature #### ## - Data available for download at the University of Stirling's DataSTORRE (http://hdl.handle.net/11667/133) ## - Upload daily max/min temperature data (Dataset C - a combination of datasets A and B) Temperature_daily_C <- read.csv("Temperature_daily_v.2018-03-13.csv") Temperature_daily_C$Date<-as.Date(Temperature_daily_C$Date) ## - Calculate mean value for each day in the time series (to create Dataset D) Temperature_daily_D<-ddply(Temperature_daily_C,.(Site,Type,Date),summarise, Temperature_c=mean(Temperature_c,na.rm=T)) # - Create full calendar of dates Temperature_daily_D<-rbind(merge(Temperature_daily_D[Temperature_daily_D$Site=="Saline",], expand.grid(Date=seq.Date(from=min(Temperature_daily_D$Date[Temperature_daily_D$Site=="Saline"],na.rm=T),to=max(Temperature_daily_D$Date[Temperature_daily_D$Site=="Saline"],na.rm=T),by="day"), Type=c("Max","Min"),Site=c("Saline")),c("Date","Type","Site"),all=T), merge(Temperature_daily_D[Temperature_daily_D$Site=="SEGC",], expand.grid(Date=seq.Date(from=min(Temperature_daily_D$Date[Temperature_daily_D$Site=="SEGC"],na.rm=T),to=max(Temperature_daily_D$Date[Temperature_daily_D$Site=="SEGC"],na.rm=T),by="day"), Type=c("Max","Min"), Site=c("SEGC")), c("Date","Type","Site"),all=T)) Temperature_daily_D$DOY<-as.integer(format(Temperature_daily_D$Date,"%j")) Temperature_daily_D$YearMonth<-as.yearmon(Temperature_daily_D$Date) Temperature_daily_D$Month<-month(Temperature_daily_D$Date) Temperature_daily_D$Year<-year(Temperature_daily_D$Date) Temperature_daily_D$Site2<-ifelse(Temperature_daily_D$Site=="Saline","Forest","Savanna") ## - Calculate monthly mean daily max/min temperature time series for both sites (Dataset E) Temperature_monthly_E<-ddply(Temperature_daily_D[Temperature_daily_D$Year<2018,],.(Site,Type,Year,Month,YearMonth),summarise, Sample_Lopé=length(Temperature_c[complete.cases(Temperature_c)]), Temperature_c=mean(Temperature_c,na.rm=T)) Temperature_monthly_E$Temperature_c[Temperature_monthly_E$Sample_Lopé<5]<-NA Temperature_monthly_E$Date<-as.Date(paste(year(Temperature_monthly_E$YearMonth),month(Temperature_monthly_E$YearMonth),"15",sep="-")) # - Complete missing values with mean value for corresponding calendar month Temperature_monthly_E<-merge(Temperature_monthly_E, ddply(Temperature_monthly_E,.(Site,Type,Month),summarise, Temperature_month_mean_c=mean(Temperature_c,na.rm=T)), c("Site","Type","Month")) Temperature_monthly_E$Temperature_interp_c<-Temperature_monthly_E$Temperature_c Temperature_monthly_E$Temperature_interp_c[is.na(Temperature_monthly_E$Temperature_interp_c)]<-Temperature_monthly_E$Temperature_month_mean_c[is.na(Temperature_monthly_E$Temperature_interp_c)] Temperature_monthly_E<-Temperature_monthly_E[order(Temperature_monthly_E$YearMonth),] Temperature_monthly_E$Temperature_standard_c<-NA Temperature_monthly_E$Temperature_standard_c[Temperature_monthly_E$Site=="Saline"&Temperature_monthly_E$Type=="Max"]<-scale(Temperature_monthly_E$Temperature_interp_c[Temperature_monthly_E$Site=="Saline"&Temperature_monthly_E$Type=="Max"]) Temperature_monthly_E$Temperature_standard_c[Temperature_monthly_E$Site=="SEGC"&Temperature_monthly_E$Type=="Max"]<-scale(Temperature_monthly_E$Temperature_interp_c[Temperature_monthly_E$Site=="SEGC"&Temperature_monthly_E$Type=="Max"]) Temperature_monthly_E$Temperature_standard_c[Temperature_monthly_E$Site=="Saline"&Temperature_monthly_E$Type=="Min"]<-scale(Temperature_monthly_E$Temperature_interp_c[Temperature_monthly_E$Site=="Saline"&Temperature_monthly_E$Type=="Min"]) Temperature_monthly_E$Temperature_standard_c[Temperature_monthly_E$Site=="SEGC"&Temperature_monthly_E$Type=="Min"]<-scale(Temperature_monthly_E$Temperature_interp_c[Temperature_monthly_E$Site=="SEGC"&Temperature_monthly_E$Type=="Min"]) ## - Calculate long-term series of minimum daily temperature combined from both sites and all equipment (Dataset F) Temperature_daily_F<-ddply(Temperature_daily_D[Temperature_daily_D$Type=="Min",],.(Year, Month=month(Date), YearMonth, DOY,Date),summarise, Temperature_c=mean(Temperature_c,na.rm=T)) Temperature_daily_F<-Temperature_daily_F[Temperature_daily_F$Year<2018,] # - Interpolate missing data Temperature_daily_F<-Temperature_daily_F Temperature_daily_F$MeanTemp_10day<-rollapply(Temperature_daily_F$Temperature_c,10,mean,na.rm=T,fill=NA) Temperature_daily_F$Temperature_interp_c<-Temperature_daily_F$Temperature_c Temperature_daily_F$Temperature_interp_c[is.na(Temperature_daily_F$Temperature_c)]<-Temperature_daily_F$MeanTemp_10day[is.na(Temperature_daily_F$Temperature_c)] Temperature_daily_F$Comments<-ifelse(is.na(Temperature_daily_F$Temperature_c),"Interpolated from 10 day running mean", "NA") Temperature_daily_F$Year_r<-scale(Temperature_daily_F$Year) Temperature_daily_F$Season<-mapvalues(Temperature_daily_F$Month,from=c(1:12),to=c("DJF","DJF","MAM","MAM","MAM","JJAS","JJAS","JJAS","JJAS","ON","ON","DJF")) ## - Calculate mean monthly minimum temperaure from combined daily dataset (to create dataset G) Temperature_monthly_G<-ddply(Temperature_daily_F,.(Year,Season,Month,YearMonth),summarise, Sample_Lopé=length(Temperature_c[complete.cases(Temperature_c)]), Temperature_c=mean(Temperature_c,na.rm=T)) Temperature_monthly_G$Temperature_c[Temperature_monthly_G$Sample_Lopé<5]<-NA Temperature_monthly_G$Date<-as.Date(paste(year(Temperature_monthly_G$YearMonth),month(Temperature_monthly_G$YearMonth),"15",sep="-")) # - Complete missing values with mean value for corresponding calendar month Temperature_monthly_G<-merge(Temperature_monthly_G, ddply(Temperature_monthly_G,.(Month),summarise, Temperature_month_mean_c=mean(Temperature_c,na.rm=T)), c("Month")) Temperature_monthly_G$Temperature_interp_c<-Temperature_monthly_G$Temperature_c Temperature_monthly_G$Temperature_interp_c[is.na(Temperature_monthly_G$Temperature_interp_c)]<-Temperature_monthly_G$Temperature_month_mean_c[is.na(Temperature_monthly_G$Temperature_interp_c)] Temperature_monthly_G<-Temperature_monthly_G[order(Temperature_monthly_G$YearMonth),] ## - Summary of temperature datasets summary(Temperature_daily_C) #Lopé daily temperature dataset including original data from each equipment at each site (Dataset C) summary(Temperature_daily_D) #Lopé mean daily temperature dataset combining data from different equipment at each site (Dataset D) summary(Temperature_monthly_E) #Lopé mean monthly temperature dataset for each site including non-interpolated and inteprolated data (Dataset E) summary(Temperature_daily_F) #Lopé mean daily minimum temperature dataset combining non-interpolated data from each site (Dataset F) summary(Temperature_monthly_G) #Lopé mean monthly minimum daily temperature dataset for both sites combined including non-interpolated and inteprolated data (Dataset G) ### Berkeley Temp ### ## - download available from http://climexp.knmi.nl/start.cgi for the grid-cell overlapping the SEGC location (0.2N, 11.6E). Berkeley_daily <- read.csv("Temperature_Berkeley_Min_daily.csv") names(Berkeley_daily)<-c("Date","Temperature_c") Berkeley_daily$Date<-as.Date(as.character(Berkeley_daily$Date),"%Y%m%d") Berkeley_daily$DOY<-format(Berkeley_daily$Date,"%j") Berkeley_daily$Year<-year(Berkeley_daily$Date) Berkeley_daily$Month<-month(Berkeley_daily$Date) Berkeley_daily$Year_r<-scale(Berkeley_daily$Year) Berkeley_daily<-Berkeley_daily[Berkeley_daily$Year>1983&Berkeley_daily$Year<2018,] ### CRU Temp ### ## - downloaded from http://climexp.knmi.nl/start.cgi for the grid-cell overlapping the SEGC location (0.2N, 11.6E). CRU_monthly <- read.table("Temperature_CRU_Min_monthly.txt") names(CRU_monthly)<-c("Year",1:12) CRU_monthly<-gather(CRU_monthly,"Month","Temperature_c",2:13) CRU_monthly$YearMonth<-as.yearmon(paste(CRU_monthly$Year,CRU_monthly$Month),"%Y %m") CRU_monthly$Year_r<-scale(CRU_monthly$Year) CRU_monthly<-CRU_monthly[CRU_monthly$Year>1983,] ### Absolute humidity ### ## - Data available for download at the University of Stirling's DataSTORRE (http://hdl.handle.net/11667/133) ## - Upload daily absolute humidity record (Dataset A) AbsHumidity_daily_A <- read.csv("Humidity_daily_v.2018-03-14.csv") AbsHumidity_daily_A$Date<-as.Date(AbsHumidity_daily_A$Date,format="%d/%m/%y") AbsHumidity_daily_A<-AbsHumidity_daily_A[complete.cases(AbsHumidity_daily_A$Date),] AbsHumidity_daily_A$Site2<-ifelse(AbsHumidity_daily_A$Site=="Saline","Forest","Savanna") AbsHumidity_daily_B<-ddply(AbsHumidity_daily_A[AbsHumidity_daily_A$Night==T,],.(Site=Site2,Date),summarise,AbsoluteHumidity_gm3=mean(AbsoluteHumidity_gm3,na.rm=T)) AbsHumidity_daily_B<-merge(expand.grid(Date=seq(from=min(AbsHumidity_daily_B$Date),to=max(AbsHumidity_daily_B$Date),"day"), Site=c("Savanna","Forest")), AbsHumidity_daily_B,c("Date","Site"),all.x=T) AbsHumidity_daily_B$Year<-year(AbsHumidity_daily_B$Date) AbsHumidity_daily_B$Month<-month(AbsHumidity_daily_B$Date) AbsHumidity_daily_B$YearMonth<-as.yearmon(AbsHumidity_daily_B$Date) AbsHumidity_daily_B$DOY<-as.integer(format(AbsHumidity_daily_B$Date,"%j")) ## - Calculate monthly mean time series for both sites (Dataset C) AbsHumidity_monthly_C<-ddply(AbsHumidity_daily_B[AbsHumidity_daily_B$Year<2018,],.(Site,Year, Month,YearMonth),summarise, Sample_Lopé=length(AbsoluteHumidity_gm3[complete.cases(AbsoluteHumidity_gm3)]), AbsoluteHumidity_gm3=mean(AbsoluteHumidity_gm3,na.rm=T)) AbsHumidity_monthly_C$Date<-as.Date(paste(AbsHumidity_monthly_C$Year,AbsHumidity_monthly_C$Month,"15",sep="-")) AbsHumidity_monthly_C$AbsoluteHumidity_gm3[AbsHumidity_monthly_C$Sample_Lopé<5]<-NA #Complete missing values with mean value for corresponding calendar month AbsHumidity_monthly_C<-merge(AbsHumidity_monthly_C, ddply(AbsHumidity_monthly_C,.(Site,Month),summarise,AbsoluteHumidity_month_mean_gm3=mean(AbsoluteHumidity_gm3,na.rm=T)), c("Site","Month")) AbsHumidity_monthly_C$AbsHumidity_interp_gm3<-AbsHumidity_monthly_C$AbsoluteHumidity_gm3 AbsHumidity_monthly_C$AbsHumidity_interp_gm3[is.na(AbsHumidity_monthly_C$AbsoluteHumidity_gm3)]<-AbsHumidity_monthly_C$AbsoluteHumidity_month_mean_gm3[is.na(AbsHumidity_monthly_C$AbsoluteHumidity_gm3)] AbsHumidity_monthly_C$AbsHumidity_standard_gm3<-scale(AbsHumidity_monthly_C$AbsHumidity_interp_gm3) AbsHumidity_monthly_C<-AbsHumidity_monthly_C[order(AbsHumidity_monthly_C$Site,AbsHumidity_monthly_C$Date),] summary(AbsHumidity_daily_A) summary(AbsHumidity_daily_B) summary(AbsHumidity_monthly_C) ### Wind ### ## - Data available for download at the University of Stirling's DataSTORRE (http://hdl.handle.net/11667/133) ## - Upload daily mean wind speed record (Dataset A) Wind_daily_A <- read.csv("Wind_daily_v.2016-07-01.csv") names(Wind_daily_A)[3]<-"Wind_m_s" Wind_daily_A$Date<-as.Date(Wind_daily_A$Date) Wind_daily_B<-ddply(Wind_daily_A,.(Date),summarise, Wind_m_s=mean(Wind_m_s,na.rm=T)) Wind_daily_B<-merge(data.frame(Date=seq.Date(from=min(Wind_daily_B$Date,na.rm=T),to=max(Wind_daily_B$Date,na.rm=T),by="day")), Wind_daily_B,"Date",all.x=T) Wind_daily_B$Year<-year(Wind_daily_B$Date) Wind_daily_B$Month<-month(Wind_daily_B$Date) Wind_daily_B$YearMonth<-as.yearmon(Wind_daily_B$Date) Wind_daily_B$DOY<-as.integer(format(Wind_daily_B$Date,"%j")) ## - Calculate monthly mean time series for both sites (Dataset C) Wind_monthly_C<-ddply(Wind_daily_B[Wind_daily_B$Year<2018,],.(Year, Month,YearMonth),summarise, Sample_Lopé=length(Wind_m_s[complete.cases(Wind_m_s)]), Wind_m_s=mean(Wind_m_s,na.rm=T)) Wind_monthly_C$Date<-as.Date(paste(Wind_monthly_C$Year,Wind_monthly_C$Month,"15",sep="-")) Wind_monthly_C$Wind_m_s[Wind_monthly_C$Sample_Lopé<5]<-NA # - Complete missing values with mean value for corresponding calendar month Wind_monthly_C<-merge(Wind_monthly_C, ddply(Wind_monthly_C,.(Month),summarise,Wind_month_mean_m_s=mean(Wind_m_s,na.rm=T)), c("Month")) Wind_monthly_C$Wind_interp_m_s<-Wind_monthly_C$Wind_m_s Wind_monthly_C$Wind_interp_m_s[is.na(Wind_monthly_C$Wind_m_s)]<-Wind_monthly_C$Wind_month_mean_m_s[is.na(Wind_monthly_C$Wind_m_s)] Wind_monthly_C$Wind_standard_m_s<-scale(Wind_monthly_C$Wind_interp_m_s) Wind_monthly_C<-Wind_monthly_C[order(Wind_monthly_C$Date),] summary(Wind_daily_A) summary(Wind_daily_B) summary(Wind_monthly_C) ### Solar ### ## - Data available for download at the University of Stirling's DataSTORRE (http://hdl.handle.net/11667/133) ## - Upload daily mean solar radiation record (Dataset A) Solar_daily_A<- read.csv("Solar_daily_v.2015-10-30.csv") names(Solar_daily_A)[3]<-"Solar_Wm2" Solar_daily_A$Date<-as.Date(Solar_daily_A$Date) Solar_daily_A$Month<-month(Solar_daily_A$Date) Solar_daily_B<-ddply(Solar_daily_A,.(Date),summarise,Solar_Wm2 = mean(Solar_Wm2,na.rm=T)) Solar_daily_B<-merge(data.frame(Date=seq(from=min(Solar_daily_B$Date),to=max(Solar_daily_B$Date),"day")), Solar_daily_B,c("Date"),all.x=T) Solar_daily_B$Year<-year(Solar_daily_B$Date) Solar_daily_B$Month<-month(Solar_daily_B$Date) Solar_daily_B$YearMonth<-as.yearmon(Solar_daily_B$Date) Solar_daily_B$DOY<-as.integer(format(Solar_daily_B$Date,"%j")) Solar_daily_B$Solar_Wm2[which(Solar_daily_B$Solar_Wm2<10)]<-NA ## - Calculate monthly mean time series for both sites (Dataset C) Solar_monthly_C<-ddply(Solar_daily_B[Solar_daily_B$Year<2018,],.(Year, Month,YearMonth),summarise, Sample_Lopé=length(Solar_Wm2[complete.cases(Solar_Wm2)]), Solar_Wm2=mean(Solar_Wm2,na.rm=T)) Solar_monthly_C$Date<-as.Date(paste(Solar_monthly_C$Year,Solar_monthly_C$Month,"15",sep="-")) Solar_monthly_C$Solar_Wm2[Solar_monthly_C$Sample_Lopé<5]<-NA # - Complete missing values with mean value for corresponding calendar month Solar_monthly_C<-merge(Solar_monthly_C, ddply(Solar_monthly_C,.(Month),summarise,Solar_month_mean_Wm2=mean(Solar_Wm2,na.rm=T)), c("Month")) Solar_monthly_C$Solar_interp_Wm2<-Solar_monthly_C$Solar_Wm2 Solar_monthly_C$Solar_interp_Wm2[is.na(Solar_monthly_C$Solar_Wm2)]<-Solar_monthly_C$Solar_month_mean_Wm2[is.na(Solar_monthly_C$Solar_Wm2)] Solar_monthly_C$Solar_standard_Wm2<-scale(Solar_monthly_C$Solar_interp_Wm2) Solar_monthly_C<-Solar_monthly_C[order(Solar_monthly_C$Date),] summary(Solar_daily_A) summary(Solar_daily_B) summary(Solar_monthly_C) ### AOT ### ## - Download available from the NASA Aerosol Robotic Network (Aeronet; https://aeronet.gsfc.nasa.gov/; Holben et al. 1998). Aerosol_daily_A <- read.csv("Aeronet_L2_140101_181231_SEGC_Lope_Gabon.lev20 2", comment.char="#") Aerosol_daily_A$Date<-Aerosol_daily_A$Date.dd.mm.yy. Aerosol_daily_A$Date<-as.Date(as.character(Aerosol_daily_A$Date),"%d:%m:%Y") Aerosol_daily_B<-merge(data.frame(Date=seq.Date(from=min(Aerosol_daily_A$Date,na.rm=T),to=max(Aerosol_daily_A$Date,na.rm=T),by="day")), Aerosol_daily_A[,c(67,3,7,13,16)],"Date",all.x=T) Aerosol_daily_B$Year<-year(Aerosol_daily_B$Date) Aerosol_daily_B$Month<-month(Aerosol_daily_B$Date) Aerosol_daily_B$YearMonth<-as.yearmon(Aerosol_daily_B$Date) Aerosol_daily_B$DOY<-as.integer(format(Aerosol_daily_B$Date,"%j")) #summary(Aerosol_daily_B$Month[Aerosol_daily_B$Data=="Missing"]) #642 days missing out of 1045 possible days between 2014-04-27 and 2017-03-06 (61% missing). #More than 70% values missing in months June-November (peaking in August, 85%). Presumably because of cloud cover data removed at quality control. #Most data available in March (35% missing vales) ## - Calculate monthly mean time series for both sites (Dataset C) Aerosol_monthly_C<-ddply(Aerosol_daily_B[Aerosol_daily_B$Year<2018,],.(Year, Month,YearMonth),summarise, Sample_Lopé=length(AOT_500[complete.cases(AOT_500)]), AOT_500=mean(AOT_500,na.rm=T)) Aerosol_monthly_C$Date<-as.Date(paste(Aerosol_monthly_C$Year,Aerosol_monthly_C$Month,"15",sep="-")) Aerosol_monthly_C$AOT_500[Aerosol_monthly_C$Sample_Lopé<5]<-NA # - Complete missing values with mean value for corresponding calendar month Aerosol_monthly_C<-merge(Aerosol_monthly_C, ddply(Aerosol_monthly_C,.(Month),summarise,AOT_500_month_mean=mean(AOT_500,na.rm=T)), c("Month")) Aerosol_monthly_C$AOT_500_interp<-Aerosol_monthly_C$AOT_500 Aerosol_monthly_C$AOT_500_interp[is.na(Aerosol_monthly_C$AOT_500)]<-Aerosol_monthly_C$AOT_500_month_mean[is.na(Aerosol_monthly_C$AOT_500)] Aerosol_monthly_C$AOT_500_standard<-scale(Aerosol_monthly_C$AOT_500_interp) Aerosol_monthly_C<-Aerosol_monthly_C[order(Aerosol_monthly_C$Date),] summary(Aerosol_daily_A) summary(Aerosol_daily_B) summary(Aerosol_monthly_C) ########################## ##### 2. Seasonality ##### ########################## ### Rainfall ### ## - Calculate rainfall means for each DOY and month Rainfall_DOY<-ddply(Rainfall_daily_D,.(DOY),summarise, MeanRain_mm=mean(Rainfall_interp_mm,na.rm=T)) Rainfall_DOY$MeanRain_10day_mm<-rollapply(Rainfall_DOY$MeanRain_mm,10,mean,partial=T) Rainfall_DOY$Dayof2019<-as.Date(as.character(Rainfall_DOY$DOY),"%j") Rainfall_DOY$Month<-month(Rainfall_DOY$Dayof2019) Rainfall_DOY<-merge(Rainfall_DOY,ddply(Rainfall_DOY,.(Month),summarise,MeanRain_month_mm=mean(MeanRain_mm,na.rm=T)),"Month") Rainfall_DOY$month_labels <- mapvalues(Rainfall_DOY$Month,from=1:12,to=c('J','F','M','A','M','J','J','A','S','O','N','D')) # - Seasonality plot (Figure 2) Rainfall_seasonal_plot<-ggplot(Rainfall_DOY,aes(x=Dayof2019))+ geom_line(data=Rainfall_DOY,aes(y=MeanRain_mm),alpha=0.7,colour="gray",size=0.4)+ geom_line(data=Rainfall_DOY, aes(y=MeanRain_10day_mm),alpha=0.7,colour="black",size=0.4)+ geom_line(data=Rainfall_DOY, aes(y=MeanRain_month_mm),alpha=1,size=0.7)+ scale_x_date(breaks=date_breaks("1 month"),labels = function(x) substring(format(as.Date(x),"%b"), 1, 1) )+ theme_classic()+ ylab("Rainfall (mm/day)")+ xlab("Month")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-03-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-06-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-10-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-12-01"))),linetype="dotted")+ #ggtitle("A.")+ theme(plot.margin=unit(c(0.5,0.5,0.5,1),"cm")) ### Temperature ### ## - Calculate temprature means for each DOY and month Temperature_DOY<-ddply(Temperature_daily_D,.(Site,Type,Month=month(Date),DOY),summarise, Temperature_c=mean(Temperature_c,na.rm=T)) Temperature_DOY<-data.frame(Temperature_DOY,Temperature_c_10days=ddply(Temperature_DOY,.(Site,Type),summarise, Temperature_c_10days=rollapply(Temperature_c,10,mean,partial=T))[,3]) Temperature_DOY<-merge(Temperature_DOY,ddply(Temperature_DOY,.(Site,Type,Month),summarise, Temperature_c_month=mean(Temperature_c,na.rm=T)), c("Site","Type","Month")) Temperature_DOY$Dayof2019<-as.Date(as.character(Temperature_DOY$DOY),"%j") Temperature_DOY$Site<-mapvalues(Temperature_DOY$Site,from=c("SEGC","Saline"),to=c("Savanna","Forest")) ## - Seasonality plot (Figure 2) MaxMin_seasonal_plot_forest<-ggplot(Temperature_DOY[Temperature_DOY$Site=="Forest",])+ geom_line(aes(x=Dayof2019,y=Temperature_c,group=Type),alpha=0.7,colour="gray",size=0.4)+ geom_line(aes(x=Dayof2019,y=Temperature_c_10days,group=Type),alpha=0.7,colour="black",size=0.4)+ geom_line(aes(x=Dayof2019,y=Temperature_c_month,group=Type),alpha=1,size=0.7)+ scale_x_date(breaks=date_breaks("1 month"),labels = function(x) substring(format(as.Date(x),"%b"), 1, 1) )+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-03-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-06-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-10-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-12-01"))),linetype="dotted")+ theme_classic()+ ylab("Temperature (c)")+ xlab("Month")+ ylim(c(20,35))+ #ggtitle("B.")+ theme(plot.margin=unit(c(0.5,0.5,0.5,1),"cm"))+ geom_label(label="Forest",x=as.Date("2019-12-01"),y=34.5) MaxMin_seasonal_plot_savanna<-ggplot(Temperature_DOY[Temperature_DOY$Site=="Savanna",])+ geom_line(aes(x=Dayof2019,y=Temperature_c,group=Type),alpha=0.7,colour="gray",size=0.4)+ geom_line(aes(x=Dayof2019,y=Temperature_c_10days,group=Type),alpha=0.7,colour="black",size=0.4)+ geom_line(aes(x=Dayof2019,y=Temperature_c_month,group=Type),alpha=1,size=0.7)+ scale_x_date(breaks=date_breaks("1 month"),labels = function(x) substring(format(as.Date(x),"%b"), 1, 1) )+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-03-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-06-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-10-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-12-01"))),linetype="dotted")+ theme_classic()+ ylab("Temperature (c)")+ xlab("Month")+ ylim(c(20,35))+ #ggtitle("B.")+ theme(plot.margin=unit(c(0.5,0.5,0.5,1),"cm"))+ geom_label(label="Savanna",x=as.Date("2019-12-01"),y=34.5) ### Absolute Humidity ### Humidity_DOY<-ddply(AbsHumidity_daily_B,.(Site,DOY),summarise, AbsoluteHumidity_gm3=mean(AbsoluteHumidity_gm3,na.rm=T)) Humidity_DOY$Date<-as.Date(as.character(Humidity_DOY$DOY),"%j") Humidity_DOY$Month<-month(Humidity_DOY$Date) Humidity_DOY<-data.frame(Humidity_DOY,AbsoluteHumidity_gm3_10days=ddply(Humidity_DOY,.(Site),summarise, AbsoluteHumidity_gm3_10days=rollapply(AbsoluteHumidity_gm3,10,mean,partial=T))[,2]) Humidity_DOY<-merge(Humidity_DOY,ddply(Humidity_DOY,.(Site,Month),summarise, AbsoluteHumidity_gm3_month=mean(AbsoluteHumidity_gm3,na.rm=T)), c("Site","Month")) #ddply(Humidity_DOY, .(Site), summarise, # Min=min(AbsoluteHumidity_gm3_month), # Max=max(AbsoluteHumidity_gm3_month)) Humidity_DOY$Dayof2019<-as.Date(as.character(Humidity_DOY$DOY),"%j") ## - Seasonality plot (Figure 2) AH_seasonal_plot_forest<-ggplot(Humidity_DOY[Humidity_DOY$Site=="Forest",])+ geom_line(aes(x=Dayof2019,y=AbsoluteHumidity_gm3),alpha=0.7,colour="gray",size=0.4)+ geom_line(aes(x=Dayof2019,y=AbsoluteHumidity_gm3_10days),alpha=0.7,colour="black",size=0.4)+ geom_line(aes(x=Dayof2019,y=AbsoluteHumidity_gm3_month),alpha=1,size=0.7)+ scale_x_date(breaks=date_breaks("1 month"),labels = function(x) substring(format(as.Date(x),"%b"), 1, 1) )+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-03-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-06-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-10-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-12-01"))),linetype="dotted")+ theme_classic()+ ylim(c(17.5,23.5))+ ylab("Absolute Humidity (g/m3)")+ xlab("Month")+ #ggtitle("B.")+ theme(plot.margin=unit(c(0.5,0.5,0.5,1),"cm"))+ geom_label(label="Forest",x=as.Date("2019-12-01"),y=23) AH_seasonal_plot_savanna<-ggplot(Humidity_DOY[Humidity_DOY$Site=="Savanna",])+ geom_line(aes(x=Dayof2019,y=AbsoluteHumidity_gm3),alpha=0.7,colour="gray",size=0.4)+ geom_line(aes(x=Dayof2019,y=AbsoluteHumidity_gm3_10days),alpha=0.7,colour="black",size=0.4)+ geom_line(aes(x=Dayof2019,y=AbsoluteHumidity_gm3_month),alpha=1,size=0.7)+ scale_x_date(breaks=date_breaks("1 month"),labels = function(x) substring(format(as.Date(x),"%b"), 1, 1) )+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-03-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-06-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-10-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-12-01"))),linetype="dotted")+ theme_classic()+ ylim(c(17.5,23.5))+ ylab("Absolute Humidity (g/m3)")+ xlab("Month")+ #ggtitle("B.")+ theme(plot.margin=unit(c(0.5,0.5,0.5,1),"cm"))+ geom_label(label="Savanna",x=as.Date("2019-12-01"),y=23) ### Solar ### Solar_DOY<-ddply(Solar_daily_B,.(Month=month(Date),DOY),summarise, Solar_Wm2=mean(Solar_Wm2,na.rm=T)) Solar_DOY$Solar_Wm2_10days<-rollapply(Solar_DOY$Solar_Wm2,10, mean,partial=T) Solar_DOY<-merge(Solar_DOY,ddply(Solar_DOY,.(Month),summarise,Solar_Wm2_month=mean(Solar_Wm2,na.rm=T)),"Month") Solar_DOY$Dayof2019<-as.Date(as.character(Solar_DOY$DOY),"%j") ## - Seasonality plot (Figure 2) Solar_seasonal_plot<-ggplot(Solar_DOY,aes(x=Dayof2019,y=Solar))+ geom_line(aes(x=Dayof2019,y=Solar_Wm2),alpha=0.7,colour="gray",size=0.4)+ geom_line(aes(x=Dayof2019,y=Solar_Wm2_10days),alpha=0.7,colour="black",size=0.4)+ geom_line(aes(x=Dayof2019,y=Solar_Wm2_month),alpha=1,size=0.7)+ scale_x_date(breaks=date_breaks("1 month"),labels = function(x) substring(format(as.Date(x),"%b"), 1, 1) )+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-03-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-06-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-10-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-12-01"))),linetype="dotted")+ theme_classic()+ ylab("Solar Radiation (W/m2)")+ xlab("Month")+ # ggtitle("D.")+ theme(legend.title=element_blank(),legend.position="none")+ theme(plot.margin=unit(c(0.5,0.5,0.5,1),"cm")) ### Wind ### Wind_DOY<-ddply(Wind_daily_B,.(Month=month(Date),DOY),summarise, Wind_m_s=mean(Wind_m_s,na.rm=T)) Wind_DOY$Wind_m_s_10days<-rollapply(Wind_DOY$Wind_m_s,10, mean,partial=T) Wind_DOY<-merge(Wind_DOY,ddply(Wind_DOY,.(Month),summarise,Wind_m_s_month=mean(Wind_m_s,na.rm=T)),"Month") Wind_DOY$Dayof2019<-as.Date(as.character(Wind_DOY$DOY),"%j") ## - Seasonality plot (Figure 2) Wind_seasonal_plot<-ggplot(Wind_DOY,aes(x=Dayof2019,y=Solar))+ geom_line(aes(x=Dayof2019,y=Wind_m_s),alpha=0.7,colour="gray",size=0.4)+ geom_line(aes(x=Dayof2019,y=Wind_m_s_10days),alpha=0.7,colour="black",size=0.4)+ geom_line(aes(x=Dayof2019,y=Wind_m_s_month),alpha=1,size=0.7)+ scale_x_date(breaks=date_breaks("1 month"),labels = function(x) substring(format(as.Date(x),"%b"), 1, 1) )+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-03-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-06-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-10-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-12-01"))),linetype="dotted")+ theme_classic()+ ylab("Wind Speed (m/s)")+ xlab("Month")+ theme(legend.title=element_blank(),legend.position="none")+ theme(plot.margin=unit(c(0.5,0.5,0.5,1),"cm"))+ #ggtitle("E.")+ guides(linetype=guide_legend(nrow=2),alpha=guide_legend(nrow=2),colour=guide_legend(nrow=2)) ### Aerosol optical thickness ### Aerosol_DOY<-ddply(Aerosol_daily_B,.(Month,DOY),summarise, MeanAOT_440=mean(AOT_440,na.rm=T), MeanAOT_500=mean(AOT_500,na.rm=T), MeanAOT_675=mean(AOT_675,na.rm=T)) Aerosol_DOY$MeanAOT_440_10days<-rollapply(Aerosol_DOY$MeanAOT_440,10, mean,na.rm=T,partial=T) Aerosol_DOY$MeanAOT_500_10days<-rollapply(Aerosol_DOY$MeanAOT_500,10, mean,na.rm=T,partial=T) Aerosol_DOY$MeanAOT_675_10days<-rollapply(Aerosol_DOY$MeanAOT_675,10, mean,na.rm=T,partial=T) Aerosol_DOY<-merge(Aerosol_DOY,ddply(Aerosol_DOY,.(Month),summarise, MeanAOT_440_month=mean(MeanAOT_440,na.rm=T), MeanAOT_500_month=mean(MeanAOT_500,na.rm=T), MeanAOT_675_month=mean(MeanAOT_675,na.rm=T)),"Month") Aerosol_DOY$Dayof2019<-as.Date(as.character(Aerosol_DOY$DOY),"%j") ## - Seasonality plot (Figure 2) AOT_seasonal_plot<-ggplot(Aerosol_DOY)+ geom_line(aes(x=Dayof2019,y=MeanAOT_500),alpha=0.7,colour="gray",size=0.4)+ geom_line(aes(x=Dayof2019,y=MeanAOT_500_10days),alpha=0.7,colour="black",size=0.4)+ geom_line(aes(x=Dayof2019,y=MeanAOT_500_month),alpha=1,size=0.7)+ scale_x_date(breaks=date_breaks("1 month"),labels = function(x) substring(format(as.Date(x),"%b"), 1, 1) )+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-03-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-06-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-10-01"))),linetype="dotted")+ geom_vline(aes(xintercept=as.numeric(as.Date("2019-12-01"))),linetype="dotted")+ theme_classic()+ ylab("Aerosol Optical Depth (500nm)")+ xlab("Month")+ theme(legend.title=element_blank(),legend.position="none")+ theme(plot.margin=unit(c(0.5,0.5,0.5,1),"cm"))+ #ggtitle("E.")+ guides(linetype=guide_legend(nrow=2),alpha=guide_legend(nrow=2),colour=guide_legend(nrow=2)) # Supplementary Figure S2 #ggplot(Aerosol_DOY)+ # geom_line(aes(x=Dayof2019,y=MeanAOT_500_month,linetype="AOT_500"),alpha=1,size=0.7)+ # geom_line(aes(x=Dayof2019,y=MeanAOT_440_month,linetype="AOT_440"),alpha=1,size=0.7)+ # geom_line(aes(x=Dayof2019,y=MeanAOT_675_month,linetype="AOT_675"),alpha=1,size=0.7)+ # scale_x_date(labels=date_format("%b"),breaks=date_breaks("1 month"))+ # geom_vline(aes(xintercept=as.numeric(as.Date("2019-03-01"))),linetype="dotted")+ # geom_vline(aes(xintercept=as.numeric(as.Date("2019-06-01"))),linetype="dotted")+ # geom_vline(aes(xintercept=as.numeric(as.Date("2019-10-01"))),linetype="dotted")+ # geom_vline(aes(xintercept=as.numeric(as.Date("2019-12-01"))),linetype="dotted")+ # theme_classic()+ # ylab("Aerosol Optical Depth")+ # theme(legend.title=element_blank(),legend.position="right")+ # xlab("Month")+ # theme(plot.margin=unit(c(0.5,0.5,0.5,1),"cm"))+ # guides(linetype=guide_legend(nrow=3),alpha=guide_legend(nrow=3),colour=guide_legend(nrow=2)) ## - Arrange all seasonal plots together - Figure 2 #ggarrange(Rainfall_seasonal_plot,Wind_seasonal_plot, # MaxMin_seasonal_plot_savanna,MaxMin_seasonal_plot_forest, # AH_seasonal_plot_savanna,AH_seasonal_plot_forest, # Solar_seasonal_plot,AOT_seasonal_plot, # ncol=2,nrow=4,align="hv", # labels=c("A","B","C","D","E","F","G","H")) ############################## ##### 3. Fourier spectra ##### ############################## ### Rainfall ### Rainfall_spec<-spectrum(Rainfall_monthly_E$Rainfall_standard_mm,demean=T,detrend=T) Rainfall_spec_df<-data.frame(Freq=Rainfall_spec$freq,Spec=Rainfall_spec$spec,Data="Rainfall",Type="",Site="Savanna") ### Temperature ### MaxTemp_Saline_spec<-spectrum(Temperature_monthly_E$Temperature_standard_c[Temperature_monthly_E$Site=="Saline"&Temperature_monthly_E$Type=="Max"],demean=T,detrend=T) MaxTemp_Saline_spec_df<-data.frame(Freq=MaxTemp_Saline_spec$freq,Spec=MaxTemp_Saline_spec$spec,Data="Temp (max)", Type="Max",Site="Forest") MaxTemp_SEGC_spec<-spectrum(Temperature_monthly_E$Temperature_standard_c[Temperature_monthly_E$Site=="SEGC"&Temperature_monthly_E$Type=="Max"],demean=T,detrend=T) MaxTemp_SEGC_spec_df<-data.frame(Freq=MaxTemp_SEGC_spec$freq,Spec=MaxTemp_SEGC_spec$spec,Data="Temp (max)", Type="Max",Site="Savanna") MinTemp_Saline_spec<-spectrum(Temperature_monthly_E$Temperature_standard_c[Temperature_monthly_E$Site=="Saline"&Temperature_monthly_E$Type=="Min"],demean=T,detrend=T) MinTemp_Saline_spec_df<-data.frame(Freq=MinTemp_Saline_spec$freq,Spec=MinTemp_Saline_spec$spec,Data="Temp (min)", Type="Min",Site="Forest") MinTemp_SEGC_spec<-spectrum(Temperature_monthly_E$Temperature_standard_c[Temperature_monthly_E$Site=="SEGC"&Temperature_monthly_E$Type=="Min"],demean=T,detrend=T) MinTemp_SEGC_spec_df<-data.frame(Freq=MinTemp_SEGC_spec$freq,Spec=MinTemp_SEGC_spec$spec,Data="Temp (min)", Type="Min",Site="Savanna") ### Humidity ### AbsHumidity_Savanna_spec<-spectrum(AbsHumidity_monthly_C$AbsHumidity_standard_gm3[AbsHumidity_monthly_C$Site=="Savanna"],demean=T,detrend=T) AbsHumidity_Savanna_spec_df<-data.frame(Freq=AbsHumidity_Savanna_spec$freq,Spec=AbsHumidity_Savanna_spec$spec,Data="Humidity", Type="",Site="Savanna") AbsHumidity_Forest_spec<-spectrum(AbsHumidity_monthly_C$AbsHumidity_standard_gm3[AbsHumidity_monthly_C$Site=="Forest"],demean=T,detrend=T) AbsHumidity_Forest_spec_df<-data.frame(Freq=AbsHumidity_Forest_spec$freq,Spec=AbsHumidity_Forest_spec$spec,Data="Humidity", Type="",Site="Forest") ### Solar ### Solar_spec<-spectrum(Solar_monthly_C$Solar_standard_Wm2,demean=T,detrend=T) Solar_spec_df<-data.frame(Freq=Solar_spec$freq,Spec=Solar_spec$spec,Data="Solar radiation", Type="",Site="Savanna") ### Wind ### Wind_spec<-spectrum(Wind_monthly_C$Wind_standard_m_s,demean=T,detrend=T) Wind_spec_df<-data.frame(Freq=Wind_spec$freq,Spec=Wind_spec$spec,Data="Wind spped", Type="",Site="Savanna") ### AOT ### AOT_spec<-spectrum(Aerosol_monthly_C$AOT_500_standard,demean=T,detrend=T) AOT_spec_df<-data.frame(Freq=AOT_spec$freq,Spec=AOT_spec$spec,Data="AOD", Type="",Site="Savanna") ## - Combine all data Spec_combined<-rbind(Rainfall_spec_df, MaxTemp_Saline_spec_df, MaxTemp_SEGC_spec_df, MinTemp_Saline_spec_df, MinTemp_SEGC_spec_df, AbsHumidity_Savanna_spec_df, AbsHumidity_Forest_spec_df, Solar_spec_df, Wind_spec_df, AOT_spec_df) ## - Make combined plot (Supplementary Figure S1) ggplot(Spec_combined[Spec_combined$Freq<0.5,],aes(x=Freq,y=Spec))+ geom_hline(yintercept=-Inf)+ geom_vline(xintercept=-Inf)+ geom_line()+ #scale_linetype_manual(values = c("dotted","solid","dashed"))+ #scale_colour_manual(values=c("gray","black","darkgray"))+ #scale_size_manual(values=c(0.8,1,0.5))+ theme_classic()+ facet_grid(Data~Site)+ geom_vline(aes(xintercept=1/12),linetype="dotted")+ geom_vline(aes(xintercept=1/6),linetype="dotted")+ theme(plot.margin=unit(c(0.5,0.5,0.5,1),"cm"))+ ylab("Spectrum")+ xlab("Frequency (1/months)") ############################### ##### 4. Long-term trends ##### ############################### ### Rainfall ### ## - Test for linear trend in total annual rainfall over time # - Run glm Rain_annual_m1<-glm(round(Rainfall_yearly_mm)~Year_r,data=Rainfall_yearly_F,family=poisson) Rain_annual_m2<-glm(round(Rainfall_yearly_mm)~1,data=Rainfall_yearly_F,family=poisson) # - Table 2 Rainfall_trend_AIC<-data.frame(Response="Rainfall", Model=c("Long-term change", "No long-term change"), Predictors=c("Year","Intercept only"), AIC=round(AIC(Rain_annual_m1,Rain_annual_m2)$AIC,digits=1), DF= AIC(Rain_annual_m1,Rain_annual_m2)$df) # - Predict trend Rainfall_yearly_F$Rainfall_yearly_glm_mm<-predict(Rain_annual_m1,Rainfall_yearly_F,type="response") # - Change in total annual rainfall (mm) per decade #(Rainfall_yearly_F$Rainfall_yearly_F_glm[Rainfall_yearly_F$Year==2017]-Rainfall_yearly_F$Rainfall_yearly_F_glm[Rainfall_yearly_F$Year==1984])/(2017-1984)*10 #Percentage change compared to mean total annual rainfall #(Rainfall_yearly_F$Rainfall_yearly_F_glm[Rainfall_yearly_F$Year==2017]-Rainfall_yearly_F$Rainfall_yearly_F_glm[Rainfall_yearly_F$Year==1984])/(2017-1984)*10/mean(Rainfall_yearly_F$Rainfall_yearly_F_mm,na.rm=T) # - Plot trend and raw data (Figure 3A) Rainfall_annual_glm_plot<-ggplot(Rainfall_yearly_F)+ geom_point(aes(x=Year,y=Rainfall_yearly_mm),colour="gray")+ geom_line(aes(x=Year,y=Rainfall_yearly_mm),colour="gray")+ geom_line(aes(x=Year,y=Rainfall_yearly_glm_mm))+ scale_x_continuous(breaks=c(1985,1990,1995,2000,2005,2010,2015))+ theme_classic()+ ylab("Rainfall (mm/year)")+ xlab("Year")+ theme(plot.margin=unit(c(0.5,0.5,0.5,0.5),"cm")) ## - Test for linear trend interacting with season # - GLM on daily data including interaction with season Rain_seasonal_m1<-glmer(round(Rainfall_mm)~Year_r*Season+(1|DOY)+(1|Year),data=Rainfall_daily_D,family=poisson) # - Check for autocorrellation #Rainfall_daily_D$Residuals<-NA #Rainfall_daily_D$Residuals[complete.cases(Rainfall_daily_D$Rainfall_mm)]<-resid(SeasonalTrend_m1) #acf_plot(Rainfall_daily_mean$Residuals,split_by=list(Rainfall_daily_mean$DOY),fun=median) #acf_n_plots(Rainfall_daily_mean$Residuals,split_by=list(Rainfall_daily_mean$DOY)) Rain_seasonal_m2<-glmer(round(Rainfall_mm)~Year_r+Season+(1|DOY)+(1|Year),data=Rainfall_daily_D,family=poisson) # - Table 3 Rainfall_seasonal_trend_AIC<-data.frame(Response="Rainfall", Model=c("Long-term change varying by season", "Long-term change not varying by season"), Predictors=c("Year * Season","Year + Season"), AIC=round(AIC(Rain_seasonal_m1,Rain_seasonal_m2)$AIC,digits=1), DF= AIC(Rain_seasonal_m1,Rain_seasonal_m2)$df) # - Remove intercept to directly compare estimates Rain_seasonal_m1_mod<-glmer(round(Rainfall_mm)~-1+Season+Year_r:Season+(1|DOY)+(1|Year),data=Rainfall_daily_D,family=poisson) # - Extract estimates Rainfall_seasonal_trend_glm<-tidy(Rain_seasonal_m1_mod,conf.int=T) Rainfall_seasonal_trend_glm$sig<-ifelse(Rainfall_seasonal_trend_glm$conf.low>0|Rainfall_seasonal_trend_glm$conf.high<0,"*","") Rainfall_seasonal_trend_glm$conf.low<-round(Rainfall_seasonal_trend_glm$conf.low,digits=2) Rainfall_seasonal_trend_glm$conf.high<-round(Rainfall_seasonal_trend_glm$conf.high,digits=2) Rainfall_seasonal_trend_glm<-unite(Rainfall_seasonal_trend_glm[Rainfall_seasonal_trend_glm$group=="fixed",],"95%CI",c("conf.low","conf.high"),sep=",") Rainfall_seasonal_trend_glm<-as.data.frame(Rainfall_seasonal_trend_glm) #Table 4 Rainfall_seasonal_trend_glm<-data.frame(Response = "Rainfall", Predictor = c("DJF","JJAS","MAM","ON","Year: DJF","Year: JJAS","Year: MAM","Year: ON"), Estimate= round(Rainfall_seasonal_trend_glm[Rainfall_seasonal_trend_glm$group=="fixed","estimate"],digits=2), SE = round(Rainfall_seasonal_trend_glm[Rainfall_seasonal_trend_glm$group=="fixed","std.error"],digits=2), "T Value" = round(Rainfall_seasonal_trend_glm[Rainfall_seasonal_trend_glm$group=="fixed","statistic"],digits=2), "95% CI" =Rainfall_seasonal_trend_glm[Rainfall_seasonal_trend_glm$group=="fixed","95%CI"], "."=Rainfall_seasonal_trend_glm[Rainfall_seasonal_trend_glm$group=="fixed","sig"]) # - Predict from the GLM #Rainfall_daily_D$Lopé_glm<-predict(Rain_seasonal_m1,Rainfall_daily_D,re.form=NA,type="response") #ddply(Rainfall_daily_D,.(Season),summarise, # Rainfall_1984=mean(Lopé_glm[Year==1984],na.rm=T), # Rainfall_2017=mean(Lopé_glm[Year==2017],na.rm=T), # Years=2017-1984, # Difference_mm=Rainfall_2017-Rainfall_1984, # Difference_mm_yrs=Difference_mm/Years, # Difference_mm_dec=Difference_mm_yrs*10, # Rainfall_mean_mm=mean(Rainfall_mm,na.rm=T), # Difference_per_dec=(Difference_mm_dec)/Rainfall_mean_mm) ### Temperature ### ## - Test for linear trend in minimum daily temperature # - Run LMM to test for linear trend over time Temp_annual_m1<-lmer(Temperature_c~Year_r+(1|DOY)+(1|Year),data=Temperature_daily_F,na.action = na.exclude) Temp_annual_m2<-lmer(Temperature_c~1+(1|DOY)+(1|Year),data=Temperature_daily_F,na.action = na.exclude) AIC(Temp_annual_m1,Temp_annual_m2) # - Check for temporal autocorrelation #Temperature_daily_F$Residuals<-NA #Temperature_daily_F$Residuals<-resid(Temp_annual_m1) #acf_plot(Temperature_daily_F$Residuals,split_by=list(Temperature_daily_F$DOY),fun=median) #acf_n_plots(Temperature_daily_F$Residuals,split_by=list(Temperature_daily_F$DOY)) # - Table 2 Temperature_trend_AIC<-data.frame(Response="Temperature", Model=c("Long-term change", "No long-term change"), Predictors=c("Year","Intercept only"), AIC=round(AIC(Temp_annual_m1,Temp_annual_m2)$AIC,digits=1), DF= AIC(Temp_annual_m1,Temp_annual_m2)$df) # - Predict from LMM Temperature_daily_F$LMM_predict<-predict(Temp_annual_m1,Temperature_daily_F,re.form=NA) Temperature_annual_LMM<-ddply(Temperature_daily_F,.(Year),summarise, Temperature_interp_mean_c=mean(Temperature_interp_c,na.rm=T), Temperature_mean_c=mean(Temperature_c,na.rm=T), Temperature_interp_mean_sample=length(Temperature_interp_c[complete.cases(Temperature_interp_c)]), Temperature_mean_sample=length(Temperature_c[complete.cases(Temperature_c)]), Temperature_LMM_c=mean(LMM_predict,na.rm=T)) Temperature_annual_LMM<-Temperature_annual_LMM[complete.cases(Temperature_annual_LMM$Year),] Temperature_annual_LMM$Temperature_mean_c[Temperature_annual_LMM$Temperature_mean_sample<5]<-NA #(Temperature_annual_LMM$Temperature_c[Temperature_annual_LMM$Year==2017]-Temperature_annual_LMM$Temperature_c[Temperature_annual_LMM$Year==1984])/(2017-1984)*10 #0.23 per decade #(Temperature_annual_LMM$Temperature_c[Temperature_annual_LMM$Year==2017]-Temperature_annual_LMM$Temperature_c[Temperature_annual_LMM$Year==1984])/(2017-1984)*10/mean(Temperature_daily_F$Temperature_c,na.rm=T) #1.1% # - Figure 3B MinTemp_annual_glm_plot<-ggplot()+ geom_point(data=Temperature_annual_LMM,aes(x=YearMonth,y=Temperature_mean_c),colour="gray",alpha=0.4)+ # geom_line(data=Temperature_ym_LMM,aes(x=YearMonth,y=Temperature_mean_c),colour="gray")+ geom_line(data=Temperature_annual_LMM[month(Temperature_ym_LMM$YearMonth)==6,],aes(x=YearMonth,y=Temperature_LMM_c))+ scale_x_continuous(breaks=c(1985,1990,1995,2000,2005,2010,2015))+ theme_classic()+ ylab("Min temperature (c)")+ xlab("Year")+ theme(plot.margin=unit(c(0.5,0.5,0.5,0.5),"cm")) ## - Test for linear trend interacting with season # - GLM on daily data including interaction with season Temp_seasonal_m1<-lmer(Temperature_c~Year_r*Season+(1|DOY)+(1|Year),data=Temperature_daily_F,na.action = na.exclude) Temp_seasonal_m2<-lmer(Temperature_c~Year_r+Season+(1|DOY)+(1|Year),data=Temperature_daily_F,na.action = na.exclude) # - Check for autocorrellation #Temperature_daily_F$Residuals<-NA #Temperature_daily_F$Residuals<-resid(Temp_seasonal_m1) #acf_plot(Temperature_daily_F$Residuals,split_by=list(Temperature_daily_F$DOY),fun=median) #acf_n_plots(Temperature_daily_F$Residuals,split_by=list(Temperature_daily_F$DOY)) # - Table 3 Temperature_seasonal_trend_AIC<-data.frame(Response="Temperature", Model=c("Long-term change varying by season", "Long-term change not varying by season"), Predictors=c("Year * Season","Year + Season"), AIC=round(AIC(Temp_seasonal_m1,Temp_seasonal_m2)$AIC,digits=1), DF= AIC(Temp_seasonal_m1,Temp_seasonal_m2)$df) Temp_seasonal_m1_mod<-lmer(Temperature_c~-1+Season+Year_r:Season+(1|DOY)+(1|Year),data=Temperature_daily_F,na.action = na.exclude) Temperature_seasonal_trend_glm<-tidy(Temp_seasonal_m1_mod,conf.int=T) Temperature_seasonal_trend_glm$sig<-ifelse(Temperature_seasonal_trend_glm$conf.low>0|Temperature_seasonal_trend_glm$conf.high<0,"*","") Temperature_seasonal_trend_glm$conf.low<-round(Temperature_seasonal_trend_glm$conf.low,digits=2) Temperature_seasonal_trend_glm$conf.high<-round(Temperature_seasonal_trend_glm$conf.high,digits=2) Temperature_seasonal_trend_glm<-unite(Temperature_seasonal_trend_glm[Temperature_seasonal_trend_glm$group=="fixed",],"95%CI",c("conf.low","conf.high"),sep=",") Temperature_seasonal_trend_glm<-as.data.frame(Temperature_seasonal_trend_glm) # - Table 4 Temperature_seasonal_trend_glm<-data.frame(Response = "Temperature", Predictor = c("DJF","JJAS","MAM","ON","Year: DJF","Year: JJAS","Year: MAM","Year: ON"), Estimate= round(Temperature_seasonal_trend_glm[Temperature_seasonal_trend_glm$group=="fixed","estimate"],digits=2), SE = round(Temperature_seasonal_trend_glm[Temperature_seasonal_trend_glm$group=="fixed","std.error"],digits=2), "T Value" = round(Temperature_seasonal_trend_glm[Temperature_seasonal_trend_glm$group=="fixed","statistic"],digits=2), "95% CI" =Temperature_seasonal_trend_glm[Temperature_seasonal_trend_glm$group=="fixed","95%CI"], "."=Temperature_seasonal_trend_glm[Temperature_seasonal_trend_glm$group=="fixed","sig"]) # - Predict from the LMM #Temperature_daily_F$Lopé_glm<-predict(Temp_seasonal_m1,Temperature_daily_F,re.form=NA,type="response") #ddply(Temperature_daily_F,.(Season),summarise, # Temp_1984=mean(Lopé_glm[Year==1984],na.rm=T), # Temp_2017=mean(Lopé_glm[Year==2017],na.rm=T), # Years=2017-1984, # Difference_c=Temp_2017-Temp_1984, # Difference_c_yrs=Difference_c/Years, # Difference_c_dec=Difference_c_yrs*10, # Temp_mean_c=mean(Temperature_c,na.rm=T), # Difference_per_dec=(Difference_c_dec)/Temp_mean_c) ### - Gridded temperature data for Lopé ### ## - Berkeley Berkeley_annual_m1<-lmer(Temperature_c~Year_r+(1|DOY)+(1|Year),data=Berkeley_daily,na.action = na.exclude) summary(Berkeley_annual_m1) # - Predict from GLM #Berkeley_daily$GLM_predict<-predict(Berkeley_annual_m1,Berkeley_daily,re.form=NA,allow.new.levels=TRUE) #Berkeley_yearly<-ddply(Berkeley_daily,.(Year),summarise, # Temperature_mean_c=mean(Temperature_c,na.rm=T),Temperature_GLM_c=mean(GLM_predict,na.rm=T)) #(Berkeley_yearly$Temperature_GLM_c[Berkeley_yearly$Year==2017]-Berkeley_yearly$Temperature_GLM_c[Berkeley_yearly$Year==1984])/(2017-1984)*10 # - 0.1607 per decade #((Berkeley_yearly$Temperature_GLM_c[Berkeley_yearly$Year==2017]-Berkeley_yearly$Temperature_GLM_c[Berkeley_yearly$Year==1984])/(2017-1984)*10)/mean(Berkeley_yearly$Temperature_mean_c,na.rm=T) # - 0.0077 ## - CRU CRU_annual_m1<-lmer(Temperature_c~Year_r+(1|Month)+(1|Year),data=CRU_monthly,na.action = na.exclude) summary(CRU_annual_m1) # - Predict from GLM #CRU_monthly$GLM_predict<-predict(CRU_annual_m1,CRU_monthly,re.form=NA,allow.new.levels=TRUE) #CRU_yearly<-ddply(CRU_monthly,.(Year),summarise,Temperature_mean_c=mean(Temperature_c,na.rm=T),Temperature_GLM_c=mean(GLM_predict,na.rm=T)) #(CRU_yearly$Temperature_GLM_c[CRU_yearly$Year==2016]-CRU_yearly$Temperature_GLM_c[CRU_yearly$Year==1984])/(2016-1984)*10 # - 0.188 per decade #((CRU_yearly$Temperature_GLM_c[CRU_yearly$Year==2016]-CRU_yearly$Temperature_GLM_c[CRU_yearly$Year==1984])/(2016-1984)*10)/mean(CRU_yearly$Temperature_mean_c,na.rm=T) # - 0.0085 #################################### ##### 5. Periodicity over time ##### #################################### ### Rainfall ### ## - Compute wavelet transform on complete (interpolated) monthly time series for rainfall Rainfall_wavelet<-wt(Rainfall_monthly_E[,c(2,9)],dj=1/12) #plot(Rainfall_wavelet,xlab="",ylab="Rainfall period (yrs)") #Figure 3C # - Extract wavelet trasnform for three components (6-month, 12-month and 2-4 years) Rainfall_wavelet_extract<-data.frame(YearMonth=Rainfall_wavelet$t, COI=Rainfall_wavelet$coi, Power_6months=Rainfall_wavelet$power.corr[which.min(abs(Rainfall_wavelet$period-0.5)),], Power_12months=Rainfall_wavelet$power.corr[which.min(abs(Rainfall_wavelet$period-1)),], Power_2_4years=colMeans(Rainfall_wavelet$power.corr[which.min(abs(Rainfall_wavelet$period-2)):which.min(abs(Rainfall_wavelet$period-4)),])) # - Remove data within the cone of influence (overly influenced by edge effects) Rainfall_wavelet_extract$Power_6months[Rainfall_wavelet_extract$COI<0.5]<-NA Rainfall_wavelet_extract$Power_12months[Rainfall_wavelet_extract$COI<1]<-NA Rainfall_wavelet_extract$Power_2_4years[Rainfall_wavelet_extract$COI<4]<-NA Rainfall_wavelet_extract$Date<-as.Date(paste(year(Rainfall_wavelet_extract$YearMonth),month(Rainfall_wavelet_extract$YearMonth),"15",sep="-")) #mean(Rainfall_wavelet_extract$Power_6months,na.rm=T)/mean(Rainfall_wavelet_extract$Power_12months,na.rm=T) #mean(Rainfall_wavelet_extract$Power_6months,na.rm=T)/mean(Rainfall_wavelet_extract$Power_2_4years,na.rm=T) # - Plot extracted wavelet components (Figure 3E) Rainfall_extract_plot<-ggplot(Rainfall_wavelet_extract)+ geom_line(aes(x=Date,y=Power_6months*0.0001,linetype="Biannual"))+ geom_line(aes(x=Date,y=Power_12months*0.0001,linetype="Annual"))+ geom_line(aes(x=Date,y=Power_2_4years*0.0001,linetype="2-4 years"))+ scale_linetype_manual(values=c("dotted","solid","longdash"))+ scale_x_date(breaks=c(as.Date("1985-01-01"),as.Date("1990-01-01"),as.Date("1995-01-01"),as.Date("2000-01-01"),as.Date("2005-01-01"),as.Date("2010-01-01"),as.Date("2015-01-01")),date_label="%Y")+ theme_classic()+ ylab("Wavelet power (/10000)")+ xlab("Year")+ theme(plot.margin=unit(c(0.5,0.5,0.5,0.5),"cm"), legend.title = element_blank(),legend.position=c(0.9,0.9)) ### Temperature ### ## - Compute wavelet transform on complete (interpolated) monthly time series for mean minimum daily temperature Temperature_wavelet<-wt(Temperature_monthly_G[,c(3,8)],dj=1/12) #plot(Temperature_wavelet,xlab="",ylab="Temperature period (yrs)") # Figure 3D # - Extract wavelet trasnform for three components (6-month, 12-month and 2-4 years) Temperature_wavelet_extract<-data.frame(YearMonth=Temperature_wavelet$t, COI=Temperature_wavelet$coi, Power_6months=Temperature_wavelet$power.corr[which.min(abs(Temperature_wavelet$period-0.5)),], Power_12months=Temperature_wavelet$power.corr[which.min(abs(Temperature_wavelet$period-1)),], Power_2_4years=colMeans(Temperature_wavelet$power.corr[which.min(abs(Temperature_wavelet$period-2)):which.min(abs(Temperature_wavelet$period-4)),])) # - Remove data within the cone of influence (overly influenced by edge effects) Temperature_wavelet_extract$Power_6months[Temperature_wavelet_extract$COI<0.5]<-NA Temperature_wavelet_extract$Power_12months[Temperature_wavelet_extract$COI<1]<-NA Temperature_wavelet_extract$Power_2_4years[Temperature_wavelet_extract$COI<4]<-NA Temperature_wavelet_extract$Date<-as.Date(paste(year(Temperature_wavelet_extract$YearMonth),month(Temperature_wavelet_extract$YearMonth),"15",sep="-")) #mean(Temperature_wavelet_extract$Power_12months,na.rm=T)/mean(Temperature_wavelet_extract$Power_6months,na.rm=T) #mean(Temperature_wavelet_extract$Power_12months,na.rm=T)/mean(Temperature_wavelet_extract$Power_2_4years,na.rm=T) # - Plot extracted wavelet components (Figure 3F) Temperature_extract_plot<-ggplot(Temperature_wavelet_extract)+ geom_line(aes(x=Date,y=Power_6months,linetype="Biannual"))+ geom_line(aes(x=Date,y=Power_12months,linetype="Annual"))+ geom_line(aes(x=Date,y=Power_2_4years,linetype="2-4 years"))+ scale_linetype_manual(values=c("dotted","solid","longdash"))+ scale_x_date(breaks=c(as.Date("1985-01-01"),as.Date("1990-01-01"),as.Date("1995-01-01"),as.Date("2000-01-01"),as.Date("2005-01-01"),as.Date("2010-01-01"),as.Date("2015-01-01")),date_label="%Y")+ theme_classic()+ ylim(c(0,100))+ ylab("Wavelet power")+ xlab("Year")+ theme(plot.margin=unit(c(0.5,0.5,0.5,0.5),"cm"), legend.title = element_blank(),legend.position=c(0.9,0.9)) ################################### ###### 6. Influence of oceans ##### ################################### ### Rainfall ### ## - Merge rainfall and teleconnection datasets Rainfall_monthly_E<-merge(Rainfall_monthly_E,Tele[,c(3,4,5,6,8,10,13)],c("Year","Month")) Rainfall_monthly_E$Year_rescaled<-scale(Rainfall_monthly_E$Year) # - Scale teleconnection data Rainfall_monthly_E$MEI_original<-Rainfall_monthly_E$MEI Rainfall_monthly_E$MEI<-as.numeric(scale(Rainfall_monthly_E$MEI)) Rainfall_monthly_E$NATL<-scale(Rainfall_monthly_E$NATL_anom) Rainfall_monthly_E$SATL<-scale(Rainfall_monthly_E$SATL_anom) Rainfall_monthly_E$IOD_original<-Rainfall_monthly_E$IOD Rainfall_monthly_E$IOD<-as.numeric(scale(Rainfall_monthly_E$IOD)) Rainfall_monthly_E<-Rainfall_monthly_E[order(Rainfall_monthly_E$YearMonth),] ## - Model the effects of ocean temperatures on rainfall at Lopé by season # - Explore collinearity #Rainfall2<-Rainfall_monthly_E[complete.cases(Rainfall_monthly_E$MEI)&complete.cases(Rainfall_monthly_E$Rainfall_month_mm),c()] #names(Rainfall2)[2]<-"Rainfall" #names(Rainfall2)[1]<-"Year" #corrplot(cor(Rainfall2),type="upper") # - GLM competition removing terms sequentially starting with full model Rainfall_Tele_m1_full<-glmer(round(Rainfall_mm)~MEI*Season+NATL*Season+SATL*Season+IOD*Season+Year_rescaled*Season+(1|Month)+(1|Year),data=Rainfall_monthly_E[complete.cases(Rainfall_monthly_E$MEI),],family="poisson",na.action=na.exclude) #Check for autocorrelation #Rainfall_monthly_E$OceansGLMResiduals<-NA #Rainfall_monthly_E$OceansGLMResiduals<-resid(Rainfall_Tele_m1_full) #acf_plot(Rainfall_monthly_E$OceansGLMResiduals,split_by=list(Rainfall_monthly_E$DOY),fun=median) #acf_n_plots(Rainfall_monthly_E$OceansGLMResiduals,split_by=list(Rainfall_monthly_E$DOY)) #acf_plot(Rainfall_monthly_E$OceansGLMResiduals,split_by=list(Rainfall_monthly_E$Year),fun=median) #acf_n_plots(Rainfall_monthly_E$OceansGLMResiduals,split_by=list(Rainfall_monthly_E$Year)) Rainfall_Tele_m2_IODonly<-glmer(round(Rainfall_mm)~MEI*Season+NATL*Season+SATL*Season+IOD+Year_rescaled*Season+(1|Month)+(1|Year),data=Rainfall_monthly_E[complete.cases(Rainfall_monthly_E$MEI),],family="poisson") Rainfall_Tele_m3_MEIonly<-glmer(round(Rainfall_mm)~MEI+NATL*Season+SATL*Season+IOD*Season+Year_rescaled*Season+(1|Month)+(1|Year),data=Rainfall_monthly_E[complete.cases(Rainfall_monthly_E$MEI),],family="poisson") Rainfall_Tele_m4_SATLonly<-glmer(round(Rainfall_mm)~MEI*Season+NATL*Season+SATL+IOD*Season+Year_rescaled*Season+(1|Month)+(1|Year),data=Rainfall_monthly_E[complete.cases(Rainfall_monthly_E$MEI),],family="poisson") Rainfall_Tele_m5_NATLonly<-glmer(round(Rainfall_mm)~MEI*Season+NATL+SATL*Season+IOD*Season+Year_rescaled*Season+(1|Month)+(1|Year),data=Rainfall_monthly_E[complete.cases(Rainfall_monthly_E$MEI),],family="poisson") # - Compare AIC values for all models (Table 5) Rainfall_oceans_AIC<-data.frame(Response="Rainfall", Predictors=c("Season + NATL:Season + SATL:Season + MEI:Season + IOD:Season + Year:Season", "Season + NATL:Season + SATL:Season + MEI:Season + IOD + Year:Season", "Season + NATL:Season + SATL:Season + MEI + IOD:Season + Year:Season", "Season + NATL:Season + SATL + MEI:Season + IOD:Season + Year:Season", "Season + NATL + SATL:Season + MEI:Season + IOD:Season + Year:Season"), AIC=round(AIC(Rainfall_Tele_m1_full,Rainfall_Tele_m2_IODonly,Rainfall_Tele_m3_MEIonly,Rainfall_Tele_m4_SATLonly,Rainfall_Tele_m5_NATLonly)$AIC,digits=2), DF=AIC(Rainfall_Tele_m1_full,Rainfall_Tele_m2_IODonly,Rainfall_Tele_m3_MEIonly,Rainfall_Tele_m4_SATLonly,Rainfall_Tele_m5_NATLonly)$df) # - Modify best model by removing intercept to allow direct comparisong of estimates Rainfall_Tele_m1_full_mod<-glmer(round(Rainfall_mm)~-1+Season+MEI:Season+NATL:Season+SATL:Season+IOD:Season+Year_rescaled:Season+(1|Month)+(1|Year),data=Rainfall_monthly_E[complete.cases(Rainfall_monthly_E$MEI),],family="poisson") Rainfall_oceans_glm<-tidy(Rainfall_Tele_m1_full_mod,conf.int=T) Rainfall_oceans_glm$sig<-ifelse(Rainfall_oceans_glm$conf.low>0|Rainfall_oceans_glm$conf.high<0,"*","") Rainfall_oceans_glm$conf.low<-round(Rainfall_oceans_glm$conf.low,digits=2) Rainfall_oceans_glm$conf.high<-round(Rainfall_oceans_glm$conf.high,digits=2) Rainfall_oceans_glm<-as.data.frame(Rainfall_oceans_glm) Rainfall_oceans_glm<-Rainfall_oceans_glm[Rainfall_oceans_glm$group=="fixed",] Rainfall_oceans_glm<-unite(Rainfall_oceans_glm,"95%CI",c("conf.low","conf.high"),sep=",") # - Table 6 Rainfall_oceans_glm<-data.frame(Predictor=c("DJF","JJAS","MAM","ON","MEI: DJF","MEI: JJAS","MEI: MAM","MEI: ON","NATL: DJF","NATL: JJAS","NATL: MAM","NATL: ON","SATL: DJF","SATL: JJAS","SATL: MAM","SATL: ON","IOD: DJF","IOD: JJAS","IOD: MAM","IOD: ON","Year: DJF","Year: JJAS","Year: MAM","Year: ON"), Estimate= round(Rainfall_oceans_glm[,"estimate"],digits=2), SE = round(Rainfall_oceans_glm[,"std.error"],digits=2), "Z Value" = round(Rainfall_oceans_glm[,"statistic"],digits=2), "95% CI" =Rainfall_oceans_glm[,"95%CI"], "."=Rainfall_oceans_glm[,"sig"]) ## - Predict from the best GLM and plot effects of NATL # - Construct new dataset Rainfall_NATL_glm<-expand.grid(MEI=0,Season=levels(as.factor(Rainfall_monthly_E$Season)),NATL_anom=unique(Rainfall_monthly_E$NATL_anom),SATL=0,IOD=0,Year_rescaled=0) Rainfall_NATL_glm$NATL<-scale(Rainfall_NATL_glm$NATL_anom) ] # - Predict from the best GLM Rainfall_NATL_glm$Rainfall_NATL_glm<-predict(Rainfall_Tele_m1_full_mod,newdata=Rainfall_NATL_glm,type="response",re.form=NA) # - Summarise rainfall data by season (to add to plot below) Rainfall_seasonal<-ddply(Rainfall_monthly_E,.(Year,Season),summarise, NATL_anom=mean(NATL_anom,na.rm=T), SATL_anom=mean(SATL_anom,na.rm=T), MEI_original=mean(MEI_original,na.rm=T), IOD_original=mean(IOD_original,na.rm=T), Rainfall_sample=length(Rainfall_mm[complete.cases(Rainfall_mm)]), Rainfall_mm=mean(Rainfall_mm,na.rm=T)) # - Plot summarised raw data and model prediction (Figure 4) Rainfall_NATL_seasonal_plot<-ggplot()+ geom_point(data=Rainfall_seasonal[Rainfall_seasonal$Season %in% c("DJF","ON", "JJAS"),],aes(x=NATL_anom,y=Rainfall_mm,shape=Season),alpha=0.6)+ geom_line(data=Rainfall_NATL_glm[Rainfall_NATL_glm$Season %in% c("DJF","ON", "JJAS"),],aes(x=NATL_anom,y=Rainfall_NATL_glm,linetype=Season))+ scale_linetype_manual(values=Linetype_lookup)+ scale_shape_manual(values=Shape_lookup)+ theme_classic()+ ylab("Rainfall (mm/day)")+ xlab("North Atlantic SST anomaly")+ ggtitle("D.")+ theme(plot.margin=unit(c(0.5,0.5,0.5,0.5),"cm"),legend.position = "none") ## - Predict from the best GLM and plot effects of SATL # - Construct new dataset Rainfall_SATL_glm<-expand.grid(MEI=0,Season=levels(as.factor(Rainfall_monthly_E$Season)),SATL_anom=Rainfall_monthly_E$SATL_anom,NATL=0,IOD=0,Year_rescaled=0) Rainfall_SATL_glm$SATL<-scale(Rainfall_SATL_glm$SATL_anom) # - Predict from the best GLM Rainfall_SATL_glm$Rainfall_SATL_glm<-predict(Rainfall_Tele_m1_full,newdata=Rainfall_SATL_glm,type="response",re.form=NA) # - Plot summarised raw data and model prediction (Figure 4) Rainfall_SATL_seasonal_plot<-ggplot()+ geom_point(data=Rainfall_seasonal,aes(x=SATL_anom,y=Rainfall_mm,shape=Season),alpha=0.6)+ geom_line(data=Rainfall_SATL_glm,aes(x=SATL_anom,y=Rainfall_SATL_glm,linetype=Season))+ scale_linetype_manual(values=Linetype_lookup)+ scale_shape_manual(values=Shape_lookup)+ theme_classic()+ ylab("Rainfall (mm/month)")+ xlab("South Atlantic SST anomaly")+ ggtitle("E.")+ theme(plot.margin=unit(c(0.5,0.5,0.5,0.5),"cm"),legend.position = "none") # - Predict from data by season Rainfall_SATL_glm$SATL_anom_round<-round((Rainfall_SATL_glm$SATL*0.328)+0.037,digits=1) Rainfall_SATL_glm$SATL_anom_round<-round(Rainfall_SATL_glm$SATL_anom,digits=1) Rainfall_SATL_diff<-ddply(Rainfall_SATL_glm,.(Season,SATL_anom_round),summarise,Rainfall=mean(Rainfall_SATL_glm,na.rm=T)) Rainfall_SATL_diff$Rainfall[Rainfall_SATL_diff$Season=="MAM"&Rainfall_SATL_diff$SATL_anom_round==0.5]-Rainfall_SATL_diff$Rainfall[Rainfall_SATL_diff$Season=="MAM"&Rainfall_SATL_diff$SATL_anom_round==-0.5] Rainfall_SATL_diff$Rainfall[Rainfall_SATL_diff$Season=="JJAS"&Rainfall_SATL_diff$SATL_anom_round==0.5]-Rainfall_SATL_diff$Rainfall[Rainfall_SATL_diff$Season=="JJAS"&Rainfall_SATL_diff$SATL_anom_round==-0.5] ## - Predict from the best GLM and plot effects of MEI # - Construct new dataset Rainfall_MEI_glm<-expand.grid(MEI_original=Rainfall_monthly_E$MEI_original,Season=levels(as.factor(Rainfall_monthly_E$Season)),NATL=0,SATL=0,IOD=0,Year_rescaled=0) Rainfall_MEI_glm$MEI<-scale(Rainfall_MEI_glm$MEI_original) # - Predict from the best GLM Rainfall_MEI_glm$Rainfall_MEI_glm<-predict(Rainfall_Tele_m1_full,newdata=Rainfall_MEI_glm,type="response",re.form=NA) # - Plot summarised raw data and model prediction (Figure 4) Rainfall_MEI_seasonal_plot<-ggplot()+ geom_point(data=Rainfall_seasonal,aes(x=MEI_original,y=Rainfall_mm,shape=Season),alpha=0.6)+ geom_line(data=Rainfall_MEI_glm,aes(x=MEI_original,y=Rainfall_MEI_glm,linetype=Season))+ scale_shape_manual(values=Shape_lookup)+ scale_linetype_manual(values=Linetype_lookup)+ theme_classic()+ ylab("Rainfall (mm/month)")+ xlab("Multivariate ENSO Index")+ ggtitle("B.")+ theme(plot.margin=unit(c(0.5,0.5,0.5,0.5),"cm"),legend.position = "none") # - Predict from the data by season Rainfall_MEI_glm$MEI_original_round<-round(Rainfall_MEI_glm$MEI_original,digits=1) Rainfall_MEI_diff<-ddply(Rainfall_MEI_glm,.(Season,MEI_original_round),summarise,Rainfall=mean(Rainfall_MEI_glm,na.rm=T)) Rainfall_MEI_diff<-Rainfall_MEI_diff[complete.cases(Rainfall_MEI_diff$MEI_original_round),] Rainfall_MEI_diff$Rainfall[Rainfall_MEI_diff$Season=="DJF"&Rainfall_MEI_diff$MEI_original_round==0.5]-Rainfall_MEI_diff$Rainfall[Rainfall_MEI_diff$Season=="DJF"&Rainfall_MEI_diff$MEI_original_round==-0.5] Rainfall_MEI_diff$Rainfall[Rainfall_MEI_diff$Season=="ON"&Rainfall_MEI_diff$MEI_original_round==0.5]-Rainfall_MEI_diff$Rainfall[Rainfall_MEI_diff$Season=="ON"&Rainfall_MEI_diff$MEI_original_round==-0.5] Rainfall_MEI_diff$Rainfall[Rainfall_MEI_diff$Season=="MAM"&Rainfall_MEI_diff$MEI_original_round==0.5]-Rainfall_MEI_diff$Rainfall[Rainfall_MEI_diff$Season=="MAM"&Rainfall_MEI_diff$MEI_original_round==-0.5] ## - Predict from the best GLM and plot effects of IOD # - Construct new dataset Rainfall_IOD_glm<-expand.grid(MEI=0,Season=levels(as.factor(Rainfall_monthly_E$Season)),SATL=0,NATL=0,IOD_original=Rainfall_monthly_E$IOD_original,Year_rescaled=0) Rainfall_IOD_glm$IOD<-scale(Rainfall_IOD_glm$IOD_original) # - Predict from the best GLM Rainfall_IOD_glm$Rainfall_IOD_glm<-predict(Rainfall_Tele_m1_full,newdata=Rainfall_IOD_glm,type="response",re.form=NA) # - Plot summarised raw data and model prediction (Figure 4) Rainfall_IOD_seasonal_plot<-ggplot()+ geom_point(data=Rainfall_seasonal,aes(x=IOD_original,y=Rainfall_mm,shape=Season),alpha=0.6)+ geom_line(data=Rainfall_IOD_glm,aes(x=IOD_original,y=Rainfall_IOD_glm,linetype=Season))+ scale_linetype_manual(values=Linetype_lookup)+ scale_shape_manual(values=Shape_lookup)+ theme_classic()+ ylab("Rainfall (mm/month)")+ xlab("Indian Ocean Dipole")+ ggtitle("F.")+ theme(plot.margin=unit(c(0.5,0.5,0.5,0.5),"cm"),legend.position = "none") # - Predict from data by season Rainfall_IOD_glm$IOD_original_round<-round(Rainfall_IOD_glm$IOD_original,digits=1) Rainfall_IOD_diff<-ddply(Rainfall_IOD_glm,.(Season,IOD_original_round),summarise,Rainfall=mean(Rainfall_IOD_glm,na.rm=T)) Rainfall_IOD_diff$Rainfall[Rainfall_IOD_diff$Season=="JJAS"&Rainfall_IOD_diff$IOD_original_round==0.5]-Rainfall_IOD_diff$Rainfall[Rainfall_IOD_diff$Season=="JJAS"&Rainfall_IOD_diff$IOD_original_round==-0.5] Rainfall_IOD_diff$Rainfall[Rainfall_IOD_diff$Season=="ON"&Rainfall_IOD_diff$IOD_original_round==0.5]-Rainfall_IOD_diff$Rainfall[Rainfall_IOD_diff$Season=="ON"&Rainfall_IOD_diff$IOD_original_round==-0.5] ### Temperature ### ## - Merge temperature and teleconnection datasets Temperature_monthly_G<-merge(Temperature_monthly_G,Tele[,c(3,4,5,6,8,10,13)],c("Year","Month")) Temperature_monthly_G$Year_rescaled<-scale(Temperature_monthly_G$Year) # - Scale teleconnection data Temperature_monthly_G$MEI_original<-Temperature_monthly_G$MEI Temperature_monthly_G$MEI<-as.numeric(scale(Temperature_monthly_G$MEI)) Temperature_monthly_G$NATL<-scale(Temperature_monthly_G$NATL_anom) Temperature_monthly_G$SATL<-scale(Temperature_monthly_G$SATL_anom) Temperature_monthly_G$IOD_original<-Temperature_monthly_G$IOD Temperature_monthly_G$IOD<-as.numeric(scale(Temperature_monthly_G$IOD)) Temperature_monthly_G<-Temperature_monthly_G[order(Temperature_monthly_G$YearMonth),] ## - Model the effects of ocean temperatures on temperature at Lopé by season # - Explore collinearity #Temperature2<-Temperature_monthly_E[complete.cases(Temperature_monthly_E$MEI)&complete.cases(Temperature_monthly_E$Temperature_month_mm),c()] #names(Temperature2)[2]<-"Temperature" #names(Temperature2)[1]<-"Year" #corrplot(cor(Temperature2),type="upper") # - GLM competition removing terms sequentially starting with full model Temperature_Tele_m1_full<-lmer(Temperature_c~MEI*Season+NATL*Season+SATL*Season+IOD*Season+Year_rescaled*Season+(1|Month)+(1|Year),data=Temperature_monthly_G[complete.cases(Temperature_monthly_G$MEI),],na.action = na.exclude) # - Check for autocorrelation #Temperature_monthly_G$OceansGLMResiduals<-NA #Temperature_monthly_G$OceansGLMResiduals<-resid(Temperature_Tele_m1_full) #acf_plot(Temperature_monthly_G$OceansGLMResiduals,split_by=list(Temperature_monthly_G$DOY),fun=median) #acf_n_plots(Temperature_monthly_G$OceansGLMResiduals,split_by=list(Temperature_monthly_G$DOY)) Temperature_Tele_m2_IODonly<-lmer(Temperature_c~MEI*Season+NATL*Season+SATL*Season+IOD+Year_rescaled*Season+(1|Month)+(1|Year),data=Temperature_monthly_G[complete.cases(Temperature_monthly_G$MEI),]) Temperature_Tele_m2_noIOD<-lmer(Temperature_c~MEI*Season+NATL*Season+SATL*Season+Year_rescaled*Season+(1|Month)+(1|Year),data=Temperature_monthly_G[complete.cases(Temperature_monthly_G$MEI),]) Temperature_Tele_m3_MEIonly<-lmer(Temperature_c~MEI+NATL*Season+SATL*Season+Year_rescaled*Season+(1|Month)+(1|Year),data=Temperature_monthly_G[complete.cases(Temperature_monthly_G$MEI),]) Temperature_Tele_m3_noMEI<-lmer(Temperature_c~NATL*Season+SATL*Season+Year_rescaled*Season+(1|Month)+(1|Year),data=Temperature_monthly_G[complete.cases(Temperature_monthly_G$MEI),]) Temperature_Tele_m4_SATLonly<-lmer(Temperature_c~MEI+NATL*Season+SATL+Year_rescaled*Season+(1|Month)+(1|Year),data=Temperature_monthly_G[complete.cases(Temperature_monthly_G$MEI),]) Temperature_Tele_m4_noSATL<-lmer(Temperature_c~MEI+NATL*Season+Year_rescaled*Season+(1|Month)+(1|Year),data=Temperature_monthly_G[complete.cases(Temperature_monthly_G$MEI),]) Temperature_Tele_m5_NATLonly<-lmer(Temperature_c~MEI+NATL+Year_rescaled*Season+(1|Month)+(1|Year),data=Temperature_monthly_G[complete.cases(Temperature_monthly_G$MEI),]) Temperature_Tele_m5_noNATL<-lmer(Temperature_c~MEI+Year_rescaled*Season+(1|Month)+(1|Year),data=Temperature_monthly_G[complete.cases(Temperature_monthly_G$MEI),]) # - Compare AIC values for all models Temp_AIC_df<-AIC(Temperature_Tele_m1_full, Temperature_Tele_m2_IODonly,Temperature_Tele_m2_noIOD, Temperature_Tele_m3_MEIonly,Temperature_Tele_m3_noMEI, Temperature_Tele_m4_SATLonly,Temperature_Tele_m4_noSATL, Temperature_Tele_m5_NATLonly,Temperature_Tele_m5_noNATL) # Table 5 Temperature_oceans_AIC<-data.frame(Response="Temperature", Predictors=c("Season + NATL:Season + SATL:Season + MEI:Season + IOD:Season + Year:Season", "Season + NATL:Season + SATL:Season + MEI:Season + IOD + Year:Season", "Season + NATL:Season + SATL:Season + MEI:Season + Year:Season", "Season + NATL:Season + SATL:Season + MEI + Year:Season", "Season + NATL:Season + SATL:Season + Year:Season", "Season + NATL:Season + SATL + MEI + Year:Season", "Season + NATL:Season + MEI + Year:Season", "Season + NATL + MEI + Year:Season", "Season + MEI + Year:Season"), AIC=round(Temp_AIC_df$AIC,digits=2), DF=Temp_AIC_df$df) # - Modify best model by removing the intercept to directly compare estimates to each other Temperature_Tele_m5_noNATL_mod<-lmer(Temperature_c~-1+Season+MEI+Year_rescaled:Season+(1|Month)+(1|Year),data=Temperature_monthly_G[complete.cases(Temperature_monthly_G$MEI),]) Temperature_oceans_glm<-tidy(Temperature_Tele_m5_noNATL_mod,conf.int=T) Temperature_oceans_glm$sig<-ifelse(Temperature_oceans_glm$conf.low>0|Temperature_oceans_glm$conf.high<0,"*","") Temperature_oceans_glm$conf.low<-round(Temperature_oceans_glm$conf.low,digits=2) Temperature_oceans_glm$conf.high<-round(Temperature_oceans_glm$conf.high,digits=2) Temperature_oceans_glm<-as.data.frame(Temperature_oceans_glm) Temperature_oceans_glm<-Temperature_oceans_glm[Temperature_oceans_glm$group=="fixed",] Temperature_oceans_glm<-unite(Temperature_oceans_glm,"95%CI",c("conf.low","conf.high"),sep=",") # - Table 7 Temperature_oceans_glm<-data.frame(Predictor=c("DJF","JJAS","MAM","ON","MEI","Year: DJF","Year: JJAS","Year: MAM","Year: ON"), Estimate= round(Temperature_oceans_glm[,"estimate"],digits=2), SE = round(Temperature_oceans_glm[,"std.error"],digits=2), "T Value" = round(Temperature_oceans_glm[,"statistic"],digits=2), "95% CI" =Temperature_oceans_glm[,"95%CI"], "."=Temperature_oceans_glm[,"sig"]) ## - Predict from the best GLM and plot effects of MEI # - Construct new dataset MinTemp_MEI_glm<-expand.grid(MEI_original=Temperature_monthly_G$MEI_original,Season=levels(as.factor(Temperature_monthly_G$Season)),NATL=0,SATL=0,IOD=0,Year_rescaled=0) MinTemp_MEI_glm$MEI<-scale(MinTemp_MEI_glm$MEI_original) # - Summarise rainfall data by season (to add to plot below) Temperature_annual<-ddply(Temperature_monthly_G,.(Year),summarise, NATL=mean(NATL,na.rm=T), SATL=mean(SATL,na.rm=T), MEI_original=mean(MEI_original,na.rm=T), MEI=mean(MEI,na.rm=T), IOD=mean(IOD,na.rm=T), Temperature_c=mean(Temperature_c,na.rm=T)) # - Plot summarised raw data and model prediction (Figure 4C) MinTemp_MEI_seasonal_plot<-ggplot()+ geom_point(data=Temperature_seasonal,aes(x=MEI_original,y=Temperature_c,shape=Season),alpha=0.6)+ geom_line(data=MinTemp_MEI_glm,aes(x=MEI_original,y=MinTemp_MEI_glm,linetype=Season))+ scale_shape_manual(values=Shape_lookup)+ scale_linetype_manual(values=Linetype_lookup)+ scale_y_continuous(breaks=c(22,23))+ theme_classic()+ ylab("Min. Temp. (c)")+ xlab("Multivariate ENSO Index")+ ggtitle("C.")+ theme(plot.margin=unit(c(0.5,0.5,0.5,0.5),"cm"),legend.position = "none") # - Predict from the best GLM MinTemp_MEI_glm$MinTemp_MEI_glm<-predict(Temperature_Tele_m5_noNATL,newdata=MinTemp_MEI_glm,type="response",re.form=NA) MinTemp_MEI_glm$MEI_original_round<-round(MinTemp_MEI_glm$MEI_original,digits=1) MinTemp_MEI_diff$MinTemp[MinTemp_MEI_diff$Season=="DJF"&MinTemp_MEI_diff$MEI_original_round==0.5]-MinTemp_MEI_diff$MinTemp[MinTemp_MEI_diff$Season=="DJF"&MinTemp_MEI_diff$MEI_original_round==-0.5] MinTemp_MEI_diff$MinTemp[MinTemp_MEI_diff$Season=="MAM"&MinTemp_MEI_diff$MEI_original_round==0.5]-MinTemp_MEI_diff$MinTemp[MinTemp_MEI_diff$Season=="MAM"&MinTemp_MEI_diff$MEI_original_round==-0.5] ### Combine all estimates ### ## - Combine data to create Figure 4A names(Rainfall_oceans_glm)[4]<-"Value" Rainfall_oceans_glm<-Rainfall_oceans_glm[c(5:20),] Rainfall_oceans_glm$Season<-c("DJF","JJAS","MAM","ON") Rainfall_oceans_glm$OceanIndex<-rep(c("MEI","NATL","SATL","IOD"),each=4) Rainfall_oceans_glm$Response<-"Rainfall" names(Temperature_oceans_glm)[4]<-"Value" Temperature_oceans_glm<-Temperature_oceans_glm[c(5,5,5,5),] Temperature_oceans_glm$Season<-c("DJF","JJAS","MAM","ON") Temperature_oceans_glm$OceanIndex<-rep(c("MEI"),each=4) Temperature_oceans_glm$Response<-"Minimum Temperature" Estimates_all_combined<-rbind(Rainfall_oceans_glm,Temperature_oceans_glm) Estimates_all_combined$Trend<-Estimates_all_combined$Estimate>0 Estimates_all_combined$Season<-factor(Estimates_all_combined$Season,levels=c("DJF","MAM","JJAS","ON","All"),ordered=T) Estimates_all_combined$OceanIndex<-as.factor(Estimates_all_combined$OceanIndex) Estimates_all_combined$OceanIndex<-factor(Estimates_all_combined$OceanIndex,levels=c("IOD","SATL","NATL","MEI"),ordered=T) Estimates_all_combined$Response<-factor(Estimates_all_combined$Response,levels=c("Rainfall","Minimum Temperature"),ordered=T) names(Estimates_all_combined)[6]<-"Sig" # - plot estimate data to create Figure 4A Summary_plot<-ggplot(Estimates_all_combined,aes(x=Season,y=OceanIndex))+ geom_point(data=Estimates_all_combined[Estimates_all_combined$Sig=="*",],aes(colour=Trend,size=abs(Estimate),alpha=abs(Value)))+ scale_size_continuous(range=c(5,10))+ scale_colour_manual(values=c("TRUE"="darkslategray3","FALSE"="coral"))+ scale_alpha_continuous(range=c(0.3,1))+ ylab("Oceanic influence")+ theme_classic()+ theme(legend.position="none",panel.border = element_rect(colour = "black", fill=NA, size=0.5), axis.line = element_line(size=0.5), panel.grid.major=element_line(size=0.2), plot.margin=unit(c(0.5,0.5,0.5,0.5),"cm"))+ ggtitle("A.")+ facet_wrap(~Response,ncol=3) ## - Combine all prediction plots to make Figure 4 NA_plot<-ggplot()+theme_minimal() legend_plot<-ggplot()+ geom_point(data=Rainfall_monthly_E,aes(x=MEI,y=Rainfall_mm,shape=Season),alpha=0.6)+ geom_line(data=Rainfall_MEI_glm,aes(x=MEI,y=Rainfall_MEI_glm,linetype=Season))+ scale_shape_manual(values=Shape_lookup)+ scale_linetype_manual(values=Linetype_lookup)+ theme_classic()+ ylab("Rainfall (mm/month)")+ xlab("Multivariate ENSO Index")+ theme(legend.position = "right") legend_plot<-get_legend(legend_plot) #grid.arrange(Summary_plot, # Rainfall_MEI_seasonal_plot,MinTemp_MEI_seasonal_plot, # Rainfall_NATL_seasonal_plot,legend_plot, # Rainfall_SATL_seasonal_plot,NA_plot, # Rainfall_IOD_seasonal_plot,NA_plot, # ncol = 2, nrow = 5,heights=c(1.2,1,1,1,1), # layout_matrix = rbind(c(1,1), c(2,3),c(4,5),c(6,7),c(8,9)))