### This script analyses pan-EU trait distributions to TRY data ###
### written by Maik Billing

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(ggmap)
library(data.table)
library(plyr)
library(ncdf4)
library(tables)
library(gridExtra)
library(grid)

path.in.TRY <- 'D:/TRY/'

DataSetID<-ObsID<-ObsDataID<-c()
lat_eu <- c(30, 70) #borders for EU          
lon_eu <- c(-10, 40)  

lat_nhem <- c(30,70)  #borders for Northern Hemisphere
lon_nhem <- c(-180,180)

#### read TRY data ####

setwd(path.in.TRY)
path.TRY<-paste(path.in.TRY,"TRY_3501_2018_final.csv",sep="")
TRY<-fread(path.TRY)

#filter by Plant growth form 3 (tree)
TRY<-TRY[TRY$V11 %in% c(3), ]

TRY.raw<-TRY #"save" raw dataset

### Get a map
map.eu <- get_map(location = c(lon = mean(lon_eu), lat = mean(lat_eu)), zoom = 3,
               maptype = "satellite", source = "google")
mapWorld <- borders("world", colour="gray50", fill="gray50") # create a layer of borders
world_map <- ggplot()+mapWorld

eu_map <- ggmap(map.eu)+
  scale_x_continuous(limits = lon_eu, expand = c(0, 0)) +
  scale_y_continuous(limits = lat_eu, expand = c(0, 0))

# E U R O P E #

##########################################################################

######## S L A ###########

### BLSG ###

TRY<-TRY.raw #load TRY data

lat<-c(45,70)
lon<-c(-10,40)

TRY<-filter.latlon(TRY,lat_eu,lon_eu)
#TRY<-filter.latlon(TRY,lat,lon)
TRY<-filter.broadleaved(TRY)
TRY<-filter.deciduous(TRY)
TRY<-filter.SLA.nopetiole(TRY)

DataSetID<-c(DataSetID,as.numeric(TRY$V1))
ObsID<-c(ObsID,as.numeric(TRY$V3))
ObsDataID<-c(ObsDataID,as.numeric(TRY$V4))

SLA.BLSG<-get.Trait(TRY)
lat.SLA.BLSG<-get.lat(TRY)
lon.SLA.BLSG<-get.lon(TRY)

df<-data.frame(lon=get.lon(TRY),lat=get.lat(TRY)) 
p<-eu_map+geom_point(aes(x = lon, y = lat), data = df,alpha=0.3,color="red")
p<-p+ggtitle("Measured SLA on broadleaved deciduous trees")
p

df<-data.frame(SLA=SLA.BLSG)
p<-ggplot(df, aes(x=SLA))+geom_histogram(binwidth = 1)
p<-p+scale_x_continuous(limits=c(0,60),breaks = pretty(c(0,60), n = 10))
p<-p+ggtitle("SLA on broadleaved deciduous trees, EU")
p<-p+ geom_vline(xintercept=mean(df$SLA), lwd=1, linetype=2, color="black") 
p<-p+geom_text(aes(label=paste("mean:",round(mean(df$SLA),2)," sd:",round(sd(df$SLA),2),sep=""),y=0,x=mean(df$SLA)),
               vjust=1,col='grey70',size=7)
p

#################################

### BLEG ###

TRY<-TRY.raw #load TRY data

lat<-c(30,50)
lon<-c(-10,40)

TRY<-filter.latlon(TRY,lat,lon)
TRY<-filter.broadleaved(TRY)
TRY<-filter.evergreen(TRY)
TRY<-filter.SLA.nopetiole(TRY)

DataSetID<-c(DataSetID,as.numeric(TRY$V1))
ObsID<-c(ObsID,as.numeric(TRY$V3))
ObsDataID<-c(ObsDataID,as.numeric(TRY$V4))

SLA.BLEG<-get.Trait(TRY)
lat.SLA.BLEG<-get.lat(TRY)
lon.SLA.BLEG<-get.lon(TRY)

df<-data.frame(lon=get.lon(TRY),lat=get.lat(TRY)) 
p<-eu_map+geom_point(aes(x = lon, y = lat), data = df,alpha=0.3,color="red")
p<-p+ggtitle("Measured SLA on broadleaved deciduous trees")
p

df<-data.frame(SLA=SLA.BLEG)
p<-ggplot(df, aes(x=SLA))+geom_histogram(binwidth = 0.1)
p<-p+scale_x_continuous(limits=c(0,60),breaks = pretty(c(0,60), n = 10))
p<-p+ggtitle("SLA on broadleaved evergreen trees, EU")
p<-p+ geom_vline(xintercept=mean(df$SLA), lwd=1, linetype=2, color="black") 
p<-p+geom_text(aes(label=paste("mean:",round(mean(df$SLA),2)," sd:",round(sd(df$SLA),2),sep=""),y=0,x=mean(df$SLA)),
               vjust=1,col='grey70',size=7)
p


###########################################

### NLEG ###

TRY<-TRY.raw #load TRY data

lat<-c(50,70)
lon<-c(-10,40)

TRY<-filter.latlon(TRY,lat_eu,lon_eu)
TRY<-filter.needleleaved(TRY)
TRY<-filter.evergreen(TRY)
TRY<-filter.SLA.undefpetiole(TRY)

DataSetID<-c(DataSetID,as.numeric(TRY$V1))
ObsID<-c(ObsID,as.numeric(TRY$V3))
ObsDataID<-c(ObsDataID,as.numeric(TRY$V4))

SLA.NLEG<-get.Trait(TRY)
lat.SLA.NLEG<-get.lat(TRY)
lon.SLA.NLEG<-get.lon(TRY)

df<-data.frame(lon=get.lon(TRY),lat=get.lat(TRY)) 
p<-eu_map+geom_point(aes(x = lon, y = lat), data = df,alpha=0.3,color="red")
p<-p+ggtitle("Measured SLA on broadleaved deciduous trees")
p

df<-data.frame(SLA=SLA.NLEG)
p<-ggplot(df, aes(x=SLA))+geom_histogram(binwidth = 0.1)
p<-p+scale_x_continuous(limits=c(0,60),breaks = pretty(c(0,60), n = 10))
p<-p+ggtitle("SLA on needlleaved evergreen trees, EU")
p<-p+ geom_vline(xintercept=mean(df$SLA), lwd=1, linetype=2, color="black") 
p<-p+geom_text(aes(label=paste("mean:",round(mean(df$SLA),2)," sd:",round(sd(df$SLA),2),sep=""),y=0,x=mean(df$SLA)),
               vjust=1,col='grey70',size=7)
p


##########################################################################

######## W O O D D E N S I T Y ###########

### BLSG ###

TRY<-TRY.raw #load TRY data

lat<-c(30,70)
lon<-c(-10,40)

TRY<-filter.latlon(TRY,lat,lon)
TRY<-filter.broadleaved(TRY)
TRY<-filter.deciduous(TRY)
TRY<-filter.WD(TRY)

DataSetID<-c(DataSetID,as.numeric(TRY$V1))
ObsID<-c(ObsID,as.numeric(TRY$V3))
ObsDataID<-c(ObsDataID,as.numeric(TRY$V4))

WD.BLSG<-get.Trait(TRY)
lat.WD.BLSG<-get.lat(TRY)
lon.WD.BLSG<-get.lon(TRY)

df<-data.frame(lon=get.lon(TRY),lat=get.lat(TRY)) 
p<-eu_map+geom_point(aes(x = lon, y = lat), data = df,alpha=0.3,color="red")
p<-p+ggtitle("Measured WD on broadleaved deciduous trees")
p

df<-data.frame(WD=WD.BLSG)
p<-ggplot(df, aes(x=WD))+geom_histogram(binwidth = 0.01)
p<-p+scale_x_continuous(limits=c(0,1.2),breaks = pretty(c(0,1.2), n = 10))
p<-p+ggtitle("WD on broadleaved deciduous trees, EU")
p<-p+ geom_vline(xintercept=mean(df$WD), lwd=1, linetype=2, color="black") 
p<-p+geom_text(aes(label=paste("mean:",round(mean(df$WD),2)," sd:",round(sd(df$WD),2),sep=""),y=0,x=mean(df$WD)),
               vjust=1,col='grey70',size=7)
p

#################################

### BLEG ###

TRY<-TRY.raw #load TRY data

lat<-c(30,50)
lon<-c(-10,40)

TRY<-filter.latlon(TRY,lat,lon)
TRY<-filter.broadleaved(TRY)
TRY<-filter.evergreen(TRY)
TRY<-filter.WD(TRY)

DataSetID<-c(DataSetID,as.numeric(TRY$V1))
ObsID<-c(ObsID,as.numeric(TRY$V3))
ObsDataID<-c(ObsDataID,as.numeric(TRY$V4))

WD.BLEG<-get.Trait(TRY)
lat.WD.BLEG<-get.lat(TRY)
lon.WD.BLEG<-get.lon(TRY)

df<-data.frame(lon=get.lon(TRY),lat=get.lat(TRY)) 
p<-eu_map+geom_point(aes(x = lon, y = lat), data = df,alpha=0.3,color="red")
p<-p+ggtitle("Measured WD on broadleaved evergreen trees")
p

df<-data.frame(WD=WD.BLEG)
p<-ggplot(df, aes(x=WD))+geom_histogram(binwidth = 0.01)
p<-p+scale_x_continuous(limits=c(0,1.2),breaks = pretty(c(0,1.2), n = 10))
p<-p+ggtitle("WD on broadleaved evergreen trees, EU")
p<-p+ geom_vline(xintercept=mean(df$WD), lwd=1, linetype=2, color="black") 
p<-p+geom_text(aes(label=paste("mean:",round(mean(df$WD),2)," sd:",round(sd(df$WD),2),sep=""),y=0,x=mean(df$WD)),
               vjust=1,col='grey70',size=7)
p

###########################################

### NLEG ###

TRY<-TRY.raw #load TRY data

lat<-c(50,70)
lon<-c(-10,40)

TRY<-filter.latlon(TRY,lat_eu,lon_eu)
TRY<-filter.needleleaved(TRY)
TRY<-filter.evergreen(TRY)
TRY<-filter.WD(TRY)

DataSetID<-c(DataSetID,as.numeric(TRY$V1))
ObsID<-c(ObsID,as.numeric(TRY$V3))
ObsDataID<-c(ObsDataID,as.numeric(TRY$V4))

WD.NLEG<-get.Trait(TRY)
lat.WD.NLEG<-get.lat(TRY)
lon.WD.NLEG<-get.lon(TRY)

df<-data.frame(lon=get.lon(TRY),lat=get.lat(TRY)) 
p<-eu_map+geom_point(aes(x = lon, y = lat), data = df,alpha=0.3,color="red")
p<-p+ggtitle("Measured WD on needlelevaed evergreen trees")
p

df<-data.frame(WD=WD.NLEG)
p<-ggplot(df, aes(x=WD))+geom_histogram(binwidth = 0.01)
p<-p+scale_x_continuous(limits=c(0,1.2),breaks = pretty(c(0,1.2), n = 10))
p<-p+ggtitle("WD on needlleaved evergreen trees, EU")
p<-p+ geom_vline(xintercept=mean(df$WD), lwd=1, linetype=2, color="black") 
p<-p+geom_text(aes(label=paste("mean:",round(mean(df$WD),2)," sd:",round(sd(df$WD),2),sep=""),y=0,x=mean(df$WD)),
               vjust=1,col='grey70',size=7)
p



##### LOAD FIT DATA #####

startyear_simulation<-1951

startyear<-2004
endyear<-2013

istart<-startyear-startyear_simulation+1
iend<-endyear-startyear_simulation+1

path.out<- "D:/Data/LPJmL-Fit/Europe/historical/"
path.in<-"D:/Data/LPJmL-Fit/Europe/historical/"

### load grid ###
variable<-"vegc"

filenames<-list.files(path.in) #get all .nc filenames in path.output.FIT
#filter filename that contain outputname
filenames<-grep(variable, filenames, value = TRUE) 
filename<-grep("1951", filenames, value = TRUE) 

ncfile<-nc_open(paste(path.in,filename,sep=""))
lats<-ncvar_get(ncfile,varid="latitude")
lons<-ncvar_get(ncfile,varid="longitude")
coord<-cbind(rep(lons,length(lats)),rep(lats,each=length(lons)))
nc_close(ncfile)

#### SLA ####

variable<-"sla"
varid<-"mass_sla"
binsname<-"sla"
zscalename<-"SLA [mm²/mg]\n"
variablenames<-"SLA"
titlename<-"Mean SLA"
scaling<-1/(455)

filenames<-list.files(path.in) #get all .nc filenames in path.output.FIT
#filter filename that contain outputname
filename<-grep(variable, filenames, value = TRUE) 
filename<-grep("51", filename, value = TRUE) 
filename<-grep("mass", filename, value = TRUE) 

data_raw<-read.mult.runs(filename, varid)
bins_SLA<-get.bins(filename,binsname)/scaling

data_pft_SLA<-apply(data_raw[,,,,istart:iend],c(1,2,3,4),mean,na.rm=T) #mean over years



#### WD ####

variable<-"wooddens"
varid<-"mass_wooddens"
binsname<-"wooddensity"
zscalename<-"WD [g/cm³]\n"
titlename<-"Mean WD"
variablenames<-"WD"
scaling<-10^6*0.455

filenames<-list.files(path.in) #get all .nc filenames in path.output.FIT
#filter filename that contain outputname
filename<-grep(variable, filenames, value = TRUE) 
filename<-grep("51", filename, value = TRUE) 
filename<-grep("mass", filename, value = TRUE) 

data_raw<-read.mult.runs(filename, varid)
bins_WD<-get.bins(filename,binsname)/scaling

data_pft_WD<-apply(data_raw[,,,,istart:iend],c(1,2,3,4),mean,na.rm=T) #mean over years




#### PLOTTING ####

col.fill.FIT<-"red"
col.fill.TRY<-"blue"
alpha<-0.4
xlim.SLA<-c(0,25)
xlim.WD<-c(0.25,1.2)

PFTs.txt<-c("BL-S","BL-E","NL")

label.try<-"TRY-Database"
label.lpj<-"LPJmL-FIT     "
xlabel<-"Specific Leaf Area [mm²/mg]"
ylabel<-"Relative Count"

### SLA ###

### BLSG ###

lat.TRY<-lat.SLA.BLSG
lon.TRY<-lon.SLA.BLSG
trait.TRY<-SLA.BLSG

dat<-data_pft_SLA[,,,1]

dat_hist<-rep(0,100)
for(i in 1:length(lon.TRY)){
  dat_hist<-dat_hist+dat[which.min(abs(lons - lon.TRY[i]))[1],which.min(abs(lats - lat.TRY[i]))[1],]
}

str(bins)
vec_hist<-c()
for(i in 1:100){
  vec_hist<-c(vec_hist,rep(bins_SLA[i],dat_hist[i]))
}

df<-data.frame(trait=c(vec_hist,trait.TRY), 
               type=c(rep(label.lpj,length(vec_hist)),rep(label.try,length(trait.TRY))))
mu <- ddply(df, "type", summarise, grp.mean=mean(trait))
sd <- ddply(df, "type", summarise, grp.sd=sd(trait))
mu
sd
p<-ggplot(df,aes(x=trait,color=type,fill=type,y=..density..))+
  geom_histogram(alpha=alpha,position="identity")
p<-p+geom_vline(data=mu, aes(xintercept=grp.mean, color=type),
           linetype="dashed",show_guide = F)
p<-p+xlim(xlim.SLA)+ylab(ylabel)+xlab(xlabel)
p<-p+ggtitle(PFTs.txt[1])
p<-p+theme(legend.title=element_blank())
p.SLA.BLSG<-p
p

### BLEG ###

lat.TRY<-lat.SLA.BLEG
lon.TRY<-lon.SLA.BLEG
trait.TRY<-SLA.BLEG

dat<-data_pft_SLA[,,,2]

dat_hist<-rep(0,100)
for(i in 1:length(lon.TRY)){
  dat_hist<-dat_hist+dat[which.min(abs(lons - lon.TRY[i]))[1],which.min(abs(lats - lat.TRY[i]))[1],]
}

str(bins)
vec_hist<-c()
for(i in 1:100){
  vec_hist<-c(vec_hist,rep(bins_SLA[i],dat_hist[i]))
}

df<-data.frame(trait=c(vec_hist,trait.TRY), 
               type=c(rep(label.lpj,length(vec_hist)),rep(label.try,length(trait.TRY))))
mu <- ddply(df, "type", summarise, grp.mean=mean(trait))
sd <- ddply(df, "type", summarise, grp.sd=sd(trait))
mu
sd
p<-ggplot(df,aes(x=trait,color=type,fill=type,y=..density..))+
  geom_histogram(alpha=alpha,position="identity")
p<-p+geom_vline(data=mu, aes(xintercept=grp.mean, color=type),
                linetype="dashed",show_guide = F)
p<-p+xlim(xlim.SLA)+ylab(ylabel)+xlab(xlabel)
p<-p+ggtitle(PFTs.txt[2])
p<-p+theme(legend.title=element_blank())
p.SLA.BLEG<-p


### NLEG ###

lat.TRY<-lat.SLA.BLEG
lon.TRY<-lon.SLA.BLEG
trait.TRY<-SLA.BLEG

dat<-apply(data_pft_SLA[,,,3:4],1:3,sum)

dat_hist<-rep(0,100)
for(i in 1:length(lon.TRY)){
  dat_hist<-dat_hist+dat[which.min(abs(lons - lon.TRY[i]))[1],which.min(abs(lats - lat.TRY[i]))[1],]
}

str(bins)
vec_hist<-c()
for(i in 1:100){
  vec_hist<-c(vec_hist,rep(bins_SLA[i],dat_hist[i]))
}

df<-data.frame(trait=c(vec_hist,trait.TRY), 
               type=c(rep(label.lpj,length(vec_hist)),rep(label.try,length(trait.TRY))))
mu <- ddply(df, "type", summarise, grp.mean=mean(trait))
mu <- ddply(df, "type", summarise, grp.mean=mean(trait))
sd <- ddply(df, "type", summarise, grp.sd=sd(trait))
mu
sd

p<-ggplot(df,aes(x=trait,color=type,fill=type,y=..density..))+
  geom_histogram(alpha=alpha,position="identity")
p<-p+geom_vline(data=mu, aes(xintercept=grp.mean, color=type),
                linetype="dashed",show_guide = F)
p<-p+xlim(xlim.SLA)+ylab(ylabel)+xlab(xlabel)
p<-p+ggtitle(PFTs.txt[3])
p<-p+theme(legend.title=element_blank())
p.SLA.NLEG<-p
#p


### WD ###

xlabel<-"Wood Density [g/cm³]"

### BLSG ###

lat.TRY<-lat.WD.BLSG
lon.TRY<-lon.WD.BLSG
trait.TRY<-WD.BLSG

dat<-data_pft_WD[,,,1]

dat_hist<-rep(0,100)
for(i in 1:length(lon.TRY)){
  dat_hist<-dat_hist+dat[which.min(abs(lons - lon.TRY[i]))[1],which.min(abs(lats - lat.TRY[i]))[1],]
}

str(bins)
vec_hist<-c()
for(i in 1:100){
  vec_hist<-c(vec_hist,rep(bins_WD[i],dat_hist[i]))
}

df<-data.frame(trait=c(vec_hist,trait.TRY), 
               type=c(rep(label.lpj,length(vec_hist)),rep(label.try,length(trait.TRY))))
mu <- ddply(df, "type", summarise, grp.mean=mean(trait))
mu <- ddply(df, "type", summarise, grp.mean=mean(trait))
sd <- ddply(df, "type", summarise, grp.sd=sd(trait))
mu
sd
p<-ggplot(df,aes(x=trait,color=type,fill=type,y=..density..))+
  geom_histogram(alpha=alpha,position="identity")
p<-p+geom_vline(data=mu, aes(xintercept=grp.mean, color=type),
                linetype="dashed",show_guide = F)
p<-p+xlim(xlim.WD)+ylab(ylabel)+xlab(xlabel)
p<-p+theme(legend.title=element_blank())
p.WD.BLSG<-p
p

### BLEG ###

lat.TRY<-lat.WD.BLEG
lon.TRY<-lon.WD.BLEG
trait.TRY<-WD.BLEG

dat<-data_pft_WD[,,,2]

dat_hist<-rep(0,100)
for(i in 1:length(lon.TRY)){
  dat_hist<-dat_hist+dat[which.min(abs(lons - lon.TRY[i]))[1],which.min(abs(lats - lat.TRY[i]))[1],]
}

str(bins)
vec_hist<-c()
for(i in 1:100){
  vec_hist<-c(vec_hist,rep(bins_WD[i],dat_hist[i]))
}

df<-data.frame(trait=c(vec_hist,trait.TRY), 
               type=c(rep(label.lpj,length(vec_hist)),rep(label.try,length(trait.TRY))))
mu <- ddply(df, "type", summarise, grp.mean=mean(trait))
mu <- ddply(df, "type", summarise, grp.mean=mean(trait))
sd <- ddply(df, "type", summarise, grp.sd=sd(trait))
mu
sd
p<-ggplot(df,aes(x=trait,color=type,fill=type,y=..density..))+
  geom_histogram(alpha=alpha,position="identity")
p<-p+geom_vline(data=mu, aes(xintercept=grp.mean, color=type),
                linetype="dashed",show_guide = F)
p<-p+xlim(xlim.WD)+ylab(ylabel)+xlab(xlabel)
#p<-p+ggtitle("Broadl. Evergreen")
p<-p+theme(legend.title=element_blank())
p.WD.BLEG<-p
#p

### NLEG ###

lat.TRY<-lat.WD.NLEG
lon.TRY<-lon.WD.NLEG
trait.TRY<-WD.NLEG

dat<-apply(data_pft_WD[,,,3:4],1:3,sum)

dat_hist<-rep(0,100)
for(i in 1:length(lon.TRY)){
  dat_hist<-dat_hist+dat[which.min(abs(lons - lon.TRY[i]))[1],which.min(abs(lats - lat.TRY[i]))[1],]
}

str(bins)
vec_hist<-c()
for(i in 1:100){
  vec_hist<-c(vec_hist,rep(bins_WD[i],dat_hist[i]))
}

df<-data.frame(trait=c(vec_hist,trait.TRY), 
               type=c(rep(label.lpj,length(vec_hist)),rep(label.try,length(trait.TRY))))
mu <- ddply(df, "type", summarise, grp.mean=mean(trait))
mu <- ddply(df, "type", summarise, grp.mean=mean(trait))
sd <- ddply(df, "type", summarise, grp.sd=sd(trait))
mu
sd
p<-ggplot(df,aes(x=trait,color=type,fill=type,y=..density..))+
  geom_histogram(alpha=alpha,position="identity")
p<-p+geom_vline(data=mu, aes(xintercept=grp.mean, color=type),
                linetype="dashed",show_guide = F)
p<-p+xlim(xlim.WD)+ylab(ylabel)+xlab(xlabel)
#p<-p+ggtitle("Needlel. Evergreen")
p<-p+theme(legend.title=element_blank())
p.WD.NLEG<-p
#p

library(ggpubr)

mytheme<- theme(legend.text=element_text(size=40),
                legend.key.size = unit(1, 'cm'),
                legend.position = "bottom")

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)}

mylegend<-g_legend(p.SLA.BLSG+mytheme)

p.commonleg <- grid.arrange(arrangeGrob(p.SLA.BLSG + theme(legend.position="none"),
                               p.SLA.BLEG + theme(legend.position="none"),
                               p.SLA.NLEG + theme(legend.position="none"),
                               p.WD.BLSG + theme(legend.position="none"),
                               p.WD.BLEG + theme(legend.position="none"),
                               p.WD.NLEG + theme(legend.position="none"),
                               nrow=2),
                   mylegend, nrow=2,heights=c(10, 1))
ggsave(paste(path.out,"validate_Traits_commonleg.jpg",sep=""),device = "jpg",dpi = 300,p.commonleg,scale=3,width=20,height=9, units="cm")


#p.commonleg<-ggarrange(p.SLA.BLSG+mytheme, p.SLA.BLEG+mytheme,p.SLA.NLEG+mytheme,
#                      p.WD.BLSG+mytheme,p.WD.BLEG+mytheme,p.WD.NLEG+mytheme,
#                      ncol=3, nrow=2, 
#                      common.legend = TRUE, legend="bottom",align="none")

#ggsave(paste(path.out,"validate_Traits_commonleg.jpg",sep=""),device = "jpg",dpi = 300,p.commonleg,scale=3,width=20,height=9, units="cm")

#ls.nolegends<-list(p.SLA.BLSG+theme(legend.position="none"),p.SLA.BLEG+theme(legend.position="none"),p.SLA.NLEG+theme(legend.position="none"),
#                   p.WD.BLSG+theme(legend.position="none"),p.WD.BLEG+theme(legend.position="none"),p.WD.NLEG+theme(legend.position="none"))


#ls<-list(p.SLA.BLSG,p.SLA.BLEG,p.SLA.NLEG,
#         p.WD.BLSG,p.WD.BLEG,p.WD.NLEG)



#plots.nolegend<-grid.arrange(grobs=ls.nolegends,
#                            ncol=3, nrow=2)
#plots<-grid.arrange(grobs=ls,ncol=3, nrow=2)
#ggsave(paste(path.out,"validate_Traits_legend.jpg",sep=""),device = "jpg",dpi = 600,plots,scale=3,width=20,height=7, units="cm")



#ggsave(paste(path.out,"validate_Traits.jpg",sep=""),device = "jpg",dpi = 300,plots.nolegend,scale=3,width=20,height=7, units="cm")



table(DataSetID)

write.table(table(DataSetID),file="D:/TRY/DataSetID.csv",dec=",",row.names = F)

#write.table(unique(ObsID),file="D:/TRY/ObsID.csv",dec=",",row.names = F)

#write.table(unique(ObsDataID),file="D:/TRY/ObsDataID.csv",dec=",",row.names = F)