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