##

#   GDD again but separeting years
#
#####

library(data.table)
library(dplyr)
library(lme4)
library(nlme)
library(ggplot2)
library(ggthemes)
library(ggpubr)
library(ggExtra)
library(imputeTS)
library(tidyr)
library(pollen)


setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/2.ObsInfo")

## load unique sites (109)
load('MothSites109.RDATA')

Moth.year = data.frame(SiteID = rep(Moth.m2$SiteID, 29), 
                       YEUREF= rep(Moth.m2$YEUREF, 29), 
                       XEUREF = rep(Moth.m2$XEUREF, 29),  
                       Year = rep(1998:2021, each=109))
head(Moth.year)



# Environmental data
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/1.EnvironmentalInfo")
Clim <- read.csv("RecbaseMoths_weather_CORRECTED.csv", header=T)
head(Clim)
dim(Clim)

Clim = Clim[, c(2:6, 8:10)]

Climnew <- read.csv("RecbaseMoths_weather_19_21.csv", header=T)
head(Climnew)
Climnew$Site = paste(Climnew$Latitude, Climnew$Longitude)

SiteMoth = Moth.year
SiteMoth$Site = paste(SiteMoth$YEUREF, SiteMoth$XEUREF)

Climnew2 = merge(Climnew, SiteMoth, by = c('Site'), all.x=TRUE)
head(Climnew2)
Climnewsub = Climnew2[, c(12, 4:7, 10,9, 8)]
head(Climnewsub)
colnames(Climnewsub) = colnames(Clim)

ClimAll = rbind(Clim, Climnewsub)


ClimAll = ClimAll[ClimAll$site %in% unique(Moth.year$SiteID),]
ClimAll = ClimAll[order(ClimAll$site, ClimAll$year, ClimAll$month, ClimAll$day),]
ClimAll = ClimAll[ClimAll$year < 2021 ,]


#-.......................... Climatic variable preparation
#-..........................
# calculating GDD
## -------------------------------------------------------------------------

Year= 1998:2020
Site = unique(ClimAll$site)

# calculating GDD

# winter time only
datesWinter = data.frame(Enddate = paste(Year,4,1, sep="-"),        # from when to end
                         Startdate = paste(Year-1,10,15, sep="-"))


datesWinter$Startdate <- as.Date(datesWinter$Startdate)     
datesWinter$Enddate   <- as.Date(datesWinter$Enddate) 

ClimAll = ClimAll %>%
  mutate(date = paste(year, month, day, sep="-"))  # 
head(ClimAll)

ClimAll$date   <- as.Date(ClimAll$date)

P.Winter.k =  vector(mode='list', length(Site))
for (k in 1:length(Site)) {
  clim.store.tmp = ClimAll[ClimAll$site == Site[k],]
  Moth.m2.i = Moth.year[Moth.year$SiteID == Site[k] ,]
  Moth.m2.i = cbind(Moth.m2.i, datesWinter)
  
  for(j in 1:nrow(Moth.m2.i)){
    P.Winter.ksite = clim.store.tmp[clim.store.tmp$date >= Moth.m2.i$Startdate[j] & clim.store.tmp$date <= Moth.m2.i$Enddate[j],]
    #
    P.Winter.k[[k]][[j]] = P.Winter.ksite
  }
}

head(P.Winter.k[[1]])
length(P.Winter.k)


#-..........................
## chilling DD calculation
## -------------------------------------------------------------------------

#sum.CDD <- function(x) sum(5-x[x < 5], na.rm = TRUE) #define function to calculate GDD5

gdd_5 = gdd0 =gdd5 = gdd10 = #gdd_5.2 = gdd0.2 =gdd5.2 = gdd10.2 = 
  cdd_5 = cdd0 = cdd5 = cdd10 = vector(mode='list', 109)

for (i in 1:length(Site)) {
  Clim.j = P.Winter.k[[i]]
  
  for (j in 1:length(Clim.j)) {
    Clim.ij = Clim.j[[j]]
    
    Clim.ij$avg = (Clim.ij$tempMinDay.tmin_ + Clim.ij$tempMaxDay.tmax_)/2
    
	## chilling degree days
    cdd5[[i]][[j]] = sum(5 - Clim.ij$avg[Clim.ij$avg < 5], na.rm = TRUE)
    
    cdd10[[i]][[j]] = sum(10 - Clim.ij$avg[Clim.ij$avg < 10], na.rm = TRUE)
    
    
    cdd0[[i]][[j]] = sum(0 - Clim.ij$avg[Clim.ij$avg < 0], na.rm = TRUE)
    
    
    cdd_5[[i]][[j]] = sum(-5 - Clim.ij$avg[Clim.ij$avg < -5], na.rm = TRUE)
    
    
	
	## growing degree days
    gdd5[[i]][[j]] = sum(Clim.ij$avg[Clim.ij$avg > 5] - 5, na.rm = TRUE)
    
    
    gdd10[[i]][[j]] = sum(Clim.ij$avg[Clim.ij$avg > 10] - 10, na.rm = TRUE)
    
     
    gdd0[[i]][[j]] = sum(Clim.ij$avg[Clim.ij$avg > 0] - 0, na.rm = TRUE)
    
    
    gdd_5[[i]][[j]] = sum(Clim.ij$avg[Clim.ij$avg > -5] - (-5), na.rm = TRUE)
    
    
  }
  
}

GDD_5 = unlist(gdd_5)
GDD0 = unlist(gdd0)
GDD5 = unlist(gdd5)
GDD10 = unlist(gdd10)

GDD_5.2 = unlist(gdd_5.2)
GDD0.2 = unlist(gdd0.2)
GDD5.2 = unlist(gdd5.2)
GDD10.2 = unlist(gdd10.2)

CDD_5 = unlist(cdd_5)
CDD0 = unlist(cdd0)
CDD5 = unlist(cdd5)
CDD10 = unlist(cdd10)


SY = data.frame(Site = rep(Site, each=length(Year)), 
                YEUREF= rep(Moth.m2$YEUREF, 29), 
                XEUREF = rep(Moth.m2$XEUREF, 29),
                Year = rep(Year, length(Site)))

Sgdd.wint = cbind(SY, CDD_5, CDD0, CDD5, CDD10, GDD_5, GDD0, GDD5, GDD10)

colnames(Sgdd.wint)[c(5:12)] = c('CDD_5w', 'CDD0w', 'CDD5w', 'CDD10w',
                        'GDD_5w', 'GDD0w', 'GDD5w', 'GDD10w')

head(Sgdd.wint)
save(Sgdd.wint, file = 'gcddwintFMI0806.RDATA')


## look at the data gdd

Sgdd.wint = as.data.frame(Sgdd.wint)
summary(Sgdd.wint)
boxplot(Sgdd.wint$CDD5w ~ Sgdd.wint$Year)




###   Snow
## -------------------------------------------------------------------------


Snow = vector(mode='list', length(Site))
for (i in 1:length(Site)) {
  Clim.j = P.Winter.k[[i]]
  
  for (j in 1:length(Clim.j)) {
    Clim.ij = Clim.j[[j]]
    
    Snow[[i]][[j]] = mean(Clim.ij$snowDay.snow_, na.rm=T)
    
  }
}

SnowAv = unlist(Snow)


#SY = data.frame(Site = rep(Site, each=29), Year = rep(Year, 109))


SnowAv = cbind(SY, SnowAv)
colnames(SnowAv)[5] = c('snow')

save(SnowAv, file = 'SnowAV0806.RDATA')
boxplot(SnowAv$snow~ Sgdd.wint$Year)


### Summer --------------------------------------------------------------------------
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## ## ---- Growing degree days


Year=1998:2020

datesSummer = data.frame(Startdate = paste(Year,04,01, sep="-"),
                         Enddate = paste(Year, 10, 15, sep="-"))
datesSummer$Startdate <- as.Date(datesSummer$Startdate)     
datesSummer$Enddate   <- as.Date(datesSummer$Enddate) 


# calculating GDD ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Year2 = 1998:2020

P.Summer.k = vector(mode='list', length(Site))
for (k in 1:length(Site)) {
  clim.store.tmp = ClimAll[ClimAll$site == Site[k],]
  Moth.m2.i = Moth.year[Moth.year$SiteID == Site[k] ,]
  Moth.m2.i = cbind(Moth.m2.i, datesSummer)
  
  for(j in 1:length(Year2)){
    P.Summer.ksite = clim.store.tmp[clim.store.tmp$date >= datesSummer$Startdate[j] & clim.store.tmp$date <= datesSummer$Enddate[j],]
    
    P.Summer.k[[k]][[j]] = P.Summer.ksite
  }
}


length(P.Summer.k[[1]]) # list element for each entry of Moth.m (Site/Year/Month/species)

gdd_5 = gdd0 =gdd5 = gdd10 = vector(mode='list', length(Site))
for (j in 1:length(Site)) {
  Clim.j = P.Summer.k[[j]]
  
  for (i in 1:length(Year)) {
    Climsum.ij = Clim.j[[i]]
    
    Climsum.ij$avg = (Climsum.ij$tempMinDay.tmin_ + Climsum.ij$tempMaxDay.tmax_)/2
    
    gdd5[[j]][[i]] = sum(Climsum.ij$avg[Climsum.ij$avg > 5] - 5, na.rm = TRUE)
    
    
    gdd10[[j]][[i]] = sum(Climsum.ij$avg[Climsum.ij$avg > 10] - 10, na.rm = TRUE)
    
    
    gdd0[[j]][[i]] = sum(Climsum.ij$avg[Climsum.ij$avg > 0] - 0, na.rm = TRUE)
    
    
    gdd_5[[j]][[i]] = sum(Climsum.ij$avg[Climsum.ij$avg > -5] - (-5), na.rm = TRUE)
    
  }
  
}
GDD_5 = unlist(gdd_5)
GDD0 = unlist(gdd0)
GDD5 = unlist(gdd5)
GDD10 = unlist(gdd10)

Sgdd.Sum = data.frame(SY, GDD_5s = GDD_5, 
                      GDD0s = GDD0, GDD5s = GDD5, GDD10s = GDD10)

head(Sgdd.Sum)
boxplot(Sgdd.Sum$GDD5s~Sgdd.Sum$Year)

save(Sgdd.Sum, file = 'Sgdd-Summer080623.RDATA')



## ---- Precipitation
## -------------------------------------------------------------------------


P.Summer.k = rep(list(list()), length(Site))
for (k in 1:length(Site)) {
  clim.store.tmp = clim.new[clim.new$site == Site[k],]
  
  for(j in 1:length(Year2)){
    P.Summer.ksite = clim.store.tmp[clim.store.tmp$date >= datesSummer$Startdate[j] & clim.store.tmp$date <= datesSummer$Enddate[j],]
    #
    P.Summer.k[[k]][[j]] = P.Summer.ksite
  }
}


length(P.Summer.k)

head(P.Summer.k[[1]])


# The list element contains the daily precipitation from beginning of Summer till sampling date
Psum_site = list()
for (k in 1:length(Site)){
  Psum_site[[k]] = lapply(P.Summer.k[[k]], function(x) {sum(x$precDay.rrday_, na.rm=T)})
}


Psum_site = vector(mode='list', length(Site))
for (j in 1:length(Site)) {
  Clim.j = P.Summer.k[[j]]
  
  for (i in 1:length(Year)) {
    Clim.ij = Clim.j[[i]]
    Psum_site[[j]][[i]] = sum(Clim.ij$precDay.rrday_, na.rm=T)
    
  }
  
}

Psum = unlist(Psum_site)
summary(Psum)
length(Psum)
head(Psum)

Pcumsum = data.frame(SY, precSum = Psum)
boxplot(Pcumsum$precSum~Pcumsum$Year)

save(Pcumsum, file = 'Pcs080623.RDATA')



# combine all clim var Info ------------------------------------------------------
load('SnowAV0806.RDATA')
load('gcddwintFMI0806.RDATA')
load('Sgdd-Summer080623.RDATA')
load('Pcs080623.RDATA')

## all have the same orders
head(SnowAv)
head(Pcumsum)
head(Sgdd.Sum)
head(Sgdd.wint)

dim(SnowAv)
dim(Pcumsum)
dim(Sgdd.Sum)
dim(Sgdd.wint)

ClimSite = cbind(Sgdd.wint[, 1:12], Sgdd.Sum$GDD_5s, Sgdd.Sum$GDD0s, Sgdd.Sum$GDD5s, Sgdd.Sum$GDD10s,
                 Pcumsum$precSum, SnowAv$snow)

head(ClimSite)
#ClimSite = ClimSite[, -c(3:4)]

colnames(ClimSite)[13:18] = c('GDD_5s', 'GDD0s', 'GDD5s', 'GDD10s', 'precip', 'snow')

save(ClimSite, file = 'ClimSite080623.RDATA')


