source: CONFIG_DEVT/ORCHIDEE_OL_TP/ENSEMBLE/make_eval_figure.R @ 6327

Last change on this file since 6327 was 5570, checked in by aclsce, 3 years ago

Created ORCHIDEE_OL_TP configuration : temporary configuration to be used during prectical session.
File size: 16.8 KB
Line 
1#!/usr/bin/env Rscript
2args = commandArgs(trailingOnly=TRUE)
3
4require(reshape2)
5require(ggplot2)
6require(grid)
7require(gridExtra)
8require(plyr)
9require(dplyr)
10require(data.table)
11
12##############################################################################################
13# initialisation variables
14#############################################################################################
15bo=c(0,0,0,17,18,6,19,8,9,10,20,21,13,14,15,16)
16names=c('','','','','P.Sylvestris','P.Pinaster','PiceaSp','','','QuercusIlex','','BetulaSp','FagusSyl.','QuercusRob.','PopulusSp','')
17leaftype=c('','','','','','','','','','','','Deciduous','Deciduous','Deciduous','Deciduous','')
18monthOrder <- c("janv.","févr.","mars","avril","mai","juin","juil.","août","sept.","oct.","nov.","déc.")
19LatOrderSite <- c('RU-Cok','FI-Kaa','FI-Sod','SE-Deg','SE-Fla','IS-Gun','FI-Hyy','RU-Zot','SE-Sk2','SE-SK1','SE-Nor','UK-Gri','DK-Fou','RU-Fyo','SE-Faj','UK-ESa','UK-EBu','UK-AMo','DK-Lva','DK-Ris','DK-Sor','RU-Ha2','RU-Ha1','RU-Ha3','NL-Lut','IE-Ca1','PL-wet','NL-Loo','NL-Hor','NL-Haa','IE-Dri','NL-Ca1','NL-Lan','UK-Her','NL-Mol','UK-PL3','BE-Bra','DE-Meh','UK-Tad','UK-Ham','DE-Obe','DE-Geb','DE-Hai','DE-Tha','DE-Gri','DE-Kli','BE-Jal','BE-Lon','DE-Wet','BE-Vie','DE-Bay','CZ-BK1','SK-Tat','CZ-wet','FR-Gri','FR-Hes','FR-Fon','HU-Mat','CH-Lae','CH-Oe2','CH-Oe1','CH-Dav','AT-Neu','HU-Bug','IT-Ren','IT-Mal','IT-MBo','IT-Lav','FR-Lq1','FR-Lq2','IT-LMa','IT-PT1','IT-Cas','FR-LBr','IT-Non','FR-Pue','IT-SRo','IT-Lec','IT-Pia','IT-Ro1','IT-Ro2','ES-VDA','IT-Amp','IT-Col','IT-Cpz','IT-Cp2','IT-BCi','ES-LMa','ES-ES1','ES-ES2','PT-Esp','PT-Mi1','PT-Mi2')
20site.sp=list(
21    P.Sylvestris=c('NL-Loo','FI-Hyy','FI-Sod'),
22    P.Pinaster=c('IT-SRo','FR-LBr'),
23    PiceaSp=c('DE-Bay','DE-Tha ','DE-Wet','IT-Ren','RU-Fyo','SE-Nor','SE-Fla'),
24    QuercusIlex=c('FR-Pue','IT-Cpz','PT-Mi1'),
25    BetulaSp=c(''),
26    FagusSyl.=c('FR-Hes','DE-Hai','IT-Col','DK-Sor'),
27    QuercusRob.=c('IT-Non','FR-Fon','IT-Ro1','IT-Ro2'),
28    PopulusSp=c('IT-PT1')
29)
30
31en.mo <- c(
32    janv.="January",
33    févr.="Febuary",
34    mars ="March",
35    avril="April",
36    mai="May",
37    juin="June",
38    juil.="July",
39    août="August",
40    sept.="September",
41    oct.="October",
42    nov.="November",
43    déc.="December")
44
45vol.thres <- list(
46    PiceaSp=list(
47        YC1=list(DE=c(644,886,26,34,645,944,30,36)),
48        YC2=list(DE=c(613,842,22,30,785,958,29,31)),
49        YC3=list(DE=c(579,795,25,32,689,1021,28,34)),
50        YC4=list(DE=c(543,747,24,31,715,1065,27,32)),
51        YC5=list(DE=c(505,698,23,30,746,1117,26,31)),
52        YC6=list(DE=c(467,648,22,29,780,1173,25,30)),
53        YC7=list(DE=c(427,598,21,28,821,1238,24,29)),
54        YC8=list(DE=c(387,547,20,26,872,1320,22,27)),
55        YC9=list(DE=c(346,496,19,25,934,1418,21,26)),       
56        YC10=list(DE=c(304,445,18,24,1011,1541,20,24)),
57        YC11=list(DE=c(262,394,16,22,1115,1701,18,23)),
58        YC12=list(DE=c(220,343,14,20,1266,1939,20,21)),
59        YC13=list(DE=c(177,291,12,18,1496,2309,15,19)),
60        YC14=list(DE=c(131,238,10,15,2361,3911,15,19))       
61    )
62)
63
64############################################################################################
65climate.final <- fread('climate.final.txt')
66climate.final2 <- fread('climate_final.txt', fill=TRUE, na.string=c("Inf","-Inf","NA"))
67print(summary(climate.final2 ))
68
69Info.Site <- fread('site_info.txt')
70
71clim.sum <- NULL; clim.sum2 <-NULL; LAI.sum <- NULL; PS.sum <- NULL; WSM.sum <- NULL
72LAI.doy <- NULL; RDI.sum=NULL; GPP.doy=NULL; tot.sum=NULL
73
74for(i in args[-c(1:6)]){
75
76    Table.Summary <- fread(paste0(args[2],'Table.Summary_',i,'.txt'))
77############################################################################################
78
79new=NULL
80for (j in LatOrderSite ) {
81        if (file.exists(paste0(args[1],j,'_ALL.nc'))==FALSE) next
82        else new <- c(new,j)
83}
84LatOrderSite <- new
85print(LatOrderSite)
86###########################################################################################
87#                    REshape data to montly value only
88############################################################################################
89
90if(args[5]=="M"){
91
92attach(Table.Summary)
93ob <-melt(tapply(GPP,list(Month,Var2,Site),quantile,probs=0.75), value.name ='GPP75')
94ob2 <-melt(tapply(GPP,list(Month,Var2,Site),quantile,probs=0.25), value.name ='GPP25')
95obf <- merge(ob,ob2,c('Var1','Var2','Var3'))
96names(obf) <- c('Month','PFT','Site','GPP75','GPP25')
97obf <- subset(obf, PFT %in% c(5,6,7,10,12,13,14,15,16,17,18,19,20,21))
98clim <- merge(obf,climate.final,c('Month','Site'))
99clim$Month <- factor(clim$Month,levels=monthOrder)
100clim$Version <- i
101detach(Table.Summary)
102clim.sum <- rbind(clim.sum,clim)
103
104Table.Sub <- subset(Table.Summary, Month %in% c("juin","juil.","août"))
105attach(Table.Sub)
106ob <-melt(tapply(NPP/GPP,list(Var2,Site),quantile,probs=0.75,na.rm=TRUE), value.name ='Ratio75')
107ob2 <-melt(tapply(NPP/GPP,list(Var2,Site),quantile,probs=0.25,na.rm=TRUE), value.name ='Ratio25')
108obf <- merge(ob,ob2,c('Var1','Var2'))
109names(obf) <- c('PFT','Site','Ratio75','Ratio25')
110obf <- subset(obf, PFT %in% c(5,6,7,10,12,13,14,15,16,17,18,19,20,21))
111clim2 <- merge(obf,climate.final,c('Site'))
112clim2 <- subset(clim2, Month=='mai')
113clim2$Version <- i
114detach(Table.Sub)
115clim.sum2 <- rbind(clim.sum2,clim2)
116
117attach(Table.Summary)
118ob <-melt(tapply(LAI,list(Month,Var2,Site),max,na.rm=TRUE), value.name ='LAI75')
119ob2 <-melt(tapply(LAI,list(Month,Var2,Site),min,na.rm=TRUE), value.name ='LAI25')
120obf <- merge(ob,ob2,c('Var1','Var2','Var3'))
121names(obf) <- c('Month','PFT','Site','LAI75','LAI25')
122obf <- subset(obf, PFT %in% c(5,6,7,10,12,13,14,15,16,17,18,19,20,21))
123LAII <- merge(obf,climate.final,c('Month','Site'))
124LAII$Month <- factor(LAII$Month,levels=monthOrder)
125LAII$Version <- i
126detach(Table.Summary)
127LAI.sum <- rbind(LAI.sum,LAII)
128
129}
130
131if(args[5]=="D"){
132
133Table.Sub <- Table.Summary
134attach(Table.Sub)
135ob <-melt(tapply(LAI,list(Var2,Site,DOY),max,na.rm=TRUE), value.name ='LAIMax')
136ob2 <-melt(tapply(LAI,list(Var2,Site,DOY),min,na.rm=TRUE), value.name ='LAIMin')
137obf <- merge(ob,ob2,c('Var1','Var2','Var3'))
138names(obf) <- c('PFT','Site','DOY','LAIMax','LAIMin')
139LAII <- subset(obf, PFT %in% c(5,6,7,10,12,13,14,15,16,17,18,19,20,21))
140LAII <- merge(obf,climate.final2,c('Site','DOY'))
141LAII <- merge(LAII,Info.Site,'Site')
142LAII$Version <- i
143detach(Table.Sub)
144LAI.doy <- rbind(LAI.doy,LAII)
145
146
147Table.Sub <- Table.Summary
148attach(Table.Sub)
149ob <-melt(tapply(GPP,list(Var2,Site,DOY),max,na.rm=TRUE), value.name ='GPPMax')
150ob2 <-melt(tapply(GPP,list(Var2,Site,DOY),min,na.rm=TRUE), value.name ='GPPMin')
151ob3 <-melt(tapply(GPP,list(Var2,Site,DOY),median,na.rm=TRUE), value.name ='GPPMed')
152obf <- merge(ob,ob2,c('Var1','Var2','Var3'))
153obf <- merge(obf,ob3,c('Var1','Var2','Var3'))
154names(obf) <- c('PFT','Site','DOY','GPPMax','GPPMin','GPPMed')
155LAII <- subset(obf, PFT %in% c(5,6,7,10,12,13,14,15,16,17,18,19,20,21))
156LAII <- merge(obf,climate.final2,c('Site','DOY'))
157LAII <- merge(LAII,Info.Site,c('Site'))
158LAII$Version <- i
159detach(Table.Sub)
160GPP.doy <- rbind(GPP.doy,LAII)
161
162
163TS <- Table.Summary[,c('Year','DOY','PLANT_STATUS','Var2','Site')]
164TS <- subset(TS, (PLANT_STATUS==5 & DOY>150) | PLANT_STATUS==3 )
165TS1 <- TS %>% group_by(Var2,Site,Year) %>% slice(which(PLANT_STATUS==3)[1])
166TS1[,'Phenology'] <- 'budburst'
167TS2 <- TS %>% group_by(Var2,Site,Year) %>% slice(which(PLANT_STATUS==5)[1])
168TS2[,'Phenology'] <- 'leaf senescence'
169TS <- rbind(TS1,TS2)
170
171
172attach(TS)
173ob <-melt(tapply(DOY,list(Phenology,Var2,Site),max,na.rm=TRUE), value.name ='PSmax')
174ob2 <-melt(tapply(DOY,list(Phenology,Var2,Site),min,na.rm=TRUE), value.name ='PSmin')
175obf <- merge(ob,ob2,c('Var1','Var2','Var3'))
176names(obf) <- c('Phenology','PFT','Site','PSmax','PSmin')
177obf <- subset(obf, PFT %in% c(5,6,7,10,12,13,14,15,16,17,18,19,20,21))
178PS <- merge(obf,climate.final,c('Site'))
179PS <- unique(PS[,c('Phenology','PFT','Site','PSmax','PSmin','Lat')])
180PS[,'Site'] <- factor(PS[,'Site'],levels=rev(LatOrderSite))
181PS$Version <- i
182detach(TS)
183PS.sum <- rbind(PS.sum,PS)
184
185tot.tab = merge(Table.Summary,Info.Site,'Site')
186tot.tab$Version <- i
187tot.sum <- rbind(tot.sum,tot.tab)
188
189}
190
191}
192print('Reformating done !')
193
194############################################################################################
195# forestry graphics
196############################################################################################
197if(args[5]=="M"){
198
199pdf(paste0('Forestry_Eval_',args[6],'.pdf'),width=16,height=10)
200for(i in 7){
201    STR.bot <- NULL
202    fm='DE'
203    ycorder <- c('YC1','YC2','YC3','YC4','YC5','YC6','YC7','YC8','YC9','YC10','YC11','YC12','YC13','YC14')
204    for(j in ycorder){
205        vol.bot.min <- vol.thres[[names[[i]]]][[j]][[fm]][1]/10000
206        vol.bot.max <- vol.thres[[names[[i]]]][[j]][[fm]][2]/10000   
207        STR.b <- subset(Table.Summary, Var2==i & WOOD_VOL > vol.bot.min & WOOD_VOL < vol.bot.max & Site %in% c('DE-Hai','DE-Tha','DE-Wet','DE-Bay'))
208        if(nrow(STR.b) == 0) next
209        STR.b[,"Yield"] <- j
210        STR.b[,"ind.min"] <- vol.thres[[names[[i]]]][[j]][[fm]][5]
211        STR.b[,"ind.max"] <- vol.thres[[names[[i]]]][[j]][[fm]][6]
212        STR.b[,"dg.min"] <- vol.thres[[names[[i]]]][[j]][[fm]][3]
213        STR.b[,"dg.max"] <- vol.thres[[names[[i]]]][[j]][[fm]][4]
214        STR.b[,"H.min"] <- vol.thres[[names[[i]]]][[j]][[fm]][7]
215        STR.b[,"H.max"] <- vol.thres[[names[[i]]]][[j]][[fm]][8]
216        STR.bot <- rbind(STR.bot,STR.b)
217    }
218    STR.bot[,"Yield"] <- factor(STR.bot[,"Yield"],levels=ycorder)
219    x2= ggplot(STR.bot,aes(as.factor(Site), IND*10000, linetype=Version))+
220        geom_boxplot(fill="grey50")+ ylim(0,14000)+
221        geom_hline(aes(yintercept = ind.max),linetype = 2, colour='red',size=1.5)+
222        geom_hline(aes(yintercept = ind.min),linetype = 2, colour='red',size=1.5)+
223        xlab("") + ylab("number of individue (tree/ha)") + coord_flip()+ facet_grid(.~Yield)
224    x3= ggplot(STR.bot,aes(as.factor(Site), DIAMETER_DOM*100, linetype=Version))+
225        geom_boxplot(fill="grey50")+ ylim(0,40) +
226        geom_hline(aes(yintercept = dg.max),linetype = 2, colour='red',size=1.5)+
227        geom_hline(aes(yintercept = dg.min),linetype = 2, colour='red',size=1.5)+
228        xlab("") + ylab("dominant diameter (cm)") + coord_flip()+ facet_grid(.~Yield)
229    x4= ggplot(STR.bot,aes(as.factor(Site), HEIGHT_DOM, linetype=Version))+
230        geom_boxplot(fill="grey50")+ ylim(0,25) +
231        geom_hline(aes(yintercept = H.max),linetype = 2, colour='red',size=1.5)+
232        geom_hline(aes(yintercept = H.min),linetype = 2, colour='red',size=1.5)+
233        xlab("") + ylab("dominant height (m)") + coord_flip()+ facet_grid(.~Yield)   
234    grid.arrange(x2,x3,x4,nrow=3)
235}
236dev.off()
237
238##############################################################################################
239# LAI WSM graphics V2
240##############################################################################################
241pdf(paste0('LAI-WSM_eval_',args[7],'.pdf'),width=16,height=10)
242for(i in c(5,7,12,13,14)){
243    subtable <- subset(LAI.sum, PFT==i  & Lat<57 & Lat>=45)
244    subtablebo <- subset(LAI.sum, PFT==bo[i] & Lat>=57 )
245    subtableso <- subset(LAI.sum, PFT==i  & Lat<45 )
246    subtableBE <- subset(subtablebo, PFT==bo[i] & Site %in% site.sp[[names[i]]])
247    subtableTE <- subset(subtable, PFT==i &  Site %in% site.sp[[names[i]]])
248    subtableME <- subset(subtableso, PFT==i & Site %in% site.sp[[names[i]]])
249    x1= ggplot(LAI.sum,aes(as.factor(Lat),y=LAI25, ymin=LAI25,ymax=LAI75,linetype=Version)) +
250        theme_bw()+ ylim(0,10)+
251        geom_crossbar(data= subtablebo, fill="deepskyblue3",alpha=0.3,fatten=1)+
252        geom_crossbar(data= subtable, fill="chartreuse4",alpha=0.3,fatten=1)+
253        geom_crossbar(data= subtableso, fill="darkorange3",alpha=0.3,fatten=1) +
254        geom_crossbar(data= subtableTE, fill="chartreuse4",fatten=1) +
255        geom_crossbar(data= subtableME, fill="darkorange3",fatten=1) +
256        geom_crossbar(data= subtableBE, fill="deepskyblue3",fatten=1)+
257        xlab("") + ylab("LAI in m2/m2") + scale_x_discrete(labels = rev(LatOrderSite)) +
258        facet_grid(.~Month,labeller = labeller(Month=en.mo)) + coord_flip()+
259        theme(panel.spacing = unit(0, "lines"),strip.text.x = element_text(face = "bold",size=18))
260    subtable <- subset(WSM.sum, PFT==i  & Lat<57 & Lat>=45)
261    subtablebo <- subset(WSM.sum, PFT==bo[i] & Lat>=57 )
262    subtableso <- subset(WSM.sum, PFT==i  & Lat<45 )
263    subtableBE <- subset(subtablebo, PFT==bo[i] & Site %in% site.sp[[names[i]]])
264    subtableTE <- subset(subtable, PFT==i &  Site %in% site.sp[[names[i]]])
265    subtableME <- subset(subtableso, PFT==i & Site %in% site.sp[[names[i]]])
266    x2= ggplot(WSM.sum,aes(as.factor(Lat),y=WSM25, ymin=WSM25,ymax=WSM75,linetype=Version)) +theme_bw()+
267        geom_crossbar(data= subtablebo, fill="deepskyblue3",alpha=0.3,fatten=1)+
268        geom_crossbar(data= subtable, fill="chartreuse4",alpha=0.3,fatten=1)+
269        geom_crossbar(data= subtableso, fill="darkorange3",alpha=0.3,fatten=1) +
270        geom_crossbar(data= subtableTE, fill="chartreuse4",fatten=1) +
271        geom_crossbar(data= subtableME, fill="darkorange3",fatten=1) +
272        geom_crossbar(data= subtableBE, fill="deepskyblue3",fatten=1)+
273        xlab("") + ylab("monthly average water stress") + ylim(0,1) + coord_flip() +
274        scale_x_discrete(labels = rev(LatOrderSite)) +
275        facet_grid(.~Month,labeller = labeller(Month=en.mo))
276    grid.arrange(x1,x2,nrow=2,top = textGrob(paste('PFT:',i,'called',names[i], "10% and 90% quantile of a 250 years simulation"),gp=gpar(fontsize=20,font=3)))
277}
278dev.off()
279###################################################################################################
280}
281
282if(args[5]=="D" & args[6]=="SBG"){
283
284######################################################################
285# Global analysis                                                    #
286######################################################################
287
288pdf(paste0('global_eval_',args[7],'.pdf'),width=16,height=10)
289for(i in c(10,13,14,5,6,7,18,19)){
290datasp=subset(GPP.doy,PFT==i & PFT==SPE,Site)
291droplevels(datasp)
292nbsite = unique(datasp)
293nofacet = dim(nbsite)[1]>1
294x1 = ggplot(subset(GPP.doy, PFT==SPE & PFT==i), aes(x=DOY,ymin=GPPMin,ymax=GPPMax,colour=Version)) +
295     geom_ribbon(alpha=0.4,size=0.2)+ xlab("DOY") + ylab("GPP") +
296     {if(nofacet)facet_wrap(.~Site, scales = "free")}+
297     geom_ribbon(aes(ymin=GPP.min,ymax=GPP.max),alpha=0.1,colour='black',size=0.2)+
298     geom_smooth(aes(y=GPP.med),colour='black',linetype=2,size=0.2,span = 0.3)+
299     geom_smooth(aes(y=GPPMed,colour=Version),linetype=2,size=0.2,span = 0.3)
300grid.arrange(x1,nrow=1,top = textGrob(paste('PFT:',i,'called',names[i],"Min and Max of 50 years simulation"),gp=gpar(fontsize=20,font=3)))
301
302x1 = ggplot(subset(LAI.doy, PFT==SPE & PFT==i), aes(x=DOY,ymin=LAIMin,ymax=LAIMax,colour=Version))+
303     geom_ribbon(alpha=0.4,size=0.2)+ xlab("DOY") + ylab("LAI") +
304     geom_line(aes(y=LAImax),linetype=2,colour='black',size=0.2)+
305     geom_line(aes(y=LAImin),linetype=2,colour='black',size=0.2)+
306     {if(nofacet)facet_wrap(.~Site, scales = "free")}
307grid.arrange(x1,nrow=1,top = textGrob(paste('PFT:',i,'called',names[i],"Interannual Min and Max"),gp=gpar(fontsize=20,font=3)))
308
309x1 = ggplot(subset(tot.sum, Var2==SPE & Var2==i), aes(x=as.Date(Date),y=LAI, colour=Version))+
310     ylab("LAI") + geom_smooth(aes(linetype=Version),size=0.2,span = 0.1)+
311     {if(nofacet)facet_wrap(.~Site, scales = "free")}+ ylim(0,10)
312grid.arrange(x1,nrow=1,top = textGrob(paste('PFT:',i,'called',names[i]),gp=gpar(fontsize=20,font=3)))
313
314x1 = ggplot(subset(tot.sum, Var2==SPE & Var2==i), aes(x=as.Date(Date),y=GPP,colour=Version))+
315     ylab("GPP") + geom_smooth(aes(linetype=Version),size=0.2,span = 0.1)+
316     {if(nofacet)facet_wrap(.~Site, scales = "free")}
317grid.arrange(x1,nrow=1,top = textGrob(paste('PFT:',i,'called',names[i]),gp=gpar(fontsize=20,font=3)))
318
319x1 = ggplot(subset(tot.sum, Var2==SPE & Var2==i), aes(x=as.Date(Date),y=IND,colour=Version))+
320     ylab("IND") + geom_line(aes(linetype=Version),size=0.2)+ ylim(0,1.4)+
321     {if(nofacet)facet_wrap(.~Site, scales = "free")}
322x2 = ggplot(subset(tot.sum, Var2==SPE & Var2==i), aes(x=as.Date(Date),y=RDI,colour=Version))+
323     ylab("RDI") + geom_line(aes(linetype=Version),size=0.2)+ ylim(0.25,1)+
324     {if(nofacet)facet_wrap(.~Site, scales = "free")}
325grid.arrange(x1,x2,nrow=2,top = textGrob(paste('PFT:',i,'called',names[i]),gp=gpar(fontsize=20,font=3)))
326
327#If(leaftype(i)=='Deciduous'){
328#   x1 = ggplot(subset(PS.sum,PFT==i),aes(Site, y=PSmin, ymin=PSmin,ymax=PSmax,colour=Version)) +
329#                  geom_crossbar(aes(group=Phenology),fill='black',alpha=0.3,fatten=1)+
330#                  geom_hline(yintercept = 90,linetype = 2)+
331#                  geom_hline(yintercept = 130,linetype = 2)+
332#                  geom_hline(yintercept = 310,linetype = 2)+
333#                  geom_hline(yintercept = 285,linetype = 2)+
334#                  xlab("") + ylab("DOY") + ylim(0,366) + coord_flip()
335#   grid.arrange(x1,nrow=1,top = textGrob(paste('PFT:',i,'called',names[i], "Min and Max of 50 years simulation"),gp=gpar(fontsize=20,font=3)))
336#}
337
338}
339dev.off()
340
341}
342
343if(args[5]=="D" & args[6]=="SFR"){
344
345}
346
347
Note: See TracBrowser for help on using the repository browser.