library(dplyr)

setwd("J:/Documents/Arbete/Fortlöpande miljöanalys/Fenologi/Skoglig fenologi vid SLU")    #define working directory

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

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

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

species<-c("Downy birch", "Silver birch", "Norway spruce","Scots pine")
speciesID<-c(1180100,1180200,2030100,2040300)

xpre<-91:250
ypre<-rep(0,160)
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,ShootProcMean=0,ShootProcSE=0)[FALSE, ]
Interpol.tsum<-data.frame(Station=0,Species=0,Tsum=0,ShootProcMean=0,ShootProcSE=0)[FALSE, ]

treedata0<- read.csv("InterpolTreeData.csv", header=TRUE, sep=",", dec=".")
treetsumdata0<- read.csv("InterpolTsumTreeData.csv", header=TRUE, sep=",", dec=".")

treedata<-subset(treedata0,treedata0$Station==selectedStation)
treetsumdata<-subset(treetsumdata0,treetsumdata0$Station==selectedStation)

specieslist<-unique(treedata$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(treedata,treedata$Species==specieslist[l])
  testspeciestsum<-subset(treetsumdata,treetsumdata$Species==specieslist[l])
  
  treedatamean<-aggregate(testspecies$ShootProc, list(testspecies$doy), mean)
  colnames(treedatamean)<-c("doy","ShootProcMean")
  treedatase<-aggregate(testspecies$ShootProc, list(testspecies$doy),st.err)
  colnames(treedatase)<-c("doy","ShootProcSE")
  treedatameanse<-merge(treedatamean,treedatase,by="doy")
  treedatameanse$Station<-selectedStation
  treedatameanse$Species<-specieslist[l]
  treedatameanse<-treedatameanse %>% relocate(Station,Species)
  Interpol.doy<-rbind(Interpol.doy,treedatameanse)
  
  
  treetsumdatamean<-aggregate(testspeciestsum$ShootDevProc, list(testspeciestsum$Tempsum), mean)
  colnames(treetsumdatamean)<-c("Tsum","ShootProcMean")
  treetsumdatase<-aggregate(testspeciestsum$ShootDevProc, list(testspeciestsum$Tempsum),st.err)
  colnames(treetsumdatase)<-c("Tsum","ShootProcSE")
  treetsumdatameanse<-merge(treetsumdatamean,treetsumdatase,by="Tsum")
  treetsumdatameanse$Station<-selectedStation
  treetsumdatameanse$Species<-specieslist[l]
  treetsumdatameanse<-treetsumdatameanse %>% relocate(Station,Species)
  Interpol.tsum<-rbind(Interpol.tsum,treetsumdatameanse)
  
  yearlist<-unique(testspecies$Year)
  
  plot(xpre,ypre, main=paste(stationName[o], species[n]), xlab="Day from Jan 1", ylab = "% shoot length",xlim=c(80,220),ylim=c(0,100),col="white")
  
  for (j in 1:length(yearlist)) {
    testyear<-subset(testspecies,testspecies$Year==yearlist[j])
    
  
    lines(testyear$doy,testyear$ShootProc, col = 6, lty="dashed", lwd=0.5)
  }
  
  lines(treedatameanse$doy,treedatameanse$ShootProcMean, col = 1, lwd=3)
  lines(treedatameanse$doy,treedatameanse$ShootProcMean-treedatameanse$ShootProcSE, col = 1, lty="dashed", lwd=2)
  lines(treedatameanse$doy,treedatameanse$ShootProcMean+treedatameanse$ShootProcSE, col = 1, lty="dashed", lwd=2)

  
  
  plot(xtspre,ytspre, main=paste(stationName[o], species[n]), xlab="Temperature sum (+5°C thresohold)", ylab = "% shoot length",xlim=c(0,1000), ylim=c(0,100),col="white")
  
  for (j in 1:length(yearlist)) {
    testyear<-subset(testspeciestsum,testspeciestsum$Year==yearlist[j])
    
    
    lines(testyear$Tempsum,testyear$ShootDevProc, col = 6, lty="dashed", lwd=0.5)
  }
  
  lines(treetsumdatameanse$Tsum,treetsumdatameanse$ShootProcMean, col = 1, lwd=3)
  lines(treetsumdatameanse$Tsum,treetsumdatameanse$ShootProcMean-treetsumdatameanse$ShootProcSE, col = 1, lty="dashed", lwd=2)
  lines(treetsumdatameanse$Tsum,treetsumdatameanse$ShootProcMean+treetsumdatameanse$ShootProcSE, col = 1, lty="dashed", lwd=2)
  
}

filename=paste0("TreeData_", 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("TreeData_", stationName[o], "_TsumMeanSE.csv")
write.csv(Interpol.tsum, filename, row.names = FALSE)             #storage of the average values per tsum, for the selected station