setwd("J:/Documents/Arbete/Fortlöpande miljöanalys/Fenologi/Skoglig fenologi vid SLU")
library(dplyr)

selectedStation=105                                    #choose a station to process, see ID:s on the next statement

stationID<-c(102,103,104,105)
stationName<-c("Svartberget", "Siljansfors", "Asa", "Tönnersjöheden")

for (o in 1:4) {
  if (selectedStation==stationID[o]) break
}

species<-c("Bilberry","Lingonberry")
speciesID<-c(2440100,2440200)

xpre<-91:320
ypre<-rep(0,230)
minxpre<-min(xpre)
maxxpre<-max(xpre)
xtspre<-seq(0,2000,by=5)
ytspre<-rep(0,401)

st.err <- function(x) {
  sd(x)/sqrt(length(x))
}

# skip the following 2 steps if storing several stations mean+SE data in the same files is desirable, except for first station
Interpol.doy<-data.frame(Station=0,Species=0,doy=0,TotProcMean=0,FlowersTotProcMean=0,UnripeTotProcMean=0,RipeTotProcMean=0,TotProcSE=0,FlowersTotProcSE=0,UnripeTotProcSE=0,RipeTotProcSE=0)[FALSE, ]
Interpol.tsum<-data.frame(Station=0,Species=0,Tsum=0,TotProcMean=0,FlowersTotProcMean=0,UnripeTotProcMean=0,RipeTotProcMean=0,TotProcSE=0,FlowersTotProcSE=0,UnripeTotProcSE=0,RipeTotProcSE=0)[FALSE, ]

berrydata0<- read.csv("InterpolBerryData.csv", header=TRUE, sep=",", dec=".")
berrydatatsum0<- read.csv("InterpolTsumBerryData.csv", header=TRUE, sep=",", dec=".")

berrydata<-subset(berrydata0,berrydata0$Station==104)
berrydatatsum<-subset(berrydatatsum0,berrydatatsum0$Station==104)

specieslist<-unique(berrydata$Species)

plot.new()
par(mfrow = c(1,2))

for (l in 1:length(specieslist)) {
  for (n in 1:4) {
    if (speciesID[n]==specieslist[l]) break
  }
  testspecies<-subset(berrydata,berrydata$Species==specieslist[l])
  
  berrydatamean<-aggregate(testspecies[, 12:15], list(testspecies$doy), mean)
  colnames(berrydatamean)<-c("doy","TotProcMean","FlowersTotProcMean","UnripeTotProcMean","RipeTotProcMean")
  berrydatase<-aggregate(testspecies[, 12:15], list(testspecies$doy), st.err)
  colnames(berrydatase)<-c("doy","TotProcSE","FlowersTotProcSE","UnripeTotProcSE","RipeTotProcSE")
  berrydatameanse<-merge(berrydatamean,berrydatase,by="doy")
  berrydatameanse$Station<-selectedStation
  berrydatameanse$Species<-specieslist[l]
  berrydatameanse<-berrydatameanse %>% relocate(Station,Species)
  Interpol.doy<-rbind(Interpol.doy,berrydatameanse)
  
  yearlist<-unique(testspecies$Year)

  plot(xpre,ypre, main=paste(stationName[o], species[n]), xlab="Day from Jan 1", ylab = "% appearence",xlim=c(120,320),ylim=c(0,100),col="white")
  
  
  # for (j in 1:length(yearlist)) {
  #   testyear<-subset(testspecies,testspecies$Year==yearlist[j])
  #   
  # 
  #   lines(testyear$doy,testyear$FlowersTotProc, col = 6, lty="dashed", lwd=0.5)
  #   lines(testyear$doy,testyear$UnripeTotProc, col = 3, lty="dashed", lwd=0.5)
  #   lines(testyear$doy,testyear$RipeTotProc, col = 4, lty="dashed", lwd=0.5)
  # }
  
  lines(berrydatameanse$doy,berrydatameanse$FlowersTotProcMean, col = 6, lwd=3)
  lines(berrydatameanse$doy,berrydatameanse$FlowersTotProcMean-berrydatameanse$FlowersTotProcSE, col = 6, lty="dashed", lwd=2)
  lines(berrydatameanse$doy,berrydatameanse$FlowersTotProcMean+berrydatameanse$FlowersTotProcSE, col = 6, lty="dashed", lwd=2)
  lines(berrydatameanse$doy,berrydatameanse$UnripeTotProcMean, col = 3, lwd=3)
  lines(berrydatameanse$doy,berrydatameanse$UnripeTotProcMean-berrydatameanse$UnripeTotProcSE, col = 3, lty="dashed", lwd=2)
  lines(berrydatameanse$doy,berrydatameanse$UnripeTotProcMean+berrydatameanse$UnripeTotProcSE, col = 3, lty="dashed", lwd=2)
  lines(berrydatameanse$doy,berrydatameanse$RipeTotProcMean, col = 4, lwd=3)
  lines(berrydatameanse$doy,berrydatameanse$RipeTotProcMean-berrydatameanse$RipeTotProcSE, col = 4, lty="dashed", lwd=2)
  lines(berrydatameanse$doy,berrydatameanse$RipeTotProcMean+berrydatameanse$RipeTotProcSE, col = 4, lty="dashed", lwd=2)
  lines(berrydatameanse$doy,berrydatameanse$TotProcMean, col = 1, lwd=2)
  lines(berrydatameanse$doy,berrydatameanse$TotProcMean-berrydatameanse$TotProcSE, col = 1, lty="dashed", lwd=2)
  lines(berrydatameanse$doy,berrydatameanse$TotProcMean+berrydatameanse$TotProcSE, col = 1, lty="dashed", lwd=2)
  
  
  testspeciestsum<-subset(berrydatatsum,berrydatatsum$Species==specieslist[l])

  berrydatatsummean<-aggregate(testspeciestsum[, 12:15], list(testspeciestsum$Tempsum), mean)
  colnames(berrydatatsummean)<-c("Tsum","TotProcMean","FlowersTotProcMean","UnripeTotProcMean","RipeTotProcMean")
  berrydatatsumse<-aggregate(testspeciestsum[, 12:15], list(testspeciestsum$Tempsum), st.err)
  colnames(berrydatatsumse)<-c("Tsum","TotProcSE","FlowersTotProcSE","UnripeTotProcSE","RipeTotProcSE")
  berrydatatsummeanse<-merge(berrydatatsummean,berrydatatsumse,by="Tsum")
  berrydatatsummeanse$Station<-selectedStation
  berrydatatsummeanse$Species<-specieslist[l]
  berrydatatsummeanse<-berrydatatsummeanse %>% relocate(Station,Species)
  Interpol.tsum<-rbind(Interpol.tsum,berrydatatsummeanse)
  
  plot(xtspre,ytspre, main=paste(stationName[o], species[n]), xlab="Temperature sum (+5°C threshold)", ylab = "% appearence",xlim=c(0,1500),ylim=c(0,100),col="white")
  
  
  # for (j in 1:length(yearlist)) {
  #   testyear<-subset(testspeciestsum,testspeciestsum$Year==yearlist[j])
  #   
  #   
  #   lines(testyear$Tempsum,testyear$FlowersTotProc, col = 6, lty="dashed", lwd=0.5)
  #   lines(testyear$Tempsum,testyear$UnripeTotProc, col = 3, lty="dashed", lwd=0.5)
  #   lines(testyear$Tempsum,testyear$RipeTotProc, col = 4, lty="dashed", lwd=0.5)
  # }
  
  lines(berrydatatsummeanse$Tsum,berrydatatsummeanse$FlowersTotProcMean, col = 6, lwd=3)
  lines(berrydatatsummeanse$Tsum,berrydatatsummeanse$FlowersTotProcMean-berrydatatsummeanse$FlowersTotProcSE, col = 6, lty="dashed", lwd=2)
  lines(berrydatatsummeanse$Tsum,berrydatatsummeanse$FlowersTotProcMean+berrydatatsummeanse$FlowersTotProcSE, col = 6, lty="dashed", lwd=2)
  lines(berrydatatsummeanse$Tsum,berrydatatsummeanse$UnripeTotProcMean, col = 3, lwd=3)
  lines(berrydatatsummeanse$Tsum,berrydatatsummeanse$UnripeTotProcMean-berrydatatsummeanse$UnripeTotProcSE, col = 3, lty="dashed", lwd=2)
  lines(berrydatatsummeanse$Tsum,berrydatatsummeanse$UnripeTotProcMean+berrydatatsummeanse$UnripeTotProcSE, col = 3, lty="dashed", lwd=2)
  lines(berrydatatsummeanse$Tsum,berrydatatsummeanse$RipeTotProcMean, col = 4, lwd=3)
  lines(berrydatatsummeanse$Tsum,berrydatatsummeanse$RipeTotProcMean-berrydatatsummeanse$RipeTotProcSE, col = 4, lty="dashed", lwd=2)
  lines(berrydatatsummeanse$Tsum,berrydatatsummeanse$RipeTotProcMean+berrydatatsummeanse$RipeTotProcSE, col = 4, lty="dashed", lwd=2)
  lines(berrydatatsummeanse$Tsum,berrydatatsummeanse$TotProcMean, col = 1, lwd=2)
  lines(berrydatatsummeanse$Tsum,berrydatatsummeanse$TotProcMean-berrydatatsummeanse$TotProcSE, col = 1, lty="dashed", lwd=2)
  lines(berrydatatsummeanse$Tsum,berrydatatsummeanse$TotProcMean+berrydatatsummeanse$TotProcSE, col = 1, lty="dashed", lwd=2)
  
  
  
  legend("topright", 
         legend = c("Flowers", "Unripe berries","Ripe berries","Totals"), 
         lwd=3,
         col = c(6,3,4,1), 
         lty= c(1,1,1,1),
         bty = "n", 
         text.col = "black", 
         horiz = F )
}

filename=paste0("BerryData_", stationName[o], "_DoyMeanSE.csv")
write.csv(Interpol.doy, filename, row.names = FALSE)              #storage of the average values per doy, for the selected station
filename=paste0("BerryData_", stationName[o], "_TsumMeanSE.csv")
write.csv(Interpol.tsum, filename, row.names = FALSE)             #storage of the average values per tsum, for the selected station