|
1 ;******************************************************** |
|
2 ; histogram normalized by rain and compute correleration |
|
3 ;******************************************************** |
|
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
|
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
|
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 = 18 |
|
20 nclass = 20 |
|
21 |
|
22 plot_type = "ps" |
|
23 plot_type_new = "png" |
|
24 |
|
25 ;************************************************ |
|
26 ; read in data: model |
|
27 ;************************************************ |
|
28 |
|
29 model_name = "b30.061n" |
|
30 model_grid = "T31" |
|
31 |
|
32 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" |
|
33 film = "b30.061n_1995-2004_MONS_climo_lnd.nc" |
|
34 ;film = "i01.03cn_1545-1569_MONS_climo.nc" |
|
35 fm = addfile(dirm+film,"r") |
|
36 |
|
37 laimod = fm->TLAI |
|
38 |
|
39 ;************************************************ |
|
40 ; read in data: observed |
|
41 ;************************************************ |
|
42 |
|
43 ob_name = "MODIS MOD 15A2 2000-2005" |
|
44 |
|
45 diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/" |
|
46 filo1 = "land_class_"+model_grid+".nc" |
|
47 filo2 = "LAI_2000-2005_MONS_"+model_grid+".nc" |
|
48 |
|
49 fo1 = addfile(diro+filo1,"r") |
|
50 fo2 = addfile(diro+filo2,"r") |
|
51 |
|
52 classob = tofloat(fo1->LAND_CLASS) |
|
53 laiob = fo2->LAI |
|
54 ;******************************************************************* |
|
55 ; Calculate "nice" bins for binning the data in equally spaced ranges |
|
56 ;******************************************************************** |
|
57 nclassn = nclass + 1 |
|
58 range = fspan(0,nclassn-1,nclassn) |
|
59 ; print (range) |
|
60 |
|
61 ; Use this range information to grab all the values in a |
|
62 ; particular range, and then take an average. |
|
63 |
|
64 nr = dimsizes(range) |
|
65 nx = nr-1 |
|
66 xvalues = new((/2,nx/),float) |
|
67 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2. |
|
68 dx = xvalues(0,1) - xvalues(0,0) ; range width |
|
69 dx4 = dx/4 ; 1/4 of the range |
|
70 xvalues(1,:) = xvalues(0,:) - dx/5. |
|
71 ;----------------------------------------------------------------- |
|
72 |
|
73 ;----------------------------------------------------------------- |
|
74 ;(B) max |
|
75 ;-------------------------------------------------------------------- |
|
76 ; get data |
|
77 |
|
78 ; observed |
|
79 laiob_max = laiob(0,:,:) |
|
80 s = laiob(:,0,0) |
|
81 laiob_max@long_name = "Leaf Area Index Max" |
|
82 |
|
83 dsizes_z = dimsizes(laiob) |
|
84 nlat = dsizes_z(1) |
|
85 nlon = dsizes_z(2) |
|
86 |
|
87 do j = 0,nlat-1 |
|
88 do i = 0,nlon-1 |
|
89 s = laiob(:,j,i) |
|
90 laiob_max(j,i) = max(s) |
|
91 end do |
|
92 end do |
|
93 |
|
94 ; print (min(y)+"/"+max(y)) |
|
95 delete (s) |
|
96 delete (dsizes_z) |
|
97 ;------------------------- |
|
98 ; model |
|
99 laimod_max = laimod(0,:,:) |
|
100 s = laimod(:,0,0) |
|
101 laimod_max@long_name = "Leaf Area Index Max" |
|
102 |
|
103 dsizes_z = dimsizes(laimod) |
|
104 nlat = dsizes_z(1) |
|
105 nlon = dsizes_z(2) |
|
106 |
|
107 do j = 0,nlat-1 |
|
108 do i = 0,nlon-1 |
|
109 s = laimod(:,j,i) |
|
110 laimod_max(j,i) = max(s) |
|
111 end do |
|
112 end do |
|
113 |
|
114 ; print (min(laimod_max)+"/"+max(laimod_max)) |
|
115 delete (s) |
|
116 delete (dsizes_z) |
|
117 ;------------------------ |
|
118 DATA11_1D = ndtooned(classob) |
|
119 DATA12_1D = ndtooned(laiob_max) |
|
120 DATA22_1D = ndtooned(laimod_max) |
|
121 |
|
122 yvalues = new((/2,nx/),float) |
|
123 mn_yvalues = new((/2,nx/),float) |
|
124 mx_yvalues = new((/2,nx/),float) |
|
125 |
|
126 do nd=0,1 |
|
127 |
|
128 ; See if we are doing model or observational data. |
|
129 |
|
130 if(nd.eq.0) then |
|
131 data_ob = DATA11_1D |
|
132 data_mod = DATA12_1D |
|
133 else |
|
134 data_ob = DATA11_1D |
|
135 data_mod = DATA22_1D |
|
136 end if |
|
137 |
|
138 ; Loop through each range and check for values. |
|
139 |
|
140 do i=0,nr-2 |
|
141 if (i.ne.(nr-2)) then |
|
142 ; print("") |
|
143 ; print("In range ["+range(i)+","+range(i+1)+")") |
|
144 idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) |
|
145 else |
|
146 ; print("") |
|
147 ; print("In range ["+range(i)+",)") |
|
148 idx = ind(range(i).le.data_ob) |
|
149 end if |
|
150 |
|
151 ; Calculate average, and get min and max. |
|
152 |
|
153 if(.not.any(ismissing(idx))) then |
|
154 yvalues(nd,i) = avg(data_mod(idx)) |
|
155 mn_yvalues(nd,i) = min(data_mod(idx)) |
|
156 mx_yvalues(nd,i) = max(data_mod(idx)) |
|
157 count = dimsizes(idx) |
|
158 else |
|
159 count = 0 |
|
160 yvalues(nd,i) = yvalues@_FillValue |
|
161 mn_yvalues(nd,i) = yvalues@_FillValue |
|
162 mx_yvalues(nd,i) = yvalues@_FillValue |
|
163 end if |
|
164 |
|
165 ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) |
|
166 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) |
|
167 |
|
168 ; Clean up for next time in loop. |
|
169 |
|
170 delete(idx) |
|
171 end do |
|
172 delete(data_ob) |
|
173 delete(data_mod) |
|
174 end do |
|
175 ;----------------------------------------------------------------- |
|
176 ; compute correlation coef and M score |
|
177 |
|
178 u = yvalues(0,:) |
|
179 v = yvalues(1,:) |
|
180 |
|
181 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
182 uu = u(good) |
|
183 vv = v(good) |
|
184 |
|
185 ccrMax = esccr(uu,vv,0) |
|
186 ; print (ccrMax) |
|
187 |
|
188 ; new eq |
|
189 bias = sum(abs(vv-uu)/(vv+uu)) |
|
190 Mmax = (1.- (bias/dimsizes(uu)))*5. |
|
191 |
|
192 print (Mmax) |
|
193 |
|
194 delete (u) |
|
195 delete (v) |
|
196 delete (uu) |
|
197 delete (vv) |
|
198 ;-------------------------------------------------------------------- |
|
199 ; histogram res |
|
200 |
|
201 resm = True |
|
202 resm@gsnMaximize = True |
|
203 resm@gsnDraw = False |
|
204 resm@gsnFrame = False |
|
205 resm@xyMarkLineMode = "Markers" |
|
206 resm@xyMarkerSizeF = 0.014 |
|
207 resm@xyMarker = 16 |
|
208 resm@xyMarkerColors = (/"Brown","Blue"/) |
|
209 ; resm@trYMinF = min(mn_yvalues) - 10. |
|
210 ; resm@trYMaxF = max(mx_yvalues) + 10. |
|
211 resm@trYMinF = min(mn_yvalues) - 2 |
|
212 resm@trYMaxF = max(mx_yvalues) + 4 |
|
213 |
|
214 resm@tiYAxisString = "Max LAI (Leaf Area Index)" |
|
215 resm@tiXAxisString = "Land Cover Type" |
|
216 |
|
217 max_bar = new((/2,nx/),graphic) |
|
218 min_bar = new((/2,nx/),graphic) |
|
219 max_cap = new((/2,nx/),graphic) |
|
220 min_cap = new((/2,nx/),graphic) |
|
221 |
|
222 lnres = True |
|
223 line_colors = (/"brown","blue"/) |
|
224 ;------------------------------------------------------------------ |
|
225 ; Start the graphics. |
|
226 |
|
227 plot_name = "histogram_max" |
|
228 title = model_name + " vs Observed" |
|
229 resm@tiMainString = title |
|
230 |
|
231 wks = gsn_open_wks (plot_type,plot_name) |
|
232 ;----------------------------- |
|
233 ; Add a boxed legend using the more simple method |
|
234 |
|
235 resm@pmLegendDisplayMode = "Always" |
|
236 ; resm@pmLegendWidthF = 0.1 |
|
237 resm@pmLegendWidthF = 0.08 |
|
238 resm@pmLegendHeightF = 0.05 |
|
239 resm@pmLegendOrthogonalPosF = -1.17 |
|
240 ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) |
|
241 ; resm@pmLegendParallelPosF = 0.18 |
|
242 resm@pmLegendParallelPosF = 0.88 ;(rightward) |
|
243 |
|
244 ; resm@lgPerimOn = False |
|
245 resm@lgLabelFontHeightF = 0.015 |
|
246 resm@xyExplicitLegendLabels = (/"observed",model_name/) |
|
247 ;----------------------------- |
|
248 tRes = True |
|
249 tRes@txFontHeightF = 0.025 |
|
250 |
|
251 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")" |
|
252 |
|
253 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) |
|
254 |
|
255 xy = gsn_csm_xy(wks,xvalues,yvalues,resm) |
|
256 ;------------------------------- |
|
257 ;Attach the vertical bar and the horizontal cap line |
|
258 |
|
259 do nd=0,1 |
|
260 lnres@gsLineColor = line_colors(nd) |
|
261 do i=0,nx-1 |
|
262 |
|
263 if(.not.ismissing(mn_yvalues(nd,i)).and. \ |
|
264 .not.ismissing(mx_yvalues(nd,i))) then |
|
265 |
|
266 ; Attach the vertical bar, both above and below the marker. |
|
267 |
|
268 x1 = xvalues(nd,i) |
|
269 y1 = yvalues(nd,i) |
|
270 y2 = mn_yvalues(nd,i) |
|
271 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
272 |
|
273 y2 = mx_yvalues(nd,i) |
|
274 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
275 |
|
276 ; Attach the horizontal cap line, both above and below the marker. |
|
277 |
|
278 x1 = xvalues(nd,i) - dx4 |
|
279 x2 = xvalues(nd,i) + dx4 |
|
280 y1 = mn_yvalues(nd,i) |
|
281 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
282 |
|
283 y1 = mx_yvalues(nd,i) |
|
284 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
285 end if |
|
286 end do |
|
287 end do |
|
288 |
|
289 draw(xy) |
|
290 frame(wks) |
|
291 |
|
292 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
|
293 ; system("rm "+plot_name+"."+plot_type) |
|
294 ; system("rm "+plot_name+"-1."+plot_type_new) |
|
295 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) |
|
296 |
|
297 clear (wks) |
|
298 |
|
299 delete (DATA11_1D) |
|
300 delete (DATA12_1D) |
|
301 delete (DATA22_1D) |
|
302 ;delete (range) |
|
303 ;delete (xvalues) |
|
304 delete (yvalues) |
|
305 delete (mn_yvalues) |
|
306 delete (mx_yvalues) |
|
307 delete (good) |
|
308 delete (max_bar) |
|
309 delete (min_bar) |
|
310 delete (max_cap) |
|
311 delete (min_cap) |
|
312 ;----------------------------------------------------------------- |
|
313 ;global res |
|
314 |
|
315 resg = True ; Use plot options |
|
316 resg@cnFillOn = True ; Turn on color fill |
|
317 resg@gsnSpreadColors = True ; use full colormap |
|
318 ; resg@cnFillMode = "RasterFill" ; Turn on raster color |
|
319 ; resg@lbLabelAutoStride = True |
|
320 resg@cnLinesOn = False ; Turn off contourn lines |
|
321 resg@mpFillOn = False ; Turn off map fill |
|
322 |
|
323 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
324 resg@cnMinLevelValF = 0. ; Min level |
|
325 resg@cnMaxLevelValF = 10. ; Max level |
|
326 resg@cnLevelSpacingF = 1. ; interval |
|
327 |
|
328 ;global contour ob |
|
329 |
|
330 delta = 0.00001 |
|
331 laiob_max = where(ismissing(laiob_max).and.(ismissing(laimod_max).or.(laimod_max.lt.delta)),0.,laiob_max) |
|
332 |
|
333 plot_name = "global_max_ob" |
|
334 title = ob_name |
|
335 resg@tiMainString = title |
|
336 |
|
337 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
338 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
339 |
|
340 plot = gsn_csm_contour_map_ce(wks,laiob_max,resg) |
|
341 frame(wks) |
|
342 |
|
343 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
|
344 ; system("rm "+plot_name+"."+plot_type) |
|
345 ; system("rm "+plot_name+"-1."+plot_type_new) |
|
346 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) |
|
347 |
|
348 clear (wks) |
|
349 ;------------------------------------------------------------------------ |
|
350 ;global contour model |
|
351 |
|
352 plot_name = "global_max_model" |
|
353 title = "Model " + model_name |
|
354 resg@tiMainString = title |
|
355 |
|
356 wks = gsn_open_wks (plot_type,plot_name) |
|
357 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
358 |
|
359 delete (plot) |
|
360 plot = gsn_csm_contour_map_ce(wks,laimod_max,resg) |
|
361 frame(wks) |
|
362 |
|
363 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
|
364 ; system("rm "+plot_name+"."+plot_type) |
|
365 ; system("rm "+plot_name+"-1."+plot_type_new) |
|
366 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) |
|
367 |
|
368 clear (wks) |
|
369 ;------------------------------------------------------------------------ |
|
370 ;global contour model vs ob |
|
371 |
|
372 plot_name = "global_max_model_vs_ob" |
|
373 |
|
374 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
375 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
376 |
|
377 delete (plot) |
|
378 plot=new(3,graphic) ; create graphic array |
|
379 |
|
380 resg@gsnFrame = False ; Do not draw plot |
|
381 resg@gsnDraw = False ; Do not advance frame |
|
382 |
|
383 ; plot correlation coef |
|
384 |
|
385 gRes = True |
|
386 gRes@txFontHeightF = 0.02 |
|
387 gRes@txAngleF = 90 |
|
388 |
|
389 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")" |
|
390 |
|
391 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
|
392 ;-------------------------------------------------------------------- |
|
393 |
|
394 ;(a) ob |
|
395 |
|
396 title = ob_name |
|
397 resg@tiMainString = title |
|
398 |
|
399 plot(0) = gsn_csm_contour_map_ce(wks,laiob_max,resg) |
|
400 |
|
401 ;(b) model |
|
402 |
|
403 title = "Model "+ model_name |
|
404 resg@tiMainString = title |
|
405 |
|
406 plot(1) = gsn_csm_contour_map_ce(wks,laimod_max,resg) |
|
407 |
|
408 ;(c) model-ob |
|
409 |
|
410 zz = laimod_max |
|
411 zz = laimod_max - laiob_max |
|
412 title = "Model_"+model_name+" - Observed" |
|
413 resg@tiMainString = title |
|
414 |
|
415 resg@cnMinLevelValF = -6. ; Min level |
|
416 resg@cnMaxLevelValF = 6. ; Max level |
|
417 resg@cnLevelSpacingF = 1. ; interval |
|
418 |
|
419 |
|
420 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) |
|
421 |
|
422 pres = True ; panel plot mods desired |
|
423 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around |
|
424 ; indiv. plots in panel |
|
425 pres@gsnMaximize = True ; fill the page |
|
426 |
|
427 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
|
428 |
|
429 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
|
430 ; system("rm "+plot_name+"."+plot_type) |
|
431 ; system("rm "+plot_name+"-1."+plot_type_new) |
|
432 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) |
|
433 |
|
434 frame (wks) |
|
435 clear (wks) |
|
436 |
|
437 delete (plot) |
|
438 ;----------------------------------------------------------------- |
|
439 ;(C) phase |
|
440 ;-------------------------------------------------------------------- |
|
441 ; get data |
|
442 |
|
443 ; observed |
|
444 laiob_phase = laiob(0,:,:) |
|
445 s = laiob(:,0,0) |
|
446 laiob_phase@long_name = "Leaf Area Index Max Month" |
|
447 |
|
448 dsizes_z = dimsizes(laiob) |
|
449 nlat = dsizes_z(1) |
|
450 nlon = dsizes_z(2) |
|
451 |
|
452 do j = 0,nlat-1 |
|
453 do i = 0,nlon-1 |
|
454 s = laiob(:,j,i) |
|
455 laiob_phase(j,i) = maxind(s) + 1 |
|
456 end do |
|
457 end do |
|
458 |
|
459 ; print (min(laiob_phase)+"/"+max(laiob_phase)) |
|
460 delete (s) |
|
461 delete (dsizes_z) |
|
462 ;------------------------- |
|
463 ; model |
|
464 laimod_phase = laimod(0,:,:) |
|
465 s = laimod(:,0,0) |
|
466 laimod_phase@long_name = "Leaf Area Index Max Month" |
|
467 |
|
468 dsizes_z = dimsizes(laimod) |
|
469 nlat = dsizes_z(1) |
|
470 nlon = dsizes_z(2) |
|
471 |
|
472 do j = 0,nlat-1 |
|
473 do i = 0,nlon-1 |
|
474 s = laimod(:,j,i) |
|
475 laimod_phase(j,i) = maxind(s) + 1 |
|
476 end do |
|
477 end do |
|
478 |
|
479 ; print (min(laimod_phase)+"/"+max(laimod_phase)) |
|
480 delete (s) |
|
481 delete (dsizes_z) |
|
482 ;------------------------ |
|
483 DATA11_1D = ndtooned(classob) |
|
484 DATA12_1D = ndtooned(laiob_phase) |
|
485 DATA22_1D = ndtooned(laimod_phase) |
|
486 |
|
487 yvalues = new((/2,nx/),float) |
|
488 mn_yvalues = new((/2,nx/),float) |
|
489 mx_yvalues = new((/2,nx/),float) |
|
490 |
|
491 do nd=0,1 |
|
492 |
|
493 ; See if we are doing model or observational data. |
|
494 |
|
495 if(nd.eq.0) then |
|
496 data_ob = DATA11_1D |
|
497 data_mod = DATA12_1D |
|
498 else |
|
499 data_ob = DATA11_1D |
|
500 data_mod = DATA22_1D |
|
501 end if |
|
502 |
|
503 ; Loop through each range and check for values. |
|
504 |
|
505 do i=0,nr-2 |
|
506 if (i.ne.(nr-2)) then |
|
507 ; print("") |
|
508 ; print("In range ["+range(i)+","+range(i+1)+")") |
|
509 idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) |
|
510 else |
|
511 ; print("") |
|
512 ; print("In range ["+range(i)+",)") |
|
513 idx = ind(range(i).le.data_ob) |
|
514 end if |
|
515 |
|
516 ; Calculate average, and get min and max. |
|
517 |
|
518 if(.not.any(ismissing(idx))) then |
|
519 yvalues(nd,i) = avg(data_mod(idx)) |
|
520 mn_yvalues(nd,i) = min(data_mod(idx)) |
|
521 mx_yvalues(nd,i) = max(data_mod(idx)) |
|
522 count = dimsizes(idx) |
|
523 else |
|
524 count = 0 |
|
525 yvalues(nd,i) = yvalues@_FillValue |
|
526 mn_yvalues(nd,i) = yvalues@_FillValue |
|
527 mx_yvalues(nd,i) = yvalues@_FillValue |
|
528 end if |
|
529 |
|
530 ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) |
|
531 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) |
|
532 |
|
533 ; Clean up for next time in loop. |
|
534 |
|
535 delete(idx) |
|
536 end do |
|
537 delete(data_ob) |
|
538 delete(data_mod) |
|
539 end do |
|
540 ;----------------------------------------------------------------- |
|
541 ; compute correlation coef and M score |
|
542 |
|
543 u = yvalues(0,:) |
|
544 v = yvalues(1,:) |
|
545 |
|
546 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
547 uu = u(good) |
|
548 vv = v(good) |
|
549 |
|
550 ccrPhase = esccr(uu,vv,0) |
|
551 ; print (ccrPhase) |
|
552 |
|
553 ; old eq |
|
554 ; bias = abs(avg(vv)-avg(uu)) |
|
555 ; new eq |
|
556 bias = avg(abs(vv-uu)) |
|
557 |
|
558 bias = where((bias.gt. 6.),12.-bias,bias) |
|
559 Mphase = ((6. - bias)/6.)*5. |
|
560 |
|
561 print (Mphase) |
|
562 |
|
563 delete (u) |
|
564 delete (v) |
|
565 delete (uu) |
|
566 delete (vv) |
|
567 ;-------------------------------------------------------------------- |
|
568 ; histogram res |
|
569 |
|
570 resm = True |
|
571 resm@gsnMaximize = True |
|
572 resm@gsnDraw = False |
|
573 resm@gsnFrame = False |
|
574 resm@xyMarkLineMode = "Markers" |
|
575 resm@xyMarkerSizeF = 0.014 |
|
576 resm@xyMarker = 16 |
|
577 resm@xyMarkerColors = (/"Brown","Blue"/) |
|
578 ; resm@trYMinF = min(mn_yvalues) - 10. |
|
579 ; resm@trYMaxF = max(mx_yvalues) + 10. |
|
580 resm@trYMinF = min(mn_yvalues) - 2 |
|
581 resm@trYMaxF = max(mx_yvalues) + 4 |
|
582 |
|
583 resm@tiYAxisString = "Max LAI (Leaf Area Index) Month" |
|
584 resm@tiXAxisString = "Land Cover Type" |
|
585 |
|
586 max_bar = new((/2,nx/),graphic) |
|
587 min_bar = new((/2,nx/),graphic) |
|
588 max_cap = new((/2,nx/),graphic) |
|
589 min_cap = new((/2,nx/),graphic) |
|
590 |
|
591 lnres = True |
|
592 line_colors = (/"brown","blue"/) |
|
593 ;------------------------------------------------------------------ |
|
594 ; Start the graphics. |
|
595 |
|
596 plot_name = "histogram_phase" |
|
597 title = model_name + " vs Observed" |
|
598 resm@tiMainString = title |
|
599 |
|
600 wks = gsn_open_wks (plot_type,plot_name) |
|
601 ;----------------------------- |
|
602 ; Add a boxed legend using the more simple method |
|
603 |
|
604 resm@pmLegendDisplayMode = "Always" |
|
605 ; resm@pmLegendWidthF = 0.1 |
|
606 resm@pmLegendWidthF = 0.08 |
|
607 resm@pmLegendHeightF = 0.05 |
|
608 resm@pmLegendOrthogonalPosF = -1.17 |
|
609 ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) |
|
610 ; resm@pmLegendParallelPosF = 0.18 |
|
611 resm@pmLegendParallelPosF = 0.88 ;(rightward) |
|
612 |
|
613 ; resm@lgPerimOn = False |
|
614 resm@lgLabelFontHeightF = 0.015 |
|
615 resm@xyExplicitLegendLabels = (/"observed",model_name/) |
|
616 ;----------------------------- |
|
617 tRes = True |
|
618 tRes@txFontHeightF = 0.025 |
|
619 |
|
620 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")" |
|
621 |
|
622 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) |
|
623 |
|
624 xy = gsn_csm_xy(wks,xvalues,yvalues,resm) |
|
625 ;------------------------------- |
|
626 ;Attach the vertical bar and the horizontal cap line |
|
627 |
|
628 do nd=0,1 |
|
629 lnres@gsLineColor = line_colors(nd) |
|
630 do i=0,nx-1 |
|
631 |
|
632 if(.not.ismissing(mn_yvalues(nd,i)).and. \ |
|
633 .not.ismissing(mx_yvalues(nd,i))) then |
|
634 |
|
635 ; Attach the vertical bar, both above and below the marker. |
|
636 |
|
637 x1 = xvalues(nd,i) |
|
638 y1 = yvalues(nd,i) |
|
639 y2 = mn_yvalues(nd,i) |
|
640 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
641 |
|
642 y2 = mx_yvalues(nd,i) |
|
643 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
644 |
|
645 ; Attach the horizontal cap line, both above and below the marker. |
|
646 |
|
647 x1 = xvalues(nd,i) - dx4 |
|
648 x2 = xvalues(nd,i) + dx4 |
|
649 y1 = mn_yvalues(nd,i) |
|
650 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
651 |
|
652 y1 = mx_yvalues(nd,i) |
|
653 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
654 end if |
|
655 end do |
|
656 end do |
|
657 |
|
658 draw(xy) |
|
659 frame(wks) |
|
660 |
|
661 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
|
662 ; system("rm "+plot_name+"."+plot_type) |
|
663 ; system("rm "+plot_name+"-1."+plot_type_new) |
|
664 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) |
|
665 |
|
666 clear (wks) |
|
667 |
|
668 delete (DATA11_1D) |
|
669 delete (DATA12_1D) |
|
670 delete (DATA22_1D) |
|
671 ;delete (range) |
|
672 ;delete (xvalues) |
|
673 delete (yvalues) |
|
674 delete (mn_yvalues) |
|
675 delete (mx_yvalues) |
|
676 delete (good) |
|
677 delete (max_bar) |
|
678 delete (min_bar) |
|
679 delete (max_cap) |
|
680 delete (min_cap) |
|
681 ;----------------------------------------------------------------- |
|
682 ;global res |
|
683 |
|
684 resg = True ; Use plot options |
|
685 resg@cnFillOn = True ; Turn on color fill |
|
686 resg@gsnSpreadColors = True ; use full colormap |
|
687 ; resg@cnFillMode = "RasterFill" ; Turn on raster color |
|
688 ; resg@lbLabelAutoStride = True |
|
689 resg@cnLinesOn = False ; Turn off contourn lines |
|
690 resg@mpFillOn = False ; Turn off map fill |
|
691 |
|
692 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
693 resg@cnMinLevelValF = 1. ; Min level |
|
694 resg@cnMaxLevelValF = 12. ; Max level |
|
695 resg@cnLevelSpacingF = 1. ; interval |
|
696 |
|
697 ;global contour ob |
|
698 |
|
699 delta = 0.00001 |
|
700 laiob_phase = where(ismissing(laiob_phase).and.(ismissing(laimod_phase).or.(laimod_phase.lt.delta)),0.,laiob_phase) |
|
701 |
|
702 plot_name = "global_phase_ob" |
|
703 title = ob_name |
|
704 resg@tiMainString = title |
|
705 |
|
706 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
707 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
708 |
|
709 plot = gsn_csm_contour_map_ce(wks,laiob_phase,resg) |
|
710 frame(wks) |
|
711 |
|
712 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
|
713 ; system("rm "+plot_name+"."+plot_type) |
|
714 ; system("rm "+plot_name+"-1."+plot_type_new) |
|
715 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) |
|
716 |
|
717 clear (wks) |
|
718 ;------------------------------------------------------------------------ |
|
719 ;global contour model |
|
720 |
|
721 plot_name = "global_phase_model" |
|
722 title = "Model " + model_name |
|
723 resg@tiMainString = title |
|
724 |
|
725 wks = gsn_open_wks (plot_type,plot_name) |
|
726 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
727 |
|
728 delete (plot) |
|
729 plot = gsn_csm_contour_map_ce(wks,laimod_phase,resg) |
|
730 frame(wks) |
|
731 |
|
732 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
|
733 ; system("rm "+plot_name+"."+plot_type) |
|
734 ; system("rm "+plot_name+"-1."+plot_type_new) |
|
735 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) |
|
736 |
|
737 clear (wks) |
|
738 ;------------------------------------------------------------------------ |
|
739 ;global contour model vs ob |
|
740 |
|
741 plot_name = "global_phase_model_vs_ob" |
|
742 |
|
743 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
744 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
745 |
|
746 delete (plot) |
|
747 plot=new(3,graphic) ; create graphic array |
|
748 |
|
749 resg@gsnFrame = False ; Do not draw plot |
|
750 resg@gsnDraw = False ; Do not advance frame |
|
751 |
|
752 ; plot correlation coef |
|
753 |
|
754 gRes = True |
|
755 gRes@txFontHeightF = 0.02 |
|
756 gRes@txAngleF = 90 |
|
757 |
|
758 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")" |
|
759 |
|
760 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
|
761 ;-------------------------------------------------------------------- |
|
762 |
|
763 ;(a) ob |
|
764 |
|
765 title = ob_name |
|
766 resg@tiMainString = title |
|
767 |
|
768 plot(0) = gsn_csm_contour_map_ce(wks,laiob_phase,resg) |
|
769 |
|
770 ;(b) model |
|
771 |
|
772 title = "Model "+ model_name |
|
773 resg@tiMainString = title |
|
774 |
|
775 plot(1) = gsn_csm_contour_map_ce(wks,laimod_phase,resg) |
|
776 |
|
777 ;(c) model-ob |
|
778 |
|
779 delete (zz) |
|
780 zz = laimod_phase |
|
781 zz = laimod_phase - laiob_phase |
|
782 title = "Model_"+model_name+" - Observed" |
|
783 resg@tiMainString = title |
|
784 |
|
785 resg@cnMinLevelValF = -6. ; Min level |
|
786 resg@cnMaxLevelValF = 6. ; Max level |
|
787 resg@cnLevelSpacingF = 1. ; interval |
|
788 |
|
789 |
|
790 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) |
|
791 |
|
792 ; pres = True ; panel plot mods desired |
|
793 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around |
|
794 ; indiv. plots in panel |
|
795 ; pres@gsnMaximize = True ; fill the page |
|
796 |
|
797 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
|
798 |
|
799 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
|
800 ; system("rm "+plot_name+"."+plot_type) |
|
801 ; system("rm "+plot_name+"-1."+plot_type_new) |
|
802 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) |
|
803 |
|
804 frame (wks) |
|
805 clear (wks) |
|
806 |
|
807 delete (plot) |
|
808 ;----------------------------------------------------------------- |
|
809 ;(D) grow day |
|
810 ;-------------------------------------------------------------------- |
|
811 ; get data |
|
812 |
|
813 day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/) |
|
814 |
|
815 ; observed |
|
816 laiob_grow = laiob(0,:,:) |
|
817 laiob_grow@long_name = "Days of Growing Season" |
|
818 |
|
819 dsizes_z = dimsizes(laiob) |
|
820 ntime = dsizes_z(0) |
|
821 nlat = dsizes_z(1) |
|
822 nlon = dsizes_z(2) |
|
823 |
|
824 do j = 0,nlat-1 |
|
825 do i = 0,nlon-1 |
|
826 nday = 0. |
|
827 do k = 0,ntime-1 |
|
828 if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then |
|
829 nday = nday + day_of_data(k) |
|
830 end if |
|
831 end do |
|
832 |
|
833 laiob_grow(j,i) = nday |
|
834 end do |
|
835 end do |
|
836 |
|
837 ; print (min(laiob_grow)+"/"+max(laiob_grow)) |
|
838 ;------------------------- |
|
839 ; model |
|
840 laimod_grow = laimod(0,:,:) |
|
841 laimod_grow@long_name = "Days of Growing Season" |
|
842 |
|
843 dsizes_z = dimsizes(laimod) |
|
844 ntime = dsizes_z(0) |
|
845 nlat = dsizes_z(1) |
|
846 nlon = dsizes_z(2) |
|
847 |
|
848 do j = 0,nlat-1 |
|
849 do i = 0,nlon-1 |
|
850 nday = 0. |
|
851 do k = 0,ntime-1 |
|
852 if (.not. ismissing(laimod(k,j,i)) .and. laimod(k,j,i) .gt. 1.0) then |
|
853 nday = nday + day_of_data(k) |
|
854 end if |
|
855 end do |
|
856 |
|
857 laimod_grow(j,i) = nday |
|
858 end do |
|
859 end do |
|
860 |
|
861 ; print (min(laimod_grow)+"/"+max(laimod_grow)) |
|
862 ;------------------------ |
|
863 DATA11_1D = ndtooned(classob) |
|
864 DATA12_1D = ndtooned(laiob_grow) |
|
865 DATA22_1D = ndtooned(laimod_grow) |
|
866 |
|
867 yvalues = new((/2,nx/),float) |
|
868 mn_yvalues = new((/2,nx/),float) |
|
869 mx_yvalues = new((/2,nx/),float) |
|
870 |
|
871 do nd=0,1 |
|
872 |
|
873 ; See if we are doing model or observational data. |
|
874 |
|
875 if(nd.eq.0) then |
|
876 data_ob = DATA11_1D |
|
877 data_mod = DATA12_1D |
|
878 else |
|
879 data_ob = DATA11_1D |
|
880 data_mod = DATA22_1D |
|
881 end if |
|
882 |
|
883 ; Loop through each range and check for values. |
|
884 |
|
885 do i=0,nr-2 |
|
886 if (i.ne.(nr-2)) then |
|
887 ; print("") |
|
888 ; print("In range ["+range(i)+","+range(i+1)+")") |
|
889 idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) |
|
890 else |
|
891 ; print("") |
|
892 ; print("In range ["+range(i)+",)") |
|
893 idx = ind(range(i).le.data_ob) |
|
894 end if |
|
895 |
|
896 ; Calculate average, and get min and max. |
|
897 |
|
898 if(.not.any(ismissing(idx))) then |
|
899 yvalues(nd,i) = avg(data_mod(idx)) |
|
900 mn_yvalues(nd,i) = min(data_mod(idx)) |
|
901 mx_yvalues(nd,i) = max(data_mod(idx)) |
|
902 count = dimsizes(idx) |
|
903 else |
|
904 count = 0 |
|
905 yvalues(nd,i) = yvalues@_FillValue |
|
906 mn_yvalues(nd,i) = yvalues@_FillValue |
|
907 mx_yvalues(nd,i) = yvalues@_FillValue |
|
908 end if |
|
909 |
|
910 ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) |
|
911 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) |
|
912 |
|
913 ; Clean up for next time in loop. |
|
914 |
|
915 delete(idx) |
|
916 end do |
|
917 delete(data_ob) |
|
918 delete(data_mod) |
|
919 end do |
|
920 ;----------------------------------------------------------------- |
|
921 ; compute correlation coef and M score |
|
922 |
|
923 u = yvalues(0,:) |
|
924 v = yvalues(1,:) |
|
925 |
|
926 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
927 uu = u(good) |
|
928 vv = v(good) |
|
929 |
|
930 ccrGrow = esccr(uu,vv,0) |
|
931 ; print (ccrGrow) |
|
932 |
|
933 ; new eq |
|
934 bias = sum(abs(vv-uu)/(vv+uu)) |
|
935 Mgrow = (1.- (bias/dimsizes(uu)))*5. |
|
936 |
|
937 print (Mgrow) |
|
938 |
|
939 delete (u) |
|
940 delete (v) |
|
941 delete (uu) |
|
942 delete (vv) |
|
943 ;-------------------------------------------------------------------- |
|
944 ; histogram res |
|
945 |
|
946 resm = True |
|
947 resm@gsnMaximize = True |
|
948 resm@gsnDraw = False |
|
949 resm@gsnFrame = False |
|
950 resm@xyMarkLineMode = "Markers" |
|
951 resm@xyMarkerSizeF = 0.014 |
|
952 resm@xyMarker = 16 |
|
953 resm@xyMarkerColors = (/"Brown","Blue"/) |
|
954 ; resm@trYMinF = min(mn_yvalues) - 10. |
|
955 ; resm@trYMaxF = max(mx_yvalues) + 10. |
|
956 resm@trYMinF = min(mn_yvalues) - 2 |
|
957 resm@trYMaxF = max(mx_yvalues) + 4 |
|
958 |
|
959 resm@tiYAxisString = "Days of Growing season" |
|
960 resm@tiXAxisString = "Land Cover Type" |
|
961 |
|
962 max_bar = new((/2,nx/),graphic) |
|
963 min_bar = new((/2,nx/),graphic) |
|
964 max_cap = new((/2,nx/),graphic) |
|
965 min_cap = new((/2,nx/),graphic) |
|
966 |
|
967 lnres = True |
|
968 line_colors = (/"brown","blue"/) |
|
969 ;------------------------------------------------------------------ |
|
970 ; Start the graphics. |
|
971 |
|
972 plot_name = "histogram_grow" |
|
973 title = model_name + " vs Observed" |
|
974 resm@tiMainString = title |
|
975 |
|
976 wks = gsn_open_wks (plot_type,plot_name) |
|
977 ;----------------------------- |
|
978 ; Add a boxed legend using the more simple method |
|
979 |
|
980 resm@pmLegendDisplayMode = "Always" |
|
981 ; resm@pmLegendWidthF = 0.1 |
|
982 resm@pmLegendWidthF = 0.08 |
|
983 resm@pmLegendHeightF = 0.05 |
|
984 resm@pmLegendOrthogonalPosF = -1.17 |
|
985 ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) |
|
986 ; resm@pmLegendParallelPosF = 0.18 |
|
987 resm@pmLegendParallelPosF = 0.88 ;(rightward) |
|
988 |
|
989 ; resm@lgPerimOn = False |
|
990 resm@lgLabelFontHeightF = 0.015 |
|
991 resm@xyExplicitLegendLabels = (/"observed",model_name/) |
|
992 ;----------------------------- |
|
993 tRes = True |
|
994 tRes@txFontHeightF = 0.025 |
|
995 |
|
996 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")" |
|
997 |
|
998 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) |
|
999 |
|
1000 xy = gsn_csm_xy(wks,xvalues,yvalues,resm) |
|
1001 ;------------------------------- |
|
1002 ;Attach the vertical bar and the horizontal cap line |
|
1003 |
|
1004 do nd=0,1 |
|
1005 lnres@gsLineColor = line_colors(nd) |
|
1006 do i=0,nx-1 |
|
1007 |
|
1008 if(.not.ismissing(mn_yvalues(nd,i)).and. \ |
|
1009 .not.ismissing(mx_yvalues(nd,i))) then |
|
1010 |
|
1011 ; Attach the vertical bar, both above and below the marker. |
|
1012 |
|
1013 x1 = xvalues(nd,i) |
|
1014 y1 = yvalues(nd,i) |
|
1015 y2 = mn_yvalues(nd,i) |
|
1016 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
1017 |
|
1018 y2 = mx_yvalues(nd,i) |
|
1019 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
1020 |
|
1021 ; Attach the horizontal cap line, both above and below the marker. |
|
1022 |
|
1023 x1 = xvalues(nd,i) - dx4 |
|
1024 x2 = xvalues(nd,i) + dx4 |
|
1025 y1 = mn_yvalues(nd,i) |
|
1026 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
1027 |
|
1028 y1 = mx_yvalues(nd,i) |
|
1029 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
1030 end if |
|
1031 end do |
|
1032 end do |
|
1033 |
|
1034 draw(xy) |
|
1035 frame(wks) |
|
1036 |
|
1037 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
|
1038 ; system("rm "+plot_name+"."+plot_type) |
|
1039 ; system("rm "+plot_name+"-1."+plot_type_new) |
|
1040 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) |
|
1041 |
|
1042 clear (wks) |
|
1043 |
|
1044 delete (DATA11_1D) |
|
1045 delete (DATA12_1D) |
|
1046 delete (DATA22_1D) |
|
1047 ;delete (range) |
|
1048 ;delete (xvalues) |
|
1049 delete (yvalues) |
|
1050 delete (mn_yvalues) |
|
1051 delete (mx_yvalues) |
|
1052 delete (good) |
|
1053 delete (max_bar) |
|
1054 delete (min_bar) |
|
1055 delete (max_cap) |
|
1056 delete (min_cap) |
|
1057 ;----------------------------------------------------------------- |
|
1058 ;global res |
|
1059 |
|
1060 resg = True ; Use plot options |
|
1061 resg@cnFillOn = True ; Turn on color fill |
|
1062 resg@gsnSpreadColors = True ; use full colormap |
|
1063 ; resg@cnFillMode = "RasterFill" ; Turn on raster color |
|
1064 ; resg@lbLabelAutoStride = True |
|
1065 resg@cnLinesOn = False ; Turn off contourn lines |
|
1066 resg@mpFillOn = False ; Turn off map fill |
|
1067 |
|
1068 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
1069 resg@cnMinLevelValF = 0. ; Min level |
|
1070 resg@cnMaxLevelValF = 360. ; Max level |
|
1071 resg@cnLevelSpacingF = 30. ; interval |
|
1072 |
|
1073 ;global contour ob |
|
1074 |
|
1075 delta = 0.00001 |
|
1076 laiob_grow = where(ismissing(laiob_grow).and.(ismissing(laimod_grow).or.(laimod_grow.lt.delta)),0.,laiob_grow) |
|
1077 |
|
1078 plot_name = "global_grow_ob" |
|
1079 title = ob_name |
|
1080 resg@tiMainString = title |
|
1081 |
|
1082 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1083 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1084 |
|
1085 plot = gsn_csm_contour_map_ce(wks,laiob_grow,resg) |
|
1086 frame(wks) |
|
1087 |
|
1088 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
|
1089 ; system("rm "+plot_name+"."+plot_type) |
|
1090 ; system("rm "+plot_name+"-1."+plot_type_new) |
|
1091 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) |
|
1092 |
|
1093 clear (wks) |
|
1094 ;------------------------------------------------------------------------ |
|
1095 ;global contour model |
|
1096 |
|
1097 plot_name = "global_grow_model" |
|
1098 title = "Model " + model_name |
|
1099 resg@tiMainString = title |
|
1100 |
|
1101 wks = gsn_open_wks (plot_type,plot_name) |
|
1102 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1103 |
|
1104 delete (plot) |
|
1105 plot = gsn_csm_contour_map_ce(wks,laimod_grow,resg) |
|
1106 frame(wks) |
|
1107 |
|
1108 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
|
1109 ; system("rm "+plot_name+"."+plot_type) |
|
1110 ; system("rm "+plot_name+"-1."+plot_type_new) |
|
1111 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) |
|
1112 |
|
1113 clear (wks) |
|
1114 ;------------------------------------------------------------------------ |
|
1115 ;global contour model vs ob |
|
1116 |
|
1117 plot_name = "global_grow_model_vs_ob" |
|
1118 |
|
1119 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1120 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1121 |
|
1122 delete (plot) |
|
1123 plot=new(3,graphic) ; create graphic array |
|
1124 |
|
1125 resg@gsnFrame = False ; Do not draw plot |
|
1126 resg@gsnDraw = False ; Do not advance frame |
|
1127 |
|
1128 ; plot correlation coef |
|
1129 |
|
1130 gRes = True |
|
1131 gRes@txFontHeightF = 0.02 |
|
1132 gRes@txAngleF = 90 |
|
1133 |
|
1134 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")" |
|
1135 |
|
1136 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
|
1137 ;-------------------------------------------------------------------- |
|
1138 |
|
1139 ;(a) ob |
|
1140 |
|
1141 title = ob_name |
|
1142 resg@tiMainString = title |
|
1143 |
|
1144 plot(0) = gsn_csm_contour_map_ce(wks,laiob_grow,resg) |
|
1145 |
|
1146 ;(b) model |
|
1147 |
|
1148 title = "Model "+ model_name |
|
1149 resg@tiMainString = title |
|
1150 |
|
1151 plot(1) = gsn_csm_contour_map_ce(wks,laimod_grow,resg) |
|
1152 |
|
1153 ;(c) model-ob |
|
1154 |
|
1155 delete (zz) |
|
1156 zz = laimod_grow |
|
1157 zz = laimod_grow - laiob_grow |
|
1158 title = "Model_"+model_name+" - Observed" |
|
1159 resg@tiMainString = title |
|
1160 |
|
1161 resg@cnMinLevelValF = -120. ; Min level |
|
1162 resg@cnMaxLevelValF = 120. ; Max level |
|
1163 resg@cnLevelSpacingF = 20. ; interval |
|
1164 |
|
1165 |
|
1166 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) |
|
1167 |
|
1168 pres = True ; panel plot mods desired |
|
1169 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around |
|
1170 ; indiv. plots in panel |
|
1171 pres@gsnMaximize = True ; fill the page |
|
1172 |
|
1173 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
|
1174 |
|
1175 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
|
1176 ; system("rm "+plot_name+"."+plot_type) |
|
1177 ; system("rm "+plot_name+"-1."+plot_type_new) |
|
1178 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) |
|
1179 |
|
1180 frame (wks) |
|
1181 clear (wks) |
|
1182 |
|
1183 delete (plot) |
|
1184 end |
|
1185 |