# Code by Alastair Potts
# 9 October 2023
# 

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


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/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)")+
  scale_x_date(date_breaks = "1 month", date_labels =  "%b %Y")+
  theme(axis.text.x=element_text(angle=60, hjust=1))
ggsave("results/FigA_PptAtSite4.jpg",width=8,height=4)
ggsave("manuscript/figures/FigA_PptAtSite4.jpg",width=8,height=4)
###

################################################################################
###



###

chirps.dat%>%
  as_tibble()%>%
  mutate(yr=year(date),
         mn=month(date))%>%
  group_by(yr,mn,Site)%>%
  summarise(ppt=sum(chirps))%>%
  mutate(dt=ymd(paste(yr,mn,"15",sep="-")))%>%
  ggplot(aes(dt,ppt))+
  geom_line()+
  facet_wrap(~Site,ncol=1)

chirps.dat%>%
  as_tibble()%>%
  mutate(date=ymd(date))%>%
  filter(date<=ymd("2021-10-13"))%>%
  filter(date>=ymd("2021-10-13")-months(3))%>%
  group_by(Site)%>%
  summarise(threemonth_ppt=sum(chirps))->month3
chirps.dat%>%
  as_tibble()%>%
  mutate(date=ymd(date))%>%
  filter(date<=ymd("2021-10-13"))%>%
  filter(date>=ymd("2021-10-13")-months(3))%>%
  group_by(Site)%>%
  summarise(threemonth_ppt=sum(chirps))->month3
chirps.dat%>%
  as_tibble()%>%
  mutate(date=ymd(date))%>%
  filter(date<=ymd("2021-10-13"))%>%
  filter(date>=ymd("2021-10-13")-months(6))%>%
  group_by(Site)%>%
  summarise(sixmonth_ppt=sum(chirps))->month6
chirps.dat%>%
  as_tibble()%>%
  mutate(date=ymd(date))%>%
  filter(date<=ymd("2021-10-13"))%>%
  filter(date>=ymd("2021-10-13")-months(9))%>%
  group_by(Site)%>%
  summarise(ninemonth_ppt=sum(chirps))->month9
chirps.dat%>%
  as_tibble()%>%
  mutate(date=ymd(date))%>%
  filter(date<=ymd("2021-10-13"))%>%
  filter(date>=ymd("2021-10-13")-months(12))%>%
  group_by(Site)%>%
  summarise(twelvemonth_ppt=sum(chirps))->month12
left_join(month3,month6)%>%
  left_join(month9)%>%
  left_join(month12)%>%
  round(0)%>%
  write.csv("results/preceding_ppt.csv")


(chirps.dat%>%
    mutate(Site=factor(Site,labels=paste("Site:",1:10)))%>%
    filter(id==1)%>%
    as_tibble()%>%
    distinct()%>%
    mutate(yr=year(date),
           mn=month(date))%>%
    group_by(yr,Site)%>%
    mutate(cumppt=cumsum(chirps))%>%
    mutate(yd=yday(date))%>%
    #filter(yr==2019,Site==1)%>%
    ungroup()->tata_cd)%>%
  ggplot(aes(yd,cumppt))+
  theme_bw()+
  geom_hline(yintercept = c(100,200),linetype=2)+
  geom_line(aes(group=yr),colour="grey",alpha=0.5)+
  geom_line(data=(
    tata_cd%>%
      group_by(Site,yd)%>%
      summarise(cumppt=mean(cumppt))
  ),aes(yd,cumppt),colour="black")+
  #geom_line(data=filter(tata_cd,yr==2020,date<ymd('2021-10-15')),aes(yd,cumppt),colour="orange",size=2)+
  geom_line(data=filter(tata_cd,yr==2021,date<ymd('2021-10-15')),aes(yd,cumppt),colour="red",size=2)+
  facet_wrap(~Site,nrow=2)+
  ylab("Cumulative precipiation (mm)")+
  xlab("Day of the year")
ggsave(paste("results/CHIRPS_yday_",format(Sys.time(),"%Y%m%d-%Hh%Mm"),".jpg",sep=""),width=14,height=7)


forGG%>%
  ggplot(aes(Date,spi))+
  geom_line()+
  geom_point(aes(colour=spi))+
  facet_wrap(~Site,ncol=1)+
  geom_hline(yintercept = 0)

forGG%>%
  ggplot(aes(as.character(year),spi))+
  geom_boxplot()+
  facet_wrap(~Site)+
  geom_hline(yintercept = 0)



### End of figure creation


# Extra analyses of interest...


# Creating a sliding window
window <- 269
chirps.dat%>%
  filter(id==1)%>%
  as_tibble()%>%
  select(date,prcp=chirps,Site)%>%
  arrange(Site,date)%>%
  group_by(Site)%>%
  mutate(window_prcp=
           slide_dbl(prcp, 
                     ~sum(.x),
                     .before=window,
                     .complete=TRUE))%>%
  drop_na(window_prcp)%>%
  mutate(date=ymd(date))%>%
  mutate(start=date-window)%>%
  select(Site,start, end=date, window_prcp)%>%
  mutate(end_month = month(end),
         end_day = day(end),
         end_year = year(end)) %>%
  group_by(Site,end_month, end_day) %>%
  mutate(threshold_05 = quantile(window_prcp, prob = 0.05),
         threshold_10 = quantile(window_prcp, prob = 0.10),
         threshold_25 = quantile(window_prcp, prob = 0.25),
         threshold_50 = quantile(window_prcp, prob = 0.5)) %>%
  ungroup()%>%
  mutate(fake_date = ymd(glue("2020-{end_month}-{end_day}")))->
  drought_data

drought_data %>%
  select(Site,fake_date,window_prcp=threshold_05) %>%
  distinct()->
  drought_line

drought_data %>%
  select(Site,fake_date,window_prcp=threshold_50) %>%
  distinct()->
  fifty_line


drought_data%>%
  #filter(end_year==2012)%>%
  # mutate(is_drought_year = end_year == 2012,
  #        end_year = fct_reorder(factor(end_year), is_drought_year)) %>%
  ggplot(aes(x= fake_date, y=window_prcp)) +
  geom_line(show.legend = FALSE,colour="grey",aes(group = end_year))+
  geom_line(data = drought_line, aes(x = fake_date, y = window_prcp),
            color = "red",size=1.5,linetype=2)+
  geom_line(data = fifty_line, aes(x = fake_date, y = window_prcp),
            color = "darkgreen",size=1.5)+
  geom_line(data=
              (drought_data%>%
                 filter(end>ymd("2021-01-01"),end<ymd("2021-12-31"))),
            aes(x=fake_date,y=window_prcp),colour="blue",size=1.5)+
  facet_wrap(~Site,ncol=2)+
  scale_x_date(date_breaks = "2 months", date_labels = "%B") +
  labs(x = NULL,
       y = "Total precipitation over previous 270 days (mm)")

# Borrowed from wes_pallette() https://github.com/karthik/wesanderson/blob/master/R/colors.R
qCols <- rev(c("#3B9AB2", "#78B7C5", "#74A089", "#E1AF00", "#F21A00"))


(chirps.dat%>%
    mutate(Month=month(date),Year=year(date))%>%
    group_by(Site,Year,Month)%>%
    summarise(PPT=sum(chirps))%>%
    ungroup()%>%
    group_by(Site,Month)%>%
    mutate(Monthq10=quantile(PPT,0.1,na.rm=T),
           Monthq33=quantile(PPT,0.33,na.rm=T),
           Monthq66=quantile(PPT,0.66,na.rm=T),
           Monthq90=quantile(PPT,0.90,na.rm=T),  
           Monthq100=quantile(PPT,1.00,na.rm=T),
           Monthq50=median(PPT,na.rm=T),
           Monthqmean=mean(PPT,na.rm=T))%>%
    mutate(MonthQuantileNo=1+as.numeric(PPT>Monthq10)+as.numeric(PPT>Monthq33)+
             as.numeric(PPT>Monthq66)+as.numeric(PPT>Monthq90))%>%
    mutate(lab=ifelse(PPT<Monthq10,"*",""))->DATall)%>%
  ungroup()%>%
  group_by(Site,Month)%>%
  summarise(mean_monthly=mean(PPT,na.rm=T),
            sd_monthly=sd(PPT,na.rm=T))%>%
  ungroup()%>%
  mutate(Year1=2022,Year2=2021,Year3=2020,Year4=2019)%>%
  pivot_longer(cols=c("Year1","Year2","Year3","Year4"),values_to = "Year")%>%
  left_join(DATall)->
  tutu

tutu%>%
  ggplot()+
  geom_line(aes(Year,Monthq50,group=Month),col="darkgrey")+
  geom_segment(aes(x=Year,xend=Year,y=Monthq50,
                   yend=PPT,
                   colour=as.factor(MonthQuantileNo)),
               lwd=1.5)+
  facet_grid(Site~Month)+
  scale_colour_manual(values=qCols,
                      labels=c("0-10","10-33","33-66","66-90","90-100"),
                      na.translate = F)+
  theme(axis.text.x = element_text(angle=90,hjust=0.5),legend.position="top")+
  guides(col=guide_legend(title="Percentile"))+
  ylab("Monthly rainfall (mm)")

tutu%>%
  mutate(YMD=ceiling_date(ym(paste(Year,Month,sep="-")),unit="month")-1)%>%
  ggplot(aes(YMD,PPT,fill=as.factor(MonthQuantileNo)))+
  geom_col()+
  scale_fill_manual(values=qCols,labels=c("0-10","10-33","33-66","66-90","90-100"))+
  geom_text(aes(YMD,PPT,label=lab),size=12,col=qCols[1])+
  facet_grid(Year~Site)

geom_segment(aes(x=YMD,xend=YMD,y=q33,yend=q66))+
  geom_point(aes(x=YMD,y=q50))+
  #facet_wrap(~Station,ncol=1)+
  coord_cartesian(xlim=c(floor_date(today(),"month")-48*30,floor_date(today(),"month")-1),
                  ylim=c(0,80))+
  guides(fill=guide_legend(title="Percentile"))+
  scale_x_continuous(breaks=seq(ceiling_date(today()-48*31,"month"),ceiling_date(today(),"month")-1,by="months")-1)+
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position="top")+
  ylab("Monthly Precipitation (mm)")+
  xlab("")