1 ; ****************************************************
2 ; combine scatter, histogram, global and zonal plots
3 ; compute all correlation coef and M score
4 ; ****************************************************
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
8 ;************************************************
13 ;************************************************
15 ;************************************************
16 ;film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
17 ;model_name = "b30.061n"
20 ;film = "newcn05_ncep_1i_ANN_climo_lnd.nc"
23 ;--------------------------------------------
24 film = "i01.10cn_1948-2004_ANN_climo.nc"
28 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
29 fm = addfile (dirm+film,"r")
35 datamod0 = data1(0,:,:)
36 datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)
37 ;-------------------------------------------
38 ;film = "i01.10casa_1948-2004_ANN_climo.nc"
39 ;model_name = "10casa"
42 ;dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
43 ;fm = addfile (dirm+film,"r")
49 ;datamod0 = data1(0,:,:)
50 ;datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:)
51 ;----------------------------------------------------
52 html_name = "table.html." + model_name
53 html_new = html_name +".new"
54 system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \
55 "mv -f "+html_new+" "+html_name)
59 ;************************************************
61 ;************************************************
62 ob_name = "LC15_Amazon_Biomass"
64 diro = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"
65 filo = "amazon_biomass_"+model_grid+".nc"
66 fo = addfile (diro+filo,"r")
71 ;************************************************
72 ; Units for these variables are:
75 ; We want to convert these to KgC/m2
76 ; ha = 100m*100m = 10,000 m2
77 ; MgC/ha*1000/10,000 = KgC/m2
79 factor_aboveground = 0.5
81 factor_unit_mod = 0.001
83 dataob = dataob * factor_aboveground * factor_unit_ob
84 datamod0 = datamod0 * factor_unit_mod
86 dataob@units = "KgC/m2"
87 datamod0@units = "KgC/m2"
89 dataob@long_name = "Amazon Biomass"
90 datamod0@long_name = "Amazon Biomass"
91 ;********************************************************
92 ; get subset of datamod0 that match dataob
97 ind_lonL = ind(xm .eq. xo(0))
98 ind_lonR = ind(xm .eq. xo(nlon-1))
99 ind_latS = ind(ym .eq. yo(0))
100 ind_latN = ind(ym .eq. yo(nlat-1))
103 datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)
107 ;----------------------------------------------------------------------
111 resg = True ; Use plot options
112 resg@cnFillOn = True ; Turn on color fill
113 resg@gsnSpreadColors = True ; use full colormap
114 ; resg@cnFillMode = "RasterFill" ; Turn on raster color
115 ; resg@lbLabelAutoStride = True
116 resg@cnLinesOn = False ; Turn off contourn lines
117 resg@mpFillOn = False ; Turn off map fill
118 resg@gsnAddCyclic = False
120 resg@gsnSpreadColors = True ; use full colormap
121 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
122 resg@cnMinLevelValF = 0. ; Min level
123 resg@cnMaxLevelValF = 30. ; Max level
124 resg@cnLevelSpacingF = 2. ; interval
126 resg@mpMinLatF = -21.1 ; range to zoom in on
127 resg@mpMaxLatF = 13.8
128 resg@mpMinLonF = 277.28
129 resg@mpMaxLonF = 326.43
130 ;------------------------------------------------------------------------
133 plot_name = "global_ob"
134 title = ob_name+" "+ model_grid
135 resg@tiMainString = title
137 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
138 gsn_define_colormap(wks,"gui_default") ; choose colormap
140 plot = gsn_csm_contour_map_ce(wks,dataob,resg)
143 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
144 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
145 "rm "+plot_name+"-*."+plot_type_new+";"+ \
146 "rm "+plot_name+"."+plot_type)
149 ;------------------------------------------------------------------------
150 ;global contour model
152 plot_name = "global_model"
153 title = "Model "+ model_name
154 resg@tiMainString = title
156 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
157 gsn_define_colormap(wks,"gui_default") ; choose colormap
159 plot = gsn_csm_contour_map_ce(wks,datamod,resg)
162 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
163 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
164 "rm "+plot_name+"-*."+plot_type_new+";"+ \
165 "rm "+plot_name+"."+plot_type)
168 ;------------------------------------------------------------------------
169 ;global contour model vs ob
171 plot_name = "global_model_vs_ob"
173 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
174 gsn_define_colormap(wks,"gui_default") ; choose colormap
177 plot=new(3,graphic) ; create graphic array
179 resg@gsnFrame = False ; Do not draw plot
180 resg@gsnDraw = False ; Do not advance frame
182 ;(d) compute correlation coef and M score
184 uu = ndtooned(datamod)
185 vv = ndtooned(dataob)
187 good = ind(.not.ismissing(uu) .and. .not.ismissing(vv))
192 ccrG = esccr(ug,vg,0)
195 MG = (ccrG*ccrG)* 5.0
197 bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
198 MG = (1. - (bias/dimsizes(ug)))*5.
200 M_biomass = sprintf("%.2f", MG)
201 system("sed s#M_biomass#"+M_biomass+"# "+html_name+" > "+html_new+";"+ \
202 "mv -f "+html_new+" "+html_name)
205 ; plot correlation coef
208 gRes@txFontHeightF = 0.02
211 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
213 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
214 ;--------------------------------------------------------------------
218 title = ob_name+" "+ model_grid
219 resg@tiMainString = title
221 plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg)
225 title = "Model "+ model_name
226 resg@tiMainString = title
228 plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg)
233 zz = datamod - dataob
234 title = "Model_"+model_name+" - Observed"
236 resg@cnMinLevelValF = -10. ; Min level
237 resg@cnMaxLevelValF = 10. ; Max level
238 resg@cnLevelSpacingF = 2. ; interval
239 resg@tiMainString = title
241 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
243 pres = True ; panel plot mods desired
244 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
245 ; indiv. plots in panel
246 pres@gsnMaximize = True ; fill the page
248 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
250 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
251 "rm "+plot_name+"."+plot_type)
255 ;-------------------------------------------------------------------
256 temp_name = "biomass." + model_name
257 system("mkdir -p " + temp_name+";"+ \
258 "cp "+ html_name + " " +temp_name+"/table.html"+";"+ \
259 "mv *.png " + temp_name +";"+ \
260 "tar cf "+ temp_name +".tar " + temp_name)
261 ;-------------------------------------------------------------------