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