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 ;**************************************************************
6 procedure set_line(lines:string,nline:integer,newlines:string)
8 ; add line to ascci/html file
10 nnewlines = dimsizes(newlines)
11 if(nline+nnewlines-1.ge.dimsizes(lines))
12 print("set_line: bad index, not setting anything.")
15 lines(nline:nline+nnewlines-1) = newlines
16 ; print ("lines = " + lines(nline:nline+nnewlines-1))
17 nline = nline + nnewlines
20 ;**************************************************************
27 ;-----------------------------------------------------
28 ; edit table.html of current model for movel1_vs_model2
30 if (isvar("compare")) then
31 html_name2 = compare+"/table.html"
32 html_new2 = html_name2 +".new"
35 ;------------------------------------------------------
36 ; edit table.html for current model
38 html_name = model_name+"/table.html"
39 html_new = html_name +".new"
41 ;------------------------------------------------------
44 fm = addfile(dirm+film10,"r")
49 dsizes = dimsizes(laimod)
54 ;-----------------------------------
57 film_l = "lnd_"+ model_grid + ".nc"
58 fm_l = addfile (dirs+film_l,"r")
59 landfrac = fm_l->landfrac
62 ;----------------------------------
63 ; read biome data: model
65 biome_name_mod = "Model PFT Class"
67 film_c = "class_pft_"+model_grid+".nc"
68 fm_c = addfile(dirs+film_c,"r")
69 classmod = fm_c->CLASS_PFT
73 ; model data has 17 land-type classes
76 ;----------------------------------------------------------
79 ;----------------------------------
80 ; read biome data: observed
82 biome_name_ob = "MODIS LandCover"
85 filo_c = "land_class_"+model_grid+".nc"
86 fo = addfile(dir_c+filo_c,"r")
87 classob = tofloat(fo->LAND_CLASS)
91 ; input observed data has 20 land-type classes
94 ;---------------------------------
95 ; read lai data: observed
97 ;ob_name = "MODIS MOD15A2 (2000-2005)"
98 ob_name = "MODIS MOD15A2 (2000-2004)"
100 dir_l = diro + "lai/"
101 ;filo = "LAI_2000-2005_MONS_"+model_grid+".nc"
102 filo = "LAI_2000-2004_MONS_"+model_grid+".nc"
103 fo = addfile(dir_l+filo,"r")
108 ;-------------------------------------------------
109 ; take into account landfrac
111 laimod = laimod * conform(laimod,landfrac,(/1,2/))
112 laiob = laiob * conform(laiob ,landfrac,(/1,2/))
116 ;************************************************
118 ;************************************************
119 resg = True ; Use plot options
120 resg@cnFillOn = True ; Turn on color fill
121 resg@gsnSpreadColors = True ; use full colormap
122 resg@cnLinesOn = False ; Turn off contourn lines
123 resg@mpFillOn = False ; Turn off map fill
124 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
126 ;************************************************
127 ; plot global biome class: (1) observed
128 ;************************************************
130 resg@cnMinLevelValF = 1. ; Min level
131 resg@cnMaxLevelValF = 19. ; Max level
132 resg@cnLevelSpacingF = 1. ; interval
134 classob@_FillValue = 1.e+36
135 classob = where(classob.eq.0,classob@_FillValue,classob)
137 plot_name = "global_class_ob"
138 title = biome_name_ob
139 resg@tiMainString = title
141 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
142 gsn_define_colormap(wks,"gui_default") ; choose colormap
144 plot = gsn_csm_contour_map_ce(wks,classob,resg)
148 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
149 "rm "+plot_name+"."+plot_type)
151 ;************************************************
152 ; plot global biome class: (2) model
153 ;************************************************
155 resg@cnMinLevelValF = 0. ; Min level
156 resg@cnMaxLevelValF = 16. ; Max level
157 resg@cnLevelSpacingF = 1. ; interval
159 plot_name = "global_class_model"
160 title = biome_name_mod
161 resg@tiMainString = title
163 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
164 gsn_define_colormap(wks,"gui_default") ; choose colormap
166 plot = gsn_csm_contour_map_ce(wks,classmod,resg)
170 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
171 "rm "+plot_name+"."+plot_type)
173 ;*******************************************************************
174 ; for html table : all 3 components (Mean, Max, Phase)
175 ;*******************************************************************
177 ; column (not including header column)
179 component = (/"Mean","Max","Phase"/)
181 col_head = (/"observed",model_name,"M_score" \
182 ,"observed",model_name,"M_score" \
183 ,"observed",model_name,"M_score" \
186 n_comp = dimsizes(component)
187 ncol = dimsizes(col_head )
189 ; row (not including header row)
191 ;----------------------------------------------------
192 ; using model biome class:
193 row_head = (/"Not Vegetated" \
194 ,"Needleleaf Evergreen Temperate Tree" \
195 ,"Needleleaf Evergreen Boreal Tree" \
196 ; ,"Needleleaf Deciduous Boreal Tree" \
197 ,"Broadleaf Evergreen Tropical Tree" \
198 ,"Broadleaf Evergreen Temperate Tree" \
199 ,"Broadleaf Deciduous Tropical Tree" \
200 ,"Broadleaf Deciduous Temperate Tree" \
201 ; ,"Broadleaf Deciduous Boreal Tree" \
202 ; ,"Broadleaf Evergreen Shrub" \
203 ,"Broadleaf Deciduous Temperate Shrub" \
204 ,"Broadleaf Deciduous Boreal Shrub" \
206 ,"C3 Non-Arctic Grass" \
213 nrow = dimsizes(row_head)
215 ; arrays to be passed to table.
216 text = new ((/nrow, ncol/),string )
221 ;********************************************************************
222 ; use land-type class to bin the data in equally spaced ranges
223 ;********************************************************************
225 ; using model biome class
228 range = fspan(0,nclass,nclass+1)
230 ; Use this range information to grab all the values in a
231 ; particular range, and then take an average.
233 nx = dimsizes(range) - 1
235 ; for 2 data: model and observed
238 ; using model biome class
240 base = ndtooned(classmod)
244 yvalues = new((/data_n,nx/),float)
245 count = new((/data_n,nx/),float)
247 ;************************************************************************
248 ; go through all components
249 ;************************************************************************
259 data_ob = dim_avg_Wrap(laiob (lat|:,lon|:,time|:))
260 data_mod = dim_avg_Wrap(laimod(lat|:,lon|:,time|:))
269 data_ob = laiob(0,:,:)
270 data_ob@long_name = "Leaf Area Index Max"
274 data_ob(j,i) = max(laiob(:,j,i))
279 data_mod = laimod(0,:,:)
280 data_mod@long_name = "Leaf Area Index Max"
284 data_mod(j,i) = max(laimod(:,j,i))
297 data_ob = laiob(0,:,:)
298 data_ob@long_name = "Leaf Area Index Max Month"
302 data_ob(j,i) = maxind(laiob(:,j,i)) + 1
307 data_mod = laimod(0,:,:)
308 data_mod@long_name = "Leaf Area Index Max Month"
312 data_mod(j,i) = maxind(laimod(:,j,i)) + 1
320 ;==============================
322 ;==============================
324 ; Loop through each range, using base
328 if (i.ne.(nx-1)) then
329 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
331 idx = ind(base.ge.range(i))
334 ; loop through each dataset
339 data = ndtooned(data_ob)
343 data = ndtooned(data_mod)
348 if (.not.any(ismissing(idx))) then
349 yvalues(n,i) = avg(data(idx))
350 count(n,i) = dimsizes(idx)
352 yvalues(n,i) = yvalues@_FillValue
356 ;#############################################################
357 ; using model biome class:
359 ; set the following 4 classes to _FillValue:
360 ; (3)Needleleaf Deciduous Boreal Tree,
361 ; (8)Broadleaf Deciduous Boreal Tree,
362 ; (9)Broadleaf Evergreen Shrub,
365 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
366 yvalues(n,i) = yvalues@_FillValue
369 ;#############################################################
377 ;=====================================
378 ; compute correlation coef and M score
379 ;=====================================
381 ; FMH: Max scores are depend on which component is being analyzed, so this
382 ; is now set in each if statement above.
390 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
394 uu_count = u_count(good)
395 vv_count = v_count(good)
397 ; compute correlation coef
401 bias = avg(abs(vv-uu))
402 bias = where((bias.gt. 6.),12.-bias,bias)
403 Mscore = ((6. - bias)/6.)*score_max
405 bias = sum(abs(vv-uu)/abs(vv+uu))
406 Mscore = (1.- (bias/dimsizes(uu)))*score_max
409 M_score = sprintf("%.2f", Mscore)
413 M_total = M_total + Mscore
415 ;================================
416 ; output M_score to score sheet
417 ;===============================
419 if (isvar("compare")) then
420 system("sed -e '1,/M_lai_"+component(m)+"/s/M_lai_"+component(m)+"/"+M_score+"/' "+html_name2+" > "+html_new2+";"+ \
421 "mv -f "+html_new2+" "+html_name2)
424 system("sed s#M_lai_"+component(m)+"#"+M_score+"# "+html_name+" > "+html_new+";"+ \
425 "mv -f "+html_new+" "+html_name)
427 ;==============================
428 ; output M_score to html table
429 ;==============================
434 text(i,n) = sprintf("%.1f",uu(i))
435 text(i,n+1) = sprintf("%.1f",vv(i))
438 text(nrow-1,n) = sprintf("%.1f",sum(uu*uu_count)/sum(uu_count))
439 text(nrow-1,n+1) = sprintf("%.1f",sum(vv*vv_count)/sum(vv_count))
440 text(nrow-1,n+2) = M_score
452 ;========================================
453 ; global res changes for each component
454 ;========================================
458 resg@cnMinLevelValF = 0.
459 resg@cnMaxLevelValF = 10.
460 resg@cnLevelSpacingF = 1.
462 data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)
466 resg@cnMinLevelValF = 0.
467 resg@cnMaxLevelValF = 10.
468 resg@cnLevelSpacingF = 1.
470 data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)
474 resg@cnMinLevelValF = 1.
475 resg@cnMaxLevelValF = 12.
476 resg@cnLevelSpacingF = 1.
478 data_mod = where(ismissing(data_ob).or.(data_mod.lt.delta),data_mod@_FillValue,data_mod)
479 data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)
483 ;=========================
484 ; global contour : ob
485 ;=========================
487 plot_name = "global_"+component(m)+"_ob"
489 resg@tiMainString = title
491 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
492 gsn_define_colormap(wks,"gui_default") ; choose colormap
494 plot = gsn_csm_contour_map_ce(wks,data_ob,resg)
499 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
500 "rm "+plot_name+"."+plot_type)
501 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
504 ;============================
505 ; global contour : model
506 ;============================
508 plot_name = "global_"+component(m)+"_model"
509 title = "Model " + model_name + " (2000-2004)"
510 resg@tiMainString = title
512 wks = gsn_open_wks (plot_type,plot_name)
513 gsn_define_colormap(wks,"gui_default")
515 plot = gsn_csm_contour_map_ce(wks,data_mod,resg)
520 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
521 "rm "+plot_name+"."+plot_type)
522 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
524 ;================================
525 ; global contour: model vs ob
526 ;================================
528 plot_name = "global_"+component(m)+"_model_vs_ob"
530 wks = gsn_open_wks (plot_type,plot_name)
531 gsn_define_colormap(wks,"gui_default")
533 plot=new(3,graphic) ; create graphic array
535 resg@gsnFrame = False ; Do not draw plot
536 resg@gsnDraw = False ; Do not advance frame
538 ; plot correlation coef
541 gRes@txFontHeightF = 0.02
544 correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")"
546 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
551 resg@tiMainString = title
553 plot(0) = gsn_csm_contour_map_ce(wks,data_ob,resg)
557 title = "Model "+ model_name + " (2000-2004)"
558 resg@tiMainString = title
560 plot(1) = gsn_csm_contour_map_ce(wks,data_mod,resg)
565 resg@cnMinLevelValF = -2.
566 resg@cnMaxLevelValF = 2.
567 resg@cnLevelSpacingF = 0.4
571 resg@cnMinLevelValF = -6.
572 resg@cnMaxLevelValF = 6.
573 resg@cnLevelSpacingF = 1.
577 resg@cnMinLevelValF = -6.
578 resg@cnMaxLevelValF = 6.
579 resg@cnLevelSpacingF = 1.
583 zz = data_mod - data_ob
584 title = "Model_"+model_name+" - MODIS MOD15A2"
585 resg@tiMainString = title
587 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
591 pres = True ; panel plot mods desired
592 pres@gsnMaximize = True ; fill the page
594 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
599 system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
600 "rm "+plot_name+"."+plot_type)
601 ;system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
606 resg@gsnFrame = True ; Do advance frame
607 resg@gsnDraw = True ; Do draw plot
611 ;**************************************************
613 ;**************************************************
614 output_html = "table_model_vs_ob.html"
616 header_text = "<H1>LAI: Model "+model_name+" vs. MODIS observations</H1>"
618 header = (/"<HTML>" \
620 ,"<TITLE>CLAMP metrics</TITLE>" \
627 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
629 ," <th bgcolor=DDDDDD rowspan=2>Biome Class</th>" \
630 ," <th bgcolor=DDDDDD colspan=3>"+component(0)+" (m2/m2)</th>" \
631 ," <th bgcolor=DDDDDD colspan=3>"+component(1)+" (m2/m2)</th>" \
632 ," <th bgcolor=DDDDDD colspan=3>"+component(2)+" (months)</th>" \
633 ," <th bgcolor=DDDDDD rowspan=2>Annual Plot</th>" \
636 ," <th bgcolor=DDDDDD >observed</th>" \
637 ," <th bgcolor=DDDDDD >"+model_name+"</th>" \
638 ," <th bgcolor=DDDDDD >M_score</th>" \
639 ," <th bgcolor=DDDDDD >observed</th>" \
640 ," <th bgcolor=DDDDDD >"+model_name+"</th>" \
641 ," <th bgcolor=DDDDDD >M_score</th>" \
642 ," <th bgcolor=DDDDDD >observed</th>" \
643 ," <th bgcolor=DDDDDD >"+model_name+"</th>" \
644 ," <th bgcolor=DDDDDD >M_score</th>" \
647 table_footer = "</table>"
651 lines = new(50000,string)
654 set_line(lines,nline,header)
655 set_line(lines,nline,table_header)
656 ;-----------------------------------------------
660 set_line(lines,nline,row_header)
672 txt11 = "<a href=./annual_biome_"+n+".png>model_vs_ob</a>"
674 set_line(lines,nline,"<th>"+txt1+"</th>")
675 set_line(lines,nline,"<th>"+txt2+"</th>")
676 set_line(lines,nline,"<th>"+txt3+"</th>")
677 set_line(lines,nline,"<th>"+txt4+"</th>")
678 set_line(lines,nline,"<th>"+txt5+"</th>")
679 set_line(lines,nline,"<th>"+txt6+"</th>")
680 set_line(lines,nline,"<th>"+txt7+"</th>")
681 set_line(lines,nline,"<th>"+txt8+"</th>")
682 set_line(lines,nline,"<th>"+txt9+"</th>")
683 set_line(lines,nline,"<th>"+txt10+"</th>")
684 set_line(lines,nline,"<th>"+txt11+"</th>")
686 set_line(lines,nline,row_footer)
688 ;-----------------------------------------------
689 set_line(lines,nline,table_footer)
690 set_line(lines,nline,footer)
692 ; Now write to an HTML file
694 idx = ind(.not.ismissing(lines))
695 if(.not.any(ismissing(idx))) then
696 asciiwrite(output_html,lines(idx))
705 ;************************************************************************
706 ; go through all ntime
707 ;************************************************************************
709 ; for 2 data: model and observed
712 ; using model biome class
714 base = ndtooned(classmod)
718 yvalues = new((/data_n,nx,ntime/),float)
720 ;==============================
722 ;==============================
724 ; Loop through each range, using base
728 if (i.ne.(nx-1)) then
729 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
731 idx = ind(base.ge.range(i))
734 ; loop through each dataset
741 data = ndtooned(laiob (m,:,:))
745 data = ndtooned(laimod(m,:,:))
750 if (.not.any(ismissing(idx))) then
751 yvalues(n,i,m) = avg(data(idx))
753 yvalues(n,i,m) = yvalues@_FillValue
756 ;#############################################################
757 ; using model biome class:
759 ; set the following 4 classes to _FillValue:
760 ; (3)Needleleaf Deciduous Boreal Tree,
761 ; (8)Broadleaf Deciduous Boreal Tree,
762 ; (9)Broadleaf Evergreen Shrub,
765 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
766 yvalues(n,i,m) = yvalues@_FillValue
768 ;#############################################################
778 good = ind(.not.ismissing(yvalues(0,:,0)))
780 n_biome = dimsizes(good)
782 ;----------------------------------------------------------------
783 ; data for tseries plot
785 yvalues_g = new((/data_n,n_biome,ntime/),float)
789 yvalues_g = yvalues(:,good,:)
793 ;*******************************************************************
795 ;*******************************************************************
796 ; for x-axis in xyplot
798 mon@long_name = "month"
800 res = True ; plot mods desired
801 res@xyLineThicknesses = (/2.0,2.0,2.0/) ; make 2nd lines thicker
802 res@xyLineColors = (/"blue","red"/) ; change line color
804 ; res@tiMainFontHeightF = 0.025 ; size of title
806 ; Add a boxed legend using the more simple method
807 res@pmLegendDisplayMode = "Always"
808 ; res@pmLegendWidthF = 0.1
809 res@pmLegendWidthF = 0.08
810 res@pmLegendHeightF = 0.06
811 ; res@pmLegendOrthogonalPosF = -1.17
812 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
813 res@pmLegendOrthogonalPosF = -0.30 ;(downward)
815 ; res@pmLegendParallelPosF = 0.18
816 res@pmLegendParallelPosF = 0.23 ;(rightward)
818 ; res@lgPerimOn = False
819 res@lgLabelFontHeightF = 0.015
820 res@xyExplicitLegendLabels = (/"observed",model_name/)
822 ;************************************************************
824 plot_data = new((/2,ntime/),float)
825 plot_data@long_name = "Leaf Area Index"
827 ;----------------------------------------------
828 ; time series plot : per biome
832 plot_name = "annual_biome_"+ m
834 wks = gsn_open_wks (plot_type,plot_name)
836 title = "LAI : "+ row_head(m)
837 res@tiMainString = title
839 plot_data(0,:) = yvalues_g(0,m,:)
840 plot_data(1,:) = yvalues_g(1,m,:)
842 plot = gsn_csm_xy(wks,mon,plot_data,res)
847 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
848 "rm "+plot_name+"."+plot_type)
851 ;--------------------------------------------
852 ; time series plot: all biome
854 plot_name = "annual_biome_"+ n_biome
856 wks = gsn_open_wks (plot_type,plot_name)
858 title = "LAI : "+ row_head(n_biome)
859 res@tiMainString = title
862 plot_data(0,k) = avg(yvalues_g(0,:,k))
863 plot_data(1,k) = avg(yvalues_g(1,:,k))
866 plot = gsn_csm_xy(wks,mon,plot_data,res)
871 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
872 "rm "+plot_name+"."+plot_type)
876 ;***************************************************************************
877 ; write total score to file
878 ;***************************************************************************
880 asciiwrite("M_save.lai", M_total)
882 ;***************************************************************************
884 ;***************************************************************************
885 output_dir = model_name+"/lai"
887 system("mv *.png *.html " + output_dir)
888 ;***************************************************************************