|
1 ;******************************************************** |
|
2 ; histogram normalized by rain and compute correleration |
|
3 ;******************************************************** |
|
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test" |
|
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test" |
|
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" |
|
7 |
|
8 procedure pminmax(data:numeric,name:string) |
|
9 begin |
|
10 print ("min/max " + name + " = " + min(data) + "/" + max(data)) |
|
11 if(isatt(data,"units")) then |
|
12 print (name + " units = " + data@units) |
|
13 end if |
|
14 end |
|
15 |
|
16 ; Main code. |
|
17 begin |
|
18 |
|
19 nclass = 20 |
|
20 |
|
21 plot_type = "ps" |
|
22 plot_type_new = "png" |
|
23 |
|
24 ;************************************************ |
|
25 ; read in data: model |
|
26 ;************************************************ |
|
27 ;film = "b30.061n_1995-2004_MONS_climo_lnd.nc" |
|
28 ;model_name = "b30.061n" |
|
29 ;model_grid = "T31" |
|
30 |
|
31 ;film = "newcn05_ncep_1i_MONS_climo_lnd.nc" |
|
32 ;model_name = "newcn" |
|
33 ;model_grid = "1.9" |
|
34 |
|
35 ;film = "i01.06cn_1798-2004_MONS_climo.nc" |
|
36 ;model_name = "06cn" |
|
37 ;model_grid = "T42" |
|
38 |
|
39 ;film = "i01.06casa_1798-2004_MONS_climo.nc" |
|
40 ;model_name = "06casa" |
|
41 ;model_grid = "T42" |
|
42 |
|
43 film = "i01.10cn_1948-2004_MONS_climo.nc" |
|
44 model_name = "10cn" |
|
45 model_grid = "T42" |
|
46 |
|
47 ;film = "i01.10casa_1948-2004_MONS_climo.nc" |
|
48 ;model_name = "10casa" |
|
49 ;model_grid = "T42" |
|
50 |
|
51 html_name = "table.html." + model_name |
|
52 html_new = html_name +".new" |
|
53 system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \ |
|
54 "mv -f "+html_new+" "+html_name) |
|
55 ;------------------------------------------------ |
|
56 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" |
|
57 fm = addfile(dirm+film,"r") |
|
58 |
|
59 laimod = fm->TLAI |
|
60 |
|
61 ;************************************************ |
|
62 ; read in data: observed |
|
63 ;************************************************ |
|
64 |
|
65 ob_name = "MODIS MOD 15A2 2000-2005" |
|
66 |
|
67 diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/" |
|
68 filo1 = "land_class_"+model_grid+".nc" |
|
69 filo2 = "LAI_2000-2005_MONS_"+model_grid+".nc" |
|
70 |
|
71 fo1 = addfile(diro+filo1,"r") |
|
72 fo2 = addfile(diro+filo2,"r") |
|
73 |
|
74 classob = tofloat(fo1->LAND_CLASS) |
|
75 laiob = fo2->LAI |
|
76 |
|
77 ;******************************************************************* |
|
78 ; for plotting table |
|
79 ;******************************************************************* |
|
80 ; table header name |
|
81 table_header_name = "LAI" |
|
82 |
|
83 ; column (not including header column) |
|
84 |
|
85 col_header1 = (/"Mean","Max","Phase","Growth"/) |
|
86 col_header2 = (/"ob","model","M" \ |
|
87 ,"ob","model","M" \ |
|
88 ,"ob","model","M" \ |
|
89 ,"ob","model","M" \ |
|
90 /) |
|
91 |
|
92 ncol1 = dimsizes(col_header1) |
|
93 ncol2 = dimsizes(col_header2) |
|
94 ncol = ncol2 |
|
95 |
|
96 ; row (not including header row) |
|
97 row_header = (/"Water Bodies" \ |
|
98 ,"Evergreen Needleleaf Forests" \ |
|
99 ,"Evergreen Broadleaf Forests" \ |
|
100 ,"Deciduous Needleleaf Forest" \ |
|
101 ,"Deciduous Broadleaf Forests" \ |
|
102 ,"Mixed Forests" \ |
|
103 ,"Closed Bushlands" \ |
|
104 ,"Open Bushlands" \ |
|
105 ,"Woody Savannas (S. Hem.)" \ |
|
106 ,"Savannas (S. Hem.)" \ |
|
107 ,"Grasslands" \ |
|
108 ,"Permanent Wetlands" \ |
|
109 ,"Croplands" \ |
|
110 ,"Urban and Built-Up" \ |
|
111 ,"Cropland/Natural Vegetation Mosaic" \ |
|
112 ,"Permanent Snow and Ice" \ |
|
113 ,"Barren or Sparsely Vegetated" \ |
|
114 ,"Unclassified" \ |
|
115 ,"Woody Savannas (N. Hem.)" \ |
|
116 ,"Savannas (N. Hem.)" \ |
|
117 ,"All biome average" \ |
|
118 /) |
|
119 nrow = dimsizes(row_header) |
|
120 |
|
121 ; arrays to be passed to table. |
|
122 value = new ((/nrow, ncol/),string ) |
|
123 |
|
124 table_length = 0.995 |
|
125 |
|
126 ; Table header |
|
127 ncr1 = (/1,1/) ; 1 row, 1 column |
|
128 xx1 = (/0.005,0.25/) ; Start and end X |
|
129 yy1 = (/0.900,0.995/) ; Start and end Y |
|
130 text1 = table_header_name |
|
131 res1 = True |
|
132 res1@txFontHeightF = 0.03 |
|
133 res1@gsFillColor = "CornFlowerBlue" |
|
134 |
|
135 ; Column header (equally space in x2) |
|
136 ncr21 = (/1,ncol1/) ; 1 rows, 4 columns |
|
137 xx21 = (/xx1(1),0.995/) ; start from end of x1 |
|
138 yy21 = (/0.9475,0.995/) ; half of y1 |
|
139 text21 = col_header1 |
|
140 res21 = True |
|
141 res21@txFontHeightF = 0.015 |
|
142 res21@gsFillColor = "Gray" |
|
143 |
|
144 ncr22 = (/1,ncol2/) ; 1 rows, 12 columns |
|
145 xx22 = xx21 ; start from end of x1 |
|
146 yy22 = (/0.900,0.9475/) ; half of y1 |
|
147 text22 = col_header2 |
|
148 res22 = True |
|
149 res22@txFontHeightF = 0.015 |
|
150 res22@gsFillColor = "Gray" |
|
151 |
|
152 ; Row header (equally space in y2) |
|
153 ncr3 = (/nrow,1/) ; 20 rows, 1 columns |
|
154 xx3 = xx1 ; same as x1 |
|
155 yy3 = (/1.0-table_length,0.900/) ; end at start of y1 |
|
156 text3 = row_header |
|
157 res3 = True |
|
158 res3@txFontHeightF = 0.01 |
|
159 res3@gsFillColor = "Gray" |
|
160 |
|
161 ; Main table body |
|
162 ncr4 = (/nrow,ncol/) ; 5 rows, 5 columns |
|
163 xx4 = xx21 ; Start and end x |
|
164 yy4 = yy3 ; Start and end Y |
|
165 text4 = new((/nrow,ncol/),string) |
|
166 |
|
167 color_fill4 = new((/nrow,ncol/),string) |
|
168 color_fill4 = "white" |
|
169 color_fill4(nrow-1,:) = "CornFlowerBlue" |
|
170 |
|
171 res4 = True ; Set up resource list |
|
172 ; res4@gsnDebug = True ; Useful to print NDC row,col values used. |
|
173 res4@txFontHeightF = 0.015 |
|
174 res4@gsFillColor = color_fill4 |
|
175 |
|
176 delete (color_fill4) |
|
177 |
|
178 ;************************************************ |
|
179 ; plot global land class: observed |
|
180 ;************************************************ |
|
181 ;global res |
|
182 |
|
183 resg = True ; Use plot options |
|
184 resg@cnFillOn = True ; Turn on color fill |
|
185 resg@gsnSpreadColors = True ; use full colormap |
|
186 ; resg@cnFillMode = "RasterFill" ; Turn on raster color |
|
187 ; resg@lbLabelAutoStride = True |
|
188 resg@cnLinesOn = False ; Turn off contourn lines |
|
189 resg@mpFillOn = False ; Turn off map fill |
|
190 |
|
191 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
192 resg@cnMinLevelValF = 1. ; Min level |
|
193 resg@cnMaxLevelValF = 19. ; Max level |
|
194 resg@cnLevelSpacingF = 1. ; interval |
|
195 |
|
196 ;global contour ob |
|
197 classob@_FillValue = 1.e+36 |
|
198 classob = where(classob.eq.0,classob@_FillValue,classob) |
|
199 |
|
200 plot_name = "global_class_ob" |
|
201 title = ob_name |
|
202 resg@tiMainString = title |
|
203 |
|
204 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
205 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
206 |
|
207 plot = gsn_csm_contour_map_ce(wks,classob,resg) |
|
208 frame(wks) |
|
209 |
|
210 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
211 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ |
|
212 "rm "+plot_name+"-*."+plot_type_new+";"+ \ |
|
213 "rm "+plot_name+"."+plot_type) |
|
214 |
|
215 clear (wks) |
|
216 ;******************************************************************* |
|
217 ; Calculate "nice" bins for binning the data in equally spaced ranges |
|
218 ;******************************************************************** |
|
219 nclassn = nclass + 1 |
|
220 range = fspan(0,nclassn-1,nclassn) |
|
221 ; print (range) |
|
222 |
|
223 ; Use this range information to grab all the values in a |
|
224 ; particular range, and then take an average. |
|
225 |
|
226 nr = dimsizes(range) |
|
227 nx = nr-1 |
|
228 xvalues = new((/2,nx/),float) |
|
229 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2. |
|
230 dx = xvalues(0,1) - xvalues(0,0) ; range width |
|
231 dx4 = dx/4 ; 1/4 of the range |
|
232 xvalues(1,:) = xvalues(0,:) - dx/5. |
|
233 ;----------------------------------------------------------------- |
|
234 ;(A) mean |
|
235 ;-------------------------------------------------------------------- |
|
236 ; get data |
|
237 |
|
238 laiob_mean = dim_avg_Wrap(laiob(lat|:,lon|:,time|:)) |
|
239 laimod_mean = dim_avg_Wrap(laimod(lat|:,lon|:,time|:)) |
|
240 |
|
241 DATA11_1D = ndtooned(classob) |
|
242 DATA12_1D = ndtooned(laiob_mean) |
|
243 DATA22_1D = ndtooned(laimod_mean) |
|
244 |
|
245 yvalues = new((/2,nx/),float) |
|
246 mn_yvalues = new((/2,nx/),float) |
|
247 mx_yvalues = new((/2,nx/),float) |
|
248 |
|
249 do nd=0,1 |
|
250 |
|
251 ; See if we are doing model or observational data. |
|
252 |
|
253 if(nd.eq.0) then |
|
254 data_ob = DATA11_1D |
|
255 data_mod = DATA12_1D |
|
256 else |
|
257 data_ob = DATA11_1D |
|
258 data_mod = DATA22_1D |
|
259 end if |
|
260 |
|
261 ; Loop through each range and check for values. |
|
262 |
|
263 do i=0,nr-2 |
|
264 if (i.ne.(nr-2)) then |
|
265 ; print("") |
|
266 ; print("In range ["+range(i)+","+range(i+1)+")") |
|
267 idx = ind((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1))) |
|
268 else |
|
269 ; print("") |
|
270 ; print("In range ["+range(i)+",)") |
|
271 idx = ind(data_ob.ge.range(i)) |
|
272 end if |
|
273 |
|
274 ; Calculate average, and get min and max. |
|
275 |
|
276 if(.not.any(ismissing(idx))) then |
|
277 yvalues(nd,i) = avg(data_mod(idx)) |
|
278 mn_yvalues(nd,i) = min(data_mod(idx)) |
|
279 mx_yvalues(nd,i) = max(data_mod(idx)) |
|
280 count = dimsizes(idx) |
|
281 else |
|
282 count = 0 |
|
283 yvalues(nd,i) = yvalues@_FillValue |
|
284 mn_yvalues(nd,i) = yvalues@_FillValue |
|
285 mx_yvalues(nd,i) = yvalues@_FillValue |
|
286 end if |
|
287 |
|
288 ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) |
|
289 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) |
|
290 |
|
291 ; Clean up for next time in loop. |
|
292 |
|
293 delete(idx) |
|
294 end do |
|
295 delete(data_ob) |
|
296 delete(data_mod) |
|
297 end do |
|
298 ;----------------------------------------------------------------- |
|
299 ; compute correlation coef and M score |
|
300 |
|
301 u = yvalues(0,:) |
|
302 v = yvalues(1,:) |
|
303 |
|
304 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
305 uu = u(good) |
|
306 vv = v(good) |
|
307 |
|
308 ccrMean = esccr(uu,vv,0) |
|
309 ; print (ccrMean) |
|
310 |
|
311 ; new eq |
|
312 bias = sum(abs(vv-uu)/(vv+uu)) |
|
313 Mmean = (1.- (bias/dimsizes(uu)))*5. |
|
314 |
|
315 M_lai_mean = sprintf("%.2f", Mmean) |
|
316 system("sed s#M_lai_mean#"+M_lai_mean+"# "+html_name+" > "+html_new+";"+ \ |
|
317 "mv -f "+html_new+" "+html_name) |
|
318 print (M_lai_mean) |
|
319 |
|
320 do i=0,nrow-2 |
|
321 text4(i,0) = sprintf("%5.2f",u(i)) |
|
322 text4(i,1) = sprintf("%5.2f",v(i)) |
|
323 text4(i,2) = "-" |
|
324 end do |
|
325 text4(nrow-1,0) = sprintf("%5.2f",avg(u)) |
|
326 text4(nrow-1,1) = sprintf("%5.2f",avg(v)) |
|
327 text4(nrow-1,2) = sprintf("%5.2f",Mmean) |
|
328 |
|
329 delete (u) |
|
330 delete (v) |
|
331 delete (uu) |
|
332 delete (vv) |
|
333 ;-------------------------------------------------------------------- |
|
334 ; histogram res |
|
335 |
|
336 resm = True |
|
337 resm@gsnMaximize = True |
|
338 resm@gsnDraw = False |
|
339 resm@gsnFrame = False |
|
340 resm@xyMarkLineMode = "Markers" |
|
341 resm@xyMarkerSizeF = 0.014 |
|
342 resm@xyMarker = 16 |
|
343 resm@xyMarkerColors = (/"Brown","Blue"/) |
|
344 ; resm@trYMinF = min(mn_yvalues) - 10. |
|
345 ; resm@trYMaxF = max(mx_yvalues) + 10. |
|
346 resm@trYMinF = min(mn_yvalues) - 2 |
|
347 resm@trYMaxF = max(mx_yvalues) + 4 |
|
348 |
|
349 resm@tiYAxisString = "Mean LAI (Leaf Area Index)" |
|
350 resm@tiXAxisString = "Land Cover Type" |
|
351 |
|
352 max_bar = new((/2,nx/),graphic) |
|
353 min_bar = new((/2,nx/),graphic) |
|
354 max_cap = new((/2,nx/),graphic) |
|
355 min_cap = new((/2,nx/),graphic) |
|
356 |
|
357 lnres = True |
|
358 line_colors = (/"brown","blue"/) |
|
359 ;------------------------------------------------------------------ |
|
360 ; Start the graphics. |
|
361 |
|
362 plot_name = "histogram_mean" |
|
363 title = model_name + " vs Observed" |
|
364 resm@tiMainString = title |
|
365 |
|
366 wks = gsn_open_wks (plot_type,plot_name) |
|
367 ;----------------------------- |
|
368 ; Add a boxed legend using the more simple method |
|
369 |
|
370 resm@pmLegendDisplayMode = "Always" |
|
371 ; resm@pmLegendWidthF = 0.1 |
|
372 resm@pmLegendWidthF = 0.08 |
|
373 resm@pmLegendHeightF = 0.05 |
|
374 resm@pmLegendOrthogonalPosF = -1.17 |
|
375 ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) |
|
376 ; resm@pmLegendParallelPosF = 0.18 |
|
377 resm@pmLegendParallelPosF = 0.88 ;(rightward) |
|
378 |
|
379 ; resm@lgPerimOn = False |
|
380 resm@lgLabelFontHeightF = 0.015 |
|
381 resm@xyExplicitLegendLabels = (/"observed",model_name/) |
|
382 ;----------------------------- |
|
383 tRes = True |
|
384 tRes@txFontHeightF = 0.025 |
|
385 |
|
386 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")" |
|
387 |
|
388 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) |
|
389 |
|
390 xy = gsn_csm_xy(wks,xvalues,yvalues,resm) |
|
391 ;------------------------------- |
|
392 ;Attach the vertical bar and the horizontal cap line |
|
393 |
|
394 do nd=0,1 |
|
395 lnres@gsLineColor = line_colors(nd) |
|
396 do i=0,nx-1 |
|
397 |
|
398 if(.not.ismissing(mn_yvalues(nd,i)).and. \ |
|
399 .not.ismissing(mx_yvalues(nd,i))) then |
|
400 |
|
401 ; Attach the vertical bar, both above and below the marker. |
|
402 |
|
403 x1 = xvalues(nd,i) |
|
404 y1 = yvalues(nd,i) |
|
405 y2 = mn_yvalues(nd,i) |
|
406 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
407 |
|
408 y2 = mx_yvalues(nd,i) |
|
409 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
410 |
|
411 ; Attach the horizontal cap line, both above and below the marker. |
|
412 |
|
413 x1 = xvalues(nd,i) - dx4 |
|
414 x2 = xvalues(nd,i) + dx4 |
|
415 y1 = mn_yvalues(nd,i) |
|
416 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
417 |
|
418 y1 = mx_yvalues(nd,i) |
|
419 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
420 end if |
|
421 end do |
|
422 end do |
|
423 |
|
424 draw(xy) |
|
425 frame(wks) |
|
426 |
|
427 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
428 "rm "+plot_name+"."+plot_type) |
|
429 |
|
430 clear (wks) |
|
431 |
|
432 delete (DATA11_1D) |
|
433 delete (DATA12_1D) |
|
434 delete (DATA22_1D) |
|
435 ;delete (range) |
|
436 ;delete (xvalues) |
|
437 delete (yvalues) |
|
438 delete (mn_yvalues) |
|
439 delete (mx_yvalues) |
|
440 delete (good) |
|
441 delete (max_bar) |
|
442 delete (min_bar) |
|
443 delete (max_cap) |
|
444 delete (min_cap) |
|
445 ;----------------------------------------------------------------- |
|
446 ;global res |
|
447 |
|
448 resg = True ; Use plot options |
|
449 resg@cnFillOn = True ; Turn on color fill |
|
450 resg@gsnSpreadColors = True ; use full colormap |
|
451 ; resg@cnFillMode = "RasterFill" ; Turn on raster color |
|
452 ; resg@lbLabelAutoStride = True |
|
453 resg@cnLinesOn = False ; Turn off contourn lines |
|
454 resg@mpFillOn = False ; Turn off map fill |
|
455 |
|
456 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
457 resg@cnMinLevelValF = 0. ; Min level |
|
458 resg@cnMaxLevelValF = 10. ; Max level |
|
459 resg@cnLevelSpacingF = 1. ; interval |
|
460 |
|
461 ;global contour ob |
|
462 |
|
463 delta = 0.00001 |
|
464 laiob_mean = where(ismissing(laiob_mean).and.(ismissing(laimod_mean).or.(laimod_mean.lt.delta)),0.,laiob_mean) |
|
465 |
|
466 plot_name = "global_mean_ob" |
|
467 title = ob_name |
|
468 resg@tiMainString = title |
|
469 |
|
470 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
471 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
472 |
|
473 plot = gsn_csm_contour_map_ce(wks,laiob_mean,resg) |
|
474 frame(wks) |
|
475 |
|
476 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
477 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ |
|
478 "rm "+plot_name+"-*."+plot_type_new+";"+ \ |
|
479 "rm "+plot_name+"."+plot_type) |
|
480 |
|
481 clear (wks) |
|
482 ;------------------------------------------------------------------------ |
|
483 ;global contour model |
|
484 |
|
485 plot_name = "global_mean_model" |
|
486 title = "Model " + model_name |
|
487 resg@tiMainString = title |
|
488 |
|
489 wks = gsn_open_wks (plot_type,plot_name) |
|
490 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
491 |
|
492 plot = gsn_csm_contour_map_ce(wks,laimod_mean,resg) |
|
493 frame(wks) |
|
494 |
|
495 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
496 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ |
|
497 "rm "+plot_name+"-*."+plot_type_new+";"+ \ |
|
498 "rm "+plot_name+"."+plot_type) |
|
499 |
|
500 clear (wks) |
|
501 ;----------------------------------------------------------------- |
|
502 ;(B) max |
|
503 ;-------------------------------------------------------------------- |
|
504 ; get data |
|
505 |
|
506 ; observed |
|
507 laiob_max = laiob(0,:,:) |
|
508 s = laiob(:,0,0) |
|
509 laiob_max@long_name = "Leaf Area Index Max" |
|
510 |
|
511 dsizes_z = dimsizes(laiob) |
|
512 nlat = dsizes_z(1) |
|
513 nlon = dsizes_z(2) |
|
514 |
|
515 do j = 0,nlat-1 |
|
516 do i = 0,nlon-1 |
|
517 s = laiob(:,j,i) |
|
518 laiob_max(j,i) = max(s) |
|
519 end do |
|
520 end do |
|
521 |
|
522 ; print (min(y)+"/"+max(y)) |
|
523 delete (s) |
|
524 delete (dsizes_z) |
|
525 ;------------------------- |
|
526 ; model |
|
527 laimod_max = laimod(0,:,:) |
|
528 s = laimod(:,0,0) |
|
529 laimod_max@long_name = "Leaf Area Index Max" |
|
530 |
|
531 dsizes_z = dimsizes(laimod) |
|
532 nlat = dsizes_z(1) |
|
533 nlon = dsizes_z(2) |
|
534 |
|
535 do j = 0,nlat-1 |
|
536 do i = 0,nlon-1 |
|
537 s = laimod(:,j,i) |
|
538 laimod_max(j,i) = max(s) |
|
539 end do |
|
540 end do |
|
541 |
|
542 ; print (min(laimod_max)+"/"+max(laimod_max)) |
|
543 delete (s) |
|
544 delete (dsizes_z) |
|
545 ;------------------------ |
|
546 DATA11_1D = ndtooned(classob) |
|
547 DATA12_1D = ndtooned(laiob_max) |
|
548 DATA22_1D = ndtooned(laimod_max) |
|
549 |
|
550 yvalues = new((/2,nx/),float) |
|
551 mn_yvalues = new((/2,nx/),float) |
|
552 mx_yvalues = new((/2,nx/),float) |
|
553 |
|
554 do nd=0,1 |
|
555 |
|
556 ; See if we are doing model or observational data. |
|
557 |
|
558 if(nd.eq.0) then |
|
559 data_ob = DATA11_1D |
|
560 data_mod = DATA12_1D |
|
561 else |
|
562 data_ob = DATA11_1D |
|
563 data_mod = DATA22_1D |
|
564 end if |
|
565 |
|
566 ; Loop through each range and check for values. |
|
567 |
|
568 do i=0,nr-2 |
|
569 if (i.ne.(nr-2)) then |
|
570 ; print("") |
|
571 ; print("In range ["+range(i)+","+range(i+1)+")") |
|
572 idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) |
|
573 else |
|
574 ; print("") |
|
575 ; print("In range ["+range(i)+",)") |
|
576 idx = ind(range(i).le.data_ob) |
|
577 end if |
|
578 |
|
579 ; Calculate average, and get min and max. |
|
580 |
|
581 if(.not.any(ismissing(idx))) then |
|
582 yvalues(nd,i) = avg(data_mod(idx)) |
|
583 mn_yvalues(nd,i) = min(data_mod(idx)) |
|
584 mx_yvalues(nd,i) = max(data_mod(idx)) |
|
585 count = dimsizes(idx) |
|
586 else |
|
587 count = 0 |
|
588 yvalues(nd,i) = yvalues@_FillValue |
|
589 mn_yvalues(nd,i) = yvalues@_FillValue |
|
590 mx_yvalues(nd,i) = yvalues@_FillValue |
|
591 end if |
|
592 |
|
593 ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) |
|
594 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) |
|
595 |
|
596 ; Clean up for next time in loop. |
|
597 |
|
598 delete(idx) |
|
599 end do |
|
600 delete(data_ob) |
|
601 delete(data_mod) |
|
602 end do |
|
603 ;----------------------------------------------------------------- |
|
604 ; compute correlation coef and M score |
|
605 |
|
606 u = yvalues(0,:) |
|
607 v = yvalues(1,:) |
|
608 |
|
609 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
610 uu = u(good) |
|
611 vv = v(good) |
|
612 |
|
613 ccrMax = esccr(uu,vv,0) |
|
614 ; print (ccrMax) |
|
615 |
|
616 ; new eq |
|
617 bias = sum(abs(vv-uu)/(vv+uu)) |
|
618 Mmax = (1.- (bias/dimsizes(uu)))*5. |
|
619 |
|
620 M_lai_max = sprintf("%.2f", Mmax) |
|
621 system("sed s#M_lai_max#"+M_lai_max+"# "+html_name+" > "+html_new+";"+ \ |
|
622 "mv -f "+html_new+" "+html_name) |
|
623 print (M_lai_max) |
|
624 |
|
625 do i=0,nrow-2 |
|
626 text4(i,3) = sprintf("%5.2f",u(i)) |
|
627 text4(i,4) = sprintf("%5.2f",v(i)) |
|
628 text4(i,5) = "-" |
|
629 end do |
|
630 text4(nrow-1,3) = sprintf("%5.2f",avg(u)) |
|
631 text4(nrow-1,4) = sprintf("%5.2f",avg(v)) |
|
632 text4(nrow-1,5) = sprintf("%5.2f",Mmax) |
|
633 |
|
634 delete (u) |
|
635 delete (v) |
|
636 delete (uu) |
|
637 delete (vv) |
|
638 ;-------------------------------------------------------------------- |
|
639 ; histogram res |
|
640 |
|
641 resm = True |
|
642 resm@gsnMaximize = True |
|
643 resm@gsnDraw = False |
|
644 resm@gsnFrame = False |
|
645 resm@xyMarkLineMode = "Markers" |
|
646 resm@xyMarkerSizeF = 0.014 |
|
647 resm@xyMarker = 16 |
|
648 resm@xyMarkerColors = (/"Brown","Blue"/) |
|
649 ; resm@trYMinF = min(mn_yvalues) - 10. |
|
650 ; resm@trYMaxF = max(mx_yvalues) + 10. |
|
651 resm@trYMinF = min(mn_yvalues) - 2 |
|
652 resm@trYMaxF = max(mx_yvalues) + 4 |
|
653 |
|
654 resm@tiYAxisString = "Max LAI (Leaf Area Index)" |
|
655 resm@tiXAxisString = "Land Cover Type" |
|
656 |
|
657 max_bar = new((/2,nx/),graphic) |
|
658 min_bar = new((/2,nx/),graphic) |
|
659 max_cap = new((/2,nx/),graphic) |
|
660 min_cap = new((/2,nx/),graphic) |
|
661 |
|
662 lnres = True |
|
663 line_colors = (/"brown","blue"/) |
|
664 ;------------------------------------------------------------------ |
|
665 ; Start the graphics. |
|
666 |
|
667 plot_name = "histogram_max" |
|
668 title = model_name + " vs Observed" |
|
669 resm@tiMainString = title |
|
670 |
|
671 wks = gsn_open_wks (plot_type,plot_name) |
|
672 ;----------------------------- |
|
673 ; Add a boxed legend using the more simple method |
|
674 |
|
675 resm@pmLegendDisplayMode = "Always" |
|
676 ; resm@pmLegendWidthF = 0.1 |
|
677 resm@pmLegendWidthF = 0.08 |
|
678 resm@pmLegendHeightF = 0.05 |
|
679 resm@pmLegendOrthogonalPosF = -1.17 |
|
680 ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) |
|
681 ; resm@pmLegendParallelPosF = 0.18 |
|
682 resm@pmLegendParallelPosF = 0.88 ;(rightward) |
|
683 |
|
684 ; resm@lgPerimOn = False |
|
685 resm@lgLabelFontHeightF = 0.015 |
|
686 resm@xyExplicitLegendLabels = (/"observed",model_name/) |
|
687 ;----------------------------- |
|
688 tRes = True |
|
689 tRes@txFontHeightF = 0.025 |
|
690 |
|
691 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")" |
|
692 |
|
693 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) |
|
694 |
|
695 xy = gsn_csm_xy(wks,xvalues,yvalues,resm) |
|
696 ;------------------------------- |
|
697 ;Attach the vertical bar and the horizontal cap line |
|
698 |
|
699 do nd=0,1 |
|
700 lnres@gsLineColor = line_colors(nd) |
|
701 do i=0,nx-1 |
|
702 |
|
703 if(.not.ismissing(mn_yvalues(nd,i)).and. \ |
|
704 .not.ismissing(mx_yvalues(nd,i))) then |
|
705 |
|
706 ; Attach the vertical bar, both above and below the marker. |
|
707 |
|
708 x1 = xvalues(nd,i) |
|
709 y1 = yvalues(nd,i) |
|
710 y2 = mn_yvalues(nd,i) |
|
711 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
712 |
|
713 y2 = mx_yvalues(nd,i) |
|
714 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
715 |
|
716 ; Attach the horizontal cap line, both above and below the marker. |
|
717 |
|
718 x1 = xvalues(nd,i) - dx4 |
|
719 x2 = xvalues(nd,i) + dx4 |
|
720 y1 = mn_yvalues(nd,i) |
|
721 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
722 |
|
723 y1 = mx_yvalues(nd,i) |
|
724 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
725 end if |
|
726 end do |
|
727 end do |
|
728 |
|
729 draw(xy) |
|
730 frame(wks) |
|
731 |
|
732 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
733 "rm "+plot_name+"."+plot_type) |
|
734 |
|
735 clear (wks) |
|
736 |
|
737 delete (DATA11_1D) |
|
738 delete (DATA12_1D) |
|
739 delete (DATA22_1D) |
|
740 ;delete (range) |
|
741 ;delete (xvalues) |
|
742 delete (yvalues) |
|
743 delete (mn_yvalues) |
|
744 delete (mx_yvalues) |
|
745 delete (good) |
|
746 delete (max_bar) |
|
747 delete (min_bar) |
|
748 delete (max_cap) |
|
749 delete (min_cap) |
|
750 ;----------------------------------------------------------------- |
|
751 ;global res |
|
752 |
|
753 resg = True ; Use plot options |
|
754 resg@cnFillOn = True ; Turn on color fill |
|
755 resg@gsnSpreadColors = True ; use full colormap |
|
756 ; resg@cnFillMode = "RasterFill" ; Turn on raster color |
|
757 ; resg@lbLabelAutoStride = True |
|
758 resg@cnLinesOn = False ; Turn off contourn lines |
|
759 resg@mpFillOn = False ; Turn off map fill |
|
760 |
|
761 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
762 resg@cnMinLevelValF = 0. ; Min level |
|
763 resg@cnMaxLevelValF = 10. ; Max level |
|
764 resg@cnLevelSpacingF = 1. ; interval |
|
765 |
|
766 ;global contour ob |
|
767 |
|
768 delta = 0.00001 |
|
769 laiob_max = where(ismissing(laiob_max).and.(ismissing(laimod_max).or.(laimod_max.lt.delta)),0.,laiob_max) |
|
770 |
|
771 plot_name = "global_max_ob" |
|
772 title = ob_name |
|
773 resg@tiMainString = title |
|
774 |
|
775 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
776 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
777 |
|
778 plot_max = gsn_csm_contour_map_ce(wks,laiob_max,resg) |
|
779 frame(wks) |
|
780 |
|
781 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
782 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ |
|
783 "rm "+plot_name+"-*."+plot_type_new+";"+ \ |
|
784 "rm "+plot_name+"."+plot_type) |
|
785 |
|
786 clear (wks) |
|
787 ;------------------------------------------------------------------------ |
|
788 ;global contour model |
|
789 |
|
790 plot_name = "global_max_model" |
|
791 title = "Model " + model_name |
|
792 resg@tiMainString = title |
|
793 |
|
794 wks = gsn_open_wks (plot_type,plot_name) |
|
795 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
796 |
|
797 plot_max = gsn_csm_contour_map_ce(wks,laimod_max,resg) |
|
798 frame(wks) |
|
799 |
|
800 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
801 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ |
|
802 "rm "+plot_name+"-*."+plot_type_new+";"+ \ |
|
803 "rm "+plot_name+"."+plot_type) |
|
804 |
|
805 clear (wks) |
|
806 ;------------------------------------------------------------------------ |
|
807 ;(C) phase |
|
808 ;-------------------------------------------------------------------- |
|
809 ; get data |
|
810 |
|
811 ; observed |
|
812 laiob_phase = laiob(0,:,:) |
|
813 s = laiob(:,0,0) |
|
814 laiob_phase@long_name = "Leaf Area Index Max Month" |
|
815 |
|
816 dsizes_z = dimsizes(laiob) |
|
817 nlat = dsizes_z(1) |
|
818 nlon = dsizes_z(2) |
|
819 |
|
820 do j = 0,nlat-1 |
|
821 do i = 0,nlon-1 |
|
822 s = laiob(:,j,i) |
|
823 laiob_phase(j,i) = maxind(s) + 1 |
|
824 end do |
|
825 end do |
|
826 |
|
827 ; print (min(laiob_phase)+"/"+max(laiob_phase)) |
|
828 delete (s) |
|
829 delete (dsizes_z) |
|
830 ;------------------------- |
|
831 ; model |
|
832 laimod_phase = laimod(0,:,:) |
|
833 s = laimod(:,0,0) |
|
834 laimod_phase@long_name = "Leaf Area Index Max Month" |
|
835 |
|
836 dsizes_z = dimsizes(laimod) |
|
837 nlat = dsizes_z(1) |
|
838 nlon = dsizes_z(2) |
|
839 |
|
840 do j = 0,nlat-1 |
|
841 do i = 0,nlon-1 |
|
842 s = laimod(:,j,i) |
|
843 laimod_phase(j,i) = maxind(s) + 1 |
|
844 end do |
|
845 end do |
|
846 |
|
847 ; print (min(laimod_phase)+"/"+max(laimod_phase)) |
|
848 delete (s) |
|
849 delete (dsizes_z) |
|
850 ;------------------------ |
|
851 DATA11_1D = ndtooned(classob) |
|
852 DATA12_1D = ndtooned(laiob_phase) |
|
853 DATA22_1D = ndtooned(laimod_phase) |
|
854 |
|
855 yvalues = new((/2,nx/),float) |
|
856 mn_yvalues = new((/2,nx/),float) |
|
857 mx_yvalues = new((/2,nx/),float) |
|
858 |
|
859 do nd=0,1 |
|
860 |
|
861 ; See if we are doing model or observational data. |
|
862 |
|
863 if(nd.eq.0) then |
|
864 data_ob = DATA11_1D |
|
865 data_mod = DATA12_1D |
|
866 else |
|
867 data_ob = DATA11_1D |
|
868 data_mod = DATA22_1D |
|
869 end if |
|
870 |
|
871 ; Loop through each range and check for values. |
|
872 |
|
873 do i=0,nr-2 |
|
874 if (i.ne.(nr-2)) then |
|
875 ; print("") |
|
876 ; print("In range ["+range(i)+","+range(i+1)+")") |
|
877 idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) |
|
878 else |
|
879 ; print("") |
|
880 ; print("In range ["+range(i)+",)") |
|
881 idx = ind(range(i).le.data_ob) |
|
882 end if |
|
883 |
|
884 ; Calculate average, and get min and max. |
|
885 |
|
886 if(.not.any(ismissing(idx))) then |
|
887 yvalues(nd,i) = avg(data_mod(idx)) |
|
888 mn_yvalues(nd,i) = min(data_mod(idx)) |
|
889 mx_yvalues(nd,i) = max(data_mod(idx)) |
|
890 count = dimsizes(idx) |
|
891 else |
|
892 count = 0 |
|
893 yvalues(nd,i) = yvalues@_FillValue |
|
894 mn_yvalues(nd,i) = yvalues@_FillValue |
|
895 mx_yvalues(nd,i) = yvalues@_FillValue |
|
896 end if |
|
897 |
|
898 ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) |
|
899 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) |
|
900 |
|
901 ; Clean up for next time in loop. |
|
902 |
|
903 delete(idx) |
|
904 end do |
|
905 delete(data_ob) |
|
906 delete(data_mod) |
|
907 end do |
|
908 ;----------------------------------------------------------------- |
|
909 ; compute correlation coef and M score |
|
910 |
|
911 u = yvalues(0,:) |
|
912 v = yvalues(1,:) |
|
913 |
|
914 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
915 uu = u(good) |
|
916 vv = v(good) |
|
917 |
|
918 ccrPhase = esccr(uu,vv,0) |
|
919 ; print (ccrPhase) |
|
920 |
|
921 ; old eq |
|
922 ; bias = abs(avg(vv)-avg(uu)) |
|
923 ; new eq |
|
924 bias = avg(abs(vv-uu)) |
|
925 |
|
926 bias = where((bias.gt. 6.),12.-bias,bias) |
|
927 Mphase = ((6. - bias)/6.)*5. |
|
928 |
|
929 M_lai_phase = sprintf("%.2f", Mphase) |
|
930 system("sed s#M_lai_phase#"+M_lai_phase+"# "+html_name+" > "+html_new+";"+ \ |
|
931 "mv -f "+html_new+" "+html_name) |
|
932 print (M_lai_phase) |
|
933 |
|
934 do i=0,nrow-2 |
|
935 text4(i,6) = sprintf("%5.2f",u(i)) |
|
936 text4(i,7) = sprintf("%5.2f",v(i)) |
|
937 text4(i,8) = "-" |
|
938 end do |
|
939 text4(nrow-1,6) = sprintf("%5.2f",avg(u)) |
|
940 text4(nrow-1,7) = sprintf("%5.2f",avg(v)) |
|
941 text4(nrow-1,8) = sprintf("%5.2f",Mphase) |
|
942 |
|
943 delete (u) |
|
944 delete (v) |
|
945 delete (uu) |
|
946 delete (vv) |
|
947 ;-------------------------------------------------------------------- |
|
948 ; histogram res |
|
949 |
|
950 resm = True |
|
951 resm@gsnMaximize = True |
|
952 resm@gsnDraw = False |
|
953 resm@gsnFrame = False |
|
954 resm@xyMarkLineMode = "Markers" |
|
955 resm@xyMarkerSizeF = 0.014 |
|
956 resm@xyMarker = 16 |
|
957 resm@xyMarkerColors = (/"Brown","Blue"/) |
|
958 ; resm@trYMinF = min(mn_yvalues) - 10. |
|
959 ; resm@trYMaxF = max(mx_yvalues) + 10. |
|
960 resm@trYMinF = min(mn_yvalues) - 2 |
|
961 resm@trYMaxF = max(mx_yvalues) + 4 |
|
962 |
|
963 resm@tiYAxisString = "Max LAI (Leaf Area Index) Month" |
|
964 resm@tiXAxisString = "Land Cover Type" |
|
965 |
|
966 max_bar = new((/2,nx/),graphic) |
|
967 min_bar = new((/2,nx/),graphic) |
|
968 max_cap = new((/2,nx/),graphic) |
|
969 min_cap = new((/2,nx/),graphic) |
|
970 |
|
971 lnres = True |
|
972 line_colors = (/"brown","blue"/) |
|
973 ;------------------------------------------------------------------ |
|
974 ; Start the graphics. |
|
975 |
|
976 plot_name = "histogram_phase" |
|
977 title = model_name + " vs Observed" |
|
978 resm@tiMainString = title |
|
979 |
|
980 wks = gsn_open_wks (plot_type,plot_name) |
|
981 ;----------------------------- |
|
982 ; Add a boxed legend using the more simple method |
|
983 |
|
984 resm@pmLegendDisplayMode = "Always" |
|
985 ; resm@pmLegendWidthF = 0.1 |
|
986 resm@pmLegendWidthF = 0.08 |
|
987 resm@pmLegendHeightF = 0.05 |
|
988 resm@pmLegendOrthogonalPosF = -1.17 |
|
989 ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) |
|
990 ; resm@pmLegendParallelPosF = 0.18 |
|
991 resm@pmLegendParallelPosF = 0.88 ;(rightward) |
|
992 |
|
993 ; resm@lgPerimOn = False |
|
994 resm@lgLabelFontHeightF = 0.015 |
|
995 resm@xyExplicitLegendLabels = (/"observed",model_name/) |
|
996 ;----------------------------- |
|
997 tRes = True |
|
998 tRes@txFontHeightF = 0.025 |
|
999 |
|
1000 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")" |
|
1001 |
|
1002 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) |
|
1003 |
|
1004 xy = gsn_csm_xy(wks,xvalues,yvalues,resm) |
|
1005 ;------------------------------- |
|
1006 ;Attach the vertical bar and the horizontal cap line |
|
1007 |
|
1008 do nd=0,1 |
|
1009 lnres@gsLineColor = line_colors(nd) |
|
1010 do i=0,nx-1 |
|
1011 |
|
1012 if(.not.ismissing(mn_yvalues(nd,i)).and. \ |
|
1013 .not.ismissing(mx_yvalues(nd,i))) then |
|
1014 |
|
1015 ; Attach the vertical bar, both above and below the marker. |
|
1016 |
|
1017 x1 = xvalues(nd,i) |
|
1018 y1 = yvalues(nd,i) |
|
1019 y2 = mn_yvalues(nd,i) |
|
1020 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
1021 |
|
1022 y2 = mx_yvalues(nd,i) |
|
1023 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
1024 |
|
1025 ; Attach the horizontal cap line, both above and below the marker. |
|
1026 |
|
1027 x1 = xvalues(nd,i) - dx4 |
|
1028 x2 = xvalues(nd,i) + dx4 |
|
1029 y1 = mn_yvalues(nd,i) |
|
1030 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
1031 |
|
1032 y1 = mx_yvalues(nd,i) |
|
1033 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
1034 end if |
|
1035 end do |
|
1036 end do |
|
1037 |
|
1038 draw(xy) |
|
1039 frame(wks) |
|
1040 |
|
1041 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1042 "rm "+plot_name+"."+plot_type) |
|
1043 |
|
1044 clear (wks) |
|
1045 |
|
1046 delete (DATA11_1D) |
|
1047 delete (DATA12_1D) |
|
1048 delete (DATA22_1D) |
|
1049 ;delete (range) |
|
1050 ;delete (xvalues) |
|
1051 delete (yvalues) |
|
1052 delete (mn_yvalues) |
|
1053 delete (mx_yvalues) |
|
1054 delete (good) |
|
1055 delete (max_bar) |
|
1056 delete (min_bar) |
|
1057 delete (max_cap) |
|
1058 delete (min_cap) |
|
1059 ;----------------------------------------------------------------- |
|
1060 ;global res |
|
1061 |
|
1062 resg = True ; Use plot options |
|
1063 resg@cnFillOn = True ; Turn on color fill |
|
1064 resg@gsnSpreadColors = True ; use full colormap |
|
1065 ; resg@cnFillMode = "RasterFill" ; Turn on raster color |
|
1066 ; resg@lbLabelAutoStride = True |
|
1067 resg@cnLinesOn = False ; Turn off contourn lines |
|
1068 resg@mpFillOn = False ; Turn off map fill |
|
1069 |
|
1070 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
1071 resg@cnMinLevelValF = 1. ; Min level |
|
1072 resg@cnMaxLevelValF = 12. ; Max level |
|
1073 resg@cnLevelSpacingF = 1. ; interval |
|
1074 |
|
1075 ;global contour ob |
|
1076 |
|
1077 delta = 0.00001 |
|
1078 laiob_phase = where(ismissing(laiob_phase).and.(ismissing(laimod_phase).or.(laimod_phase.lt.delta)),0.,laiob_phase) |
|
1079 |
|
1080 plot_name = "global_phase_ob" |
|
1081 title = ob_name |
|
1082 resg@tiMainString = title |
|
1083 |
|
1084 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1085 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1086 |
|
1087 plot = gsn_csm_contour_map_ce(wks,laiob_phase,resg) |
|
1088 frame(wks) |
|
1089 |
|
1090 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1091 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1092 "rm "+plot_name+"-*."+plot_type_new+";"+ \ |
|
1093 "rm "+plot_name+"."+plot_type) |
|
1094 |
|
1095 clear (wks) |
|
1096 ;------------------------------------------------------------------------ |
|
1097 ;global contour model |
|
1098 |
|
1099 plot_name = "global_phase_model" |
|
1100 title = "Model " + model_name |
|
1101 resg@tiMainString = title |
|
1102 |
|
1103 wks = gsn_open_wks (plot_type,plot_name) |
|
1104 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1105 |
|
1106 delete (plot) |
|
1107 plot = gsn_csm_contour_map_ce(wks,laimod_phase,resg) |
|
1108 frame(wks) |
|
1109 |
|
1110 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1111 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1112 "rm "+plot_name+"-*."+plot_type_new+";"+ \ |
|
1113 "rm "+plot_name+"."+plot_type) |
|
1114 |
|
1115 clear (wks) |
|
1116 ;----------------------------------------------------------------- |
|
1117 ;(D) grow day |
|
1118 ;-------------------------------------------------------------------- |
|
1119 ; get data |
|
1120 |
|
1121 day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/) |
|
1122 |
|
1123 ; observed |
|
1124 laiob_grow = laiob(0,:,:) |
|
1125 laiob_grow@long_name = "Days of Growing Season" |
|
1126 |
|
1127 dsizes_z = dimsizes(laiob) |
|
1128 ntime = dsizes_z(0) |
|
1129 nlat = dsizes_z(1) |
|
1130 nlon = dsizes_z(2) |
|
1131 |
|
1132 do j = 0,nlat-1 |
|
1133 do i = 0,nlon-1 |
|
1134 nday = 0. |
|
1135 do k = 0,ntime-1 |
|
1136 if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then |
|
1137 nday = nday + day_of_data(k) |
|
1138 end if |
|
1139 end do |
|
1140 |
|
1141 laiob_grow(j,i) = nday |
|
1142 end do |
|
1143 end do |
|
1144 |
|
1145 ; print (min(laiob_grow)+"/"+max(laiob_grow)) |
|
1146 ;------------------------- |
|
1147 ; model |
|
1148 laimod_grow = laimod(0,:,:) |
|
1149 laimod_grow@long_name = "Days of Growing Season" |
|
1150 |
|
1151 dsizes_z = dimsizes(laimod) |
|
1152 ntime = dsizes_z(0) |
|
1153 nlat = dsizes_z(1) |
|
1154 nlon = dsizes_z(2) |
|
1155 |
|
1156 do j = 0,nlat-1 |
|
1157 do i = 0,nlon-1 |
|
1158 nday = 0. |
|
1159 do k = 0,ntime-1 |
|
1160 if (.not. ismissing(laimod(k,j,i)) .and. laimod(k,j,i) .gt. 1.0) then |
|
1161 nday = nday + day_of_data(k) |
|
1162 end if |
|
1163 end do |
|
1164 |
|
1165 laimod_grow(j,i) = nday |
|
1166 end do |
|
1167 end do |
|
1168 |
|
1169 ; print (min(laimod_grow)+"/"+max(laimod_grow)) |
|
1170 ;------------------------ |
|
1171 DATA11_1D = ndtooned(classob) |
|
1172 DATA12_1D = ndtooned(laiob_grow) |
|
1173 DATA22_1D = ndtooned(laimod_grow) |
|
1174 |
|
1175 yvalues = new((/2,nx/),float) |
|
1176 mn_yvalues = new((/2,nx/),float) |
|
1177 mx_yvalues = new((/2,nx/),float) |
|
1178 |
|
1179 do nd=0,1 |
|
1180 |
|
1181 ; See if we are doing model or observational data. |
|
1182 |
|
1183 if(nd.eq.0) then |
|
1184 data_ob = DATA11_1D |
|
1185 data_mod = DATA12_1D |
|
1186 else |
|
1187 data_ob = DATA11_1D |
|
1188 data_mod = DATA22_1D |
|
1189 end if |
|
1190 |
|
1191 ; Loop through each range and check for values. |
|
1192 |
|
1193 do i=0,nr-2 |
|
1194 if (i.ne.(nr-2)) then |
|
1195 ; print("") |
|
1196 ; print("In range ["+range(i)+","+range(i+1)+")") |
|
1197 idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) |
|
1198 else |
|
1199 ; print("") |
|
1200 ; print("In range ["+range(i)+",)") |
|
1201 idx = ind(range(i).le.data_ob) |
|
1202 end if |
|
1203 |
|
1204 ; Calculate average, and get min and max. |
|
1205 |
|
1206 if(.not.any(ismissing(idx))) then |
|
1207 yvalues(nd,i) = avg(data_mod(idx)) |
|
1208 mn_yvalues(nd,i) = min(data_mod(idx)) |
|
1209 mx_yvalues(nd,i) = max(data_mod(idx)) |
|
1210 count = dimsizes(idx) |
|
1211 else |
|
1212 count = 0 |
|
1213 yvalues(nd,i) = yvalues@_FillValue |
|
1214 mn_yvalues(nd,i) = yvalues@_FillValue |
|
1215 mx_yvalues(nd,i) = yvalues@_FillValue |
|
1216 end if |
|
1217 |
|
1218 ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) |
|
1219 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) |
|
1220 |
|
1221 ; Clean up for next time in loop. |
|
1222 |
|
1223 delete(idx) |
|
1224 end do |
|
1225 delete(data_ob) |
|
1226 delete(data_mod) |
|
1227 end do |
|
1228 ;----------------------------------------------------------------- |
|
1229 ; compute correlation coef and M score |
|
1230 |
|
1231 u = yvalues(0,:) |
|
1232 v = yvalues(1,:) |
|
1233 |
|
1234 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
1235 uu = u(good) |
|
1236 vv = v(good) |
|
1237 |
|
1238 ccrGrow = esccr(uu,vv,0) |
|
1239 ; print (ccrGrow) |
|
1240 |
|
1241 ; new eq |
|
1242 bias = sum(abs(vv-uu)/(vv+uu)) |
|
1243 Mgrow = (1.- (bias/dimsizes(uu)))*5. |
|
1244 |
|
1245 M_lai_grow = sprintf("%.2f", Mgrow) |
|
1246 system("sed s#M_lai_grow#"+M_lai_grow+"# "+html_name+" > "+html_new+";"+ \ |
|
1247 "mv -f "+html_new+" "+html_name) |
|
1248 print (M_lai_grow) |
|
1249 |
|
1250 do i=0,nrow-2 |
|
1251 text4(i,9) = sprintf("%5.2f",u(i)) |
|
1252 text4(i,10) = sprintf("%5.2f",v(i)) |
|
1253 text4(i,11) = "-" |
|
1254 end do |
|
1255 text4(nrow-1,9) = sprintf("%5.2f",avg(u)) |
|
1256 text4(nrow-1,10) = sprintf("%5.2f",avg(v)) |
|
1257 text4(nrow-1,11) = sprintf("%5.2f",Mgrow) |
|
1258 |
|
1259 delete (u) |
|
1260 delete (v) |
|
1261 delete (uu) |
|
1262 delete (vv) |
|
1263 ;-------------------------------------------------------------------- |
|
1264 ; histogram res |
|
1265 |
|
1266 resm = True |
|
1267 resm@gsnMaximize = True |
|
1268 resm@gsnDraw = False |
|
1269 resm@gsnFrame = False |
|
1270 resm@xyMarkLineMode = "Markers" |
|
1271 resm@xyMarkerSizeF = 0.014 |
|
1272 resm@xyMarker = 16 |
|
1273 resm@xyMarkerColors = (/"Brown","Blue"/) |
|
1274 ; resm@trYMinF = min(mn_yvalues) - 10. |
|
1275 ; resm@trYMaxF = max(mx_yvalues) + 10. |
|
1276 resm@trYMinF = min(mn_yvalues) - 2 |
|
1277 resm@trYMaxF = max(mx_yvalues) + 4 |
|
1278 |
|
1279 resm@tiYAxisString = "Days of Growing season" |
|
1280 resm@tiXAxisString = "Land Cover Type" |
|
1281 |
|
1282 max_bar = new((/2,nx/),graphic) |
|
1283 min_bar = new((/2,nx/),graphic) |
|
1284 max_cap = new((/2,nx/),graphic) |
|
1285 min_cap = new((/2,nx/),graphic) |
|
1286 |
|
1287 lnres = True |
|
1288 line_colors = (/"brown","blue"/) |
|
1289 ;------------------------------------------------------------------ |
|
1290 ; Start the graphics. |
|
1291 |
|
1292 plot_name = "histogram_grow" |
|
1293 title = model_name + " vs Observed" |
|
1294 resm@tiMainString = title |
|
1295 |
|
1296 wks = gsn_open_wks (plot_type,plot_name) |
|
1297 ;----------------------------- |
|
1298 ; Add a boxed legend using the more simple method |
|
1299 |
|
1300 resm@pmLegendDisplayMode = "Always" |
|
1301 ; resm@pmLegendWidthF = 0.1 |
|
1302 resm@pmLegendWidthF = 0.08 |
|
1303 resm@pmLegendHeightF = 0.05 |
|
1304 resm@pmLegendOrthogonalPosF = -1.17 |
|
1305 ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) |
|
1306 ; resm@pmLegendParallelPosF = 0.18 |
|
1307 resm@pmLegendParallelPosF = 0.88 ;(rightward) |
|
1308 |
|
1309 ; resm@lgPerimOn = False |
|
1310 resm@lgLabelFontHeightF = 0.015 |
|
1311 resm@xyExplicitLegendLabels = (/"observed",model_name/) |
|
1312 ;----------------------------- |
|
1313 tRes = True |
|
1314 tRes@txFontHeightF = 0.025 |
|
1315 |
|
1316 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")" |
|
1317 |
|
1318 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) |
|
1319 |
|
1320 xy = gsn_csm_xy(wks,xvalues,yvalues,resm) |
|
1321 ;------------------------------- |
|
1322 ;Attach the vertical bar and the horizontal cap line |
|
1323 |
|
1324 do nd=0,1 |
|
1325 lnres@gsLineColor = line_colors(nd) |
|
1326 do i=0,nx-1 |
|
1327 |
|
1328 if(.not.ismissing(mn_yvalues(nd,i)).and. \ |
|
1329 .not.ismissing(mx_yvalues(nd,i))) then |
|
1330 |
|
1331 ; Attach the vertical bar, both above and below the marker. |
|
1332 |
|
1333 x1 = xvalues(nd,i) |
|
1334 y1 = yvalues(nd,i) |
|
1335 y2 = mn_yvalues(nd,i) |
|
1336 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
1337 |
|
1338 y2 = mx_yvalues(nd,i) |
|
1339 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
1340 |
|
1341 ; Attach the horizontal cap line, both above and below the marker. |
|
1342 |
|
1343 x1 = xvalues(nd,i) - dx4 |
|
1344 x2 = xvalues(nd,i) + dx4 |
|
1345 y1 = mn_yvalues(nd,i) |
|
1346 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
1347 |
|
1348 y1 = mx_yvalues(nd,i) |
|
1349 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
1350 end if |
|
1351 end do |
|
1352 end do |
|
1353 |
|
1354 draw(xy) |
|
1355 frame(wks) |
|
1356 |
|
1357 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1358 "rm "+plot_name+"."+plot_type) |
|
1359 |
|
1360 clear (wks) |
|
1361 |
|
1362 delete (DATA11_1D) |
|
1363 delete (DATA12_1D) |
|
1364 delete (DATA22_1D) |
|
1365 ;delete (range) |
|
1366 ;delete (xvalues) |
|
1367 delete (yvalues) |
|
1368 delete (mn_yvalues) |
|
1369 delete (mx_yvalues) |
|
1370 delete (good) |
|
1371 delete (max_bar) |
|
1372 delete (min_bar) |
|
1373 delete (max_cap) |
|
1374 delete (min_cap) |
|
1375 ;----------------------------------------------------------------- |
|
1376 ;global res |
|
1377 |
|
1378 resg = True ; Use plot options |
|
1379 resg@cnFillOn = True ; Turn on color fill |
|
1380 resg@gsnSpreadColors = True ; use full colormap |
|
1381 ; resg@cnFillMode = "RasterFill" ; Turn on raster color |
|
1382 ; resg@lbLabelAutoStride = True |
|
1383 resg@cnLinesOn = False ; Turn off contourn lines |
|
1384 resg@mpFillOn = False ; Turn off map fill |
|
1385 |
|
1386 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
1387 resg@cnMinLevelValF = 60. ; Min level |
|
1388 resg@cnMaxLevelValF = 360. ; Max level |
|
1389 resg@cnLevelSpacingF = 20. ; interval |
|
1390 |
|
1391 ;global contour ob |
|
1392 |
|
1393 laiob_grow@_FillValue = 1.e+36 |
|
1394 laiob_grow = where(laiob_grow .lt. 10.,laiob_grow@_FillValue,laiob_grow) |
|
1395 |
|
1396 plot_name = "global_grow_ob" |
|
1397 title = ob_name |
|
1398 resg@tiMainString = title |
|
1399 |
|
1400 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1401 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1402 |
|
1403 plot = gsn_csm_contour_map_ce(wks,laiob_grow,resg) |
|
1404 frame(wks) |
|
1405 |
|
1406 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1407 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1408 "rm "+plot_name+"-*."+plot_type_new+";"+ \ |
|
1409 "rm "+plot_name+"."+plot_type) |
|
1410 |
|
1411 clear (wks) |
|
1412 ;------------------------------------------------------------------------ |
|
1413 ;global contour model |
|
1414 |
|
1415 laimod_grow@_FillValue = 1.e+36 |
|
1416 laimod_grow = where(laimod_grow .lt. 10.,laimod_grow@_FillValue,laimod_grow) |
|
1417 |
|
1418 plot_name = "global_grow_model" |
|
1419 title = "Model " + model_name |
|
1420 resg@tiMainString = title |
|
1421 |
|
1422 wks = gsn_open_wks (plot_type,plot_name) |
|
1423 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1424 |
|
1425 delete (plot) |
|
1426 plot = gsn_csm_contour_map_ce(wks,laimod_grow,resg) |
|
1427 frame(wks) |
|
1428 |
|
1429 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1430 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1431 "rm "+plot_name+"-*."+plot_type_new+";"+ \ |
|
1432 "rm "+plot_name+"."+plot_type) |
|
1433 |
|
1434 clear (wks) |
|
1435 ;------------------------------------------------------------------------ |
|
1436 ;global contour model vs ob |
|
1437 |
|
1438 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
1439 resg@cnMinLevelValF = 0. ; Min level |
|
1440 resg@cnMaxLevelValF = 10. ; Max level |
|
1441 resg@cnLevelSpacingF = 1. ; interval |
|
1442 |
|
1443 plot_name = "global_mean_model_vs_ob" |
|
1444 |
|
1445 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1446 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1447 |
|
1448 delete (plot) |
|
1449 plot=new(3,graphic) ; create graphic array |
|
1450 |
|
1451 resg@gsnFrame = False ; Do not draw plot |
|
1452 resg@gsnDraw = False ; Do not advance frame |
|
1453 |
|
1454 ; plot correlation coef |
|
1455 |
|
1456 gRes = True |
|
1457 gRes@txFontHeightF = 0.02 |
|
1458 gRes@txAngleF = 90 |
|
1459 |
|
1460 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")" |
|
1461 |
|
1462 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
|
1463 ;-------------------------------------------------------------------- |
|
1464 |
|
1465 ;(a) ob |
|
1466 |
|
1467 title = ob_name |
|
1468 resg@tiMainString = title |
|
1469 |
|
1470 plot(0) = gsn_csm_contour_map_ce(wks,laiob_mean,resg) |
|
1471 |
|
1472 ;(b) model |
|
1473 |
|
1474 title = "Model "+ model_name |
|
1475 resg@tiMainString = title |
|
1476 |
|
1477 plot(1) = gsn_csm_contour_map_ce(wks,laimod_mean,resg) |
|
1478 |
|
1479 ;(c) model-ob |
|
1480 |
|
1481 zz = laimod_mean |
|
1482 zz = laimod_mean - laiob_mean |
|
1483 title = "Model_"+model_name+" - Observed" |
|
1484 resg@tiMainString = title |
|
1485 |
|
1486 resg@cnMinLevelValF = -2. ; Min level |
|
1487 resg@cnMaxLevelValF = 2. ; Max level |
|
1488 resg@cnLevelSpacingF = 0.4 ; interval |
|
1489 |
|
1490 |
|
1491 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) |
|
1492 |
|
1493 pres = True ; panel plot mods desired |
|
1494 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around |
|
1495 ; indiv. plots in panel |
|
1496 pres@gsnMaximize = True ; fill the page |
|
1497 |
|
1498 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
|
1499 |
|
1500 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1501 "rm "+plot_name+"."+plot_type) |
|
1502 |
|
1503 clear (wks) |
|
1504 delete (plot) |
|
1505 ;----------------------------------------------------------------- |
|
1506 ;global contour model vs ob |
|
1507 |
|
1508 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
1509 resg@cnMinLevelValF = 0. ; Min level |
|
1510 resg@cnMaxLevelValF = 10. ; Max level |
|
1511 resg@cnLevelSpacingF = 1. ; interval |
|
1512 |
|
1513 plot_name = "global_max_model_vs_ob" |
|
1514 |
|
1515 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1516 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1517 |
|
1518 ; delete(plot) |
|
1519 plot=new(3,graphic) ; create graphic array |
|
1520 |
|
1521 resg@gsnFrame = False ; Do not draw plot |
|
1522 resg@gsnDraw = False ; Do not advance frame |
|
1523 |
|
1524 ; plot correlation coef |
|
1525 |
|
1526 gRes = True |
|
1527 gRes@txFontHeightF = 0.02 |
|
1528 gRes@txAngleF = 90 |
|
1529 |
|
1530 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")" |
|
1531 |
|
1532 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
|
1533 ;-------------------------------------------------------------------- |
|
1534 ;(a) ob |
|
1535 |
|
1536 title = ob_name |
|
1537 resg@tiMainString = title |
|
1538 |
|
1539 plot(0) = gsn_csm_contour_map_ce(wks,laiob_max,resg) |
|
1540 |
|
1541 ;(b) model |
|
1542 |
|
1543 title = "Model "+ model_name |
|
1544 resg@tiMainString = title |
|
1545 |
|
1546 plot(1) = gsn_csm_contour_map_ce(wks,laimod_max,resg) |
|
1547 |
|
1548 ;(c) model-ob |
|
1549 |
|
1550 delete (zz) |
|
1551 zz = laimod_max |
|
1552 zz = laimod_max - laiob_max |
|
1553 title = "Model_"+model_name+" - Observed" |
|
1554 resg@tiMainString = title |
|
1555 |
|
1556 resg@cnMinLevelValF = -6. ; Min level |
|
1557 resg@cnMaxLevelValF = 6. ; Max level |
|
1558 resg@cnLevelSpacingF = 1. ; interval |
|
1559 |
|
1560 |
|
1561 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) |
|
1562 |
|
1563 pres = True ; panel plot mods desired |
|
1564 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around |
|
1565 ; indiv. plots in panel |
|
1566 pres@gsnMaximize = True ; fill the page |
|
1567 |
|
1568 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
|
1569 |
|
1570 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1571 "rm "+plot_name+"."+plot_type) |
|
1572 |
|
1573 clear (wks) |
|
1574 delete (plot) |
|
1575 ;----------------------------------------------------------------- |
|
1576 ;global contour model vs ob |
|
1577 |
|
1578 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
1579 resg@cnMinLevelValF = 1. ; Min level |
|
1580 resg@cnMaxLevelValF = 12. ; Max level |
|
1581 resg@cnLevelSpacingF = 1. ; interval |
|
1582 |
|
1583 plot_name = "global_phase_model_vs_ob" |
|
1584 |
|
1585 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1586 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1587 |
|
1588 ; delete (plot) |
|
1589 plot=new(3,graphic) ; create graphic array |
|
1590 |
|
1591 resg@gsnFrame = False ; Do not draw plot |
|
1592 resg@gsnDraw = False ; Do not advance frame |
|
1593 |
|
1594 ; plot correlation coef |
|
1595 |
|
1596 gRes = True |
|
1597 gRes@txFontHeightF = 0.02 |
|
1598 gRes@txAngleF = 90 |
|
1599 |
|
1600 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")" |
|
1601 |
|
1602 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
|
1603 ;-------------------------------------------------------------------- |
|
1604 ;(a) ob |
|
1605 |
|
1606 title = ob_name |
|
1607 resg@tiMainString = title |
|
1608 |
|
1609 plot(0) = gsn_csm_contour_map_ce(wks,laiob_phase,resg) |
|
1610 |
|
1611 ;(b) model |
|
1612 |
|
1613 title = "Model "+ model_name |
|
1614 resg@tiMainString = title |
|
1615 |
|
1616 plot(1) = gsn_csm_contour_map_ce(wks,laimod_phase,resg) |
|
1617 |
|
1618 ;(c) model-ob |
|
1619 |
|
1620 delete (zz) |
|
1621 zz = laimod_phase |
|
1622 zz = laimod_phase - laiob_phase |
|
1623 title = "Model_"+model_name+" - Observed" |
|
1624 resg@tiMainString = title |
|
1625 |
|
1626 resg@cnMinLevelValF = -6. ; Min level |
|
1627 resg@cnMaxLevelValF = 6. ; Max level |
|
1628 resg@cnLevelSpacingF = 1. ; interval |
|
1629 |
|
1630 |
|
1631 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) |
|
1632 |
|
1633 pres = True ; panel plot mods desired |
|
1634 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around |
|
1635 ; indiv. plots in panel |
|
1636 pres@gsnMaximize = True ; fill the page |
|
1637 |
|
1638 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
|
1639 |
|
1640 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1641 "rm "+plot_name+"."+plot_type) |
|
1642 |
|
1643 clear (wks) |
|
1644 delete (plot) |
|
1645 ;------------------------------------------------------------------ |
|
1646 ;global contour model vs ob |
|
1647 |
|
1648 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
1649 resg@cnMinLevelValF = 60. ; Min level |
|
1650 resg@cnMaxLevelValF = 360. ; Max level |
|
1651 resg@cnLevelSpacingF = 20. ; interval |
|
1652 |
|
1653 plot_name = "global_grow_model_vs_ob" |
|
1654 |
|
1655 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1656 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1657 |
|
1658 ; delete (plot) |
|
1659 plot=new(3,graphic) ; create graphic array |
|
1660 |
|
1661 resg@gsnFrame = False ; Do not draw plot |
|
1662 resg@gsnDraw = False ; Do not advance frame |
|
1663 |
|
1664 ; plot correlation coef |
|
1665 |
|
1666 gRes = True |
|
1667 gRes@txFontHeightF = 0.02 |
|
1668 gRes@txAngleF = 90 |
|
1669 |
|
1670 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")" |
|
1671 |
|
1672 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
|
1673 ;-------------------------------------------------------------------- |
|
1674 ;(a) ob |
|
1675 |
|
1676 title = ob_name |
|
1677 resg@tiMainString = title |
|
1678 |
|
1679 plot(0) = gsn_csm_contour_map_ce(wks,laiob_grow,resg) |
|
1680 |
|
1681 ;(b) model |
|
1682 |
|
1683 title = "Model "+ model_name |
|
1684 resg@tiMainString = title |
|
1685 |
|
1686 plot(1) = gsn_csm_contour_map_ce(wks,laimod_grow,resg) |
|
1687 |
|
1688 ;(c) model-ob |
|
1689 |
|
1690 delete (zz) |
|
1691 zz = laimod_grow |
|
1692 zz = laimod_grow - laiob_grow |
|
1693 title = "Model_"+model_name+" - Observed" |
|
1694 resg@tiMainString = title |
|
1695 |
|
1696 resg@cnMinLevelValF = -100. ; Min level |
|
1697 resg@cnMaxLevelValF = 100. ; Max level |
|
1698 resg@cnLevelSpacingF = 20. ; interval |
|
1699 |
|
1700 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) |
|
1701 |
|
1702 pres = True ; panel plot mods desired |
|
1703 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around |
|
1704 ; indiv. plots in panel |
|
1705 pres@gsnMaximize = True ; fill the page |
|
1706 |
|
1707 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
|
1708 |
|
1709 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1710 "rm "+plot_name+"."+plot_type) |
|
1711 |
|
1712 clear (wks) |
|
1713 delete (plot) |
|
1714 ;************************************************** |
|
1715 ; plot lai table |
|
1716 ;************************************************** |
|
1717 |
|
1718 plot_name = "table_lai" |
|
1719 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1720 |
|
1721 gsn_table(wks,ncr1,xx1,yy1,text1,res1) |
|
1722 gsn_table(wks,ncr21,xx21,yy21,text21,res21) |
|
1723 gsn_table(wks,ncr22,xx22,yy22,text22,res22) |
|
1724 gsn_table(wks,ncr3,xx3,yy3,text3,res3) |
|
1725 gsn_table(wks,ncr4,xx4,yy4,text4,res4) |
|
1726 |
|
1727 frame(wks) |
|
1728 |
|
1729 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1730 "rm "+plot_name+"."+plot_type) |
|
1731 ;------------------------------------------------------------------- |
|
1732 temp_name = "lai." + model_name |
|
1733 system("mkdir -p " + temp_name+";"+ \ |
|
1734 "cp "+ html_name + " " +temp_name+"/table.html"+";"+ \ |
|
1735 "mv *.png " + temp_name +";"+ \ |
|
1736 "tar cf "+ temp_name +".tar " + temp_name) |
|
1737 ;------------------------------------------------------------------- |
|
1738 end |
|
1739 |