# This script compares site data to Fluxnet and Modis data #

dev.off()
rm(list=ls(all=TRUE))

#import user-defined functions and graphics parameters
source("D:/R/Github/graphics_param.R")
source("D:/R/Github/functions.R")

library(tables)
library(ggplot2)
library(gridExtra)
library(ggpubr)
library(ncdf4)
library(lubridate)

path.FIT.data<- 'D:/Data/LPJmL-Fit/Validation/'

runs<-c("run305","run302","run300","run299","run303","run301","run304")
runs<-c("seitseminen","kalkalpen","bialowieza","hainich","laegeren","dundo","peneda")

site.names<-c("Seitseminen","Kalkalpen","Bialowieza","Hainich","Laegeren","Dundo","PenedaGeres")
plot.names<-c("Seitseminen, Tmean = 3.5°C","Kalkalpen, Tmean = 5.6°C","Bialowieza, Tmean = 7.2°C","Hainich, Tmean = 7.8°C","Lägern, Tmean = 8.7°C","Dundo, Tmean = 11.5°C","Peneda-Gerês, Tmean = 11.5°C")

site.names<-c("Seitseminen","Kalkalpen","Bialowieza","Hainich","Laegeren","Dundo")
plot.names<-c("Seitseminen, Tmean = 3.5°C","Kalkalpen, Tmean = 5.6°C","Bialowieza, Tmean = 7.2°C","Hainich, Tmean = 7.8°C","Lägern, Tmean = 8.7°C","Dundo, Tmean = 11.5°C")

site.names<-c("Seitseminen","Kalkalpen","Bialowieza","Hainich","Laegeren","Dundo")
plot.names<-c("NP Seitseminen","NP Kalkalpen","NP Bialowieza","NP Hainich","Nature Reserve Lägern","Forest Reserve Dundo")

error.names<-c("ME","MAE","RMSE","MPE","MAPE","RMSPE","cor","NMSE","NME")
error.names<-c("ME","cor","NMSE")

years<-2004:2013
theme_set(theme_bw(base_size = 30))
theme_update(axis.text=element_text(size=28))
theme_update(plot.title=element_text(size=33,hjust = 0.5,face="bold"))
linesize<-1
color.FIT<-"coral1"
color.MODIS<-"dodgerblue3"
color.FLUX<-"gold2"
alpha<-0.2



### GPP ###
years<-2004:2013
months<-1:12
ls.gpp.annual<-ls.gpp.monthly<-ls.gpp.monthly.long<-ls.gpp.monthly_nolegend<-ls.gpp.monthly.long_nolegend<-list()
array.error<-array(NA,c(length(site.names),length(error.names)*3))
for(isite in 1:length(site.names)){
  print(site.names[isite])
  FIT.data<-get.FIT.GPP(runs[isite])
  MODIS.data<-get.MODIS.GPP(site.names[isite])
  
  for(ierror in 1:length(error.names)){
    array.error[isite,ierror]<-cal.error(error.name = error.names[ierror],
                                         sim = FIT.data$vec_monthly_mean,obs=MODIS.data$vec_monthly_mean)
    array.error[isite,ierror+3]<-cal.error(error.name = error.names[ierror],
                                           sim = FIT.data$mean,obs=MODIS.data$mean)
    array.error[isite,ierror+6]<-cal.error(error.name = error.names[ierror],
                                           sim = FIT.data$annual_sum,obs=MODIS.data$annual_sum)
  }
  
  if(site.names[isite]=="Laegeren" | site.names[isite]=="Hainich"){
    

    FLUX.data<-get.FLUX(site.names[isite])
    
    if(site.names[isite]=="Hainich"){
      for(ierror in 1:length(error.names)){
        print(cal.error(error.name = error.names[ierror],
                  sim = FIT.data$vec_monthly_mean[1:48],obs=FLUX.data$gpp_vec_mean[49:96]))
        print(cal.error(error.name = error.names[ierror],
                  sim = FIT.data$mean,obs=FLUX.data$mean_gpp))
        print(cal.error(error.name = error.names[ierror],
                  sim = FIT.data$annual_sum[1:4],obs=FLUX.data$sum_gpp[5:8]))
      }
    }
    
    if(site.names[isite]=="Laegeren"){
      for(ierror in 1:length(error.names)){
        print(cal.error(error.name = error.names[ierror],
                        sim = FIT.data$vec_monthly_mean[1:84],obs=FLUX.data$gpp_vec_mean[1:84]))
        print(cal.error(error.name = error.names[ierror],
                        sim = FIT.data$mean,obs=FLUX.data$mean_gpp))
        print(cal.error(error.name = error.names[ierror],
                        sim = FIT.data$annual_sum[1:7],obs=FLUX.data$sum_gpp[1:7]))
      }
    }
    
    
    df<-data.frame(data=c(FIT.data$annual_sum,MODIS.data$annual_sum,FLUX.data$sum_gpp),
                   sd_high=c(FIT.data$annual_sum+FIT.data$sd_annual_sum,MODIS.data$annual_sum+MODIS.data$sd_annual_sum,FLUX.data$sum_gpp+0),
                   sd_low=c(FIT.data$annual_sum-FIT.data$sd_annual_sum,MODIS.data$annual_sum-MODIS.data$sd_annual_sum,FLUX.data$sum_gpp-0),
                   year=c(rep(years,2),FLUX.data$years),
                   type=c(rep(c("LPJmL-FIT","MODIS"),each=length(years)),rep("Fluxnet",length(FLUX.data$years))),
                   col=c(rep(c(color.FIT, color.MODIS),each=length(years)),rep(color.FLUX,length(FLUX.data$years))))
    p<-ggplot(df)
    p<-p+geom_ribbon(data=df,alpha=alpha,aes(x=year,ymin=sd_low,ymax=sd_high,fill=type,linetype=NA))
    p<-p+geom_line(aes(x=year,y=data,colour=type),size=linesize)
    cols <- c("LPJmL-FIT" = color.FIT, "MODIS" = color.MODIS, "Fluxnet" = color.FLUX)
    p<-p+scale_color_manual(name='',breaks = c("LPJmL-FIT", "MODIS","Fluxnet"),values=cols)
    p<-p+scale_fill_manual(name='',breaks = c("LPJmL-FIT", "MODIS","Fluxnet"),values=cols,guide = FALSE)
    p<-p+ylab("GPP [gC/m²]")+xlab("Year")
    p<-p+scale_x_continuous(breaks = pretty(df$year, n = 5))
    p<-p+xlim(2004,2013)
    p<-p+ggtitle(plot.names[isite])
    p
  }else{
    df<-data.frame(data=c(FIT.data$annual_sum,MODIS.data$annual_sum),
                   sd_high=c(FIT.data$annual_sum+FIT.data$sd_annual_sum,MODIS.data$annual_sum+MODIS.data$sd_annual_sum),
                   sd_low=c(FIT.data$annual_sum-FIT.data$sd_annual_sum,MODIS.data$annual_sum-MODIS.data$sd_annual_sum),
                   year=rep(years,2),
                   type=rep(c("LPJmL-FIT","MODIS"),each=length(years)),
                   col=rep(c(color.FIT, color.MODIS),each=length(years)))
    
    
    p<-ggplot(df)
    p<-p+geom_ribbon(data=df,alpha=alpha,aes(x=year,ymin=sd_low,ymax=sd_high,fill=type,linetype=NA))
    p<-p+geom_line(aes(x=year,y=data,colour=type),size=linesize)
    p<-p+scale_color_manual(name='',breaks = c("LPJmL-FIT", "MODIS"),values=c(color.FIT, color.MODIS))
    p<-p+scale_fill_manual(name='',breaks = c("LPJmL-FIT", "MODIS"),values=c(color.FIT, color.MODIS),guide = FALSE)
    p<-p+ylab("GPP [gC/m²]")+xlab("Year")
    p<-p+ggtitle(plot.names[isite])
    p<-p+scale_x_continuous(breaks = pretty(df$year, n = 5))
    p<-p+xlim(2004,2013)
    p
  }
  
  ls.gpp.annual[[isite]] <- p
  
  if(site.names[isite]=="Laegeren" | site.names[isite]=="Hainich"){
    FLUX.data<-get.FLUX(site.names[isite])
    df<-data.frame(data=c(FIT.data$mean,MODIS.data$mean,FLUX.data$mean_gpp),
                  sd_high=c(FIT.data$mean+FIT.data$sd,MODIS.data$mean+MODIS.data$sd,FLUX.data$mean_gpp+FLUX.data$sd_gpp),
                  sd_low=c(FIT.data$mean-FIT.data$sd,MODIS.data$mean-MODIS.data$sd,FLUX.data$mean_gpp-FLUX.data$sd_gpp),
                  year=rep(months,3),
                  type=c(rep(c("LPJmL-FIT","MODIS"),each=12),rep("Fluxnet",12)),
                  col=c(rep(c(color.FIT, color.MODIS),each=12),rep(color.FLUX,12)))
  
    p<-ggplot(df)
    p<-p+geom_ribbon(data=df,alpha=alpha,aes(x=year,ymin=sd_low,ymax=sd_high,fill=type,linetype=NA))
    p<-p+geom_line(aes(x=year,y=data,colour=type),size=linesize)
    cols <- c("LPJmL-FIT" = color.FIT, "MODIS" = color.MODIS, "Fluxnet" = color.FLUX)
    p<-p+scale_color_manual(name='',breaks = c("LPJmL-FIT", "MODIS","Fluxnet"),values=cols)
    p<-p+scale_fill_manual(name='',breaks = c("LPJmL-FIT", "MODIS","Fluxnet"),values=cols,guide = FALSE)
    p<-p+ylab("GPP [gC/m²]")+xlab("Month")
    p<-p+ggtitle(plot.names[isite])
    p<-p+scale_x_continuous(breaks = pretty(df$year, n = 5))
    p
  }else{
    df<-data.frame(data=c(FIT.data$mean,MODIS.data$mean),
                   sd_high=c(FIT.data$mean+FIT.data$sd,MODIS.data$mean+MODIS.data$sd),
                   sd_low=c(FIT.data$mean-FIT.data$sd,MODIS.data$mean-MODIS.data$sd),
                   year=rep(months,2),
                   type=rep(c("LPJmL-FIT","MODIS"),each=length(months)),
                   col=rep(c(color.FIT, color.MODIS),each=length(months)))
    
    
    p<-ggplot(df)
    p<-p+geom_ribbon(data=df,alpha=alpha,aes(x=year,ymin=sd_low,ymax=sd_high,fill=type,linetype=NA))
    p<-p+geom_line(aes(x=year,y=data,colour=type),size=linesize)
    p<-p+scale_color_manual(name='',breaks = c("LPJmL-FIT", "MODIS"),values=c(color.FIT, color.MODIS))
    p<-p+scale_fill_manual(name='',breaks = c("LPJmL-FIT", "MODIS"),values=c(color.FIT, color.MODIS),guide = FALSE)
    p<-p+ylab("GPP [gC/m²]")+xlab("Month")
    p<-p+ggtitle(plot.names[isite])
    p<-p+scale_x_continuous(breaks = pretty(df$year, n = 5))
    p
  }
  
  
  ls.gpp.monthly[[isite]] <- p
  ls.gpp.monthly_nolegend[[isite]] <- p+theme(legend.position="none")
  
  if(site.names[isite]=="Laegeren" | site.names[isite]=="Hainich"){
    df<-data.frame(data=c(FIT.data$vec_monthly_mean,MODIS.data$vec_monthly_mean,FLUX.data$gpp_vec_mean),
                   sd_high=c(FIT.data$vec_monthly_mean+FIT.data$vec_monthly_sd,MODIS.data$vec_monthly_mean+MODIS.data$vec_monthly_sd,FLUX.data$gpp_vec_mean+0),
                   sd_low=c(FIT.data$vec_monthly_mean-FIT.data$vec_monthly_sd,MODIS.data$vec_monthly_mean-MODIS.data$vec_monthly_sd,FLUX.data$gpp_vec_mean-0),
                   year=c(rep(rep(years,each=12)+rep((1:12)/12-1/24,length(years)),2),rep(FLUX.data$years,each=12)+rep((1:12)/12-1/24,length(FLUX.data$years))),
                   type=c(rep(c("LPJmL-FIT","MODIS"),each=12*length(years)),rep("Fluxnet",12*length(FLUX.data$years))),
                   col=c(rep(c(color.FIT, color.MODIS),each=12*length(years)),rep(color.FLUX,12*length(FLUX.data$years))))
    
    p<-ggplot(df)
    p<-p+geom_ribbon(data=df,alpha=alpha,aes(x=year,ymin=sd_low,ymax=sd_high,fill=type,linetype=NA))
    p<-p+geom_line(aes(x=year,y=data,colour=type),size=linesize)
    cols <- c("LPJmL-FIT" = color.FIT, "MODIS" = color.MODIS, "Fluxnet" = color.FLUX)
    p<-p+scale_color_manual(name='',breaks = c("LPJmL-FIT", "MODIS","Fluxnet"),values=cols)
    p<-p+scale_fill_manual(name='',breaks = c("LPJmL-FIT", "MODIS","Fluxnet"),values=cols,guide = FALSE)
    p<-p+ylab("GPP [gC/m²]")+xlab("Year")
    p<-p+ggtitle(plot.names[isite])
    p<-p+scale_x_continuous(breaks = pretty(df$year, n = 5))
    p<-p+xlim(2004,2013)
    p
  }else{
    df<-data.frame(data=c(FIT.data$mean,MODIS.data$mean),
                   sd_high=c(FIT.data$mean+FIT.data$sd,MODIS.data$mean+MODIS.data$sd),
                   sd_low=c(FIT.data$mean-FIT.data$sd,MODIS.data$mean-MODIS.data$sd),
                   year=rep(months,2),
                   type=rep(c("LPJmL-FIT","MODIS"),each=length(months)),
                   col=rep(c(color.FIT, color.MODIS),each=length(months)))
    df<-data.frame(data=c(FIT.data$vec_monthly_mean,MODIS.data$vec_monthly_mean),
                   sd_high=c(FIT.data$vec_monthly_mean+FIT.data$vec_monthly_sd,MODIS.data$vec_monthly_mean+MODIS.data$vec_monthly_sd),
                   sd_low=c(FIT.data$vec_monthly_mean-FIT.data$vec_monthly_sd,MODIS.data$vec_monthly_mean-MODIS.data$vec_monthly_sd),
                   year=rep(rep(years,each=12)+rep((1:12)/12-1/24,length(years)),2),
                   type=rep(c("LPJmL-FIT","MODIS"),each=12*length(years)),
                   col=rep(c(color.FIT, color.MODIS),each=12*length(years)))
    
    p<-ggplot(df)
    p<-p+geom_ribbon(data=df,alpha=alpha,aes(x=year,ymin=sd_low,ymax=sd_high,fill=type,linetype=NA))
    p<-p+geom_line(aes(x=year,y=data,colour=type),size=linesize)
    p<-p+scale_color_manual(name='',breaks = c("LPJmL-FIT", "MODIS"),values=c(color.FIT, color.MODIS))
    p<-p+scale_fill_manual(name='',breaks = c("LPJmL-FIT", "MODIS"),values=c(color.FIT, color.MODIS),guide = FALSE)
    p<-p+ylab("GPP [gC/m²]")+xlab("Year")
    p<-p+ggtitle(plot.names[isite])
    p<-p+scale_x_continuous(breaks = pretty(df$year, n = 5))
    p
  }
  
  
  ls.gpp.monthly.long[[isite]] <- p
  ls.gpp.monthly.long_nolegend[[isite]] <- p+theme(legend.position="none")
  
}

colnames(array.error) <- rep(error.names,3)
rownames(array.error) <- site.names

print(array.error)

df.error <- rbind(
  data.frame(data.frame(array.error[,1:3]),plot = 'GPP', what = factor(site.names, levels = site.names), 
             row.names= NULL, check.names = FALSE), 
  data.frame(data.frame(array.error[,4:6]),plot = 'GPP seasonal dyn.',what = factor(site.names, levels = site.names), 
             row.names = NULL,check.names = FALSE),
  data.frame(data.frame(array.error[,7:9]),plot = 'GPP interannual', what = factor(site.names, levels = site.names), 
             row.names= NULL, check.names = FALSE)
)

mytable <- tabular(Heading()*what ~ plot*(`ME` +`cor`+`NMSE`)*Heading()*(identity),data=df.error)
mytable   

grid.table(mytable)

#plot.all.gpp.monthly<-grid.arrange(grobs=ls.gpp.monthly,
                                   #ncol=3, nrow=2)

plot.all.gpp.monthly_nolegend<-grid.arrange(grobs=ls.gpp.monthly_nolegend,
                                   ncol=3, nrow=2)
plot.all.gpp.monthly.long_nolegend<-grid.arrange(grobs=ls.gpp.monthly.long_nolegend,
                                            ncol=3, nrow=2)

path.out<- "D:/Data/LPJmL-Fit/validation/"
#ggsave(paste(path.out,"GPP.jpg",sep=""),device = "jpg",dpi = 600,plot.all.gpp.monthly,scale=3,width=20,height=7, units="cm")

ggsave(paste(path.out,"GPP_nolegend.jpg",sep=""),device = "jpg",dpi = 600,plot.all.gpp.monthly_nolegend,scale=3,width=20,height=7, units="cm")



for(isite in 1:length(site.names)){
  tab<-array.error[isite,,]
  colnames(tab) <- error.names
  rownames(tab) <- c("monthly mean","annual sum","time series")
  print(site.names[isite])
  print(tab)
}

test<-array.error
dim(test)<-c(3*length(error.names),length(site.names))

array.error

test<-t(test)

test

######### combine everything


plot.all.gpp.annual<-grid.arrange(grobs=ls.gpp.annual,
                       ncol=3, nrow=3)
plot.all.gpp.annual


plot.all.gpp.monthly<-grid.arrange(grobs=ls.gpp.monthly,
                       ncol=3, nrow=3)
plot.all.gpp.monthly

plot.all.gpp.monthly.long<-grid.arrange(grobs=ls.gpp.monthly.long,
                                   ncol=3, nrow=3)
plot.all.gpp.monthly.long



p.all.gpp<-ggarrange(plot.all.gpp.monthly_nolegend,plot.all.gpp.monthly.long_nolegend,labels=c("a)","b)"),
                     align="none",font.label = list(size = 40),hjust=10,
                     ncol=1,nrow=2)
p.all.gpp

ggsave(paste(path.out,"GPP_all.jpg",sep=""),device = "jpg",dpi = 300,p.all.gpp,scale=4,width=18,height=12, units="cm")


#### OLD ####

#### NPP ####

ls.npp<-list()
for(isite in 1:length(site.names)){
  
  FIT.data<-get.FIT.NPP(runs[isite])
  MODIS.data<-get.MODIS.NPP(site.names[isite])
  
  df<-data.frame(data=c(FIT.data$annual_sum,MODIS.data$mean),
                 sd_high=c(FIT.data$annual_sum+FIT.data$sd_annual_sum,MODIS.data$mean+MODIS.data$sd),
                 sd_low=c(FIT.data$annual_sum-FIT.data$sd_annual_sum,MODIS.data$mean-MODIS.data$sd),
                 year=rep(years,2),
                 type=rep(c("LPJmL-FIT","MODIS"),each=length(years)),
                 col=rep(c(color.FIT, color.MODIS),each=length(years)))
  
  
  p<-ggplot(df)
  p<-p+geom_ribbon(data=df,alpha=alpha,aes(x=year,ymin=sd_low,ymax=sd_high,fill=type,linetype=NA))
  p<-p+geom_line(aes(x=year,y=data,colour=type),size=linesize)
  p<-p+scale_color_manual(name='',breaks = c("LPJmL-FIT", "MODIS"),values=c(color.FIT, color.MODIS))
  p<-p+scale_fill_manual(name='',breaks = c("LPJmL-FIT", "MODIS"),values=c(color.FIT, color.MODIS),guide = FALSE)
  p<-p+ylab("NPP [gC/m²]")+xlab("Year")
  p<-p+ggtitle(plot.names[isite])
  p<-p+scale_x_continuous(breaks = pretty(df$year, n = 5))
  #p
  
  ls.npp[[isite]] <- p
}


######### combine everything


legend <- g_legend(p)

plot.all.npp<-grid.arrange(grobs=ls.npp,
                           ncol=3, nrow=3)

#plot.all.npp



### fpar ###

ls.fpar<-list()
for(isite in 1:length(site.names)){
  
  FIT.data<-get.FIT.fpar(runs[isite])
  MODIS.data<-get.MODIS.fpar(site.names[isite])
  
  df<-data.frame(data=c(FIT.data$mean,MODIS.data$mean),
                 sd_high=c(FIT.data$mean+FIT.data$sd,MODIS.data$mean+MODIS.data$sd),
                 sd_low=c(FIT.data$mean-FIT.data$sd,MODIS.data$mean-MODIS.data$sd),
                 year=rep(months,2),
                 type=rep(c("LPJmL-FIT","MODIS"),each=12),
                 col=rep(c(color.FIT, color.MODIS),each=12))
  
  p<-ggplot(df)
  p<-p+geom_ribbon(data=df,alpha=alpha,aes(x=year,ymin=sd_low,ymax=sd_high,fill=type,linetype=NA))
  p<-p+geom_line(aes(x=year,y=data,colour=type),size=linesize)
  p<-p+scale_color_manual(name='',breaks = c("LPJmL-FIT", "MODIS"),values=c(color.FIT, color.MODIS))
  p<-p+scale_fill_manual(name='',breaks = c("LPJmL-FIT", "MODIS"),values=c(color.FIT, color.MODIS),guide = FALSE)
  p<-p+ylab("fpar")+xlab("Month")
  p<-p+ggtitle(plot.names[isite])
  p<-p+scale_x_continuous(breaks = pretty(df$year, n = 5))
  p
  
  ls.fpar[[isite]] <- p
}


######### combine everything


legend <- g_legend(p)

plot.all.fpar<-grid.arrange(grobs=ls.fpar,
                       ncol=3, nrow=3)
plot.all.fpar

#dev.off()

#pdf("D:/Maik/Documents/R/PlotValidation/plot_Fluxes_2.pdf",height=10,width = 14, paper="a4r")


#plot.all.npp<-grid.arrange(grobs=ls.npp,
                            #ncol=3, nrow=3)
#plot.all.gpp.annual<-grid.arrange(grobs=ls.gpp.annual,
                            #ncol=3, nrow=3)
plot.all.gpp.monthly<-grid.arrange(grobs=ls.gpp.monthly,
                            ncol=2, nrow=4)
#plot.all.fpar<-grid.arrange(grobs=ls.fpar,
                            #ncol=3, nrow=3)

#dev.off()

ggsave(paste("D:/R/PlotValidation/plot_Fluxes_","NPP","_neu.pdf",sep=""), arrangeGrob(grobs = ls.npp,nrow=2,ncol=3),scale=3,width=20,height=7, units="cm")
ggsave(paste("D:/R/PlotValidation/plot_Fluxes_","GPP_annual","_neu.pdf",sep=""), arrangeGrob(grobs = ls.gpp.annual,nrow=2,ncol=3),scale=3,width=20,height=7, units="cm")
ggsave(paste("D:/R/PlotValidation/plot_Fluxes_","GPP_monthly","_neu.pdf",sep=""), arrangeGrob(grobs = ls.gpp.monthly,nrow=2,ncol=3),scale=3,width=20,height=7, units="cm")
ggsave(paste("D:/R/PlotValidation/plot_Fluxes_","GPP_monthly_long","_neu.pdf",sep=""), arrangeGrob(grobs = ls.gpp.monthly.long,nrow=2,ncol=3),scale=3,width=23,height=5, units="cm")
ggsave(paste("D:/R/PlotValidation/plot_Fluxes_","FPAR","_neu.pdf",sep=""), arrangeGrob(grobs = ls.fpar,nrow=2,ncol=3),scale=3,width=20,height=7, units="cm")


#widths=unit(c(300,300,300,300), c("mm", "mm","mm","mm")