|
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 ;************************************************ |
|
9 begin |
|
10 |
|
11 plot_type = "ps" |
|
12 plot_type_new = "png" |
|
13 ;************************************************ |
|
14 ; read in model data |
|
15 ;************************************************ |
|
16 ;film = "b30.061n_1995-2004_ANN_climo_lnd.nc" |
|
17 ;model_name = "b30.061n" |
|
18 ;model_grid = "T31" |
|
19 |
|
20 ;film = "newcn05_ncep_1i_ANN_climo_lnd.nc" |
|
21 ;model_name = "newcn" |
|
22 ;model_grid = "1.9" |
|
23 ;-------------------------------------------- |
|
24 film = "i01.10cn_1948-2004_ANN_climo.nc" |
|
25 model_name = "10cn" |
|
26 model_grid = "T42" |
|
27 |
|
28 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" |
|
29 fm = addfile (dirm+film,"r") |
|
30 |
|
31 ;unit: gC/m2 |
|
32 data1 = fm->LIVESTEMC |
|
33 data2 = fm->DEADSTEMC |
|
34 data3 = fm->LEAFC |
|
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" |
|
40 ;model_grid = "T42" |
|
41 |
|
42 ;dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" |
|
43 ;fm = addfile (dirm+film,"r") |
|
44 |
|
45 ;unit: gC/m2 |
|
46 ;factor_WOODC = 0.7 |
|
47 ;data1 = fm->WOODC |
|
48 ;data2 = fm->LEAFC |
|
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) |
|
56 |
|
57 xm = fm->lon |
|
58 ym = fm->lat |
|
59 ;************************************************ |
|
60 ; read in ob data |
|
61 ;************************************************ |
|
62 ob_name = "LC15_Amazon_Biomass" |
|
63 |
|
64 diro = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/" |
|
65 filo = "amazon_biomass_"+model_grid+".nc" |
|
66 fo = addfile (diro+filo,"r") |
|
67 |
|
68 dataob = fo->BIOMASS |
|
69 xo = fo->lon |
|
70 yo = fo->lat |
|
71 ;************************************************ |
|
72 ; Units for these variables are: |
|
73 ; dataob : MgC/ha |
|
74 ; datamod0 : gC/m2 |
|
75 ; We want to convert these to KgC/m2 |
|
76 ; ha = 100m*100m = 10,000 m2 |
|
77 ; MgC/ha*1000/10,000 = KgC/m2 |
|
78 |
|
79 factor_aboveground = 0.5 |
|
80 factor_unit_ob = 0.1 |
|
81 factor_unit_mod = 0.001 |
|
82 |
|
83 dataob = dataob * factor_aboveground * factor_unit_ob |
|
84 datamod0 = datamod0 * factor_unit_mod |
|
85 |
|
86 dataob@units = "KgC/m2" |
|
87 datamod0@units = "KgC/m2" |
|
88 |
|
89 dataob@long_name = "Amazon Biomass" |
|
90 datamod0@long_name = "Amazon Biomass" |
|
91 ;******************************************************** |
|
92 ; get subset of datamod0 that match dataob |
|
93 |
|
94 nlon = dimsizes(xo) |
|
95 nlat = dimsizes(yo) |
|
96 |
|
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)) |
|
101 |
|
102 datamod = dataob |
|
103 datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR) |
|
104 |
|
105 ; print (datamod) |
|
106 |
|
107 ;---------------------------------------------------------------------- |
|
108 ; global contour |
|
109 |
|
110 ;res |
|
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 |
|
119 |
|
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 |
|
125 |
|
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 ;------------------------------------------------------------------------ |
|
131 ;global contour ob |
|
132 |
|
133 plot_name = "global_ob" |
|
134 title = ob_name+" "+ model_grid |
|
135 resg@tiMainString = title |
|
136 |
|
137 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
138 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
139 |
|
140 plot = gsn_csm_contour_map_ce(wks,dataob,resg) |
|
141 frame(wks) |
|
142 |
|
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) |
|
147 |
|
148 clear (wks) |
|
149 ;------------------------------------------------------------------------ |
|
150 ;global contour model |
|
151 |
|
152 plot_name = "global_model" |
|
153 title = "Model "+ model_name |
|
154 resg@tiMainString = title |
|
155 |
|
156 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
157 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
158 |
|
159 plot = gsn_csm_contour_map_ce(wks,datamod,resg) |
|
160 frame(wks) |
|
161 |
|
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) |
|
166 |
|
167 clear (wks) |
|
168 ;------------------------------------------------------------------------ |
|
169 ;global contour model vs ob |
|
170 |
|
171 plot_name = "global_model_vs_ob" |
|
172 |
|
173 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
174 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
175 |
|
176 delete (plot) |
|
177 plot=new(3,graphic) ; create graphic array |
|
178 |
|
179 resg@gsnFrame = False ; Do not draw plot |
|
180 resg@gsnDraw = False ; Do not advance frame |
|
181 |
|
182 ;(d) compute correlation coef and M score |
|
183 |
|
184 uu = ndtooned(datamod) |
|
185 vv = ndtooned(dataob) |
|
186 |
|
187 good = ind(.not.ismissing(uu) .and. .not.ismissing(vv)) |
|
188 |
|
189 ug = uu(good) |
|
190 vg = vv(good) |
|
191 |
|
192 ccrG = esccr(ug,vg,0) |
|
193 ; print (ccrG) |
|
194 |
|
195 MG = (ccrG*ccrG)* 5.0 |
|
196 ; new eq |
|
197 bias = sum(abs(ug-vg)/(abs(ug)+abs(vg))) |
|
198 MG = (1. - (bias/dimsizes(ug)))*5. |
|
199 |
|
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) |
|
203 print (M_biomass) |
|
204 |
|
205 ; plot correlation coef |
|
206 |
|
207 gRes = True |
|
208 gRes@txFontHeightF = 0.02 |
|
209 gRes@txAngleF = 90 |
|
210 |
|
211 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")" |
|
212 |
|
213 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
|
214 ;-------------------------------------------------------------------- |
|
215 |
|
216 ;(a) ob |
|
217 |
|
218 title = ob_name+" "+ model_grid |
|
219 resg@tiMainString = title |
|
220 |
|
221 plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg) |
|
222 |
|
223 ;(b) model |
|
224 |
|
225 title = "Model "+ model_name |
|
226 resg@tiMainString = title |
|
227 |
|
228 plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg) |
|
229 |
|
230 ;(c) model-ob |
|
231 |
|
232 zz = datamod |
|
233 zz = datamod - dataob |
|
234 title = "Model_"+model_name+" - Observed" |
|
235 |
|
236 resg@cnMinLevelValF = -10. ; Min level |
|
237 resg@cnMaxLevelValF = 10. ; Max level |
|
238 resg@cnLevelSpacingF = 2. ; interval |
|
239 resg@tiMainString = title |
|
240 |
|
241 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) |
|
242 |
|
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 |
|
247 |
|
248 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
|
249 |
|
250 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
251 "rm "+plot_name+"."+plot_type) |
|
252 |
|
253 clear (wks) |
|
254 delete (plot) |
|
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 ;------------------------------------------------------------------- |
|
262 end |