##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") library(cplm) library(biwavelet) ########################################## ##### 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 cpglm Rain_trend_m1.cpglmm<-cpglmm(Rainfall_mm~Year_r+(1|DOY), data=Rainfall_daily_D,na.action=na.exclude) Rain_trend_m2.cpglmm<-cpglmm(Rainfall_mm~1+(1|DOY), data=Rainfall_daily_D,na.action=na.exclude) # - Table 2 Rainfall_trend_AIC<-data.frame(Response="Rainfall", Model=c("Long-term change", "No long-term change"), Predictors=c("Year","Intercept only"), DF= AIC(Rain_trend_m1.cpglmm,Rain_trend_m2.cpglmm)$df, AIC=round(AIC(Rain_trend_m1.cpglmm,Rain_trend_m2.cpglmm)$AIC,digits=1)) Rainfall_trend_AIC$Delta_AIC<-Rainfall_trend_AIC$AIC-min(Rainfall_trend_AIC$AIC) # - 95%CI c(round(summary(Rain_trend_m1.cpglmm)$coefs[2,1] - 1.96* summary(Rain_trend_m1.cpglmm)$coefs[2,2],digits=2), round(summary(Rain_trend_m1.cpglmm)$coefs[2,1] + 1.96* summary(Rain_trend_m1.cpglmm)$coefs[2,2],digits=2)) # - Predict trend Rainfall_daily_D$Rainfall_daily_mm_predict<-predict(Rain_trend_m1.cpglmm,Rainfall_daily_D,type="response") Rainfall_yearly_F$Rainfall_annual_mm_predict<-ddply(Rainfall_daily_D,.(Year),summarise, sum(Rainfall_daily_mm_predict))[,2] # - Change in total annual rainfall (mm) per decade = -75mm per decade #(Rainfall_yearly_F$Rainfall_annual_mm_predict[Rainfall_yearly_F$Year==2017]-Rainfall_yearly_F$Rainfall_annual_mm_predict[Rainfall_yearly_F$Year==1984])/(2017-1984)*10 #Percentage change compared to mean total annual rainfall = -0.0547 #(Rainfall_yearly_F$Rainfall_annual_mm_predict[Rainfall_yearly_F$Year==2017]-Rainfall_yearly_F$Rainfall_annual_mm_predict[Rainfall_yearly_F$Year==1984])/(2017-1984)*10/mean(Rainfall_yearly_F$Rainfall_annual_mm_predict,na.rm=T) # - Plot trend and raw data (Figure 3A) Rainfall_annual_glm_plot<-ggplot()+geom_line(data=Rainfall_yearly_F[c(1,34),],aes(x=Year,y=Rainfall_annual_mm_predict))+ geom_point(data=Rainfall_yearly_F,aes(x=Year,y=Rainfall_yearly_mm),colour="gray")+ geom_line(data=Rainfall_yearly_F,aes(x=Year,y=Rainfall_yearly_mm),colour="gray")+ 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 # - CPGLM on daily data including interaction with season Rain_seasonal_m1.cpglmm<-cpglmm(Rainfall_mm~Season*Year_r+(1|DOY), data=Rainfall_daily_D,na.action=na.exclude) Rain_seasonal_m2.cpglmm<-cpglmm(Rainfall_mm~Season+Year_r+(1|DOY), data=Rainfall_daily_D,na.action=na.exclude) # - Table 3 Rainfall_seasonal_trend_AIC<-data.frame(Response="Rainfall", Model=c("Long-term change by season", "Long-term change not by season"), Predictors=c("Year * Season","Year + Season"), DF= AIC(Rain_seasonal_m1.cpglmm,Rain_seasonal_m2.cpglmm)$df, AIC=round(AIC(Rain_seasonal_m1.cpglmm,Rain_seasonal_m2.cpglmm)$AIC,digits=1)) Rainfall_seasonal_trend_AIC$Delta_AIC<-Rainfall_seasonal_trend_AIC$AIC-min(Rainfall_seasonal_trend_AIC$AIC) #Remove intercept to directly compare estimates Rain_seasonal_m1.cpglmm_mod<-cpglmm(Rainfall_mm~-1+Season+Season:Year_r+(1|DOY), data=Rainfall_daily_D,na.action=na.exclude) # - Extract estimates Rain_seasonal_m1.cpglmm_mod_coefs <- data.frame(Response="Rainfall", Predictor=c("DJF","JJAS","MAM","ON", "Year: DJF", "Year: JJAS","Year: MAM","Year: ON"), Est = round(Rain_seasonal_m1.cpglmm_mod$fixef,digits=2) , S.E. = round(summary(Rain_seasonal_m1.cpglmm_mod)$coefs[,2],digits=2), T.Value= round(summary(Rain_seasonal_m1.cpglmm_mod)$coefs[,3],digits=2)) Rain_seasonal_m1.cpglmm_mod_coefs$LL = round(Rain_seasonal_m1.cpglmm_mod_coefs$Est - 1.96 * Rain_seasonal_m1.cpglmm_mod_coefs$S.E., digits=2) Rain_seasonal_m1.cpglmm_mod_coefs$UL = round(Rain_seasonal_m1.cpglmm_mod_coefs$Est + 1.96 * Rain_seasonal_m1.cpglmm_mod_coefs$S.E., digits=2) row.names(Rain_seasonal_m1.cpglmm_mod_coefs)<-NULL #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_seasonal<-ddply(Rainfall_daily_D,.(Year,Season),summarise, Rainfall_daily_mean=mean(Rainfall_mm,na.rm=T), Rainfall_daily_predict_mm=mean(Rainfall_daily_mm_predict)) # ddply(Rainfall_seasonal,.(Season),summarise, # Rainfall_1984=Rainfall_daily_predict_mm[Year==1984], # Rainfall_2017=Rainfall_daily_predict_mm[Year==2017], # 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_daily_mean,na.rm=T), # Difference_per_dec=round((Difference_mm_dec/Rainfall_mean_mm)*100,digits=2)) ### 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),data=Temperature_daily_F,na.action = na.exclude) Temp_annual_m2<-lmer(Temperature_c~1+(1|DOY),data=Temperature_daily_F,na.action = na.exclude) # - 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"), DF= AIC(Temp_annual_m1,Temp_annual_m2)$df, AIC=round(AIC(Temp_annual_m1,Temp_annual_m2)$AIC,digits=1)) Temperature_trend_AIC$Delta_AIC<-Temperature_trend_AIC$AIC-min(Temperature_trend_AIC$AIC) # - 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_LMM_c[Temperature_annual_LMM$Year==2017]-Temperature_annual_LMM$Temperature_LMM_c[Temperature_annual_LMM$Year==1984])/(2017-1984)*10 #0.25 per decade #(Temperature_annual_LMM$Temperature_LMM_c[Temperature_annual_LMM$Year==2017]-Temperature_annual_LMM$Temperature_LMM_c[Temperature_annual_LMM$Year==1984])/(2017-1984)*10/mean(Temperature_daily_F$Temperature_c,na.rm=T) #1.1% # - Figure 3B Temperature_ym_LMM<-ddply(Temperature_daily_F,.(YearMonth),summarise, Temperature_mean_c=mean(Temperature_c,na.rm=T), Temperature_mean_sample=length(Temperature_c[complete.cases(Temperature_c)]), Temperature_LMM_c=mean(LMM_predict,na.rm=T)) Temperature_ym_LMM$Temperature_mean_c[Temperature_ym_LMM$Temperature_mean_sample<5]<-NA MinTemp_annual_glm_plot<-ggplot()+ geom_point(data=Temperature_ym_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_ym_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),data=Temperature_daily_F,na.action = na.exclude) Temp_seasonal_m2<-lmer(Temperature_c~Year_r+Season+(1|DOY),data=Temperature_daily_F,na.action = na.exclude) # - Table 3 Temperature_seasonal_trend_AIC<-data.frame(Response="Temperature", Model=c("Long-term change by season", "Long-term change not by season"), Predictors=c("Year * Season","Year + Season"), DF= AIC(Temp_seasonal_m1,Temp_seasonal_m2)$df, AIC=round(AIC(Temp_seasonal_m1,Temp_seasonal_m2)$AIC,digits=1)) Temperature_seasonal_trend_AIC$Delta_AIC<-Temperature_seasonal_trend_AIC$AIC-min(Temperature_seasonal_trend_AIC$AIC) # - Tabke 4 - Modify best model by removing intercept Temp_seasonal_m1_mod<-lmer(Temperature_c~-1+Season+Year_r:Season+(1|DOY),data=Temperature_daily_F,na.action = na.exclude) summary(Temp_seasonal_m1_mod)$coefficients Temp_seasonal_m1_mod_coefs <- data.frame(Response="Temperature", Predictor=c("DJF","JJAS","MAM","ON","Year: DJF", "Year: JJAS","Year: MAM","Year: ON"), Est = round(summary(Temp_seasonal_m1_mod)$coefficients[,1],digits=2), S.E. = round(summary(Temp_seasonal_m1_mod)$coefficients[,2],digits=2), T.Value= round(summary(Temp_seasonal_m1_mod)$coefficients[,3],digits=2)) Temp_seasonal_m1_mod_coefs$LL = round(Temp_seasonal_m1_mod_coefs$Est - 1.96 * Temp_seasonal_m1_mod_coefs$S.E., digits=2) Temp_seasonal_m1_mod_coefs$UL = round(Temp_seasonal_m1_mod_coefs$Est + 1.96 * Temp_seasonal_m1_mod_coefs$S.E., digits=2) row.names(Temp_seasonal_m1_mod_coefs)<-NULL # - 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),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.16 per decade ## - CRU CRU_annual_m1<-lmer(Temperature_c~Year_r+(1|Month),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,type="response") #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.19 per decade #################################### ##### 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 ##### ################################### ## - Load teleconnection data # - Multivariate ENSO Index (MEI; Wolter & Timlin 1993; Wolter & Timlin 1998) sourced from the NOAA website (https://www.esrl.noaa.gov/psd/enso/mei/index.html) # - Indian Ocean Dipole (IOD) Dipole Mode Index (Saji & Yamagata 2003) sourced from the NOAA website (https://www.esrl.noaa.gov/psd/gcos_wgsp/Timeseries/DMI/) # - SST anomalies for the tropical north Atlantic (NATL, 5°–20°N, 60°–30°W) and the south equatorial Atlantic (SATL ,0°–20°S, 30°W–10°E) sourced from the NOAA National Weather Service Climate Prediction Center (http://www.cpc.ncep.noaa.gov/data/indices/). Tele <- read.csv("../Data/Teleconnections_combined") Tele$YearMonth<-as.yearmon(as.Date(paste(Tele$Year,Tele$Month,"01",sep="-"))) ### Rainfall ### ## - Merge rainfall and teleconnection datasets Rainfall_Tele_monthly<-merge(Rainfall_monthly_E[,c("Month","Year","YearMonth","Season","Rainfall_interp_mm")],Tele[,c("Year","Month","MEI","NATL","NATL_anom","SATL","SATL_anom","IOD","NAO")],c("Year","Month")) Rainfall_Tele_monthly<-Rainfall_Tele_monthly[order(Rainfall_Tele_monthly$YearMonth),] Rainfall_Tele_monthly<-Rainfall_Tele_monthly[complete.cases(Rainfall_Tele_monthly$MEI),] #MEI #plot(cbind(Rainfall_Tele_monthly[,"YearMonth"],Rainfall_Tele_monthly[,"MEI"]),type="l") #MEI_wavelet<-wt(cbind(Rainfall_Tele_monthly[,"YearMonth"],Rainfall_Tele_monthly[,"MEI"]),dj=1/12) #plot(MEI_wavelet,xlab="",ylab="MEI (yrs)") Rain_MEI_wtc<-wtc(cbind(Rainfall_Tele_monthly[,"YearMonth"], scale(Rainfall_Tele_monthly[,"Rainfall_interp_mm"])), cbind(Rainfall_Tele_monthly[,"YearMonth"], scale(Rainfall_Tele_monthly[,"MEI"])), pad=T,mother="morlet",dj=1/12,nrands=1000) #Remove wavelet coherence values within the cone of influence test<-NA for( i in 1:length(Rain_MEI_wtc$coi)) test<-cbind(test,Rain_MEI_wtc$period>Rain_MEI_wtc$coi[i]) test<-test[,-1] Rain_MEI_wtc$rsq[test]<-NA Rain_wtc_global<-data.frame(period=Rain_MEI_wtc$period, global.coherence=rowMeans(Rain_MEI_wtc$rsq,na.rm=T), x="Rainfall", y="MEI") #NATL_anom_anom #plot(cbind(Rainfall_Tele_monthly[,"YearMonth"], Rainfall_Tele_monthly[,"NATL_anom"]),type="l") #NATL_anom_wavelet<-wt(cbind(Rainfall_Tele_monthly[,"YearMonth"],Rainfall_Tele_monthly[,"NATL_anom"]), dj=1/12) #plot(NATL_anom_wavelet,xlab="",ylab="NATL_anom (yrs)") Rain_NATL_anom_wtc<-wtc(cbind(Rainfall_Tele_monthly[,"YearMonth"], scale(Rainfall_Tele_monthly[,"Rainfall_interp_mm"])), cbind(Rainfall_Tele_monthly[,"YearMonth"], scale(Rainfall_Tele_monthly[,"NATL_anom"])), pad=T,mother="morlet",dj=1/12,nrands=1000) #Remove wavelet coherence values within the cone of influence test<-NA for( i in 1:length(Rain_NATL_anom_wtc$coi)) test<-cbind(test,Rain_NATL_anom_wtc$period>Rain_NATL_anom_wtc$coi[i]) test<-test[,-1] Rain_NATL_anom_wtc$rsq[test]<-NA Rain_wtc_global<-rbind(Rain_wtc_global, data.frame(period=Rain_NATL_anom_wtc$period, global.coherence=rowMeans(Rain_NATL_anom_wtc$rsq,na.rm=T), x="Rainfall", y="NATL_anom")) #SATL_anom #plot(cbind(Rainfall_Tele_monthly[,"YearMonth"], Rainfall_Tele_monthly[,"SATL_anom"]),type="l") #SATL_anom_wavelet<-wt(cbind(Rainfall_Tele_monthly[,"YearMonth"],Rainfall_Tele_monthly[,"SATL_anom"]),dj=1/12) #plot(SATL_anom_wavelet,xlab="",ylab="SATL_anom (yrs)") Rain_SATL_anom_wtc<-wtc(cbind(Rainfall_Tele_monthly[,"YearMonth"], scale(Rainfall_Tele_monthly[,"Rainfall_interp_mm"])), cbind(Rainfall_Tele_monthly[,"YearMonth"], scale(Rainfall_Tele_monthly[,"SATL_anom"])), pad=T,mother="morlet",dj=1/12,nrands=1000) #Remove wavelet coherence values within the cone of influence test<-NA for( i in 1:length(Rain_SATL_anom_wtc$coi)) test<-cbind(test,Rain_SATL_anom_wtc$period>Rain_SATL_anom_wtc$coi[i]) test<-test[,-1] Rain_SATL_anom_wtc$rsq[test]<-NA Rain_wtc_global<-rbind(Rain_wtc_global, data.frame(period=Rain_SATL_anom_wtc$period, global.coherence=rowMeans(Rain_SATL_anom_wtc$rsq,na.rm=T), x="Rainfall", y="SATL_anom")) #IOD #plot(cbind(Rainfall_Tele_monthly[,"YearMonth"],Rainfall_Tele_monthly[,"IOD"]),type="l") #IOD_wavelet<-wt(cbind(Rainfall_Tele_monthly[,"YearMonth"],Rainfall_Tele_monthly[,"IOD"]),dj=1/12) #plot(IOD_wavelet,xlab="",ylab="IOD (yrs)") Rain_IOD_wtc<-wtc(cbind(Rainfall_Tele_monthly[,"YearMonth"], scale(Rainfall_Tele_monthly[,"Rainfall_interp_mm"])), cbind(Rainfall_Tele_monthly[,"YearMonth"], scale(Rainfall_Tele_monthly[,"IOD"])), pad=T,mother="morlet",dj=1/12,nrands=1000) #Remove wavelet coherence values within the cone of influence test<-NA for( i in 1:length(Rain_IOD_wtc$coi)) test<-cbind(test,Rain_IOD_wtc$period>Rain_IOD_wtc$coi[i]) test<-test[,-1] Rain_IOD_wtc$rsq[test]<-NA Rain_wtc_global<-rbind(Rain_wtc_global, data.frame(period=Rain_IOD_wtc$period, global.coherence=rowMeans(Rain_IOD_wtc$rsq,na.rm=T), x="Rainfall", y="IOD")) ### Temperature ### ## - Merge temperature and teleconnection datasets Temperature_Tele_monthly<-merge(Temperature_monthly_G[,c("Month","Year","YearMonth","Season","Temperature_interp_c"),], Tele[,c("Year","Month","MEI","NATL","NATL_anom","SATL","SATL_anom","IOD","NAO")],c("Year","Month")) Temperature_Tele_monthly<-Temperature_Tele_monthly[order(Temperature_Tele_monthly$YearMonth),] Temperature_Tele_monthly<-Temperature_Tele_monthly[complete.cases(Temperature_Tele_monthly$MEI),] #MEI #plot(cbind(Temperature_Tele_monthly[,"YearMonth"],Temperature_Tele_monthly[,"MEI"]),type="l") #MEI_wavelet<-wt(cbind(Temperature_Tele_monthly[,"YearMonth"],Temperature_Tele_monthly[,"MEI"]),dj=1/12) #plot(MEI_wavelet,xlab="",ylab="MEI (yrs)") Temp_MEI_wtc<-wtc(cbind(Temperature_Tele_monthly[,"YearMonth"], scale(Temperature_Tele_monthly[,"Temperature_interp_c"])), cbind(Temperature_Tele_monthly[,"YearMonth"], scale(Temperature_Tele_monthly[,"MEI"])), pad=T,mother="morlet",dj=1/12,nrands=1000) #Remove wavelet coherence values within the cone of influence test<-NA for( i in 1:length(Temp_MEI_wtc$coi)) test<-cbind(test,Temp_MEI_wtc$period>Temp_MEI_wtc$coi[i]) test<-test[,-1] Temp_MEI_wtc$rsq[test]<-NA Temp_wtc_global<-data.frame(period=Temp_MEI_wtc$period, global.coherence=rowMeans(Temp_MEI_wtc$rsq,na.rm=T), x="Temperature", y="MEI") #NATL_anom # plot(cbind(Temperature_Tele_monthly[,"YearMonth"], Temperature_Tele_monthly[,"NATL_anom"]),type="l") # NATL_anom_wavelet<-wt(cbind(Temperature_Tele_monthly[,"YearMonth"], Temperature_Tele_monthly[,"NATL_anom"]),dj=1/12) Temp_NATL_anom_wtc<-wtc(cbind(Temperature_Tele_monthly[,"YearMonth"], scale(Temperature_Tele_monthly[,"Temperature_interp_c"])), cbind(Temperature_Tele_monthly[,"YearMonth"], scale(Temperature_Tele_monthly[,"NATL_anom"])), pad=T,mother="morlet",dj=1/12,nrands=1000) #Remove wavelet coherence values within the cone of influence test<-NA for( i in 1:length(Temp_NATL_anom_wtc$coi)) test<-cbind(test,Temp_NATL_anom_wtc$period>Temp_NATL_anom_wtc$coi[i]) test<-test[,-1] Temp_NATL_anom_wtc$rsq[test]<-NA Temp_wtc_global<-rbind(Temp_wtc_global, data.frame(period=Temp_NATL_anom_wtc$period, global.coherence=rowMeans(Temp_NATL_anom_wtc$rsq,na.rm=T), x="Temperature", y="NATL_anom")) #SATL_anom #plot(cbind(Temperature_Tele_monthly[,"YearMonth"],Temperature_Tele_monthly[,"SATL_anom"]),type="l") #SATL_anom_wavelet<-wt(cbind(Temperature_Tele_monthly[,"YearMonth"], Temperature_Tele_monthly[,"SATL_anom"]),dj=1/12) Temp_SATL_anom_wtc<-wtc(cbind(Temperature_Tele_monthly[,"YearMonth"], scale(Temperature_Tele_monthly[,"Temperature_interp_c"])), cbind(Temperature_Tele_monthly[,"YearMonth"], scale(Temperature_Tele_monthly[,"SATL_anom"])), pad=T,mother="morlet",dj=1/12,nrands=1000) #Remove wavelet coherence values within the cone of influence test<-NA for( i in 1:length(Temp_SATL_anom_wtc$coi)) test<-cbind(test,Temp_SATL_anom_wtc$period>Temp_SATL_anom_wtc$coi[i]) test<-test[,-1] Temp_SATL_anom_wtc$rsq[test]<-NA Temp_wtc_global<-rbind(Temp_wtc_global, data.frame(period=Temp_SATL_anom_wtc$period, global.coherence=rowMeans(Temp_SATL_anom_wtc$rsq,na.rm=T), x="Temperature", y="SATL_anom")) #IOD #plot(cbind(Temperature_Tele_monthly[,"YearMonth"],Temperature_Tele_monthly[,"IOD"]),type="l") #IOD_wavelet<-wt(cbind(Temperature_Tele_monthly[,"YearMonth"],Temperature_Tele_monthly[,"IOD"]),dj=1/12) #plot(IOD_wavelet,xlab="",ylab="IOD (yrs)") Temp_IOD_wtc<-wtc(cbind(Temperature_Tele_monthly[,"YearMonth"], scale(Temperature_Tele_monthly[,"Temperature_interp_c"])), cbind(Temperature_Tele_monthly[,"YearMonth"], scale(Temperature_Tele_monthly[,"IOD"])), pad=T,mother="morlet",dj=1/12,nrands=1000) #Remove wavelet coherence values within the cone of influence test<-NA for( i in 1:length(Temp_IOD_wtc$coi)) test<-cbind(test,Temp_IOD_wtc$period>Temp_IOD_wtc$coi[i]) test<-test[,-1] Temp_IOD_wtc$rsq[test]<-NA Temp_wtc_global<-rbind(Temp_wtc_global, data.frame(period=Temp_IOD_wtc$period, global.coherence=rowMeans(Temp_IOD_wtc$rsq,na.rm=T), x="Temperature", y="IOD")) ## - Figure 4 - Plot all wavelet coherency plots together #Note on phase arrows on wavelet coherence plots # Arrows pointing to the right mean that x and y are in phase. # Arrows pointing to the left mean that x and y are in anti-phase. # Arrows pointing up mean that y leads x by π/2. # Arrows pointing down mean that x leads y by π/2. layout(matrix(c(1,2,3,4,5,6,7,8), 4, 2, byrow = TRUE)) par(oma = c(0, 0, 0, 0), mar = c(2,4,2,2) ) #layout.show() plot(Rain_MEI_wtc, plot.cb = FALSE, plot.coi=TRUE,lwd.sig=1.5,plot.phase = TRUE,ylab="Period (yrs)") title(main="A. Rainfall-MEI",adj=0) plot(Temp_MEI_wtc, plot.cb = FALSE, plot.coi=TRUE,lwd.sig=1.5,plot.phase = TRUE,ylab="Period (yrs)") title(main="B. Temp-MEI",adj=0) plot(Rain_NATL_anom_wtc, plot.cb = FALSE, plot.coi=TRUE,lwd.sig=1.5,plot.phase = TRUE,ylab="Period (yrs)") title(main="C. Rain-NATL",adj=0) plot(Temp_NATL_anom_wtc, plot.cb = FALSE, plot.coi=TRUE,lwd.sig=1.5,plot.phase = TRUE,ylab="Period (yrs)") title(main="D. Temp-NATL",adj=0) plot(Rain_SATL_anom_wtc, plot.cb = FALSE, plot.coi=TRUE,lwd.sig=1.5,plot.phase = TRUE,ylab="Period (yrs)") title(main="E. Rain-SATL",adj=0) plot(Temp_SATL_anom_wtc, plot.cb = FALSE, plot.coi=TRUE,lwd.sig=1.5,plot.phase = TRUE,ylab="Period (yrs)") title(main="F. Temp-SATL",adj=0) plot(Rain_IOD_wtc, plot.cb = FALSE, plot.coi=TRUE,lwd.sig=1.5,plot.phase = TRUE,ylab="Period (yrs)") title(main="G. Rain-IOD",adj=0) plot(Temp_IOD_wtc, plot.cb = FALSE, plot.coi=TRUE,lwd.sig=1.5,plot.phase = TRUE,ylab="Period (yrs)") title(main="H. Temp-IOD",adj=0) ## - Supplemental Figure 3 - Plot global wavelet coherency together Temp_Rain_wtc_global<-rbind(Rain_wtc_global,Temp_wtc_global) ggplot(Temp_Rain_wtc_global,aes(y=global.coherence,x=period))+ geom_line(aes(linetype=y,colour=y))+ scale_color_brewer(type="div")+ scale_y_continuous(limits=c(0,1))+ scale_x_log10(breaks=c(0.25,0.5,1,2,4,8))+ #scale_x_continuous(trans=reverselog_trans(10),breaks=c(0.25,0.5,1,2,4,8))+ theme_classic()+ theme(legend.title=element_blank())+ xlab("Period (years)")+ ylab("Wavelet Coherence")+ #coord_flip()+ facet_wrap(~x,ncol=1)