Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
1 ; ****************************************************
2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
5 ;************************************************
11 ;------------------------------------------------------
12 ; edit table.html of current model for movel1_vs_model2
14 if (isvar("compare")) then
15 html_name2 = compare+"/table.html"
16 html_new2 = html_name2 +".new"
19 ;------------------------------------------------------
20 ; edit table.html for current model
22 html_name = model_name+"/table.html"
23 html_new = html_name +".new"
25 ;---------------------------------
26 ; get model data: landfrac and area
28 film_l = "lnd_"+ model_grid + ".nc"
29 fm_l = addfile (dirs+film_l,"r")
31 landfrac = fm_l->landfrac
36 ; change area from km**2 to m**2
39 ;-----------------------------
40 ; take into account landfrac
42 area0 = area0 * landfrac
44 ;--------------------------------------------
47 fm = addfile (dirm+film11,"r")
49 if (BGC .eq. "cn") then
53 datamod0 = data1(0,:,:)
54 datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)
57 if (BGC .eq. "casa") then
61 datamod0 = data1(0,:,:)
62 datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:)
72 ;------------------------------------------------
73 ; read amazon mask data
75 dir_m = diro + "biomass/"
76 fil_m = "amazon_mask_"+ model_grid + ".nc"
77 fm = addfile (dir_m+fil_m,"r")
79 mask_amazon0 = fm->mask_amazon
83 ;------------------------------------------------
86 ob_name = "LC15 Amazon Biomass"
88 dir_b = diro + "biomass/"
89 fil_b = "amazon_biomass_"+model_grid+".nc"
90 fo = addfile (dir_b+fil_b,"r")
98 ;************************************************
99 ; Units for these variables are:
102 ; We want to convert these to kgC/m2
103 ; ha = 100m*100m = 10,000 m2
104 ; MgC/ha*1000/10,000 = kgC/m2
106 factor_aboveground = 0.5
108 factor_unit_mod = 0.001
110 dataob = dataob * factor_aboveground * factor_unit_ob
111 datamod0 = datamod0 * factor_unit_mod
113 dataob@units = "kg C m~S~-2~N~"
114 datamod0@units = "kg C m~S~-2~N~"
116 dataob@long_name = "Amazon Biomass"
117 datamod0@long_name = "Amazon Biomass"
118 ;********************************************************
119 ; get subset of datamod0 that match dataob
124 ind_lonL = ind(xm .eq. xo(0))
125 ind_lonR = ind(xm .eq. xo(nlon-1))
126 ind_latS = ind(ym .eq. yo(0))
127 ind_latN = ind(ym .eq. yo(nlat-1))
130 datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)
133 area(:,:) = area0(ind_latS:ind_latN,ind_lonL:ind_lonR)
136 mask_amazon(:,:) = mask_amazon0(ind_latS:ind_latN,ind_lonL:ind_lonR)
138 mask_amazon@units = ""
139 ;********************************************************
140 ; sum over amazom_mask area:
142 ; Peta g = 1.e15 g = 1.e12 kg
145 ; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.)
147 Sum_area = sum(area*mask_amazon)*factor_unit
149 Sum_ob = sum(dataob*area*mask_amazon)*factor_unit
150 Sum_mod = sum(datamod*area*mask_amazon)*factor_unit
152 avg_ob = Sum_ob /Sum_area
153 avg_mod = Sum_mod/Sum_area
155 Sum_biomass_ob = sprintf("%.2f",Sum_ob )
156 Sum_biomass_mod = sprintf("%.2f",Sum_mod)
158 ;----------------------------------------------------------------------
161 resg = True ; Use plot options
162 resg@cnFillOn = True ; Turn on color fill
163 resg@gsnSpreadColors = True ; use full colormap
164 resg@cnLinesOn = False ; Turn off contourn lines
165 resg@mpFillOn = False ; Turn off map fill
166 resg@gsnAddCyclic = False
168 resg@gsnSpreadColors = True ; use full colormap
169 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
171 resg@mpMinLatF = -21.1 ; range to zoom in on
172 resg@mpMaxLatF = 13.8
173 resg@mpMinLonF = 277.28
174 resg@mpMaxLonF = 326.43
175 ;------------------------------------------------------------------------
178 plot_name = "mask_ob"
180 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
181 gsn_define_colormap(wks,"gui_default") ; choose colormap
183 ;-----------------------------------------
187 gRes@txFontHeightF = 0.02
190 area_sum_text = "(mask area = "+sprintf("%.2f", Sum_area)+"e12 m2)"
192 gsn_text_ndc(wks,area_sum_text,0.50,0.80,gRes)
193 ;-----------------------------------------
195 resg@cnMinLevelValF = 0. ; Min level
196 resg@cnMaxLevelValF = 1. ; Max level
197 resg@cnLevelSpacingF = 0.1 ; interval
199 resg@tiMainString = "Amazon Mask: grid = "+ model_grid
201 plot = gsn_csm_contour_map_ce(wks,mask_amazon,resg)
205 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
206 "rm "+plot_name+"."+plot_type)
207 ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
209 ;------------------------------------------------------------------------
212 resg@cnMinLevelValF = 0. ; Min level
213 resg@cnMaxLevelValF = 30. ; Max level
214 resg@cnLevelSpacingF = 2. ; interval
216 plot_name = "global_ob"
217 title = ob_name+" "+ model_grid
218 resg@tiMainString = title
220 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
221 gsn_define_colormap(wks,"gui_default") ; choose colormap
223 plot = gsn_csm_contour_map_ce(wks,dataob,resg)
227 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
228 "rm "+plot_name+"."+plot_type)
229 ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
231 ;------------------------------------------------------------------------
234 resg@cnMinLevelValF = 0. ; Min level
235 resg@cnMaxLevelValF = 30. ; Max level
236 resg@cnLevelSpacingF = 2. ; interval
238 plot_name = "global_model"
239 title = "Model "+ model_name + " (2000)"
240 resg@tiMainString = title
242 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
243 gsn_define_colormap(wks,"gui_default") ; choose colormap
245 plot = gsn_csm_contour_map_ce(wks,datamod,resg)
249 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
250 "rm "+plot_name+"."+plot_type)
251 ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
253 ;------------------------------------------------------------------------
254 ; contour model vs ob
256 plot_name = "global_model_vs_ob"
258 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
259 gsn_define_colormap(wks,"gui_default") ; choose colormap
262 plot=new(3,graphic) ; create graphic array
264 resg@gsnFrame = False ; Do not draw plot
265 resg@gsnDraw = False ; Do not advance frame
267 ;(d) compute correlation coef and M score
269 uu = ndtooned(datamod)
270 vv = ndtooned(dataob)
272 good = ind(.not.ismissing(uu) .and. .not.ismissing(vv))
277 ccrG = esccr(ug,vg,0)
281 ; Miomass = (ccrG*ccrG)* score_max
283 bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
284 print ("bias=" + bias + "dimsizes(ug)=" + dimsizes(ug))
285 print ("other bias=" + sum(abs(ug-vg)/(ug+vg)))
286 Mbiomass = (1. - (bias/dimsizes(ug)))*score_max
287 print ("Mbiomass=" + Mbiomass)
288 M_biomass = sprintf("%.2f", Mbiomass)
289 print ("M_biomass=" + M_biomass)
291 if (isvar("compare")) then
292 system("sed -e '1,/M_biomass/s/M_biomass/"+M_biomass+"/' "+html_name2+" > "+html_new2+";"+ \
293 "mv -f "+html_new2+" "+html_name2)
296 system("sed s#M_biomass#"+M_biomass+"# "+html_name+" > "+html_new+";"+ \
297 "mv -f "+html_new+" "+html_name)
299 ; plot correlation coef
302 gRes@txFontHeightF = 0.02
305 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
307 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
308 ;--------------------------------------------------------------------
312 title = ob_name+" "+ model_grid
313 resg@tiMainString = title
315 plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg)
319 title = "Model "+ model_name + " (2000)"
320 resg@tiMainString = title
322 plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg)
327 zz = datamod - dataob
328 title = "Model "+model_name+" (2000) - Observed"
330 resg@cnMinLevelValF = -10. ; Min level
331 resg@cnMaxLevelValF = 10. ; Max level
332 resg@cnLevelSpacingF = 2. ; interval
333 resg@tiMainString = title
335 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
337 pres = True ; panel plot mods desired
338 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
339 ; indiv. plots in panel
340 pres@gsnMaximize = True ; fill the page
342 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
348 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
349 "rm "+plot_name+"."+plot_type)
350 ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
352 resg@gsnFrame = True ; draw plot
353 resg@gsnDraw = True ; advance frame
354 ;------------------------------------------------------------------------
355 ; contour ob : masked
357 resg@cnMinLevelValF = 0. ; Min level
358 resg@cnMaxLevelValF = 30. ; Max level
359 resg@cnLevelSpacingF = 2. ; interval
361 plot_name = "global_mask_ob"
362 title = ob_name+" "+ model_grid
363 resg@tiMainString = title
365 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
366 gsn_define_colormap(wks,"gui_default") ; choose colormap
367 ;-----------------------------------------
368 ; plot average over mask region
371 gRes@txFontHeightF = 0.02
374 area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_ob)+" kg C m~S~-2~N~)"
376 gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)
377 ;-----------------------------------------
379 zo = dataob*mask_amazon
380 zo = where((mask_amazon .le. 0.01), zo@_FillValue, zo)
381 plot = gsn_csm_contour_map_ce(wks,zo,resg)
386 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
387 "rm "+plot_name+"."+plot_type)
388 ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
390 ;------------------------------------------------------------------------
391 ; contour model: masked
393 resg@cnMinLevelValF = 0. ; Min level
394 resg@cnMaxLevelValF = 30. ; Max level
395 resg@cnLevelSpacingF = 2. ; interval
397 plot_name = "global_mask_model"
398 title = "Model "+ model_name + " (2000)"
399 resg@tiMainString = title
401 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
402 gsn_define_colormap(wks,"gui_default") ; choose colormap
403 ;-----------------------------------------
404 ; plot average over mask region
407 gRes@txFontHeightF = 0.02
410 area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_mod)+" kg C m~S~-2~N~)"
412 gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)
413 ;-----------------------------------------
415 zm = datamod*mask_amazon
416 zm = where((mask_amazon .le. 0.01), zm@_FillValue, zm)
417 plot = gsn_csm_contour_map_ce(wks,zm,resg)
422 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
423 "rm "+plot_name+"."+plot_type)
424 ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
426 ;------------------------------------------------------------------------
427 ; contour model vs ob: masked
429 plot_name = "global_mask_vs_ob"
431 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
432 gsn_define_colormap(wks,"gui_default") ; choose colormap
434 plot=new(3,graphic) ; create graphic array
436 resg@gsnFrame = False ; Do not draw plot
437 resg@gsnDraw = False ; Do not advance frame
439 ;(d) compute correlation coef and M score
452 good = ind((uu .gt. 0.) .and. (vv .gt. 0.))
457 ccrG = esccr(ug,vg,0)
459 ; Miomass = (ccrG*ccrG)*score_max
461 bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
462 Mbiomass2 = (1. - (bias/dimsizes(ug)))*score_max
463 M_biomask = sprintf("%.2f", Mbiomass2)
465 if (isvar("compare")) then
466 system("sed -e '1,/M_biomask/s/M_biomask/"+M_biomask+"/' "+html_name2+" > "+html_new2+";"+ \
467 "mv -f "+html_new2+" "+html_name2)
468 system("sed -e '1,/Sum_biomass_ob/s/Sum_biomass_ob/"+Sum_biomass_ob+"/' "+html_name2+" > "+html_new2+";"+ \
469 "mv -f "+html_new2+" "+html_name2)
470 system("sed -e '1,/Sum_biomass_mod/s/Sum_biomass_mod/"+Sum_biomass_mod+"/' "+html_name2+" > "+html_new2+";"+ \
471 "mv -f "+html_new2+" "+html_name2)
474 system("sed s#M_biomask#"+M_biomask+"# "+html_name+" > "+html_new+";"+ \
475 "mv -f "+html_new+" "+html_name)
476 system("sed s#Sum_biomass_ob#"+Sum_biomass_ob+"# "+html_name+" > "+html_new+";"+ \
477 "mv -f "+html_new+" "+html_name)
478 system("sed s#Sum_biomass_mod#"+Sum_biomass_mod+"# "+html_name+" > "+html_new+";"+ \
479 "mv -f "+html_new+" "+html_name)
480 ;--------------------------------------------------------------------
481 ; plot correlation coef
484 gRes@txFontHeightF = 0.02
487 correlation_text = "(correlation coef = "+sprintf("%.2f", ccrG)+")"
489 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
490 ;--------------------------------------------------------------------
494 title = ob_name+" "+ model_grid
495 resg@tiMainString = title
497 plot(0) = gsn_csm_contour_map_ce(wks,zo,resg)
501 title = "Model "+ model_name + " (2000)"
502 resg@tiMainString = title
504 plot(1) = gsn_csm_contour_map_ce(wks,zm,resg)
510 title = "Model "+model_name+" (2000) - Observed"
512 resg@cnMinLevelValF = -10. ; Min level
513 resg@cnMaxLevelValF = 10. ; Max level
514 resg@cnLevelSpacingF = 2. ; interval
515 resg@tiMainString = title
517 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
519 pres = True ; panel plot mods desired
520 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
521 ; indiv. plots in panel
522 pres@gsnMaximize = True ; fill the page
524 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
529 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
530 "rm "+plot_name+"."+plot_type)
531 ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
533 ;***************************************************************************
534 ; add total score and write to file
535 ;***************************************************************************
536 M_total = Mbiomass + Mbiomass2
538 asciiwrite("M_save.biomass", M_total)
540 ;***************************************************************************
542 ;***************************************************************************
543 output_dir = model_name+"/biomass"
545 system("mv *.png " + output_dir)
546 ;***************************************************************************