library(abind)

rescale.data<-function(df_test,variables,plot=F){
  library(rcompanion)
  
  data<-as.data.frame(df_test)
  if(nrow(data)>5000){
    irow<-sample(1:nrow(data), 5000, replace=F)
    data<-data[irow,]
  }
  
  
  lambdas<-rep(NA,length(variables))
  
  for(ivar in 1:length(variables)){
    if(plot){
      plotNormalHistogram(data[,variables[ivar]],main = paste("Histogram of" , variables[ivar]),xlab=variables[ivar])
      qqnorm(data[,variables[ivar]],main = paste("Normal Q-Q Plot of" , variables[ivar]))
      qqline(data[,variables[ivar]],col="red")
    }
    
    lambda<-transformTukey(data[,variables[ivar]],
                           plotit=F,quiet=T,returnLambda = T)
    lambdas[ivar]<-lambda
    
    if(lambda==0){
      data[,paste(variables[ivar],"_trans",sep="")]<-log(data[,variables[ivar]])
    }else{
      if(lambda>0){
        data[,paste(variables[ivar],"_trans",sep="")]<-data[,variables[ivar]]^lambda
      }else{
        if(lambda<0){
          data[,paste(variables[ivar],"_trans",sep="")]<--1*data[,variables[ivar]]^lambda
        }
      }
    }
    if(plot){
      plotNormalHistogram(data[,paste(variables[ivar],"_trans",sep="")],main = paste("Histogram of" , paste(variables[ivar],"_trans, lambda=",lambda,sep="")),xlab=paste(variables[ivar],"_trans",sep=""))
      qqnorm(data[,paste(variables[ivar],"_trans",sep="")],main = paste("Normal Q-Q Plot of" , paste(variables[ivar],"_trans, lambda=",lambda,sep="")))
      qqline(data[,paste(variables[ivar],"_trans",sep="")],col="red")
    }
  }
  
  
  ### rescale data ###
  df_test_trans<-df_test
  for(ivar in 1:length(variables)){
    lambda<-lambdas[ivar]
    
    if(lambda==0){
      df_test_trans[,variables[ivar]]<-log(df_test_trans[,variables[ivar]])
    }else{
      if(lambda>0){
        df_test_trans[,variables[ivar]]<-df_test_trans[,variables[ivar]]^lambda
      }else{
        if(lambda<0){
          df_test_trans[,variables[ivar]]<--1*df_test_trans[,variables[ivar]]^lambda
        }
      }
    }
  }
  return(df_test_trans)
}


read.mult.runs<-function(filenames, varid){
  nruns<-length(filenames)
  for(irun in 1:nruns){
    ncfile<-nc_open(paste(path.in,filenames[irun],sep=""))
    if(irun==1){
      data<-ncvar_get(ncfile, varid=varid)
    }else{
      data<-data+ncvar_get(ncfile, varid=varid)
    }
    nc_close(ncfile)
  }
  return(data/nruns)
}

read.mult.runs.path<-function(filenames,path.in, varid){
  nruns<-length(filenames)
  for(irun in 1:nruns){
    ncfile<-nc_open(paste(path.in,filenames[irun],sep=""))
    data.append<-ncvar_get(ncfile, varid=varid)
    if(irun==1){
      data<-data.append
    }else{
      if(any(!dim(data.append)==dim(data))){
        data.append<-abind(data.append,array(NA,c(dim(data.append)[1:4])),rev.along = 1)
      }
      data<-data+data.append
    }
    nc_close(ncfile)
  }
  return(data/nruns)
}

read.mult.runs.2<-function(filenames, varid){
  nruns<-length(filenames)
  ncfile<-nc_open(paste(path.in,filenames[1],sep="")) #"irun = 1"
  data.all<-ncvar_get(ncfile, varid=varid)
  nc_close(ncfile)
  for(irun in 2:nruns){
    ncfile<-nc_open(paste(path.in,filenames[irun],sep=""))
    if(irun==2){
      data.all<-abind(data.all,ncvar_get(ncfile, varid=varid),rev.along = 0)
    }else{
      data.all<-abind(data.all,ncvar_get(ncfile, varid=varid),rev.along = 1)
    }
    nc_close(ncfile)}
  return(data.all)
}



read.mult.runs.2.trait<-function(filenames, varid,binsname){
  print(1/nruns)
  nruns<-length(filenames)
  ncfile<-nc_open(paste(path.in,filenames[1],sep="")) #"irun = 1"
  data.all<-ncvar_get(ncfile, varid=varid)
  bins<-ncvar_get(ncfile,varid=binsname)
  nc_close(ncfile)

  data_pft<-apply(data.all,c(1,2,3,5),sum,na.rm=T) #sum over all pfts
  data.all<-apply(data_pft,c(1,2,4),function(x) sum(x*bins/sum(x,na.rm=T)))
  
  for(irun in 2:nruns){
    print(irun/nruns)
    ncfile<-nc_open(paste(path.in,filenames[irun],sep=""))
    data.dummy<-ncvar_get(ncfile, varid=varid)
    nc_close(ncfile)
    data_pft<-apply(data.dummy,c(1,2,3,5),sum,na.rm=T) #sum over all pfts
    var_mean<-apply(data_pft,c(1,2,4),function(x) sum(x*bins/sum(x,na.rm=T)))
    if(irun==2){
      data.all<-abind(data.all,var_mean,rev.along = 0)
    }else{
      data.all<-abind(data.all,var_mean,rev.along = 1)
    }
    }
  return(data.all)
}


get.bins<-function(filenames,binsname){
  ncfile<-nc_open(paste(path.in,filenames[1],sep=""))
  bins<-ncvar_get(ncfile, varid=binsname)
  nc_close(ncfile)
  return(bins)
}

get.bins.path<-function(filenames,path.in,binsname){
  ncfile<-nc_open(paste(path.in,filenames[1],sep=""))
  bins<-ncvar_get(ncfile, varid=binsname)
  nc_close(ncfile)
  return(bins)
}

shannon_ind<-function(x){
  return(-sum((x/sum(x))*log(x/sum(x)),na.rm=T))
}

remove_no_data<-function(var){
  df_nodata<- data.frame(data=as.vector(apply(!nodata,c(1,2),any,na.rm=F)),lon=coord[,1],lat=coord[,2])
  df_nodata$data<-!df_nodata$data
  df_nodata<-na.omit(df_nodata)
  df_nodata<-df_nodata[which(df_nodata$data),]
  
  df <- data.frame(data=as.vector(var),lon=coord[,1],lat=coord[,2])
  df <-subset(df, !is.na(data)) #remove water tiles
  matched<-row.match(df_nodata[,c("lat","lon")], df[,c("lat","lon")])
  df[matched[!is.na(matched)],"data"]<-NA
  df<-na.omit(df)
  df<-merge(df,df.clm)
  return(df)
}

remove_no_data_pft<-function(var,ipft){
  df_nodata<- data.frame(data=as.vector(nodata[,,ipft]),lon=coord[,1],lat=coord[,2])
  df_nodata<-na.omit(df_nodata)
  df_nodata<-df_nodata[which(df_nodata$data),]
  
  df <- data.frame(data=as.vector(var),lon=coord[,1],lat=coord[,2])
  df <-subset(df, !is.na(data)) #remove water tiles
  matched<-row.match(df_nodata[,c("lat","lon")], df[,c("lat","lon")])
  df[matched[!is.na(matched)],"data"]<-NA
  #df<-na.omit(df)
  df<-merge(df,df.clm)
  return(df)
}

sd_function<-function(x,bins){
  mu<-sum(x*bins/sum(x))
  p<-bins*x/(sum(bins*x))
  return(sqrt(sum(p*(bins-mu)^2)))
}

ceiling_dec <- function(x, level=1) round(x + 5*10^(-level-1), level)
floor_dec <- function(x, level=1) round(x - 5*10^(-level-1), level)

shannon_ind<-function(x){
  return(-sum((x/sum(x))*log(x/sum(x)),na.rm=T))
}

RaosQ<-function(x){
  len<-length(x)
  x<-x/sum(x)
  Q<-0
  for(i in 1:(len-1)){
    for(j in (i+1):len){
      Q<-Q+abs(i-j)*x[i]*x[j]
    }
  }
  return(Q)
}


plot_trait<-function(limits.mean=NA,limits.Rao=NA,limits.shannon=NA,limits.pfts=NA){
  ### for all PFTs
  
  col.zscale<-rev(col.zscale)
  
  point.size<-1.7
  mytheme<-theme(axis.text=element_text(size=30),
                 axis.title=element_text(size=30))
  
  data_var<-apply(data_pft[,,,],c(1,2,3),sum,na.rm=T) #sum over all pfts
  var_mean<-apply(data_var,c(1,2),function(x) sum(x*bins/sum(x,na.rm=T)))
  var_sd<-apply(data_var,c(1,2),function(x) sd_function(x,bins))
  var_shannon<-apply(data_var,c(1,2),function(x) shannon_ind(x))
  var_Rao<-apply(data_var,c(1,2), function(x) RaosQ(x))
  
  ### mean variable ###
  #mask deserts and water tiles
  
  df<-remove_no_data(var_mean)
  
  p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
  p<-p+ggtitle(paste(titlename,", all PFTs",sep=""))+ theme(panel.background = element_rect(fill = col.plot.bg))
  if(is.na(limits.mean)){
    p<-p+labs(subtitle=paste("max: ",ceiling_dec(max(df$data,na.rm=T),2)," min: ",floor_dec(min(df$data,na.rm=T),2),sep=""))
    limits.mean<-c(min(df$data),max(df$data))
  }
  p<-p+scale_fill_gradientn(limits=limits.mean,colors=col.zscale, name=zscalename,na.value=col.nodata)
  p<-p+geom_tile(data=df_nodata,aes(y=lat,x=lon),color=col.nodata,fill=col.nodata,ylim=c())
  p<-p+guides(fill=guide_colourbar(barwidth = 2,barheight = 20))
  p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
  p<-p+ylim(latmin,latmax)
  print(p) #plot
  
  ### plot climate scatterplots ### MCWD
  
  p<-ggplot(df,aes(x=temp,y=MCWD,color=data))+geom_point(size=point.size)
  p<-p<-p+scale_color_gradientn(limits=limits.mean,colors=col.zscale,name=zscalename)
  p<-p+ggtitle(paste(titlename,", all PFTs",sep=""))
  p<-p+mytheme
  p<-p+xlab("T [°C]")+xlim(limits.temp)+ylim(limits.MCWD)
  p<-p+theme()
  print(p)
  
  ### Rao's Q ###
  
  df<-remove_no_data(var_Rao)
  
  p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
  p<-p+ggtitle(paste("Rao's Q",", all PFTs",sep=""))+ theme(panel.background = element_rect(fill = col.plot.bg))
  if(is.na(limits.Rao)){
    p<-p+labs(subtitle=paste("max: ",ceiling_dec(max(df$data,na.rm=T),2)," min: ",floor_dec(min(df$data,na.rm=T),2),sep=""))
    limits.Rao<-c(min(df$data),max(df$data))
  }
  p<-p+scale_fill_gradientn(limits=limits.Rao,colors=col.zscale2,name="Rao's Q",na.value=col.nodata)
  p<-p+geom_tile(data=df_nodata,aes(y=lat,x=lon),color=col.nodata,fill=col.nodata)
  
  p<-p+guides(fill=guide_colourbar(barwidth = 2,barheight = 20))
  p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
  p<-p+ylim(latmin,latmax)
  print(p) #plot
  
  ### plot climate scatterplots ### MCWD
  p<-ggplot(df,aes(x=temp,y=MCWD,color=data))+geom_point(size=point.size)
  p<-p<-p+scale_color_gradientn(colors=col.zscale2,name="Rao's Q")
  p<-p+ggtitle(paste("Rao's Q",", all PFTs",sep=""))
  p<-p+mytheme
  p<-p+xlab("T [°C]")+xlim(limits.temp)+ylim(limits.MCWD)
  print(p)
  
  ### shannon index ###
  df<-remove_no_data(var_shannon)
  
  p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
  p<-p+ggtitle(paste("Mean Shannon Index",", all PFTs",sep=""))+ theme(panel.background = element_rect(fill = col.plot.bg))
  if(is.na(limits.shannon)){
    p<-p+labs(subtitle=paste("max: ",ceiling_dec(max(df$data,na.rm=T),2)," min: ",floor_dec(min(df$data,na.rm=T),2),sep=""))
    limits.shannon<-c(min(df$data),max(df$data))
  }
  p<-p+scale_fill_gradientn(limits=limits.shannon,colors=col.zscale2,name="Shannon Index",na.value=col.nodata)
  p<-p+geom_tile(data=df_nodata,aes(y=lat,x=lon),color=col.nodata,fill=col.nodata)
  
  p<-p+guides(fill=guide_colourbar(barwidth = 2,barheight = 20))
  p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
  p<-p+ylim(latmin,latmax)
  print(p) #plot
  
  ### plot climate scatterplots ### MCWD
  p<-ggplot(df,aes(x=temp,y=MCWD,color=data))+geom_point(size=point.size)
  p<-p<-p+scale_color_gradientn(limits=limits.shannon,colors=col.zscale2,name="Shannon Index")
  p<-p+ggtitle(paste("Mean Shannon Index",", all PFTs",sep=""))
  p<-p+mytheme
  p<-p+xlab("T [°C]")+xlim(limits.temp)+ylim(limits.MCWD)
  print(p)
  col.zscale<-rev(col.zscale)
  
  for(ipft in 1:(length(namepft))){ #for every pft
    data_var<-data_pft[,,,ipft]
    var_mean<-apply(data_var,c(1,2),function(x) sum(x*bins/sum(x)))
    var_sd<-apply(data_var,c(1,2),function(x) sd_function(x,bins))
    var_shannon<-apply(data_var,c(1,2),function(x) shannon_ind(x))
    var_Rao<-apply(data_var,c(1,2),function(x) RaosQ(x))
    
    df<-remove_no_data_pft(var_mean,ipft)
    
    p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
    p<-p+ggtitle(paste(titlename,", ",namepft[ipft],sep=""))+ theme(panel.background = element_rect(fill = col.plot.bg))
    if(is.na(limits.pfts)){
      p<-p+scale_fill_gradientn(colors=col.zscale,name=zscalename,na.value=col.nodata)
      p<-p+labs(subtitle=paste("max: ",ceiling_dec(max(df$data,na.rm=T),2)," min: ",floor_dec(min(df$data,na.rm=T),2),sep=""))
    }else{
      p<-p+scale_fill_gradientn(limits=c(limits.pfts[ipft,1,1],limits.pfts[ipft,1,2]),colors=col.zscale,name=zscalename,na.value=col.nodata)
    }
    p<-p+geom_tile(data=df_nodata,aes(y=lat,x=lon),color=col.nodata,fill=col.nodata)
    p<-p+guides(fill=guide_colourbar(barwidth = 2,barheight = 20))
    p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
    p<-p+ylim(latmin,latmax)
    print(p) #plot
    
    ### plot climate scatterplots ### MCWD
    
    p<-ggplot(df,aes(x=temp,y=MCWD,color=data))+geom_point(size=point.size)
    p<-p<-p+scale_color_gradientn(colors=col.zscale,name=zscalename,na.value=col.nodata)
    p<-p+ggtitle(paste(titlename,", ",namepft[ipft],sep=""))
    p<-p+mytheme
    p<-p+xlab("T [°C]")+xlim(limits.temp)+ylim(limits.MCWD)
    print(p)
    
    ### Rao's Q ###
    
    df<-remove_no_data_pft(var_Rao,ipft)
    
    p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
    p<-p+ggtitle(paste("Rao's Q",", ",namepft[ipft],sep=""))+ theme(panel.background = element_rect(fill = col.plot.bg))
    if(is.na(limits.pfts)){
      p<-p+scale_fill_gradientn(colors=col.zscale2,name="Rao's Q",na.value=col.nodata)
      p<-p+labs(subtitle=paste("max: ",ceiling_dec(max(df$data,na.rm=T),2)," min: ",floor_dec(min(df$data,na.rm=T),2),sep=""))
    }else{
      p<-p+scale_fill_gradientn(limits=c(limits.pfts[ipft,2,1],limits.pfts[ipft,2,2]),colors=col.zscale2,name="Rao's Q",na.value=col.nodata)}
    p<-p+geom_tile(data=df_nodata,aes(y=lat,x=lon),color=col.nodata,fill=col.nodata)
    p<-p+guides(fill=guide_colourbar(barwidth = 2,barheight = 20))
    p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
    p<-p+ylim(latmin,latmax)
    print(p) #plot
    
    ### plot climate scatterplots ### MCWD
    p<-ggplot(df,aes(x=temp,y=MCWD,color=data))+geom_point(size=point.size)
    p<-p<-p+scale_color_gradientn(colors=col.zscale2,name="Rao's Q",na.value=col.nodata)
    p<-p+ggtitle(paste("Rao's Q",", ",namepft[ipft],sep=""))
    p<-p+mytheme
    p<-p+xlab("T [°C]")+xlim(limits.temp)+ylim(limits.MCWD)
    print(p)
    
    ### shannon index ###
    df<-remove_no_data_pft(var_shannon,ipft)
    
    p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
    p<-p+ggtitle(paste("Mean Shannon Index",", ",namepft[ipft],sep=""))+ theme(panel.background = element_rect(fill = col.plot.bg))
    if(is.na(limits.pfts)){
      p<-p+scale_fill_gradientn(colors=col.zscale,name="Shannon Index",na.value=col.nodata)
      p<-p+labs(subtitle=paste("max: ",ceiling_dec(max(df$data,na.rm=T),2)," min: ",floor_dec(min(df$data,na.rm=T),2),sep=""))
    }else{
      p<-p+scale_fill_gradientn(limits=c(limits.pfts[ipft,3,1],limits.pfts[ipft,3,2]),colors=col.zscale2,name="Shannon Index",na.value=col.nodata)
    }
    p<-p+geom_tile(data=df_nodata,aes(y=lat,x=lon),color=col.nodata,fill=col.nodata)
    p<-p+guides(fill=guide_colourbar(barwidth = 2,barheight = 20))
    p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
    p<-p+ylim(latmin,latmax)
    print(p) #plot
    
    ### plot climate scatterplots ### MCWD
    p<-ggplot(df,aes(x=temp,y=MCWD,color=data))+geom_point(size=point.size)
    p<-p<-p+scale_color_gradientn(colors=col.zscale2,name="Shannon Index",na.value=col.nodata)
    p<-p+ggtitle(paste("Mean Shannon Index",", ",namepft[ipft],sep=""))
    p<-p+xlab("T [°C]")+xlim(limits.temp)+ylim(limits.MCWD)
    p<-p+mytheme
    print(p)
    
    
    ### same scale for all ###############
    
    df<-remove_no_data_pft(var_mean,ipft)
    
    p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
    p<-p+scale_fill_gradientn(limits=limits.mean,colors=col.zscale,name=zscalename,na.value=col.nodata)
    p<-p+geom_tile(data=df_nodata,aes(y=lat,x=lon),color=col.nodata,fill=col.nodata)
    p<-p+ggtitle(paste(titlename,", ",namepft[ipft],", same scale",sep=""))+ theme(panel.background = element_rect(fill = col.plot.bg))
    p<-p+guides(fill=guide_colourbar(barwidth = 2,barheight = 20))
    p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
    p<-p+ylim(latmin,latmax)
    print(p) #plot
    
    ### plot climate scatterplots ### MCWD
    
    p<-ggplot(df,aes(x=temp,y=MCWD,color=data))+geom_point(size=point.size)
    p<-p<-p+scale_color_gradientn(limits=limits.mean,colors=col.zscale,name=zscalename,na.value=col.nodata)
    p<-p+ggtitle(paste(titlename,", ",namepft[ipft],", same scale",sep=""))
    p<-p+xlab("T [°C]")+xlim(limits.temp)+ylim(limits.MCWD)
    p<-p+mytheme
    print(p)
    
    ### Rao's Q ###
    
    df<-remove_no_data_pft(var_Rao,ipft)
    
    p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
    p<-p+scale_fill_gradientn(limits=limits.Rao,colors=col.zscale2,name="Rao's Q",na.value=col.nodata)
    p<-p+geom_tile(data=df_nodata,aes(y=lat,x=lon),color=col.nodata,fill=col.nodata)
    p<-p+ggtitle(paste("Rao's Q",", ",namepft[ipft],", same scale",sep=""))+ theme(panel.background = element_rect(fill = col.plot.bg))
    p<-p+guides(fill=guide_colourbar(barwidth = 2,barheight = 20))
    p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
    p<-p+ylim(latmin,latmax)
    print(p) #plot
    
    ### plot climate scatterplots ### MCWD
    p<-ggplot(df,aes(x=temp,y=MCWD,color=data))+geom_point(size=point.size)
    p<-p<-p+scale_color_gradientn(limits=limits.Rao,colors=col.zscale2,name="Rao's Q",na.value=col.nodata)
    p<-p+ggtitle(paste("Rao's Q",", ",namepft[ipft],", same scale",sep=""))
    p<-p+xlab("T [°C]")+xlim(limits.temp)+ylim(limits.MCWD)
    p<-p+mytheme
    print(p)
    
    ### shannon index ###
    df<-remove_no_data_pft(var_shannon,ipft)
    
    p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
    p<-p+scale_fill_gradientn(limits=limits.shannon,colors=col.zscale2,name="Shannon Index",na.value=col.nodata)
    p<-p+geom_tile(data=df_nodata,aes(y=lat,x=lon),color=col.nodata,fill=col.nodata)
    p<-p+ggtitle(paste("Mean Shannon Index",", ",namepft[ipft],", same scale",sep=""))+ theme(panel.background = element_rect(fill = col.plot.bg))
    p<-p+guides(fill=guide_colourbar(barwidth = 2,barheight = 20))
    p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
    p<-p+ylim(latmin,latmax)
    print(p) #plot
    
    ### plot climate scatterplots ### MCWD
    p<-ggplot(df,aes(x=temp,y=MCWD,color=data))+geom_point(size=point.size)
    p<-p<-p+scale_color_gradientn(limits=limits.shannon,colors=col.zscale2,name="Shannon Index",na.value=col.nodata)
    p<-p+ggtitle(paste("Mean Shannon Index",", ",namepft[ipft],", same scale",sep=""))
    p<-p+xlab("T [°C]")+xlim(limits.temp)+ylim(limits.MCWD)
    p<-p+mytheme
    print(p)
    
  }
  
}

filter.latlon<-function(TRY,lat,lon){
  TRY<-TRY[complete.cases(TRY[ , 11:12]),]
  TRY<-TRY[TRY$V15 > lat[1], ]
  TRY<-TRY[TRY$V15 < lat[2], ]
  TRY<-TRY[TRY$V16 > lon[1], ]
  TRY<-TRY[TRY$V16 < lon[2], ]
  return(TRY)
}
filter.broadleaved<-function(TRY){
  return(TRY[TRY$V12 %in% 1,])
}
filter.needleleaved<-function(TRY){
  return(TRY[TRY$V12 %in% 2,])
}
filter.deciduous<-function(TRY){
  return(TRY[TRY$V13 %in% 1,])
}
filter.evergreen<-function(TRY){
  return(TRY[TRY$V13 %in% 2,])
}


filter.WD<-function(TRY){
  return(TRY[TRY$V6 %in% 4,])
}

filter.SLA.petiole<-function(TRY){
  return(TRY[TRY$V5 %in% 3116,])
}

filter.SLA.nopetiole<-function(TRY){
  return(TRY[TRY$V5 %in% 3115,])
}

filter.SLA.undefpetiole<-function(TRY){
  return(TRY[TRY$V5 %in% 3117,])
}

filter.SLA.notclearpetiole<-function(TRY){
  return(TRY[TRY$V5 %in% 12,])
}

get.Trait<-function(TRY){
  return(as.numeric(TRY$V8))
}
get.lat<-function(TRY){
  return(as.numeric(TRY$V15))
}
get.lon<-function(TRY){
  return(as.numeric(TRY$V16))
}


g_legend<-function(a.gplot){ 
  tmp <- ggplot_gtable(ggplot_build(a.gplot)) 
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box") 
  legend <- tmp$grobs[[leg]] 
  return(legend)} 

get.pheno <- function(site){
  
  path<-"D:/Maik/Documents/Data/UFZ/PhenoMetrics/"
  #method<-"Method_DLogistics" #double logistic function for interpolation
  method<-"Method_LinIP" #linear interpolation
  modistype<-"terra"
  thresholdmethod<-"LT=0.65" #(GT global threshold, or LT local theshold)
  
  dir<-paste(path,method,"/",site, sep="")
  filenames<-grep(modistype, list.files(dir,pattern=".nc"), value = TRUE)
  
  #get dimensions of files (lat, lon)
  ncfile<-nc_open(paste(dir,"/",filenames[1],sep=""))
  dims<-dim(ncvar_get(ncfile, varid=paste("GU-",thresholdmethod,sep="")))
  nc_close(ncfile)
  greenup<-sen<-array(NA,c(dims,15))
  
  for(ifile in 1:length(filenames)){
    ncfile<-nc_open(paste(dir,"/",filenames[ifile],sep=""))
    names(ncfile$var)
    greenup[,,ifile]<-ncvar_get(ncfile, varid=paste("GU-",thresholdmethod,sep=""))
    sen[,,ifile]<-ncvar_get(ncfile, varid=paste("SEN-",thresholdmethod,sep=""))
    nc_close(ncfile)
  }
  return(list("greenup"=greenup,"sen"=sen))
}

#site<-"NorthernLimestone"
#site<-"PenedaGeres"

#x<-get.pheno(site="PenedaGeres")
#plot(apply(x$greenup,3,mean,na.rm=T))


get.MODIS.GPP <- function(site){
  
  dayspermonth<-c(31,28,31,30,31,30,31,31,30,31,30,31)
  
  path<-paste("D:/Data/GIS/MODIS/",site,"/MOD17A2H.006_aid0001.nc",sep="")
  
  #get dimensions of files (lat, lon)
  ncfile<-nc_open(path)
  data<-ncvar_get(ncfile, varid="Gpp_500m")
  time<-ncvar_get(ncfile, varid="time")
  lon<-ncvar_get(ncfile, varid="lon")
  lat<-ncvar_get(ncfile, varid="lat")
  nc_close(ncfile)
  
  #load modis pixels that have enough tree cover
  ind_modis<-read.csv(paste("D:/Data/GIS/MODIS/",site,"/",site,"_Modispixels_80_500m.csv",sep=""))$x
  mask<-rep(NA,length(lon)*length(lat)) #empty vector
  mask[ind_modis]<-1 #1 if cell has enough tree cover
  dim(mask)<-c(length(lon),length(lat)) #convert to matrix
  #multiply for every time dimension mask matrix and lat lon dimensions of data
  for(i in 1:length(time)){
    data[,,i]<-mask*data[,,i]
  }
  
  time<-as.Date("2000-01-01")+time #convert time 
  years<-year(time) #get all years
  months<-month(time) #get all months
  data_monthly<-array(NA,c(length(lon),length(lat),14,12))
  for(ilon in 1:length(lon)){ #over all lons and lats
    for(ilat in 1:length(lat)){
      df<-data.frame(data=data[ilon,ilat,],year=years,month=months) #create dataframe
      df$data<-df$data*1000 #kgC/m² to gC/m²
      df$data<-df$data/8 #gC/m² to gC/m²*day
      for(imonth in 1:12){ #multiply GPP rate by days of month (mean gC/m²*month)
        df$data[df$month==imonth]<-df$data[df$month==imonth]*dayspermonth[imonth]
      }
      mean_month<-as.vector(aggregate(df$data,by=list(years,months),FUN=mean, na.rm=TRUE)$x) #aggregate months and years
      mean_month<-c(NA,mean_month) #add first missing month in 2000
      dim(mean_month)<-c(length(unique(years)),12) #set vector to array of years and months
      data_monthly[ilon,ilat,,]<-mean_month #add to grid cell
    }
  }
  
  data_monthly<-data_monthly[,,5:14,] #select 2004-2013
  
  vec_monthly_mean<-as.vector(t(apply(data_monthly,c(3,4),mean,na.rm=T)))
  vec_monthly_sd<-as.vector(t(apply(data_monthly,c(3,4),sd,na.rm=T)))
  
  mean<-apply(data_monthly,4,mean,na.rm=T) #mean over all years and cells
  sd<-apply(data_monthly,4,sd,na.rm=T) #sd over all years and cells
  
  annual_sum<-apply(apply(data_monthly,c(1,2,3),sum),3,mean,na.rm=T)
  sd_annual_sum<-apply(apply(data_monthly,c(1,2,3),sum),3,sd,na.rm=T)
  
  return(list("mean"=mean,"sd"=sd,"annual_sum"=annual_sum,"sd_annual_sum"=sd_annual_sum,"vec_monthly_mean"=vec_monthly_mean,"vec_monthly_sd"=vec_monthly_sd))
}

get.MODIS.NPP <- function(site){
  
  path<-paste("D:/Data/GIS/MODIS/",site,"/MOD17A3.055_aid0001.nc",sep="")
  
  #get dimensions of files (lat, lon)
  ncfile<-nc_open(path)
  data<-ncvar_get(ncfile, varid="Npp_1km")
  time<-ncvar_get(ncfile, varid="time")
  lon<-ncvar_get(ncfile, varid="lon")
  lat<-ncvar_get(ncfile, varid="lat")
  nc_close(ncfile)
  
  #load modis pixels that have enough tree cover
  ind_modis<-read.csv(paste("D:/Data/GIS/MODIS/",site,"/",site,"_Modispixels_80_1000m.csv",sep=""))$x
  mask<-rep(NA,length(lon)*length(lat)) #empty vector
  mask[ind_modis]<-1 #1 if cell has enough tree cover
  dim(mask)<-c(length(lon),length(lat)) #convert to matrix
  #multiply for every time dimension mask matrix and lat lon dimensions of data
  
  data<-data*1000 #kgC/m² to gC/m²
  
  if(length(dim(data))>1){
    for(i in 1:length(time)){
      data[,,i]<-mask*data[,,i]
    }
    data<-data[,,5:14]#select years 2004-2013
    
    mean<-apply(data,3,mean,na.rm=T) #mean over all years
    sd<-apply(data,3,sd,na.rm=T) #sd over all years
  }else{#only 1 pixel
    data<-data[5:14]#select years 2004-2013
    
    mean<-data
    sd<-rep(0,length(data))
  }
  
  
  return(list("mean"=mean,"sd"=sd))
}

get.MODIS.fpar <- function(site){
  
  path<-paste("D:/Maik/Documents/Data/MODIS/data/nc/",site,"/MOD15A2H.006_aid0001.nc",sep="")
  
  #get dimensions of files (lat, lon)
  ncfile<-nc_open(path)
  data<-ncvar_get(ncfile, varid="Fpar_500m")
  time<-ncvar_get(ncfile, varid="time")
  lon<-ncvar_get(ncfile, varid="lon")
  lat<-ncvar_get(ncfile, varid="lat")
  nc_close(ncfile)
  
  #load modis pixels that have enough tree cover
  ind_modis<-read.csv(paste("D:/Maik/Documents/Data/MODIS/data/nc/",site,"_Modispixels_80.csv",sep=""))$x
  mask<-rep(NA,length(lon)*length(lat)) #empty vector
  mask[ind_modis]<-1 #1 if cell has enough tree cover
  dim(mask)<-c(length(lon),length(lat)) #convert to matrix
  #multiply for every time dimension mask matrix and lat lon dimensions of data
  for(i in 1:length(time)){
    data[,,i]<-mask*data[,,i]
  }
  
  time<-as.Date("2000-01-01")+time #convert time 
  years<-year(time) #get all years
  months<-month(time) #get all months
  df<-data.frame(data=apply(data,3,mean,na.rm=T),year=years,month=months) #create dataframe
  mean_month<-as.vector(aggregate(df$data,by=list(years,months),FUN=mean, na.rm=TRUE)$x) #aggregate months and years
  mean_month<-c(NA,mean_month) #add first missing month in 2000
  
  dim(mean_month)<-c(length(unique(years)),12) #set vector to array
  
  mean_month<-mean_month[5:14,]#select years 2004-2013
  
  mean<-apply(mean_month,2,mean,na.rm=T) #mean over all years
  sd<-apply(mean_month,2,sd,na.rm=T) #sd over all years
  return(list("mean"=mean,"sd"=sd))
}


get.FIT.NPP<-function(run){
  
  nyears<-615
  temp.resolution<-12
  npatches<-16
  size<-4
  
  filenames<-gsub("","", list.files(paste(path.FIT.data,sep=""),pattern="\\.bin$")) #get all .nc filenames in path.output.FIT
  filenames<-grep(run, filenames, value = TRUE) 
  filename<-grep("mnpp", filenames, value = TRUE) 
  
  binfile  <- file(paste(path.FIT.data,"/",filename,sep=""),"rb") # read binary mode
  data.bin<-readBin(binfile,double(),n=npatches*temp.resolution*nyears, size=size)
  close(binfile)
  
  dim(data.bin)<-c(npatches,temp.resolution,nyears)
  data.sum<-apply(data.bin[,,(nyears-9):nyears],c(1,3),sum) #sum over all patches, last years only 
  
  data.sum<-data.sum*1 #scale output data
  
  annual_sum<-apply(data.sum,2,mean,na.rm=T) #mean over all years
  sd_annual_sum<-apply(data.sum,2,sd,na.rm=T) #sd over all years
  
  mean<-apply(data.bin[,,(nyears-9):nyears],c(2),mean) #mean over all patches and last years only
  sd<-apply(data.bin[,,(nyears-9):nyears],c(2),sd) #sd over all patches and last years only 
  
  return(list("mean"=mean,"sd"=sd,"annual_sum"=annual_sum,"sd_annual_sum"=sd_annual_sum))
}

get.FIT.GPP<-function(run){
  
  nyears<-113
  temp.resolution<-12
  npatches<-50
  size<-4
  
  filenames<-gsub("","", list.files(paste(path.FIT.data,sep=""),pattern="\\.bin$")) #get all .nc filenames in path.output.FIT
  filenames<-grep(run, filenames, value = TRUE) 
  filename<-grep("mgpp", filenames, value = TRUE) 
  
  binfile  <- file(paste(path.FIT.data,"/",filename,sep=""),"rb") # read binary mode
  data.bin<-readBin(binfile,double(),n=npatches*temp.resolution*nyears, size=size)
  close(binfile)
  
  dim(data.bin)<-c(npatches,temp.resolution,nyears)
  
  data.mean.patches<-apply(data.bin[,,(nyears-9):nyears],c(2,3),mean) #mean over all patches, last years only 
  data.sd.patches<-apply(data.bin[,,(nyears-9):nyears],c(2,3),sd)
  
  annual_sum<-apply(data.mean.patches,2,sum,na.rm=T) #sum over all months
  sd_annual_sum<-apply(data.sd.patches,2,sum,na.rm=T) #sd
  
  mean<-apply(data.bin[,,(nyears-9):nyears],c(2),mean) #mean over all patches and last years only
  sd<-apply(data.bin[,,(nyears-9):nyears],c(2),sd) #sd over all patches and last years only 
  
  return(list("mean"=mean,"sd"=sd,"annual_sum"=annual_sum,"sd_annual_sum"=sd_annual_sum,"vec_monthly_mean"=as.vector(data.mean.patches),"vec_monthly_sd"=as.vector(data.sd.patches)))
}


get.FIT.fpar<-function(run){
  
  nyears<-113
  temp.resolution<-12
  npatches<-200
  size<-4
  
  filenames<-gsub("","", list.files(paste(path.FIT.data,run,sep=""),pattern="\\.bin$")) #get all .nc filenames in path.output.FIT
  filename<-grep("fpar", filenames, value = TRUE) 
  
  binfile  <- file(paste(path.FIT.data,run,"/",filename,sep=""),"rb") # read binary mode
  data.bin<-readBin(binfile,double(),n=npatches*temp.resolution*nyears, size=size)
  close(binfile)
  
  dim(data.bin)<-c(npatches,temp.resolution,nyears)
  
  mean<-apply(data.bin[,,(nyears-9):nyears],c(2),mean) #mean over all patches and last years only
  sd<-apply(data.bin[,,(nyears-9):nyears],c(2),sd) #sd over all patches and last years only 
  
  return(list("mean"=mean,"sd"=sd))
}

get.FLUX<-function(site){
  path.in.flux<-paste('D:/Data/Fluxnet/',site,"/",sep="")
  
  filenames<-gsub("","", list.files(path.in.flux,pattern="\\.txt$")) #get all .nc filenames in path.output.FIT
  filenames<-grep("L4_m", filenames, value = TRUE) #Level 4 data, Monthly aggregated
  #print(filenames)
  
  resp<-nee<-gpp<-matrix(NA,12,length(filenames)) #create empty matrices
  
  for(ifile in 1:length(filenames)){
    tab<-read.csv(paste(path.in.flux,filenames[ifile],sep=""),header=T) # read table for one year
    
    resp[,ifile]<-tab$Reco_st*tab$n_days #respiration converted from gC/m²*day gC/m²*month
    nee[,ifile]<-tab$NEE_st_fMDS*tab$n_days #Net ecosystem exchange converted from gC/m²*day gC/m²*month
    gpp[,ifile]<-tab$GPP_st_MDS*tab$n_days #gross primary production converted from gC/m²*day gC/m²*month
  }
  
  gpp[gpp<0]<-0
  gpp_vec_mean<-as.vector(gpp)
  resp[gpp<0]<-NA
  nee[gpp<0]<-NA
  
  mean.resp<-apply(resp,1,mean,na.rm=T)
  mean.nee<-apply(nee,1,mean,na.rm=T)
  mean.gpp<-apply(gpp,1,mean,na.rm=T)
  
  sd.resp<-apply(resp,1,sd,na.rm=T)
  sd.nee<-apply(nee,1,sd,na.rm=T)
  sd.gpp<-apply(gpp,1,sd,na.rm=T)
  
  sum.gpp<-apply(gpp,2,sum,na.rm=T)
  years<-if(site=="Hainich"){2000:2007}else{if(site=="Laegeren"){2004:2010}}
  mean.nee<-apply(nee,1,mean,na.rm=T)
  mean.gpp<-apply(gpp,1,mean,na.rm=T)
  
  return(list(mean_resp=mean.resp,mean_nee=mean.nee,mean_gpp=mean.gpp,
              sd_resp=sd.resp,sd_nee=sd.nee,sd_gpp=sd.gpp,
              sum_gpp=sum.gpp,years=years,gpp_vec_mean=gpp_vec_mean))
}



get.FIT.SLA<-function(run){
  
  ivariable<-1
  nyears<-113
  
  variables<-c("sla","longevity","wooddens","height")
  weights<-c("ind","ind","ind","ind")
  weights<-c("mass","mass","mass","mass")
  units<-c("in mm²/mg","in a","in g/cm³","in m")
  scalings<-c(1000*0.455,1,(10^-6)/0.455,1)
  temp.resolution<-1
  bins.minmax<- matrix(c(0.005,0.07,0.2,4,0.7*10^5,6.5*10^5,0,50), nrow=length(variables), ncol=2, byrow = TRUE)
  size<-4
  nbins<-100
  npatches<-200
  nruns<-1
  
  mean.trait<-array(NA,c(nruns,length(variables)))
  
  df.data<-array(NA,c(nruns,nyears))
  df.data.hov<-array(NA,c(nruns,nyears*nbins))
  
  filenames<-gsub("","", list.files(paste(path.FIT.data,run,"/",sep=""),pattern="\\.bin$")) #get all .nc filenames in path.output.FIT
  
  filename<-grep(variables[ivariable], filenames, value = TRUE) 
  filename<-grep(weights[ivariable], filename, value = TRUE) 
  
  file  <- file(paste(path.FIT.data,run,"/",filename,sep=""),"rb")
  data.bin<-readBin(file,numeric(),n =npatches*nbins*nyears, size=size)
  close(file)
  
  dim(data.bin)<-c(npatches,nbins,nyears)
  data.bin.sum.ind<-apply(data.bin,c(1,3),sum) #sum up all indiviuals for each patch and year
  
  bins<-seq(bins.minmax[ivariable,1], bins.minmax[ivariable,2],length.out = 100)*scalings[ivariable]# get bins
  
  data.bin.sum.variable<-apply(data.bin,c(1,3),function(x) t(x)%*%bins)# sum up all trait values (bin*number of indiviuals)
  
  data.bin.mean<-apply(data.bin.sum.variable/data.bin.sum.ind,2,mean) #mean trait values in each year
  
  data.bin.sum.ind.patches<-apply(data.bin,c(2,3),sum)# sum of individuals over all cells/patches per year and bin
  
  sum.ind.peryear<-apply(data.bin.sum.ind.patches,2,sum) #sum of individuals over all cells/patches and bins per year
  
  for(i in 1:nyears){data.bin.sum.ind.patches[,i]<-data.bin.sum.ind.patches[,i]/sum.ind.peryear[i]}
  
  #data.bin.sum.ind.patches... normalized counts
  
  return(list("mean"=data.bin.mean,"counts"=data.bin.sum.ind.patches,"bins"=bins))
}

get.FIT.WD<-function(run){
  
  ivariable<-3
  nyears<-113
  
  variables<-c("sla","longevity","wooddens","height")
  weights<-c("ind","ind","ind","ind")
  weights<-c("mass","mass","mass","mass")
  units<-c("in mm²/mg","in a","in g/cm³","in m")
  scalings<-c(1000*0.455,1,(10^-6)/0.455,1)
  temp.resolution<-1
  bins.minmax<- matrix(c(0.005,0.07,0.2,4,0.7*10^5,6.5*10^5,0,50), nrow=length(variables), ncol=2, byrow = TRUE)
  size<-4
  nbins<-100
  npatches<-200
  nruns<-1
  
  mean.trait<-array(NA,c(nruns,length(variables)))
  
  df.data<-array(NA,c(nruns,nyears))
  df.data.hov<-array(NA,c(nruns,nyears*nbins))
  
  filenames<-gsub("","", list.files(paste(path.FIT.data,run,"/",sep=""),pattern="\\.bin$")) #get all .nc filenames in path.output.FIT
  
  filename<-grep(variables[ivariable], filenames, value = TRUE) 
  filename<-grep(weights[ivariable], filename, value = TRUE) 
  
  file  <- file(paste(path.FIT.data,run,"/",filename,sep=""),"rb")
  data.bin<-readBin(file,numeric(),n =npatches*nbins*nyears, size=size)
  close(file)
  
  dim(data.bin)<-c(npatches,nbins,nyears)
  data.bin.sum.ind<-apply(data.bin,c(1,3),sum) #sum up all indiviuals for each patch and year
  
  bins<-seq(bins.minmax[ivariable,1], bins.minmax[ivariable,2],length.out = 100)*scalings[ivariable]# get bins
  
  data.bin.sum.variable<-apply(data.bin,c(1,3),function(x) t(x)%*%bins)# sum up all trait values (bin*number of indiviuals)
  
  data.bin.mean<-apply(data.bin.sum.variable/data.bin.sum.ind,2,mean) #mean trait values in each year
  
  data.bin.sum.ind.patches<-apply(data.bin,c(2,3),sum)# sum of individuals over all cells/patches per year and bin
  
  sum.ind.peryear<-apply(data.bin.sum.ind.patches,2,sum) #sum of individuals over all cells/patches and bins per year
  
  for(i in 1:nyears){data.bin.sum.ind.patches[,i]<-data.bin.sum.ind.patches[,i]/sum.ind.peryear[i]}
  
  #data.bin.sum.ind.patches... normalized counts
  
  return(list("mean"=data.bin.mean,"counts"=data.bin.sum.ind.patches,"bins"=bins))
}

get.FIT.LL<-function(run){
  
  ivariable<-2
  nyears<-113
  
  variables<-c("sla","longevity","wooddens","height")
  weights<-c("ind","ind","ind","ind")
  weights<-c("mass","mass","mass","mass")
  units<-c("in mm²/mg","in a","in g/cm³","in m")
  scalings<-c(1000*0.455,1,(10^-6)/0.455,1)
  temp.resolution<-1
  bins.minmax<- matrix(c(0.005,0.07,0.2,4,0.7*10^5,6.5*10^5,0,50), nrow=length(variables), ncol=2, byrow = TRUE)
  size<-4
  nbins<-100
  npatches<-200
  nruns<-1
  
  mean.trait<-array(NA,c(nruns,length(variables)))
  
  df.data<-array(NA,c(nruns,nyears))
  df.data.hov<-array(NA,c(nruns,nyears*nbins))
  
  filenames<-gsub("","", list.files(paste(path.FIT.data,run,"/",sep=""),pattern="\\.bin$")) #get all .nc filenames in path.output.FIT
  
  filename<-grep(variables[ivariable], filenames, value = TRUE) 
  filename<-grep(weights[ivariable], filename, value = TRUE) 
  
  file  <- file(paste(path.FIT.data,run,"/",filename,sep=""),"rb")
  data.bin<-readBin(file,numeric(),n =npatches*nbins*nyears, size=size)
  close(file)
  
  dim(data.bin)<-c(npatches,nbins,nyears)
  data.bin.sum.ind<-apply(data.bin,c(1,3),sum) #sum up all indiviuals for each patch and year
  
  bins<-seq(bins.minmax[ivariable,1], bins.minmax[ivariable,2],length.out = 100)*scalings[ivariable]# get bins
  
  data.bin.sum.variable<-apply(data.bin,c(1,3),function(x) t(x)%*%bins)# sum up all trait values (bin*number of indiviuals)
  
  data.bin.mean<-apply(data.bin.sum.variable/data.bin.sum.ind,2,mean) #mean trait values in each year
  
  data.bin.sum.ind.patches<-apply(data.bin,c(2,3),sum)# sum of individuals over all cells/patches per year and bin
  
  sum.ind.peryear<-apply(data.bin.sum.ind.patches,2,sum) #sum of individuals over all cells/patches and bins per year
  
  for(i in 1:nyears){data.bin.sum.ind.patches[,i]<-data.bin.sum.ind.patches[,i]/sum.ind.peryear[i]}
  
  #data.bin.sum.ind.patches... normalized counts
  
  return(list("mean"=data.bin.mean,"counts"=data.bin.sum.ind.patches,"bins"=bins))
}

get.FIT.Height<-function(run){
  
  ivariable<-4
  nyears<-113
  
  variables<-c("sla","longevity","wooddens","height")
  weights<-c("ind","ind","ind","ind")
  weights<-c("mass","mass","mass","mass")
  units<-c("in mm²/mg","in a","in g/cm³","in m")
  scalings<-c(1000*0.455,1,(10^-6)/0.455,1)
  temp.resolution<-1
  bins.minmax<- matrix(c(0.005,0.07,0.2,4,0.7*10^5,6.5*10^5,0,50), nrow=length(variables), ncol=2, byrow = TRUE)
  size<-4
  nbins<-100
  npatches<-200
  nruns<-1
  
  mean.trait<-array(NA,c(nruns,length(variables)))
  
  df.data<-array(NA,c(nruns,nyears))
  df.data.hov<-array(NA,c(nruns,nyears*nbins))
  
  filenames<-gsub("","", list.files(paste(path.FIT.data,run,"/",sep=""),pattern="\\.bin$")) #get all .nc filenames in path.output.FIT
  
  filename<-grep(variables[ivariable], filenames, value = TRUE) 
  filename<-grep(weights[ivariable], filename, value = TRUE) 
  
  file  <- file(paste(path.FIT.data,run,"/",filename,sep=""),"rb")
  data.bin<-readBin(file,numeric(),n =npatches*nbins*nyears, size=size)
  close(file)
  
  dim(data.bin)<-c(npatches,nbins,nyears)
  data.bin.sum.ind<-apply(data.bin,c(1,3),sum) #sum up all indiviuals for each patch and year
  
  bins<-seq(bins.minmax[ivariable,1], bins.minmax[ivariable,2],length.out = 100)*scalings[ivariable]# get bins
  
  data.bin.sum.variable<-apply(data.bin,c(1,3),function(x) t(x)%*%bins)# sum up all trait values (bin*number of indiviuals)
  
  data.bin.mean<-apply(data.bin.sum.variable/data.bin.sum.ind,2,mean) #mean trait values in each year
  
  data.bin.sum.ind.patches<-apply(data.bin,c(2,3),sum)# sum of individuals over all cells/patches per year and bin
  
  sum.ind.peryear<-apply(data.bin.sum.ind.patches,2,sum) #sum of individuals over all cells/patches and bins per year
  
  for(i in 1:nyears){data.bin.sum.ind.patches[,i]<-data.bin.sum.ind.patches[,i]/sum.ind.peryear[i]}
  
  #data.bin.sum.ind.patches... normalized counts
  
  return(list("mean"=data.bin.mean,"counts"=data.bin.sum.ind.patches,"bins"=bins))
}

get.TRY.SLA<-function(){
  path.in.TRY <- 'D:/TRY/'
  path.TRY<-paste(path.in.TRY,"TRY_Maik_latlon.RData",sep="")
  path.lookup<-paste(path.in.TRY,"TRY_Categorical_Traits_Lookup_Table_2011_10_12_JK_csv.csv",sep="")
  
  #TRY<-read.table(path.TRY,sep=",")
  load(path.TRY)
  
  lookup<-read.csv2(path.lookup)
  
  #filter by Plant growth form 3 (tree)
  TRY<-TRY[TRY$V7 %in% c(3), ]
  
  TRY.raw<-TRY
  TRY<-TRY.raw
  
  #BLSG
  lat<-c(45,70)
  lon<-c(-10,40)
  ### filter by lat, lon
  #exclude NAs in lat lon
  TRY<-TRY[complete.cases(TRY[ , 11:12]),]
  TRY<-TRY[TRY$V11 > lat[1], ]
  TRY<-TRY[TRY$V11 < lat[2], ]
  TRY<-TRY[TRY$V12 > lon[1], ]
  TRY<-TRY[TRY$V12 < lon[2], ]
  
  #filter by Plant leaf type 1 (broadleafed)
  TRY<-TRY[TRY$V8 %in% c(1), ]
  #filter by Phenology 1 (deciduous)
  TRY<-TRY[TRY$V9 %in% c(1), ]
  SLA<-TRY$V5[which(TRY$V6==12)]
  SLA.BLSG<-hist(SLA,breaks=100,plot=F)
  
  #BLEG
  TRY<-TRY.raw
  
  lat<-c(30,45)
  lon<-c(-10,40)
  ### filter by lat, lon
  #exclude NAs in lat lon
  TRY<-TRY[complete.cases(TRY[ , 11:12]),]
  TRY<-TRY[TRY$V11 > lat[1], ]
  TRY<-TRY[TRY$V11 < lat[2], ]
  TRY<-TRY[TRY$V12 > lon[1], ]
  TRY<-TRY[TRY$V12 < lon[2], ]
  #filter by Plant leaf type 1 (broadleafed)
  TRY<-TRY[TRY$V8 %in% c(1), ]
  #filter by Phenology 1 (evergreen)
  TRY<-TRY[TRY$V9 %in% c(2), ]
  SLA<-TRY$V5[which(TRY$V6==12)]
  SLA.BLEG<-hist(SLA,breaks=100,plot=F)
  
  TRY<-TRY.raw
  
  lat<-c(45,70)
  lon<-c(-10,40)
  ### filter by lat, lon
  #exclude NAs in lat lon
  TRY<-TRY[complete.cases(TRY[ , 11:12]),]
  TRY<-TRY[TRY$V11 > lat[1], ]
  TRY<-TRY[TRY$V11 < lat[2], ]
  TRY<-TRY[TRY$V12 > lon[1], ]
  TRY<-TRY[TRY$V12 < lon[2], ]
  
  #filter by Plant leaf type 1 (needleleaved)
  TRY<-TRY[TRY$V8 %in% c(2), ]
  #filter by Phenology 1 (evergreen)
  TRY<-TRY[TRY$V9 %in% c(2), ]
  SLA<-TRY$V5[which(TRY$V6==12)]
  SLA.NLEG<-hist(SLA,breaks=100,plot=F)
  
  return(list("BLSG_counts"=SLA.BLSG$counts/sum(SLA.BLSG$counts),
              "BLSG_breaks"=SLA.BLSG$mids,
              "BLEG_counts"=SLA.BLEG$counts/sum(SLA.BLEG$counts),
              "BLEG_breaks"=SLA.BLEG$mids,
              "NLEG_counts"=SLA.NLEG$counts/sum(SLA.NLEG$counts),
              "NLEG_breaks"=SLA.NLEG$mids
  ))
}

get.TRY.WD<-function(){
  path.in.TRY <- 'D:/Maik/Documents/TRY/'
  path.TRY<-paste(path.in.TRY,"TRY_Maik_latlon.RData",sep="")
  path.lookup<-paste(path.in.TRY,"TRY_Categorical_Traits_Lookup_Table_2011_10_12_JK_csv.csv",sep="")
  
  #TRY<-read.table(path.TRY,sep=",")
  load(path.TRY)
  
  lookup<-read.csv2(path.lookup)
  #filter by Plant growth form 3 (tree)
  lat<-c(30,70)
  lon<-c(-10,40)
  ### filter by lat, lon
  #exclude NAs in lat lon
  TRY<-TRY[complete.cases(TRY[ , 11:12]),]
  TRY<-TRY[TRY$V11 > lat[1], ]
  TRY<-TRY[TRY$V11 < lat[2], ]
  TRY<-TRY[TRY$V12 > lon[1], ]
  TRY<-TRY[TRY$V12 < lon[2], ]
  
  TRY<-TRY[TRY$V7 %in% c(3), ]
  
  TRY.raw<-TRY
  
  TRY<-TRY.raw
  #filter by Plant leaf type 1 (broadleafed)
  TRY<-TRY[TRY$V8 %in% c(1), ]
  #filter by Phenology 1 (deciduous)
  TRY<-TRY[TRY$V9 %in% c(1), ]
  WD<-TRY$V5[which(TRY$V6==4)]
  WD.BLSG<-hist(WD,breaks=80,plot=F)
  
  TRY<-TRY.raw
  #filter by Plant leaf type 1 (broadleafed)
  TRY<-TRY[TRY$V8 %in% c(1), ]
  #filter by Phenology 1 (evergreen)
  TRY<-TRY[TRY$V9 %in% c(2), ]
  WD<-TRY$V5[which(TRY$V6==4)]
  WD.BLEG<-hist(WD,breaks=100,plot=F)
  
  TRY<-TRY.raw
  #filter by Plant leaf type 1 (needleleaved)
  TRY<-TRY[TRY$V8 %in% c(2), ]
  #filter by Phenology 1 (evergreen)
  TRY<-TRY[TRY$V9 %in% c(2), ]
  WD<-TRY$V5[which(TRY$V6==4)]
  WD.NLEG<-hist(WD,breaks=50,plot=F)
  
  return(list("BLSG_counts"=WD.BLSG$counts/sum(WD.BLSG$counts),
              "BLSG_breaks"=WD.BLSG$mids,
              "BLEG_counts"=WD.BLEG$counts/sum(WD.BLEG$counts),
              "BLEG_breaks"=WD.BLEG$mids,
              "NLEG_counts"=WD.NLEG$counts/sum(WD.NLEG$counts),
              "NLEG_breaks"=WD.NLEG$mids)
  )
}

cal.error<-function(error.name,obs,sim){
  error<-NA
  if(error.name=="ME"){
    error<-sum(obs-sim)/length(obs)
    error<-round(error,digits = 1)
  }
  if(error.name=="MAE"){
    error<-sum(abs(obs-sim))/length(obs)
  }
  if(error.name=="RMSE"){
    error<-sqrt(sum((obs-sim)^2)/length(obs))
  }
  if(error.name=="MPE"){
    er<-(obs-sim)/obs
    er[!is.finite(er)] <- NA
    error<-100*sum(er,na.rm=T)/length(obs)
  }
  if(error.name=="MAPE"){
    er<-(obs-sim)/obs
    er[!is.finite(er)] <- NA
    error<-100*sum(abs(er),na.rm=T)/length(obs)
  }
  if(error.name=="RMSPE"){
    er<-(obs-sim)/obs
    er[!is.finite(er)] <- NA
    error<-100*sqrt(sum((er)^2,na.rm=T)/length(obs))
  }
  if(error.name=="cor"){
    error<-cor(obs,sim)
    error<-round(error,digits = 2)
  }
  if(error.name=="NMSE"){ #normalized MSE, Kelly et al. 2013
    error<-sum((sim-obs)^2)/sum((obs-mean(obs,na.rm=T))^2)
    error<-round(error,digits = 2)
  }
  if(error.name=="NME"){ #normalized ME, Kelly et al. 2013
    error<-sum(abs(sim-obs))/sum(abs(obs-mean(obs,na.rm=T)))
  }
  return(error)
  
}
#ls<-get.FIT.Height("run267")

#plot(ls$bins,ls$counts[,100])
#plot(ls$mean)

#flux<-get.FLUX("Hainich")

#plot(flux$mean_gpp)
#plot(flux$sd_gpp)

#plot(flux$years,flux$sum_gpp)