#Initial config
nsample = 0 #Maximum number of non-NA cells to sample. Set to 0 if you want all cells.
stepx_plot=50 #Step of the X plot axis label

xlabel="Quantile value (mm/day)"
ylabel="Quantile"
title="Daily precipitation @ 5km resolution (1980-2008)"
plotpath="plots/quantile_5+12+50+80_remap5.pdf"
qprobs=c(900:1000/1000)

#Zoom setup
zoom <- c(10, 90, 0.99, 1) #Left, right, bottom, top margins for zoomed area
box  <- c(30, 270, 0.905, 0.971) #Left, right, bottom, top margins for the are where to plot the zoom

#Definition of file variables
filename <- "/home/clima-archive3/afantini/tesi/learning/obs/EURO4M-APGD/EURO4M_1980-2008.nc"
filename[2] <- "/home/netapp-clima-scratch/afantini/ERAINT_12KM/output/daily/masked/masked_EURO4M_5_PR+TAS_Med_CORDEX_STS.1980-2008_remap_12to5.nc"
filename[3] <- "/home/netapp-clima-scratch/afantini/ERAINT_50KM/output/daily/masked/masked_EURO4M_5_PR+TAS_Med_CORDEX_STS.1980-2008_remap_50to12to5.nc"
filename[4] <- "/home/netapp-clima-scratch/afantini/ERAINT_50KM/output/daily/masked/masked_EURO4M_5_PR+TAS_Med_CORDEX_STS.1980-2008_remap_50to5.nc"
filename[5] <- "/home/clima-archive/afantini/ERAIN/cat/masked/cropped/ERAIN_pr_1980-2008_daily_masked_remap_80to50to12to5_cropped.nc"
filename[6] <- "/home/clima-archive/afantini/ERAIN/cat/masked/cropped/ERAIN_pr_1980-2008_daily_masked_remap_80to5_cropped.nc"

name="EURO4M-APGD \n 5km"
name[2]="RegCM4.4 \n 12->5km"
name[3]="RegCM4.4 \n 50->12->5km"
name[4]="RegCM4.4 \n 50->5km"
name[5]="ERA-Interim \n 80->50->12->5km"
name[6]="ERA-Interim \n 80->5km"

var="pr"
var[2]="pr"
var[3]="pr"
var[4]="pr"
var[5]="TP"
var[6]="TP"

gain=1
gain[2]=86400
gain[3]=86400
gain[4]=86400
gain[5]=86400
gain[6]=86400

colors=c("#006600", "#001F5C", "#B26B00", "#FFCC99", "#500000", "#FF3399")

# #Load the necessary libraries
library(raster)
m <- NULL
for (i in 1:length(filename) ) {
  cat(   "\n",format(Sys.time(), "%H:%M:%S" ), "    "   )
  cat("\t\tReading file", i, "of", length(filename),":\n", filename[i], "\n ")
pr <- brick(filename[i], varname=var[i], na.rm=T)


  if (nsample > ncell(pr) || nsample == 0) {
    cat("\t\t\tFile read. Calculating quantiles over all", ncell(pr), "cells...")
    qs <- quantile(getValues(pr)*gain[i], probs=qprobs, na.rm=T, type=8, names=F)
  } else {
    cnsample = nsample
    cat("\t\t\tFile read. Calculating quantiles over", cnsample, "random cells (out of", ncell(pr),"cells) for each timestep...")
    qs <- quantile(sampleRandom(pr, cnsample)*gain[i], probs=qprobs, na.rm=T, type=8, names=F)
  }

cat(" Done!\n")

assign(  paste("df", i, sep=""), data.frame(qprobs, qs, rep(name[i], length(qs)))  )
m[i] <- max(qs)
gc()
}

df_tot <- eval(parse(text=paste("rbind(", paste("df", 1:length(filename), sep="", collapse=", "), ")")))
colnames(df_tot) <- c('p','v','Dataset')

library(ggplot2)
library(grid)#For unit()



plotbig <- ggplot(data=df_tot, aes(x=v, y=p, color=Dataset)) +
geom_line(alpha = 0.7, size=.5) +
geom_point(size=.8) +
coord_cartesian(ylim=c(0.9, 1.001), xlim=c(0,max(m)+2)) +
scale_color_manual(values=colors) +
xlab(xlabel) + ylab(ylabel) +
ggtitle(title) +
theme_bw(base_size=12) +
theme(panel.grid.major = element_line(size = .5, color = "grey"), axis.line = element_line(size=.4, color = "black"), legend.position = c(.85,.5),text = element_text(size=14), legend.key.height=unit(3,"line"))

plotsmall <-  ggplotGrob(
plotbig + coord_cartesian(xlim=c(zoom[1],zoom[2]), ylim=c(zoom[3], zoom[4]+0.0001)) + theme_bw(base_size=10) + theme(legend.position="none", axis.title.x = element_blank(), axis.title.y = element_blank(), title = element_blank(), plot.margin=unit(c(0,1,0,0),"mm"), plot.background = element_rect(colour = "black"))
)




pdf(plotpath, useDingbats=FALSE)
print(
plotbig + annotation_custom(grob=plotsmall, xmin = box[1], xmax = box[2], ymin = box[3], ymax= box[4]) +
annotate("rect", xmin = zoom[1], xmax = zoom[2], ymin = zoom[3], ymax = zoom[4], alpha = 0, colour = "black") +
annotate("segment", x = box[1], xend = zoom[1], y = box[4], yend = zoom[3],colour="black") +
annotate("segment", x = box[2], xend = zoom[2], y = box[4], yend = zoom[4],colour="black")
)
dev.off()