Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
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 = "i01.03cn_1545-1569_ANN_climo.nc"
17 ;model_name = "i01.03cn"
20 ;film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
21 ;model_name = "b30.061n"
24 ;film = "newcn05_ncep_1i_ANN_climo_lnd.nc"
28 ;film = "i01.06cn_1798-2004_ANN_climo.nc"
32 ;film = "i01.06casa_1798-2004_ANN_climo.nc"
33 ;model_name = "06casa"
37 film = "i01.10cn_1948-2004_ANN_climo.nc"
41 ;film = "i01.10casa_1948-2004_ANN_climo.nc"
42 ;model_name = "10casa"
44 ;----------------------------------------------------
45 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
46 fm = addfile (dirm+film,"r")
52 datamod0 = data1(0,:,:)
53 datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)
57 ;************************************************
59 ;************************************************
60 ob_name = "LC15_Amazon_Biomass"
62 diro = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"
63 filo = "amazon_biomass_"+model_grid+".nc"
64 fo = addfile (diro+filo,"r")
69 ;************************************************
70 ; Units for these variables are:
73 ; We want to convert these to KgC/m2
74 ; ha = 100m*100m = 10,000 m2
75 ; MgC/ha*1000/10,000 = KgC/m2
77 factor_aboveground = 0.5
79 factor_unit_mod = 0.001
81 dataob = dataob * factor_aboveground * factor_unit_ob
82 datamod0 = datamod0 * factor_unit_mod
84 dataob@units = "KgC/m2"
85 datamod0@units = "KgC/m2"
87 dataob@long_name = "Amazon Biomass"
88 datamod0@long_name = "Amazon Biomass"
89 ;********************************************************
90 ; get subset of datamod0 that match dataob
95 ind_lonL = ind(xm .eq. xo(0))
96 ind_lonR = ind(xm .eq. xo(nlon-1))
97 ind_latS = ind(ym .eq. yo(0))
98 ind_latN = ind(ym .eq. yo(nlat-1))
101 datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)
105 ;----------------------------------------------------------------------
109 resg = True ; Use plot options
110 resg@cnFillOn = True ; Turn on color fill
111 resg@gsnSpreadColors = True ; use full colormap
112 ; resg@cnFillMode = "RasterFill" ; Turn on raster color
113 ; resg@lbLabelAutoStride = True
114 resg@cnLinesOn = False ; Turn off contourn lines
115 resg@mpFillOn = False ; Turn off map fill
116 resg@gsnAddCyclic = False
118 resg@gsnSpreadColors = True ; use full colormap
119 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
120 resg@cnMinLevelValF = 0. ; Min level
121 resg@cnMaxLevelValF = 30. ; Max level
122 resg@cnLevelSpacingF = 2. ; interval
124 resg@mpMinLatF = -21.1 ; range to zoom in on
125 resg@mpMaxLatF = 13.8
126 resg@mpMinLonF = 277.28
127 resg@mpMaxLonF = 326.43
128 ;------------------------------------------------------------------------
131 plot_name = "global_ob"
132 title = ob_name+" "+ model_grid
133 resg@tiMainString = title
135 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
136 gsn_define_colormap(wks,"gui_default") ; choose colormap
138 plot = gsn_csm_contour_map_ce(wks,dataob,resg)
141 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
142 ; system("rm "+plot_name+"."+plot_type)
143 system("rm "+plot_name+"-1."+plot_type_new)
144 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
147 ;------------------------------------------------------------------------
148 ;global contour model
150 plot_name = "global_model"
151 title = "Model "+ model_name
152 resg@tiMainString = title
154 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
155 gsn_define_colormap(wks,"gui_default") ; choose colormap
157 plot = gsn_csm_contour_map_ce(wks,datamod,resg)
160 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
161 ; system("rm "+plot_name+"."+plot_type)
162 system("rm "+plot_name+"-1."+plot_type_new)
163 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
166 ;------------------------------------------------------------------------
167 ;global contour model vs ob
169 plot_name = "global_model_vs_ob"
171 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
172 gsn_define_colormap(wks,"gui_default") ; choose colormap
175 plot=new(3,graphic) ; create graphic array
177 resg@gsnFrame = False ; Do not draw plot
178 resg@gsnDraw = False ; Do not advance frame
180 ;(d) compute correlation coef and M score
182 uu = ndtooned(datamod)
183 vv = ndtooned(dataob)
185 good = ind(.not.ismissing(uu) .and. .not.ismissing(vv))
190 ccrG = esccr(ug,vg,0)
193 MG = (ccrG*ccrG)* 5.0
195 bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
196 MG = (1. - (bias/dimsizes(ug)))*5.
199 ; plot correlation coef
202 gRes@txFontHeightF = 0.02
205 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
207 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
208 ;--------------------------------------------------------------------
212 title = ob_name+" "+ model_grid
213 resg@tiMainString = title
215 plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg)
219 title = "Model "+ model_name
220 resg@tiMainString = title
222 plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg)
227 zz = datamod - dataob
228 title = "Model_"+model_name+" - Observed"
230 resg@cnMinLevelValF = -10. ; Min level
231 resg@cnMaxLevelValF = 10. ; Max level
232 resg@cnLevelSpacingF = 2. ; interval
233 resg@tiMainString = title
235 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
237 pres = True ; panel plot mods desired
238 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
239 ; indiv. plots in panel
240 pres@gsnMaximize = True ; fill the page
242 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
244 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
245 ; system("rm "+plot_name+"."+plot_type)
246 ; system("rm "+plot_name+"-1."+plot_type_new)
247 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
252 ;------------------------------------------------------------------------
253 ; temp_name = "temp." + model_name
254 ; system("mkdir -p " + temp_name)
255 ; system("mv *.png " + temp_name)
256 ; system("tar cf "+ temp_name +".tar " + temp_name)
257 ;------------------------------------------------------------------------