### This script analyses Biomass, Height, FPC, GPP, Traits (SLA, WD) and Climate gradients of a pan-EU run ###
### written by Maik Billing

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

library(ncdf4)
library(ggplot2)
library(gridExtra)
library(ggpubr)
library(cowplot)
library(scales)

#import user-defined functions and graphics parameters
source("C:/Users/maikb/Documents/Dissertation/R/Github/graphics_param.R")
source("C:/Users/maikb/Documents/Dissertation/R/Github/functions.R")


#set startyear of simulation data and timeframe for evaluation
startyear_simulation<-1951
startyear<-2004
endyear<-2013

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

path.in<-"C:/Users/maikb/Documents/Dissertation/Data/LPJmL-Fit/Europe/historical/" #set path of data
path.out<- "C:/Users/maikb/Documents/Dissertation/Data/LPJmL-Fit/Europe/historical/" #set path of graphics output

mk.sites<-c("Seitseminen","Laegern","Dundo")
mk.lats<-c(61.75,47.25,44.75)
mk.lons<-c(23.25,8.25,14.75)

#### Load climate data ####
path.input<-"C:/Users/maikb/Documents/Dissertation/Data/LPJmL-Fit/Input/"
load(file=paste(path.input,"climate_EU.Rdata",sep=""))

#### Load climate grid ####
headersize<-51
file  <- file(paste(path.input,"grid_europe.clm",sep=""),"rb") # read binary mode
seek(file, where = headersize, origin="start")
x<-readBin(file,integer(),n = 2*5095, size=2)/100
lon.clm <- x[c(1:5095)*2-1]
lat.clm <- x[c(1:5095)*2]
close(file)

intersect(which(lat.clm==52.75),which(lon.clm==25.25))

#set everything above latmax to NA
temp.annual.mean[which(lat.clm>latmax),]<-NA
prec.sum[which(lat.clm>latmax),]<-NA
MCWD[which(lat.clm>latmax),]<-NA
swdown.annual.mean[which(lat.clm>latmax),]<-NA
pet.annual.sum[which(lat.clm>latmax),]<-NA

#create dataframe with climate data
df.clm<-data.frame(MCWD=apply(MCWD[,istart:iend],1,mean),temp=apply(temp.annual.mean[,istart:iend],1,mean),lon=lon.clm,lat=lat.clm,
                   prec=apply(prec.sum[,istart:iend],1,mean),
                   PET=apply(pet.annual.sum[,istart:iend],1,mean))



#### Get coordinates and landcells of simulation data ####
filenames<-list.files(path.in) #get all .nc filenames in path.output.FIT
#filter filename that contain outputname
filenames<-grep("vegc", filenames, value = TRUE) 
filename<-grep("1951", filenames, value = TRUE) 
filename<-grep("_1", filenames, value = TRUE) 

ncfile<-nc_open(paste(path.in,filename,sep=""))
lats<-ncvar_get(ncfile,varid="latitude")
lons<-ncvar_get(ncfile,varid="longitude")
data_vegc<-ncvar_get(ncfile, varid="VegC")
nc_close(ncfile)

coord<-cbind(rep(lons,length(lats)),rep(lats,each=length(lons)))
landcells<-which(!is.na(data_vegc),arr.ind = T)
landmask<-!is.na(data_vegc)
threslat<-length(lats[which(lats<latmax)]) #indice of latitude which corresponds to latmax

#### Plot VegC and extract cells with Vegc<vegc_thres (deserts) ####

variable<-"vegc"

filenames<-list.files(path.in) #get all .nc filenames in path.output.FIT

filenames<-grep(variable, filenames, value = TRUE) #filter filename that contain outputname
filenames<-grep(startyear_simulation, filenames, value = TRUE) 

data_vegc<-read.mult.runs(filenames,varid = "VegC") #read all VegC data of all runs
data_vegc<-apply(data_vegc[,,istart:iend],c(1,2),mean,na.rm=T) #mean over evaluation periode
data_vegc<-data_vegc/1000 #g to kg 

#mask water tiles
df <- data.frame(data=as.vector(data_vegc),lon=coord[,1],lat=coord[,2]) #create dataframe
df <-subset(df, !is.na(data_vegc)) #remove water tiles

p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
p<-p+scale_fill_gradientn(limits=c(0, 23),colors=col.zscale,name="Mean Vegetation carbon [kg/m²]",na.value = col.zscale[length(col.zscale)])
p<-p+ theme(panel.background = element_rect(fill = col.plot.bg))
p<-p+ylab("Latitude [°]")+xlab("Longitude [°]")
p<-p+guides(fill=guide_colourbar(barwidth = 2,barheight = 20))
p<-p+ylim(latmin,latmax)
p #plot
p.vegc<-p

#extract "desert" cells
desert<-data_vegc<vegc_thres

#### Plot Height and extract cells with height < height_thes ####
variable<-"height"
varid<-"mass_height"
binsname<-"height"
scaling<-1

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(startyear_simulation, filename, value = TRUE) 
filenames<-grep("mass", filename, value = TRUE) 

bins<-get.bins(filenames,binsname)
data_height<-read.mult.runs(filenames,varid = varid)

data_pft<-apply(data_height[,,,,istart:iend],c(1,2,3,4),mean,na.rm=T) #mean over years
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))) #calculate mean

df <- data.frame(data=as.vector(var_mean),lon=coord[,1],lat=coord[,2]) #create dataframe
df <-subset(df, landmask) #remove water tiles

p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
p<-p+scale_fill_gradientn(limits=c(0, 30),colors=col.zscale,name="Mean Height [m]",na.value = col.zscale[1])
p<-p+ theme(panel.background = element_rect(fill = col.plot.bg))
p<-p+ylab("Latitude [°]")+xlab("Longitude [°]")
p<-p+guides(fill=guide_colourbar(barwidth = 2,barheight = 20))
p<-p+ylim(latmin,latmax)
p #plot
p.height<-p

data_bins<-apply(data_pft,c(1,2,4),mean,na.rm=T) #mean over bins
noheight_pft<-data_bins<max(which(bins<height_thres)) #extract no_height cells by pft

#### Plot VegC & Height together ####

mytheme2<-theme(legend.text = element_text(size=30),
                        legend.position="top",
                legend.title = element_text(size=30))
fillbar<-guides(fill=guide_colourbar(barwidth = 20,barheight = 2,title.position="top", title.vjust= -20,direction = "horizontal"))

plot.vegc.height<-grid.arrange(grobs=list(p.vegc+mytheme2+fillbar,p.height+mytheme2+fillbar), ncol=2, nrow=1)

p.all<-ggarrange(p.vegc+mytheme2+fillbar, p.height+mytheme2+fillbar, ncol=2, nrow=1, 
                 common.legend = F,labels=c("a)","b)"),align="none",font.label = list(size = 35))
ggsave(paste(path.out,"vegc_height.jpg",sep=""),device = "jpg",dpi = 300,p.all,scale=3,width=18,height=7, units="cm")


#### Plot FPC and extract cells with nopft (FPC < fpc_thres)####

variable<-"fpc"
#zscalename<-"Cover of Plant Functional Type"
zscalename<-""
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=""))
data_fpc<-ncvar_get(ncfile, varid="FPC")
nc_close(ncfile)

data_pft<-apply(data_fpc[,,2:5,istart:iend],c(1,2,3),mean,na.rm=T)

nopft_pft<-data_pft<fpc_thres #extract cells with nopft by pft

ls.fpc<-ls.fpc_nolegend<-list()

for(ipft in 1:5){
  data<-apply(data_fpc[,,ipft+1,istart:iend],c(1,2),mean,na.rm=T)
  
  #mask deserts and water tiles
  df <- data.frame(data=as.vector(data),lon=coord[,1],lat=coord[,2])
  df <-subset(df, !is.na(data)) #remove water tiles
  
  if(ipft<5){df_nodata<-noheight_pft[,,ipft]|desert}else{df_nodata<-desert}
  df_nodata<- data.frame(data=as.vector(df_nodata),lon=coord[,1],lat=coord[,2])
  df_nodata<-df_nodata[which(df_nodata$data),]
  
  p<-ggplot(data=df, aes(y=lat, x=lon)) 
  p<-p+geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
  p<-p+scale_fill_gradientn(limits=c(0,1),colors=col.zscale3,name='Cover of Plant Funtional Type',na.value = "green")
  p<-p+geom_tile(data=df_nodata,aes(y=lat,x=lon),color=col.zscale[1],fill=col.zscale[1])
  p<-p+ggtitle(paste(pfts.txt[ipft],sep=""))+ theme(panel.background = element_rect(fill = col.plot.bg))
  p<-p+ylab("Latitude [°]")+xlab("Longitude [°]")
  p<-p+guides(fill=guide_colourbar(barwidth = 10,barheight = 1))
  p<-p+ylim(latmin,latmax)
  p
  #ls.fpc[[ipft]]<-p+colorbar+mytheme2
  ls.fpc_nolegend[[ipft]]<-p+theme(legend.position="none")
  leg<-get_legend(p+guides(fill=guide_colourbar(barwidth = 25,barheight = 3,title.position="top",title.vjust=-10))+
                             theme(legend.position="top")+theme(legend.text = element_text(size=25),
                                                                legend.title = element_text(size=30),
                                                                legend.title.align=0.5))
}

hlay <- rbind(c(1,2,NA),c(3,4,5))

plot.all.fpc_nolegend<-grid.arrange(grobs=ls.fpc_nolegend,ncol=3, nrow=2, layout_matrix=hlay)
plot.all.fpc<-grid.arrange(grobs=ls.fpc,ncol=3, nrow=2)

ls.fpc_onelegend<-ls.fpc_nolegend
ls.fpc_onelegend[[6]]<-ls.fpc_nolegend[[5]]
ls.fpc_onelegend[[5]]<-ls.fpc_nolegend[[4]]
ls.fpc_onelegend[[4]]<-ls.fpc_nolegend[[3]]
ls.fpc_onelegend[[3]]<-leg

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

ggsave(paste(path.out,"fpc_nolegend1.jpg",sep=""),device = "jpg",dpi = 600,plot.all.fpc_nolegend,scale=3,width=10,height=5, units="cm")
ggsave(paste(path.out,"fpc.jpg",sep=""),device = "jpg",dpi = 600,plot.all.fpc,scale=3,width=10,height=5, units="cm")
ggsave(paste(path.out,"fpc_onelegend.jpg",sep=""),device = "jpg",dpi = 300,plot.all.fpc_onelegend,scale=5,width=10,height=5, units="cm")


#### Combine nopft, nobiomass and noheight to nodata ####
nodata<-nopft_pft|noheight_pft
for(ipft in 1:4){
  nodata[,,ipft]<-nodata[,,ipft]|desert
}

#### Plot dominant PFT ####

fpc_max<-apply(data_fpc[,,2:7,istart:iend],c(1,2,3),mean,na.rm=T)
#fpc_max<-maskdeserts(fpc_max,desert)

dim(fpc_max)<-c(dim(fpc_max)[1]*dim(fpc_max)[2],dim(fpc_max)[3])

domPFT<-pfts.txt[max.col(fpc_max)]
#remove no data
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
domPFT[which(df_nodata$data)]<-"No Data"

df_nodata<-na.omit(df_nodata)
df_nodata<-df_nodata[which(df_nodata$data),]

matched<-row.match(df_nodata[,c("lat","lon")], df[,c("lat","lon")])
df$data[matched]<-factor("No Data")

#mask deserts and water tiles
df <- data.frame(data=domPFT,lon=coord[,1],lat=coord[,2])
nonwater_tiles<-df
df <-subset(df, !is.na(data)) #remove water tiles

breaks.txt<-c("Broadl. Summergreen","Broadl. Evergreen","Temp. Needlel. Evergreen","Boreal Needlel. Evergreen","C3 Grass","No Data")

p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
p<-p+scale_fill_manual(name="dominant PFT",breaks = pfts.txt,values=colpfts.txt)
p<-p+ggtitle(paste("Dominant PFT, ",startyear,"-",endyear,sep=""))+ theme(panel.background = element_rect(fill = col.plot.bg))
p<-p+ylab("Latitude [°]")+xlab("Longitude [°]")
print(p) #plot

# climate scatterplot
df<-merge(df,df.clm) #merge with climate data

p<-ggplot(df,aes(x=temp,y=MCWD,color=data))+geom_point(size=2)
p<-p+ggtitle("Dominant PFT")
p<-p+scale_color_manual(name="Dominant PFT",breaks = breaks.txt,values=colpfts.txt)
p<-p+xlab("Mean Temperature [°C]")+xlim(limits.temp)+ylim(limits.MCWD)+ylab("MCWD")
p<-p+ guides(colour = guide_legend(override.aes = list(size=7))) 
p

#### Plot FireC ####

variable<-"firec"

filenames<-list.files(path.in) #get all .nc filenames in path.output.FIT

filenames<-grep(variable, filenames, value = TRUE) #filter filename that contain outputname
filenames<-grep(startyear_simulation, filenames, value = TRUE) 

data_vegc<-read.mult.runs(filenames,varid = "FireC") #read all VegC data of all runs
data_vegc<-apply(data_vegc[,,istart:iend],c(1,2),mean,na.rm=T) #mean over evaluation periode
data_vegc<-data_vegc #g to kg 

#mask water tiles
df <- data.frame(data=as.vector(data_vegc),lon=coord[,1],lat=coord[,2]) #create dataframe
df <-subset(df, !is.na(data_vegc)) #remove water tiles

p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
p<-p+scale_fill_gradientn(colors=col.zscale,name="Mean Fire carbon [gC/m²]",na.value = col.zscale[length(col.zscale)])
p<-p+ theme(panel.background = element_rect(fill = col.plot.bg))
p<-p+ylab("Latitude [°]")+xlab("Longitude [°]")
p<-p+guides(fill=guide_colourbar(barwidth = 2,barheight = 20))
p<-p+ylim(latmin,latmax)
p #plot
p.firec<-p

ggsave(paste(path.out,"FireC.jpg",sep=""),device = "jpg",dpi = 300,p,scale=5,width=7,height=7, units="cm")


#### Plot GPP ####

variable<-"gpp"
varid<-"GPP"
scaling<-1

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

data_raw<-read.mult.runs(filenames,varid = varid)

data_gpp<-array(NA,c(dim(data_raw)[1:2],12,dim(data_raw)[3]/12))

for(i in 1:dim(data_raw)[1]){
  for(j in 1:dim(data_raw)[2]){
    gpp_cell<-data_raw[i,j,]
    dim(gpp_cell)<-c(12,dim(data_raw)[3]/12)
    data_gpp[i,j,,]<-gpp_cell
  }
}

data_gpp_annual<-apply(data_gpp[,,,istart:iend],c(1,2),sum)

data_gpp_annual<-data_gpp_annual/1000/length(c(istart:iend))

df<-data.frame(data=as.vector(data_gpp_annual),lon=coord[,1],lat=coord[,2])
df<-df[!is.na(df$data),]

theme_set(theme_bw(base_size = 30))
theme_update(axis.text=element_text(size=30))
theme_update(plot.title=element_text(size=30,hjust = 0.5))

mytheme2<-theme(legend.text = element_text(size=30),
                legend.position="top",
                legend.title = element_text(size=30))
fillbar<-guides(fill=guide_colourbar(barwidth = 20,barheight = 2,title.position="top", title.vjust= -20,direction = "horizontal"))


p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
p<-p+ theme(panel.background = element_rect(fill = col.plot.bg))
p<-p+scale_fill_gradientn(colors=col.zscale,name="Annual GPP [kgC/m²]",na.value=col.zscale[1])
#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)
#p<-p+ scale_fill_distiller(limits=c(limits.pfts[1],limits.pfts[2]),palette = "Spectral",name=zscalename,na.value=col.nodata)
print(p) #plot

ggsave(paste(path.out,"gpp.jpg",sep=""),device = "jpg",dpi = 300,p+mytheme2+fillbar,scale=3,width=7,height=7, units="cm")



#### Load and plot 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) 
filenames<-grep("mass", filename, value = TRUE) 

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

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

#plot_trait()

col.scale2<-rainbow(10)

point.size<-1.7
mytheme<-theme(axis.text=element_text(size=30),
               axis.title=element_text(size=30))

data_pft<-data_pft_SLA

lim.min<-lim.max<-lim.min.sd<-lim.max.sd<-c()
for(ipft in 1:4){
  data_var<-data_pft[,,,ipft]
  var_mean<-apply(data_var,c(1,2),function(x) sum(x*bins_SLA/sum(x)))
  var_sd<-apply(data_var,c(1,2),function(x) sd_function(x,bins_SLA))
  #coefficient of variation
  var_sd <- var_sd/var_mean
  
  df<-remove_no_data_pft(var_mean,ipft)
  df_sd<-remove_no_data_pft(var_sd,ipft)
  lim.min[ipft]<-min(df$data,na.rm=T)
  lim.max[ipft]<-max(df$data,na.rm=T)
  lim.min.sd[ipft]<-min(df_sd$data,na.rm=T)
  lim.max.sd[ipft]<-max(df_sd$data,na.rm=T)
}
limits.pfts<-c(min(lim.min),max(lim.max))
limits.pfts.sd<-c(min(lim.min.sd),max(lim.max.sd))

limits.pfts<-c(min(lim.min),22)

p.all<-list()
p.all.sd<-list()

for(ipft in 1:4){ #for every pft
  data_var<-data_pft[,,,ipft]
  var_mean<-apply(data_var,c(1,2),function(x) sum(x*bins_SLA/sum(x)))
  var_sd<-apply(data_var,c(1,2),function(x) sd_function(x,bins_SLA))
  #coefficient of variation
  var_sd <- var_sd/var_mean
  
  df<-remove_no_data_pft(var_mean,ipft)
  df_sd<-remove_no_data_pft(var_sd,ipft)
  df$data[which(df$data>limits.pfts[2])]<-limits.pfts[2]
  
  p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
  p<-p+ggtitle(pfts.txt[ipft])+ theme(panel.background = element_rect(fill = col.plot.bg))
  #p<-p+scale_fill_gradientn(limits=c(limits.pfts[1],limits.pfts[2]),colors=col.scale2,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 = 20,barheight = 2))
  p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
  p<-p+ylim(latmin,latmax)
  p<-p+ scale_fill_distiller(limits=c(limits.pfts[1],limits.pfts[2]),palette = "Spectral",name=zscalename,na.value=col.nodata,oob=squish)
  print(p) #plot
  p.all[[ipft]]<-p
  
  p<-ggplot(data=df_sd, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
  p<-p+ggtitle(pfts.txt[ipft])+ theme(panel.background = element_rect(fill = col.plot.bg))
  #p<-p+scale_fill_gradientn(limits=c(limits.pfts[1],limits.pfts[2]),colors=col.scale2,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 = 20,barheight = 2))
  p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
  p<-p+ylim(latmin,latmax)
  p<-p+ scale_fill_distiller(limits=c(limits.pfts.sd[1],limits.pfts.sd[2]),palette = "Spectral",name=paste("COV of ",zscalename,sep=""),na.value=col.nodata,oob=squish)
  print(p) #plot
  p.all.sd[[ipft]]<-p
}

p.all.mean.trait.sla<-ggarrange(p.all[[1]],p.all[[2]],p.all[[3]],p.all[[4]], ncol=2, nrow=2, 
                     common.legend = T,legend="bottom",labels=c("a)","b)","c)","d)"),align="none",font.label = list(size = 30))
p.all.mean.trait.sla

ggsave(paste(path.out,"SLAbyPFT.jpg",sep=""),device = "jpg",dpi = 300,p.all.mean.trait.sla,scale=5,width=8,height=7, units="cm")

p.all.sd.trait.sla<-ggarrange(p.all.sd[[1]],p.all.sd[[2]],p.all.sd[[3]],p.all.sd[[4]], ncol=2, nrow=2, 
                                common.legend = T,legend="bottom",labels=c("a)","b)","c)","d)"),align="none",font.label = list(size = 30))
p.all.sd.trait.sla

ggsave(paste(path.out,"SLAbyPFT_sd.jpg",sep=""),device = "jpg",dpi = 300,p.all.sd.trait.sla,scale=5,width=8,height=7, units="cm")




## all PFTs ##
data_SLA<-apply(data_pft_SLA[,,,],c(1,2,3),sum,na.rm=T) #mean over years
var_mean<-apply(data_SLA,c(1,2),function(x) sum(x*bins_SLA/sum(x)))

df<-remove_no_data(var_mean)
df$data[which(df$data>limits.pfts[2])]<-limits.pfts[2]

p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
p<-p+theme(panel.background = element_rect(fill = col.plot.bg))
#p<-p+scale_fill_gradientn(limits=c(limits.pfts[1],limits.pfts[2]),colors=col.scale2,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 = 20,barheight = 2))
p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
p<-p+ylim(latmin,latmax)
p<-p+ scale_fill_distiller(limits=c(limits.pfts[1],20),palette = "Spectral",name="Specific Leaf Area [mm²/mg]",na.value=col.nodata,oob=squish)
p<-p+theme(legend.position="top")+theme(legend.text = element_text(size=15))
p<-p+guides(fill=guide_colourbar(barwidth = 20,barheight = 1.5,title.position="top", title.vjust= 0))
print(p) #plot

ggsave(paste(path.out,"SLA.jpg",sep=""),device = "jpg",dpi = 300,p,scale=4,width=7,height=7, units="cm")

p<-p+ylim(latmin,latmax)
p<-p+ scale_fill_distiller(limits=c(limits.pfts[1],20),palette = "Spectral",name="Trait",na.value=col.nodata)
p<-p+xlab("")+ylab("")
p<-p+theme(axis.text=element_blank(),
           axis.ticks=element_blank(),
           legend.position = "right",
           legend.direction = "vertical",
           legend.text = element_blank())
p<-p+guides(fill=guide_colourbar(barwidth = 2,barheight = 20))
p<-p+guides(fill=guide_colorbar(title.hjust = 0.5,barwidth=2, barheight=20,
                                title.position = "right",
                                title.theme = element_text(angle=90,size=40)))
p


#### Load Butler Daten ####

#butler.data<-read.csv2("D:/Data/Butler/sla_PNAS_2017.csv")

#str(butler.data)


#### Plot and load 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)
filenames<-grep("mass", filename, value = TRUE) 

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

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

#plot_trait()

data_pft<-data_pft_WD

lim.min<-lim.max<-lim.min.sd<-lim.max.sd<-c()
for(ipft in 1:4){
  data_var<-data_pft[,,,ipft]
  var_mean<-apply(data_var,c(1,2),function(x) sum(x*bins_WD/sum(x)))
  var_sd<-apply(data_var,c(1,2),function(x) sd_function(x,bins_WD))
  #coefficient of variation
  var_sd <- var_sd/var_mean
  df<-remove_no_data_pft(var_mean,ipft)
  df_sd<-remove_no_data_pft(var_sd,ipft)
  lim.min[ipft]<-min(df$data,na.rm=T)
  lim.max[ipft]<-max(df$data,na.rm=T)
  lim.min.sd[ipft]<-min(df_sd$data,na.rm=T)
  lim.max.sd[ipft]<-max(df_sd$data,na.rm=T)
}
limits.pfts<-c(min(lim.min),max(lim.max))
limits.pfts.sd<-c(min(lim.min.sd),max(lim.max.sd))

#limits.pfts<-c(min(lim.min),22)

p.all<-list()
p.all.sd<-list()

for(ipft in 1:4){ #for every pft
  data_var<-data_pft[,,,ipft]
  var_mean<-apply(data_var,c(1,2),function(x) sum(x*bins_WD/sum(x)))
  var_sd<-apply(data_var,c(1,2),function(x) sd_function(x,bins_WD))
  #coefficient of variation
  var_sd <- var_sd/var_mean
  
  df<-remove_no_data_pft(var_mean,ipft)
  df_sd<-remove_no_data_pft(var_sd,ipft)
  df$data[which(df$data>limits.pfts[2])]<-limits.pfts[2]
  
  p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
  p<-p+ggtitle(pfts.txt[ipft])+ theme(panel.background = element_rect(fill = col.plot.bg))
  #p<-p+scale_fill_gradientn(limits=c(limits.pfts[1],limits.pfts[2]),colors=col.scale2,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 = 20,barheight = 2))
  p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
  p<-p+ylim(latmin,latmax)
  p<-p+ scale_fill_distiller(limits=c(limits.pfts[1],limits.pfts[2]),palette = "Spectral",name=zscalename,na.value=col.nodata,oob=squish)
  print(p) #plot
  p.all[[ipft]]<-p
  
  p<-ggplot(data=df_sd, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
  p<-p+ggtitle(pfts.txt[ipft])+ theme(panel.background = element_rect(fill = col.plot.bg))
  #p<-p+scale_fill_gradientn(limits=c(limits.pfts[1],limits.pfts[2]),colors=col.scale2,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 = 25,barheight = 2))
  p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
  p<-p+ylim(latmin,latmax)
  p<-p+ scale_fill_distiller(limits=c(limits.pfts.sd[1],limits.pfts.sd[2]),palette = "Spectral",name=paste("COV of ",zscalename,sep=""),na.value=col.nodata,oob=squish)
  print(p) #plot
  p.all.sd[[ipft]]<-p
}

p.all.mean.trait.wd<-ggarrange(p.all[[1]],p.all[[2]],p.all[[3]],p.all[[4]], ncol=2, nrow=2, 
                                common.legend = T,legend="bottom",labels=c("a)","b)","c)","d)"),align="none",font.label = list(size = 30))

p.all.mean.trait.wd

ggsave(paste(path.out,"WDbyPFT.jpg",sep=""),device = "jpg",dpi = 300,p.all.mean.trait.wd,scale=5,width=8,height=7, units="cm")

p.all.sd.trait.wd<-ggarrange(p.all.sd[[1]],p.all.sd[[2]],p.all.sd[[3]],p.all.sd[[4]], ncol=2, nrow=2, 
                              common.legend = T,legend="bottom",labels=c("a)","b)","c)","d)"),align="none",font.label = list(size = 30))
p.all.sd.trait.wd

ggsave(paste(path.out,"WDbyPFT_sd.jpg",sep=""),device = "jpg",dpi = 300,p.all.sd.trait.wd,scale=5,width=8,height=7, units="cm")


## all PFTs ##
data_WD<-apply(data_pft_WD[,,,],c(1,2,3),sum,na.rm=T) #mean over years
var_mean<-apply(data_WD,c(1,2),function(x) sum(x*bins_WD/sum(x)))

df<-remove_no_data(var_mean)
df$data[which(df$data>limits.pfts[2])]<-limits.pfts[2]

p<-ggplot(data=df, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
p<-p+theme(panel.background = element_rect(fill = col.plot.bg))
#p<-p+scale_fill_gradientn(limits=c(limits.pfts[1],limits.pfts[2]),colors=col.scale2,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 = 20,barheight = 2))
p<-p+xlab("Longitude [°]")+ylab("Latitude [°]")
p<-p+ylim(latmin,latmax)
p<-p+ scale_fill_distiller(limits=c(limits.pfts[1],limits.pfts[2]),palette = "Spectral",name="Wood Density [g/cm²]",na.value=col.nodata,oob=squish)
p<-p+theme(legend.position="top")+theme(legend.text = element_text(size=15))
p<-p+guides(fill=guide_colourbar(barwidth = 20,barheight = 1.5,title.position="top", title.vjust= 0))
print(p) #plot

ggsave(paste(path.out,"WD.jpg",sep=""),device = "jpg",dpi = 300,p,scale=3,width=7,height=7, units="cm")


#### Plot climate-trait-pft plots ####

limits.temp<-c(-1,22)

### climate scatter plots by PFT and precipitation ###

data_pft<-data_pft_SLA

for(ipft in 1:4){
  data_var<-data_pft[,,,ipft]
  var_mean<-apply(data_var,c(1,2),function(x) sum(x*bins_SLA/sum(x)))
  df<-remove_no_data_pft(var_mean,ipft)
  df$PFT<-pfts.txt[ipft]
  if(ipft==1){
    df_final<-df
  }else{
    df_final<-rbind(df_final,df)
  }
}

df_final<-df_final[which(!is.na(df_final$data)),]
df_final<-df_final[which(!is.na(df_final$temp)),]

#df_plot<-df_final[df_final$PFT==pfts.txt[1],]

df_final$PFT <- as.factor(df_final$PFT)
df_final$PFT <- factor(df_final$PFT,levels=unique(df_final$PFT))

p<-ggplot(df_final,aes(x=temp,y=prec,color=data,shape=PFT))+geom_point(size=4)
p<-p+scale_color_distiller(palette = "Spectral",name="SLA [mm²/mg]")
#p<-p+scale_color_gradientn(colors=col.scale1,name="SLA [mm²/mg]",na.value = col.scale1[length(col.scale1)])
p<-p+scale_shape_manual(values=c(15,17,4,1))
#p<-p+scale_color_distiller(palette = "Spectral")
p<-p+mytheme
p<-p+xlab("MAT [°C]")+xlim(limits.temp)+ylab("MAP [mm]")
p<-p+theme(legend.title=element_text(size=30),
           legend.text = element_text(size=25),
           legend.position = "top")
p<-p+guides(color=guide_colourbar(barwidth = 20,barheight = 2,title.hjust = -0.4,title.vjust=1),
            shape=F)
print(p)
p.clm.prec.temp.sla<-p

data_pft<-data_pft_WD

for(ipft in 1:4){
  data_var<-data_pft[,,,ipft]
  var_mean<-apply(data_var,c(1,2),function(x) sum(x*bins_WD/sum(x)))
  df<-remove_no_data_pft(var_mean,ipft)
  df$PFT<-pfts.txt[ipft]
  if(ipft==1){
    df_final<-df
  }else{
    df_final<-rbind(df_final,df)
  }
}

df_final<-df_final[which(!is.na(df_final$data)),]
df_final$PFT <- as.factor(df_final$PFT)
df_final$PFT <- factor(df_final$PFT,levels=unique(df_final$PFT))

#df_plot<-df_final[df_final$PFT==pfts.txt[1],]

p<-ggplot(df_final,aes(x=temp,y=prec,color=data,shape=PFT))+geom_point(size=4)
p<-p+scale_color_distiller(palette = "Spectral",name="WD [g/cm³]")
#p<-p+scale_color_gradientn(colors=col.scale1,name="WD [g/cm³]",na.value = col.scale1[length(col.scale1)])
p<-p+scale_shape_manual(values=c(15,17,4,1))
#p<-p+scale_color_distiller(palette = "Spectral")
p<-p+mytheme
p<-p+xlab("MAT [°C]")+xlim(limits.temp)+ylab("MAP [mm]")
p<-p+theme(legend.title=element_text(size=30),
      legend.text = element_text(size=25),
      legend.position = "top")

p
p_pft_leg<-p+guides(color=F,shape = guide_legend(title="PFT:",override.aes = list(size=7)))
p_pft_leg<-p_pft_leg+theme(legend.title=element_text(size=30),
                           legend.text = element_text(size=30),
                           legend.key.size = unit(2,units="cm"),
                           legend.position = "bottom")
leg_pft <- cowplot::get_legend(p_pft_leg)

p<-p+guides(color=guide_colourbar(barwidth = 20,barheight = 2,title.hjust = -0.4,title.vjust=1),
            shape = F)

print(p)
p.clm.prec.temp.wd<-p

#### Annotate Sites ####
anno.size<-8
#seitseminen
isite<-1
index<-which(df.clm$lat %in% mk.lats[isite] & df.clm$lon %in% mk.lons[isite])
mat.mk<-df.clm$temp[index]
map.mk<-df.clm$prec[index]
p.clm.prec.temp.sla.anno<-p.clm.prec.temp.sla+annotate("segment", x = mat.mk, xend = mat.mk, y = map.mk, yend = map.mk-400,
         colour = "black")+
              annotate("text", x = mat.mk, y = map.mk-500, label = mk.sites[isite],size=anno.size)
p.clm.prec.temp.sla.anno
#laegern
isite<-2
index<-which(df.clm$lat %in% mk.lats[isite] & df.clm$lon %in% mk.lons[isite])
mat.mk<-df.clm$temp[index]
map.mk<-df.clm$prec[index]
p.clm.prec.temp.sla.anno<-p.clm.prec.temp.sla.anno+ annotate("segment", x = mat.mk, xend = mat.mk+2, y = map.mk, yend = map.mk+850,colour = "black")+
                                                    annotate("text", x = mat.mk+2, y = map.mk+950, label = mk.sites[isite],size=anno.size)
p.clm.prec.temp.sla.anno
#dundo
isite<-3
index<-which(df.clm$lat %in% mk.lats[isite] & df.clm$lon %in% mk.lons[isite])
mat.mk<-df.clm$temp[index]
map.mk<-df.clm$prec[index]
p.clm.prec.temp.sla.anno<-p.clm.prec.temp.sla.anno+ annotate("segment", x = mat.mk, xend = mat.mk, y = map.mk, yend = map.mk+450,colour = "black")+
                                                    annotate("text", x = mat.mk, y = map.mk+550, label = mk.sites[isite],size=anno.size)
p.clm.prec.temp.sla.anno


# add the legend to the row we made earlier. Give it one-third of the width
# of one plot (via rel_widths).
p.all <- plot_grid( p.clm.prec.temp.sla.anno,p.clm.prec.temp.wd,
                labels = c("a)", "b)", ""),label_size=30,hjust = 0,
                ncol = 2)
p.all<-plot_grid( p.all, leg_pft, nrow = 2,rel_heights = c(1,0.2))
p.all

ggsave(paste(path.out,"clm_prec_temp_all.jpg",sep=""),device = "jpg",dpi = 1000,p.all,scale=3,width=20,height=10, units="cm")



### climate scatter plots by PFT and precipitation FOR EACH PFT ###

data_pft<-data_pft_SLA

for(ipft in 1:4){
  data_var<-data_pft[,,,ipft]
  var_mean<-apply(data_var,c(1,2),function(x) sum(x*bins_SLA/sum(x)))
  df<-remove_no_data_pft(var_mean,ipft)
  df$PFT<-pfts.txt[ipft]
  if(ipft==1){
    df_final<-df
  }else{
    df_final<-rbind(df_final,df)
  }
}

df_final<-df_final[which(!is.na(df_final$data)),]
df_final<-df_final[which(!is.na(df_final$temp)),]

#df_plot<-df_final[df_final$PFT==pfts.txt[1],]

df_final$PFT <- as.factor(df_final$PFT)
df_final$PFT <- factor(df_final$PFT,levels=unique(df_final$PFT))

p<-ggplot(df_final,aes(x=temp,y=prec,color=data))+geom_point(size=2)
p<-p+scale_color_distiller(palette = "Spectral",name="SLA [mm²/mg]")
p<-p+facet_wrap(df_final$PFT)
#p<-p+scale_color_gradientn(colors=col.scale1,name="SLA [mm²/mg]",na.value = col.scale1[length(col.scale1)])
#p<-p+scale_shape_manual(values=c(15,17,4,1))
#p<-p+scale_color_distiller(palette = "Spectral")
p<-p+mytheme
p<-p+xlab("MAT [°C]")+xlim(limits.temp)+ylab("MAP [mm]")
p<-p+theme(legend.title=element_text(size=30),
           legend.text = element_text(size=25))
p<-p+guides(color=guide_colourbar(barwidth = 20,barheight = 2))
print(p)

p.clm.prec.temp.sla<-p

data_pft<-data_pft_WD

for(ipft in 1:4){
  data_var<-data_pft[,,,ipft]
  var_mean<-apply(data_var,c(1,2),function(x) sum(x*bins_WD/sum(x)))
  df<-remove_no_data_pft(var_mean,ipft)
  df$PFT<-pfts.txt[ipft]
  if(ipft==1){
    df_final<-df
  }else{
    df_final<-rbind(df_final,df)
  }
}

df_final<-df_final[which(!is.na(df_final$data)),]
df_final$PFT <- as.factor(df_final$PFT)
df_final$PFT <- factor(df_final$PFT,levels=unique(df_final$PFT))

#df_plot<-df_final[df_final$PFT==pfts.txt[1],]

p<-ggplot(df_final,aes(x=temp,y=prec,color=data))+geom_point(size=2)
p<-p+scale_color_distiller(palette = "Spectral",name="WD [g/cm³]")
p<-p+facet_wrap(df_final$PFT)
#p<-p+scale_color_gradientn(colors=col.scale1,name="SLA [mm²/mg]",na.value = col.scale1[length(col.scale1)])
#p<-p+scale_shape_manual(values=c(15,17,4,1))
#p<-p+scale_color_distiller(palette = "Spectral")
p<-p+mytheme
p<-p+xlab("MAT [°C]")+xlim(limits.temp)+ylab("MAP [mm]")
p<-p+theme(legend.title=element_text(size=30),
           legend.text = element_text(size=25))
p<-p+guides(color=guide_colourbar(barwidth = 20,barheight = 2))
print(p)

p.clm.prec.temp.wd<-p


ggsave(paste(path.out,"clm_prec_temp_sla.jpg",sep=""),device = "jpg",dpi = 1000,p.clm.prec.temp.sla,scale=2,width=20,height=10, units="cm")
ggsave(paste(path.out,"clm_prec_temp_wd.jpg",sep=""),device = "jpg",dpi = 1000,p.clm.prec.temp.wd,scale=2,width=20,height=10, units="cm")

p.all <- plot_grid( p.clm.prec.temp.sla+mytheme2,p.clm.prec.temp.wd+mytheme2,
                    labels = c("a)", "b)"),label_size=30,hjust = 0,
                    ncol = 2)
#p.all<-plot_grid( p.all, leg_pft, nrow = 2,rel_heights = c(1,0.2))
p.all

ggsave(paste(path.out,"clm_prec_temp_all.jpg",sep=""),device = "jpg",dpi = 1000,p.all,scale=3,width=20,height=10, units="cm")



#### Old climate scatter plots ####

theme_set(theme_bw(base_size = 35))
theme_update(axis.text=element_text(size=33))
theme_update(plot.title=element_text(size=33,hjust = 0.5))
limits.temp<-c(-2,22)

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

data_var_SLA<-apply(data_pft_SLA[,,,],c(1,2,3),sum,na.rm=T) #sum over all pfts
var_mean_SLA<-apply(data_var_SLA,c(1,2),function(x) sum(x*bins_SLA/sum(x,na.rm=T)))
var_shannon_SLA<-apply(data_var_SLA,c(1,2),function(x) shannon_ind(x))

### mean variable ###
#mask deserts and water tiles

df<-remove_no_data(var_mean_SLA)

### 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=c(4.5,20),colors=col.scale1,name=zscalename,na.value = col.scale1[length(col.scale1)])
#p<-p+scale_color_distiller(palette = "Spectral")
p<-p+mytheme
p<-p+xlab("Mean Annual Temperature [°C]")+xlim(limits.temp)+ylim(limits.MCWD)+ylab("Minimum Cumulative\n Water Deficit")
p<-p+guides(color=guide_colourbar(barwidth = 2,barheight = 20))
print(p)
p.clm.SLA<-p

p<-ggplot(df,aes(x=temp,y=PET,color=data))+geom_point(size=point.size)
p<-p<-p+scale_color_gradientn(limits=c(4.5,20),colors=col.scale1,name=zscalename,na.value = col.scale1[length(col.scale1)])
#p<-p+scale_color_distiller(palette = "Spectral")
p<-p+mytheme
p<-p+xlab("Mean Annual Temperature [°C]")+xlim(limits.temp)+ylim(limits.MCWD)+ylab("Annual PET [mm]")
p<-p+guides(color=guide_colourbar(barwidth = 2,barheight = 20))
print(p)
p.clm.SLA<-p

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

data_var_WD<-apply(data_pft_WD[,,,],c(1,2,3),sum,na.rm=T) #sum over all pfts
var_mean_WD<-apply(data_var_WD,c(1,2),function(x) sum(x*bins_WD/sum(x,na.rm=T)))
var_shannon_WD<-apply(data_var_WD,c(1,2),function(x) shannon_ind(x))

### mean variable ###
#mask deserts and water tiles

df<-remove_no_data(var_mean_WD)

### 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=c(0.5,1.3),colors=col.scale1,name=zscalename,na.value = col.scale1[length(col.scale1)])
p<-p+mytheme
p<-p+xlab("Mean Annual Temperature [°C]")+xlim(limits.temp)+ylim(limits.MCWD)+ylab("Minimum Cumulative\n Water Deficit")
p<-p+guides(color=guide_colourbar(barwidth = 2,barheight = 20))
print(p)
p.clm.WD<-p

mytheme2<-theme(legend.position="top")+theme(legend.text = element_text(size=25))
colorbar<-guides(color=guide_colourbar(barwidth = 20,barheight = 2,title.position="top", title.vjust= 10))

plot.all.clm<-grid.arrange(grobs=list(p.clm.SLA+mytheme2+colorbar,p.clm.WD+mytheme2+colorbar), ncol=2, nrow=1)

p.all.clm<-ggarrange(p.clm.SLA+mytheme2+colorbar, p.clm.WD+mytheme2+colorbar, ncol=2, nrow=1, 
                 common.legend = F,labels=c("a)","b)"),align="none",font.label = list(size = 30))

plot.all.clm_nolegend<-grid.arrange(grobs=list(p.clm.SLA+theme(legend.position="none"),p.clm.WD+theme(legend.position="none")),ncol=2, nrow=1)

ggsave(paste(path.out,"CLM_nolegend.jpg",sep=""),device = "jpg",dpi = 600,plot.all.clm_nolegend,scale=3,width=20,height=7, units="cm")
ggsave(paste(path.out,"CLM.jpg",sep=""),device = "jpg",dpi = 300,p.all.clm,scale=3,width=20,height=7, units="cm")



### Shannon Index ###


limits.shannon<-c(min(c(remove_no_data(var_shannon_SLA)$data,remove_no_data(var_shannon_WD)$data)),
                  max(c(remove_no_data(var_shannon_SLA)$data,remove_no_data(var_shannon_WD)$data)))

df<-remove_no_data(var_shannon_SLA)
df_nodata[which(!df_nodata$data),]<-NA
df_nodata<-na.omit(df_nodata)

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.zscale3,name="Shannon Index",na.value=col.nodata)
p<-p+geom_tile(data=df_nodata,aes(y=lat,x=lon),color="grey85",fill="grey85")
p<-p+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+scale_colour_identity(name = 'the fill', guide = 'legend',labels = c('m1'))
p<-p+ylim(latmin,latmax)
print(p) #plot
p.shannon.SLA<-p


df<-remove_no_data(var_shannon_WD)
df_nodata[which(!df_nodata$data),]<-NA
df_nodata<-na.omit(df_nodata)

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.zscale3,name="Shannon Index",na.value=col.nodata)
p<-p+geom_tile(data=df_nodata,aes(y=lat,x=lon),color="grey85",fill="grey85")
p<-p+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+scale_colour_identity(name = 'the fill', guide = 'legend',labels = c('m1'))
p<-p+ylim(latmin,latmax)
print(p) #plot
p.shannon.WD<-p


plot.all.shannon<-grid.arrange(grobs=list(p.shannon.SLA,p.shannon.WD), ncol=2, nrow=1)

plot.all.shannon_nolegend<-grid.arrange(grobs=list(p.shannon.SLA+theme(legend.position="none")+ggtitle("Specific Leaf Area"),
                                                   p.shannon.WD+theme(legend.position="none")+ggtitle("Wood Density")),ncol=2, nrow=1)

plot.all.shannon_onelegend<-grid.arrange(grobs=list(p.shannon.SLA+theme(legend.position="none")+ggtitle("Specific Leaf Area"),
                                                   p.shannon.WD+theme(legend.position="none")+ggtitle("Wood Density"),
                                                   get_legend(p)),ncol=3, nrow=1,
                                         widths=c(0.4,0.4,0.2))


ggsave(paste(path.out,"shannon_nolegend.jpg",sep=""),device = "jpg",dpi = 600,plot.all.shannon_nolegend,scale=3,width=20,height=7, units="cm")
ggsave(paste(path.out,"shannon.jpg",sep=""),device = "jpg",dpi = 600,plot.all.shannon,scale=3,width=20,height=7, units="cm")
ggsave(paste(path.out,"shannon_onelegend.jpg",sep=""),device = "jpg",dpi = 300,plot.all.shannon_onelegend,scale=3,width=20,height=7, units="cm")