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(raster)
library(ncdf4)
library(rgdal)
library(ggplot2)
library(gridExtra)
library(ggpubr)
library(grid)
library(cowplot)

#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<-"D:/Data/LPJmL-Fit/Europe/historical/" #set path of data
path.out<- "D:/Data/LPJmL-Fit/Europe/historical/" #set path of graphics output

#### load grid ####
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)
data_vegc<-apply(data_vegc[,,istart:iend],c(1,2),mean,na.rm=T)
data_vegc<-data_vegc/1000

#create desert
desert<-data_vegc<vegc_thres

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)

#### load vegc ####

filenames<-list.files(path.in) #get all .nc filenames in path.output.FIT
#filter filename that contain outputname
filenames<-grep("vegc", filenames, value = TRUE) 
filenames<-grep("1951", filenames, value = TRUE) 
data_vegc<-read.mult.runs(filenames,varid = "VegC")

data_vegc<-apply(data_vegc[,,istart:iend],c(1,2),mean,na.rm=T)
data_vegc<-data_vegc/1000
data_vegc<-data_vegc*(2/3)




#### load validation height data ####

spdf <- readOGR(dsn = "D:/Data/Height/All_countries_hmax_hlorey", layer = "All_countries_hmax_hlorey")
df.spdf <- as.data.frame(spdf)
df.spdf <- df.spdf[, c("d_lat", "d_lon", "new_hlorey")]
df.spdf$lat<-round((df.spdf$d_lat + 0.25) * 2) / 2 - 0.25
df.spdf$lon<-round((df.spdf$d_lon + 0.25) * 2) / 2 - 0.25

df.height.val <-aggregate(df.spdf$new_hlorey, by=list(df.spdf$lat,df.spdf$lon),
                   FUN=mean, na.rm=TRUE)
names(df.height.val)<-c("lat","lon","val_data")
df.val_sd <-aggregate(df.spdf$new_hlorey, by=list(df.spdf$lat,df.spdf$lon),
                      FUN=sd, na.rm=TRUE)
plot(rasterFromXYZ(df.height.val))

df.landcells<-data.frame(land=as.vector(landmask[,]),lon=coord[,1],lat=coord[,2])
df.val.land<-merge(df.landcells,df.height.val,all.x=T)
df.val.land<-df.val.land[df.val.land$land,]

#if na, try to build mean of surrounding pixels
i.na<-which(is.na(df.val.land$val_data))
#if or all pixels
i.na<-which(df.val.land$land)
df.val.land.agg<-df.val.land
for(i in i.na){
  lat<-df.val.land$lat[i]
  lon<-df.val.land$lon[i]
  neighbors<-data.frame(lat=c(lat+0.5,lat+0.5,lat,lat-0.5,lat-0.5,lat-0.5,lat,lat+0.5),
                        lon=c(lon,lon+0.5,lon+0.5,lon+0.5,lon,lon-0.5,lon-0.5,lon-0.5))
  df.val.land.agg$val_data[i]<-mean(merge(df.val.land,neighbors)$val_data,na.rm=T)
}




#### load simulated height ####

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
var_mean[which(landmask[,]&is.na(var_mean),arr.ind=T)]<-0
df.height.sim <- data.frame(data=as.vector(var_mean),lon=coord[,1],lat=coord[,2]) #create dataframe
df.height.sim<-df.height.sim[complete.cases(df.height.sim),]

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




#create nopft
variable<-"fpc"
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

nodata<-nopft_pft|noheight_pft
for(ipft in 1:4){
  nodata[,,ipft]<-nodata[,,ipft]|desert
}

df_nodata<-data.frame(nodata=as.vector(!apply(!nodata[,,],c(1,2),any)),lon=coord[,1],lat=coord[,2])

#ggplot(df_nodata,aes(x=lon,y=lat))+geom_raster(aes(fill=nodata))

#### Load Biomass observed ####

load("D:/Data/Biomass/Forest_Aboveground_Biomassv3.rdata")
#r.avi<-t(r.avi)
plot(r.avi)
#r.avi.eu<-crop(r.avi, extent(-11, 36,30, 65))

df.vegc.obs<-as.data.frame(rasterToPoints(r.avi))

df.vegc.obs$lon<-round(df.vegc.obs$x * 4) / 4
df.vegc.obs$lat<-round(df.vegc.obs$y * 4) / 4
df.vegc.obs$avi<-df.vegc.obs$layer/10000*0.455*1000
#df.vegc.obs$avi<-df.vegc.obs$avi*1.3
df.vegc.obs<-df.vegc.obs[,c("lat","lon","avi")]

df.vegc.sim <- data.frame(data=as.vector(data_vegc),lon=coord[,1],lat=coord[,2],
                      nodata_vegc=as.vector(!apply(!nodata,c(1,2),any)))
df.vegc.sim <-subset(df.vegc.sim, !is.na(data_vegc)) #remove water tiles

df.vegc.diff<-merge(df.vegc.sim,df.vegc.obs,all.x=T)
df.vegc.diff$diff<-df.vegc.diff$data-df.vegc.diff$avi

df.vegc.diff$nodata_avi<-(is.na(df.vegc.diff$avi))

df.height.diff<-merge(df.height.sim,df.val.land.agg)
df.height.diff$diff<--df.height.diff$val_data+df.height.diff$data
str(df.height.diff)

df.height.diff<-merge(df.height.diff,df_nodata)
df.vegc.diff<-merge(df.vegc.diff,df_nodata)

df.height.diff[df.height.diff$nodata&!is.na(df.height.diff$data),c("data","diff","val_data")]<-NA
df.vegc.diff[df.vegc.diff$nodata&!is.na(df.vegc.diff$data),c("data","diff","avi")]<-NA

#### Plotting ####

max.vegc<-max(c(df.vegc.diff$data,df.vegc.diff$avi),na.rm=T)
max.height<-max(c(df.height.diff$data,df.height.diff$val_data),na.rm=T)

p<-ggplot(data=df.height.diff, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
p<-p+scale_fill_gradientn(limits=c(0, max.height),colors=col.zscale,name="Mean height simulated [m]",na.value = col.nodata)
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.sim<-p

p<-ggplot(data=df.height.diff, aes(y=lat, x=lon)) + geom_raster(aes(fill=val_data)) +coord_fixed(ratio=1)
p<-p+scale_fill_gradientn(limits=c(0, max.height),colors=col.zscale,name="Mean height, Healey et al. (2015) [m]",na.value = col.nodata)
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.obs<-p

p<-ggplot(data=df.height.diff, aes(y=lat, x=lon)) + geom_raster(aes(fill=diff)) +coord_fixed(ratio=1)
p<-p+scale_fill_gradientn(colours = c("red", "white", "green"), values = scales::rescale(c(min(df.height.diff$diff,na.rm=T),0,max(df.height.diff$diff,na.rm=T))),name="Mean height difference [m]",na.value = col.nodata)
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.diff<-p

p<-ggplot(data=df.vegc.diff, aes(y=lat, x=lon)) + geom_raster(aes(fill=data)) +coord_fixed(ratio=1)
p<-p+scale_fill_gradientn(limits=c(0, max.vegc),colors=col.zscale,name="Mean AGB simulated [kgC/mē]",na.value = col.nodata)
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.sim<-p

p<-ggplot(data=df.vegc.diff, aes(y=lat, x=lon)) + geom_raster(aes(fill=avi)) +coord_fixed(ratio=1)
p<-p+scale_fill_gradientn(limits=c(0, max.vegc),colors=col.zscale,name="Mean AGB, Thurner et al. (2014) [kgC/mē]",na.value = col.nodata)
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.obs<-p

p<-ggplot(data=df.vegc.diff, aes(y=lat, x=lon)) + geom_raster(aes(fill=diff)) +coord_fixed(ratio=1)
p<-p+scale_fill_gradientn(colours = c("red", "white", "green"), values = scales::rescale(c(min(df.vegc.diff$diff,na.rm=T),0,max(df.vegc.diff$diff,na.rm=T))),name="Mean AGB difference [kgC/mē]",na.value = col.nodata)
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.diff<-p


my_guide<-guides(fill = guide_colourbar(title.position="top", title.hjust = 0.5,
                                        barwidth = 20,barheight = 2))
my_theme<-theme(legend.title=element_text(size=30),
                legend.text = element_text(size=25),
                legend.position = "top")


p.all <- plot_grid( p.height.sim+my_theme+my_guide,
                    p.height.obs+my_theme+my_guide,
                    p.height.diff+my_theme+my_guide,
                    p.vegc.sim+my_theme+my_guide,
                    p.vegc.obs+my_theme+my_guide,
                    p.vegc.diff+my_theme+my_guide,
                    
                    labels = c("a)", "b)", "c)","d)","e)","f)"),label_size=30,hjust = 0,
                    ncol = 3)
ggsave(paste(path.out,"Biomass_Height_validation.jpg",sep=""),device = "jpg",dpi = 200,p.all,scale=4,width=18,height=10, units="cm")