# Code by Alastair Potts
# 9 October 2023
# 

setwd("D:\\Google Drive\\0_SpekboomResearchGroup2.0\\3_Experiments\\RegionalTransect21(10)_ParentPlantMaterialExperiment\\manuscript/")


library(groundhog)
set.groundhog.folder("C:\\Groundhog_R\\R4.2.3_2023-04-01")
GroundhogDay <- '2023-04-01'

#groundhog.library("tidyverse",GroundhogDay)
groundhog.library("slider",GroundhogDay)
groundhog.library("glue",GroundhogDay)
groundhog.library("lubridate",GroundhogDay)
groundhog.library("drought",GroundhogDay)
groundhog.library("dplyr",GroundhogDay)
groundhog.library("tidyr",GroundhogDay)
groundhog.library("ggplot2",GroundhogDay)
groundhog.library("patchwork",GroundhogDay)

qCols <- rev(c("#3B9AB2", "#78B7C5", "#74A089", "#E1AF00", "#F21A00"))

read.csv("data_manuscript/chirps_data.csv")->chirps.dat

### Converting daily values into monthly values
chirps.dat%>%
  as_tibble()%>%
  select(Site,date,chirps)%>%
  mutate(year=year(date),
         month=month(date))%>%
  group_by(Site,year,month)%>%
  summarise(ppt=sum(chirps))%>%
  ungroup()->
  chirps.dat_monthly

# Generating 12 month spi values
chirps.dat_monthly%>%
  group_by(Site)%>%
  #filter(year<2022)%>%
  group_by(Site)%>%
  mutate(spi=
           SDI(ppt,12)%>%as.data.frame()%>%
           pivot_longer(cols=1:12,names_to = "Temp",values_to = "spi3")%>%
           select("spi3")%>%unlist())%>%
  mutate(Date=(ym(paste(year,month))))->
  forGG

forGG%>%
  group_by(Site,month)%>%
  summarise(u_ppt=mean(ppt))->
  meanPPT
###

### Figure 1
forGG%>%
  ungroup()%>%
  left_join(meanPPT)%>%
  group_by(month)%>%
  mutate(posneg=ifelse(ppt>u_ppt,"darkgreen","red"))%>%
  ungroup()%>%
  filter(year>=2020,year<2022)%>%
  mutate(Date=ymd(paste(year,month,14)))%>%
  filter(Site==7)%>%
  ggplot()+
  #geom_point(aes(Date,u_ppt))+
  geom_line(aes(Date,u_ppt),linewidth=1.5,alpha=0.5)+
  #geom_line(aes(Date,q33))+
  geom_segment(aes(x=Date,xend=Date,y=u_ppt,
                   yend=ppt,
                   colour=as.factor(posneg)),
               linewidth=3.5)+
  scale_colour_manual(values=qCols[c(5,1)],labels=c("+","-"))+
  geom_point(aes(Date,u_ppt),size=3)+
  #geom_line(aes(Date,ppt))+
  geom_vline(xintercept=ymd(c("2021-01-01","2022-01-1")))+
  geom_vline(xintercept=ymd("2021-10-13"),linetype='dashed',lwd=2)+
  theme_bw()+
  theme(legend.position="none")+
  xlab("")+
  ylab("Monthly precipitation (mm)")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())-> B
B
#ggsave("results/FigA_PptAtSite4.jpg",width=8,height=4)
#ggsave("manuscript/figures/FigA_PptAtSite4.jpg",width=8,height=4)
###


forGG%>%
  filter(Site==7)%>%
  filter(year>=2020,year<2022)%>%
  #mutate(ym=ceiling_date(Date,unit="month")-1)%>%
  mutate(ym=ymd(paste(year,month,14)))%>%
  ggplot(aes(ym,spi))+
  geom_hline(yintercept = c(-2,-1,1,2),linetype="dashed")+
  geom_segment(aes(x=ym,y=0,xend=ym,yend=spi,color=spi),size=3)+
  scale_color_gradient2(midpoint=0,high="blue",mid=("green"),low=("red"))+
  geom_line(size=1.5,alpha=0.5)+
  theme_bw()+
  geom_vline(xintercept=ymd(c("2021-01-01","2022-01-1")))+
  scale_y_continuous(breaks=seq(-3,3),limits=c(-2,2))+
  # scale_x_date(date_breaks = "year",date_labels = "%Y",
  #              limits = c(ymd("1984-01-01"), ymd("2022-12-31")),
  #              expand=c(0,0))+
  ylab("12-month SPI")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
        legend.position="none")+
  xlab("")+
  geom_vline(xintercept=ymd("2021-10-13"),linetype='dashed',lwd=2)+
  geom_hline(yintercept = c(0),size=1.3)+
  scale_x_date(date_breaks = "1 month", date_labels =  "%b %Y")+
  theme(axis.text.x=element_text(angle=60, hjust=1))->
  A


(B+ggtitle("A)"))/(A+ggtitle("B)"))
ggsave("FinalFiguresV2/Fig_1.jpg",width=8,height=6.5)