|
1 ;************************************************************** |
|
2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
|
3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
|
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" |
|
5 ;************************************************************** |
|
6 procedure set_line(lines:string,nline:integer,newlines:string) |
|
7 begin |
|
8 ; add line to ascci/html file |
|
9 |
|
10 nnewlines = dimsizes(newlines) |
|
11 if(nline+nnewlines-1.ge.dimsizes(lines)) |
|
12 print("set_line: bad index, not setting anything.") |
|
13 return |
|
14 end if |
|
15 lines(nline:nline+nnewlines-1) = newlines |
|
16 ; print ("lines = " + lines(nline:nline+nnewlines-1)) |
|
17 nline = nline + nnewlines |
|
18 return |
|
19 end |
|
20 ;************************************************************** |
|
21 ; Main code. |
|
22 begin |
|
23 |
|
24 plot_type = "ps" |
|
25 plot_type_new = "png" |
|
26 |
|
27 ;----------------------------------------------------- |
|
28 ; edit table.html of current model for movel1_vs_model2 |
|
29 |
|
30 if (isvar("compare")) then |
|
31 html_name2 = compare+"/table.html" |
|
32 html_new2 = html_name2 +".new" |
|
33 end if |
|
34 |
|
35 ;------------------------------------------------------ |
|
36 ; edit table.html for current model |
|
37 |
|
38 html_name = model_name+"/table.html" |
|
39 html_new = html_name +".new" |
|
40 |
|
41 ;------------------------------------------------------ |
|
42 ; read data: model |
|
43 |
|
44 fm = addfile(dirm+film10,"r") |
|
45 laimod = fm->TLAI |
|
46 |
|
47 delete (fm) |
|
48 |
|
49 dsizes = dimsizes(laimod) |
|
50 ntime = dsizes(0) |
|
51 nlat = dsizes(1) |
|
52 nlon = dsizes(2) |
|
53 |
|
54 ;----------------------------------- |
|
55 ; get landfrac data |
|
56 |
|
57 film_l = "lnd_"+ model_grid + ".nc" |
|
58 fm_l = addfile (dirs+film_l,"r") |
|
59 landfrac = fm_l->landfrac |
|
60 |
|
61 delete (fm_l) |
|
62 ;---------------------------------- |
|
63 ; read biome data: model |
|
64 |
|
65 biome_name_mod = "Model PFT Class" |
|
66 |
|
67 film_c = "class_pft_"+model_grid+".nc" |
|
68 fm_c = addfile(dirs+film_c,"r") |
|
69 classmod = fm_c->CLASS_PFT |
|
70 |
|
71 delete (fm_c) |
|
72 |
|
73 ; model data has 17 land-type classes |
|
74 nclass_mod = 17 |
|
75 |
|
76 ;---------------------------------------------------------- |
|
77 ; read data: ob |
|
78 |
|
79 ;---------------------------------- |
|
80 ; read biome data: observed |
|
81 |
|
82 biome_name_ob = "MODIS LandCover" |
|
83 |
|
84 dir_c = diro + "lai/" |
|
85 filo_c = "land_class_"+model_grid+".nc" |
|
86 fo = addfile(dir_c+filo_c,"r") |
|
87 classob = tofloat(fo->LAND_CLASS) |
|
88 |
|
89 delete (fo) |
|
90 |
|
91 ; input observed data has 20 land-type classes |
|
92 nclass_ob = 20 |
|
93 |
|
94 ;--------------------------------- |
|
95 ; read lai data: observed |
|
96 |
|
97 ;ob_name = "MODIS MOD 15A2 2000-2005" |
|
98 ob_name = "MODIS MOD 15A2 2000-2004" |
|
99 |
|
100 dir_l = diro + "lai/" |
|
101 ;filo = "LAI_2000-2005_MONS_"+model_grid+".nc" |
|
102 filo = "LAI_2000-2004_MONS_"+model_grid+".nc" |
|
103 fo = addfile(dir_l+filo,"r") |
|
104 laiob = fo->LAI |
|
105 |
|
106 delete (fo) |
|
107 |
|
108 ;------------------------------------------------- |
|
109 ; take into account landfrac |
|
110 |
|
111 laimod = laimod * conform(laimod,landfrac,(/1,2/)) |
|
112 laiob = laiob * conform(laiob ,landfrac,(/1,2/)) |
|
113 |
|
114 delete (landfrac) |
|
115 |
|
116 ;************************************************ |
|
117 ; global res |
|
118 ;************************************************ |
|
119 resg = True ; Use plot options |
|
120 resg@cnFillOn = True ; Turn on color fill |
|
121 resg@gsnSpreadColors = True ; use full colormap |
|
122 resg@cnLinesOn = False ; Turn off contourn lines |
|
123 resg@mpFillOn = False ; Turn off map fill |
|
124 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
125 |
|
126 ;************************************************ |
|
127 ; plot global biome class: (1) observed |
|
128 ;************************************************ |
|
129 |
|
130 resg@cnMinLevelValF = 1. ; Min level |
|
131 resg@cnMaxLevelValF = 19. ; Max level |
|
132 resg@cnLevelSpacingF = 1. ; interval |
|
133 |
|
134 classob@_FillValue = 1.e+36 |
|
135 classob = where(classob.eq.0,classob@_FillValue,classob) |
|
136 |
|
137 plot_name = "global_class_ob" |
|
138 title = biome_name_ob |
|
139 resg@tiMainString = title |
|
140 |
|
141 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
142 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
143 |
|
144 plot = gsn_csm_contour_map_ce(wks,classob,resg) |
|
145 |
|
146 delete (wks) |
|
147 |
|
148 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
149 "rm "+plot_name+"."+plot_type) |
|
150 |
|
151 ;************************************************ |
|
152 ; plot global biome class: (2) model |
|
153 ;************************************************ |
|
154 |
|
155 resg@cnMinLevelValF = 0. ; Min level |
|
156 resg@cnMaxLevelValF = 16. ; Max level |
|
157 resg@cnLevelSpacingF = 1. ; interval |
|
158 |
|
159 plot_name = "global_class_model" |
|
160 title = biome_name_mod |
|
161 resg@tiMainString = title |
|
162 |
|
163 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
164 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
165 |
|
166 plot = gsn_csm_contour_map_ce(wks,classmod,resg) |
|
167 |
|
168 delete (wks) |
|
169 |
|
170 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
171 "rm "+plot_name+"."+plot_type) |
|
172 |
|
173 ;******************************************************************* |
|
174 ; for html table : all 3 components (Mean, Max, Phase) |
|
175 ;******************************************************************* |
|
176 |
|
177 ; column (not including header column) |
|
178 |
|
179 component = (/"Mean","Max","Phase"/) |
|
180 |
|
181 col_head = (/"observed",model_name,"M_score" \ |
|
182 ,"observed",model_name,"M_score" \ |
|
183 ,"observed",model_name,"M_score" \ |
|
184 /) |
|
185 |
|
186 n_comp = dimsizes(component) |
|
187 ncol = dimsizes(col_head ) |
|
188 |
|
189 ; row (not including header row) |
|
190 |
|
191 ;---------------------------------------------------- |
|
192 ; using model biome class: |
|
193 row_head = (/"Not Vegetated" \ |
|
194 ,"Needleleaf Evergreen Temperate Tree" \ |
|
195 ,"Needleleaf Evergreen Boreal Tree" \ |
|
196 ; ,"Needleleaf Deciduous Boreal Tree" \ |
|
197 ,"Broadleaf Evergreen Tropical Tree" \ |
|
198 ,"Broadleaf Evergreen Temperate Tree" \ |
|
199 ,"Broadleaf Deciduous Tropical Tree" \ |
|
200 ,"Broadleaf Deciduous Temperate Tree" \ |
|
201 ; ,"Broadleaf Deciduous Boreal Tree" \ |
|
202 ; ,"Broadleaf Evergreen Shrub" \ |
|
203 ,"Broadleaf Deciduous Temperate Shrub" \ |
|
204 ,"Broadleaf Deciduous Boreal Shrub" \ |
|
205 ,"C3 Arctic Grass" \ |
|
206 ,"C3 Non-Arctic Grass" \ |
|
207 ,"C4 Grass" \ |
|
208 ,"Corn" \ |
|
209 ; ,"Wheat" \ |
|
210 ,"All Biomes" \ |
|
211 /) |
|
212 |
|
213 nrow = dimsizes(row_head) |
|
214 |
|
215 ; arrays to be passed to table. |
|
216 text = new ((/nrow, ncol/),string ) |
|
217 |
|
218 ; total M_score |
|
219 M_total = 0. |
|
220 |
|
221 ;******************************************************************** |
|
222 ; use land-type class to bin the data in equally spaced ranges |
|
223 ;******************************************************************** |
|
224 |
|
225 ; using model biome class |
|
226 nclass = nclass_mod |
|
227 |
|
228 range = fspan(0,nclass,nclass+1) |
|
229 |
|
230 ; Use this range information to grab all the values in a |
|
231 ; particular range, and then take an average. |
|
232 |
|
233 nx = dimsizes(range) - 1 |
|
234 |
|
235 ; for 2 data: model and observed |
|
236 data_n = 2 |
|
237 |
|
238 ; using model biome class |
|
239 |
|
240 base = ndtooned(classmod) |
|
241 |
|
242 ; output |
|
243 |
|
244 yvalues = new((/data_n,nx/),float) |
|
245 count = new((/data_n,nx/),float) |
|
246 |
|
247 ;************************************************************************ |
|
248 ; go through all components |
|
249 ;************************************************************************ |
|
250 |
|
251 do m = 0,n_comp-1 |
|
252 |
|
253 ;=================== |
|
254 ; get data: |
|
255 ;=================== |
|
256 ; (A) Mean |
|
257 |
|
258 if (m .eq. 0) then |
|
259 data_ob = dim_avg_Wrap(laiob (lat|:,lon|:,time|:)) |
|
260 data_mod = dim_avg_Wrap(laimod(lat|:,lon|:,time|:)) |
|
261 end if |
|
262 |
|
263 ; (B) Max |
|
264 |
|
265 if (m .eq. 1) then |
|
266 |
|
267 ; observed |
|
268 data_ob = laiob(0,:,:) |
|
269 data_ob@long_name = "Leaf Area Index Max" |
|
270 |
|
271 do j = 0,nlat-1 |
|
272 do i = 0,nlon-1 |
|
273 data_ob(j,i) = max(laiob(:,j,i)) |
|
274 end do |
|
275 end do |
|
276 |
|
277 ; model |
|
278 data_mod = laimod(0,:,:) |
|
279 data_mod@long_name = "Leaf Area Index Max" |
|
280 |
|
281 do j = 0,nlat-1 |
|
282 do i = 0,nlon-1 |
|
283 data_mod(j,i) = max(laimod(:,j,i)) |
|
284 end do |
|
285 end do |
|
286 |
|
287 end if |
|
288 |
|
289 ; (C) phase |
|
290 |
|
291 if (m .eq. 2) then |
|
292 |
|
293 ; observed |
|
294 data_ob = laiob(0,:,:) |
|
295 data_ob@long_name = "Leaf Area Index Max Month" |
|
296 |
|
297 do j = 0,nlat-1 |
|
298 do i = 0,nlon-1 |
|
299 data_ob(j,i) = maxind(laiob(:,j,i)) + 1 |
|
300 end do |
|
301 end do |
|
302 |
|
303 ; model |
|
304 data_mod = laimod(0,:,:) |
|
305 data_mod@long_name = "Leaf Area Index Max Month" |
|
306 |
|
307 do j = 0,nlat-1 |
|
308 do i = 0,nlon-1 |
|
309 data_mod(j,i) = maxind(laimod(:,j,i)) + 1 |
|
310 end do |
|
311 end do |
|
312 |
|
313 end if |
|
314 |
|
315 ;============================== |
|
316 ; put data into bins |
|
317 ;============================== |
|
318 |
|
319 ; Loop through each range, using base |
|
320 |
|
321 do i=0,nx-1 |
|
322 |
|
323 if (i.ne.(nx-1)) then |
|
324 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1))) |
|
325 else |
|
326 idx = ind(base.ge.range(i)) |
|
327 end if |
|
328 |
|
329 ; loop through each dataset |
|
330 |
|
331 do n = 0,data_n-1 |
|
332 |
|
333 if (n .eq. 0) then |
|
334 data = ndtooned(data_ob) |
|
335 end if |
|
336 |
|
337 if (n .eq. 1) then |
|
338 data = ndtooned(data_mod) |
|
339 end if |
|
340 |
|
341 ; Calculate average |
|
342 |
|
343 if (.not.any(ismissing(idx))) then |
|
344 yvalues(n,i) = avg(data(idx)) |
|
345 count(n,i) = dimsizes(idx) |
|
346 else |
|
347 yvalues(n,i) = yvalues@_FillValue |
|
348 count(n,i) = 0 |
|
349 end if |
|
350 |
|
351 ;############################################################# |
|
352 ; using model biome class: |
|
353 ; |
|
354 ; set the following 4 classes to _FillValue: |
|
355 ; (3)Needleleaf Deciduous Boreal Tree, |
|
356 ; (8)Broadleaf Deciduous Boreal Tree, |
|
357 ; (9)Broadleaf Evergreen Shrub, |
|
358 ; (16)Wheat |
|
359 |
|
360 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then |
|
361 yvalues(n,i) = yvalues@_FillValue |
|
362 count(n,i) = 0 |
|
363 end if |
|
364 ;############################################################# |
|
365 |
|
366 delete(data) |
|
367 end do ; n-loop |
|
368 |
|
369 delete(idx) |
|
370 end do ; i-loop |
|
371 |
|
372 ;===================================== |
|
373 ; compute correlation coef and M score |
|
374 ;===================================== |
|
375 |
|
376 score_max = 5.0 |
|
377 |
|
378 u = yvalues(0,:) |
|
379 v = yvalues(1,:) |
|
380 u_count = count(0,:) |
|
381 v_count = count(1,:) |
|
382 |
|
383 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
384 |
|
385 uu = u(good) |
|
386 vv = v(good) |
|
387 uu_count = u_count(good) |
|
388 vv_count = v_count(good) |
|
389 |
|
390 ; compute correlation coef |
|
391 cc = esccr(uu,vv,0) |
|
392 |
|
393 if (n .eq. 2) then |
|
394 bias = avg(abs(vv-uu)) |
|
395 bias = where((bias.gt. 6.),12.-bias,bias) |
|
396 Mscore = ((6. - bias)/6.)*score_max |
|
397 else |
|
398 bias = sum(abs(vv-uu)/abs(vv+uu)) |
|
399 Mscore = (1.- (bias/dimsizes(uu)))*score_max |
|
400 end if |
|
401 |
|
402 M_score = sprintf("%.2f", Mscore) |
|
403 |
|
404 ; compute M_total |
|
405 |
|
406 M_total = M_total + Mscore |
|
407 |
|
408 ;================================ |
|
409 ; output M_score to score sheet |
|
410 ;=============================== |
|
411 |
|
412 if (isvar("compare")) then |
|
413 system("sed -e '1,/M_lai_"+component(m)+"/s/M_lai_"+component(m)+"/"+M_score+"/' "+html_name2+" > "+html_new2+";"+ \ |
|
414 "mv -f "+html_new2+" "+html_name2) |
|
415 end if |
|
416 |
|
417 system("sed s#M_lai_"+component(m)+"#"+M_score+"# "+html_name+" > "+html_new+";"+ \ |
|
418 "mv -f "+html_new+" "+html_name) |
|
419 |
|
420 ;============================== |
|
421 ; output M_score to html table |
|
422 ;============================== |
|
423 |
|
424 n = m*3 |
|
425 |
|
426 do i=0,nrow-2 |
|
427 text(i,n) = sprintf("%.1f",uu(i)) |
|
428 text(i,n+1) = sprintf("%.1f",vv(i)) |
|
429 text(i,n+2) = "-" |
|
430 end do |
|
431 text(nrow-1,n) = sprintf("%.1f",sum(uu*uu_count)/sum(uu_count)) |
|
432 text(nrow-1,n+1) = sprintf("%.1f",sum(vv*vv_count)/sum(vv_count)) |
|
433 text(nrow-1,n+2) = M_score |
|
434 |
|
435 delete (u) |
|
436 delete (v) |
|
437 delete (uu) |
|
438 delete (vv) |
|
439 delete (u_count) |
|
440 delete (v_count) |
|
441 delete (uu_count) |
|
442 delete (vv_count) |
|
443 delete (good) |
|
444 |
|
445 ;======================================== |
|
446 ; global res changes for each component |
|
447 ;======================================== |
|
448 delta = 0.00001 |
|
449 |
|
450 if (m .eq. 0) then |
|
451 resg@cnMinLevelValF = 0. |
|
452 resg@cnMaxLevelValF = 10. |
|
453 resg@cnLevelSpacingF = 1. |
|
454 |
|
455 data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob) |
|
456 end if |
|
457 |
|
458 if (m .eq. 1) then |
|
459 resg@cnMinLevelValF = 0. |
|
460 resg@cnMaxLevelValF = 10. |
|
461 resg@cnLevelSpacingF = 1. |
|
462 |
|
463 data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob) |
|
464 end if |
|
465 |
|
466 if (m .eq. 2) then |
|
467 resg@cnMinLevelValF = 1. |
|
468 resg@cnMaxLevelValF = 12. |
|
469 resg@cnLevelSpacingF = 1. |
|
470 |
|
471 data_mod = where(ismissing(data_ob).or.(data_mod.lt.delta),data_mod@_FillValue,data_mod) |
|
472 data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob) |
|
473 |
|
474 end if |
|
475 |
|
476 ;========================= |
|
477 ; global contour : ob |
|
478 ;========================= |
|
479 |
|
480 plot_name = "global_"+component(m)+"_ob" |
|
481 title = ob_name |
|
482 resg@tiMainString = title |
|
483 |
|
484 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
485 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
486 |
|
487 plot = gsn_csm_contour_map_ce(wks,data_ob,resg) |
|
488 |
|
489 delete (wks) |
|
490 delete (plot) |
|
491 |
|
492 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
493 "rm "+plot_name+"."+plot_type) |
|
494 |
|
495 ;============================ |
|
496 ; global contour : model |
|
497 ;============================ |
|
498 |
|
499 plot_name = "global_"+component(m)+"_model" |
|
500 title = "Model " + model_name |
|
501 resg@tiMainString = title |
|
502 |
|
503 wks = gsn_open_wks (plot_type,plot_name) |
|
504 gsn_define_colormap(wks,"gui_default") |
|
505 |
|
506 plot = gsn_csm_contour_map_ce(wks,data_mod,resg) |
|
507 |
|
508 delete (wks) |
|
509 delete (plot) |
|
510 |
|
511 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
512 "rm "+plot_name+"."+plot_type) |
|
513 |
|
514 ;================================ |
|
515 ; global contour: model vs ob |
|
516 ;================================ |
|
517 |
|
518 plot_name = "global_"+component(m)+"_model_vs_ob" |
|
519 |
|
520 wks = gsn_open_wks (plot_type,plot_name) |
|
521 gsn_define_colormap(wks,"gui_default") |
|
522 |
|
523 plot=new(3,graphic) ; create graphic array |
|
524 |
|
525 resg@gsnFrame = False ; Do not draw plot |
|
526 resg@gsnDraw = False ; Do not advance frame |
|
527 |
|
528 ; plot correlation coef |
|
529 |
|
530 gRes = True |
|
531 gRes@txFontHeightF = 0.02 |
|
532 gRes@txAngleF = 90 |
|
533 |
|
534 correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")" |
|
535 |
|
536 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
|
537 |
|
538 ; plot ob |
|
539 |
|
540 title = ob_name |
|
541 resg@tiMainString = title |
|
542 |
|
543 plot(0) = gsn_csm_contour_map_ce(wks,data_ob,resg) |
|
544 |
|
545 ; plot model |
|
546 |
|
547 title = "Model "+ model_name |
|
548 resg@tiMainString = title |
|
549 |
|
550 plot(1) = gsn_csm_contour_map_ce(wks,data_mod,resg) |
|
551 |
|
552 ; plot model-ob |
|
553 |
|
554 if (m .eq. 0) then |
|
555 resg@cnMinLevelValF = -2. |
|
556 resg@cnMaxLevelValF = 2. |
|
557 resg@cnLevelSpacingF = 0.4 |
|
558 end if |
|
559 |
|
560 if (m .eq. 1) then |
|
561 resg@cnMinLevelValF = -6. |
|
562 resg@cnMaxLevelValF = 6. |
|
563 resg@cnLevelSpacingF = 1. |
|
564 end if |
|
565 |
|
566 if (m .eq. 2) then |
|
567 resg@cnMinLevelValF = -6. |
|
568 resg@cnMaxLevelValF = 6. |
|
569 resg@cnLevelSpacingF = 1. |
|
570 end if |
|
571 |
|
572 zz = data_mod |
|
573 zz = data_mod - data_ob |
|
574 title = "Model_"+model_name+" - Observed" |
|
575 resg@tiMainString = title |
|
576 |
|
577 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) |
|
578 |
|
579 ; plot panel |
|
580 |
|
581 pres = True ; panel plot mods desired |
|
582 pres@gsnMaximize = True ; fill the page |
|
583 |
|
584 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
|
585 |
|
586 delete (wks) |
|
587 delete (plot) |
|
588 |
|
589 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
590 "rm "+plot_name+"."+plot_type) |
|
591 |
|
592 delete (data_ob) |
|
593 delete (data_mod) |
|
594 |
|
595 resg@gsnFrame = True ; Do advance frame |
|
596 resg@gsnDraw = True ; Do draw plot |
|
597 |
|
598 end do ; m-loop |
|
599 |
|
600 ;************************************************** |
|
601 ; html table |
|
602 ;************************************************** |
|
603 output_html = "table_model_vs_ob.html" |
|
604 |
|
605 header_text = "<H1>LAI: Model "+model_name+" vs. MODIS observations</H1>" |
|
606 |
|
607 header = (/"<HTML>" \ |
|
608 ,"<HEAD>" \ |
|
609 ,"<TITLE>CLAMP metrics</TITLE>" \ |
|
610 ,"</HEAD>" \ |
|
611 ,header_text \ |
|
612 /) |
|
613 footer = "</HTML>" |
|
614 |
|
615 table_header = (/ \ |
|
616 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \ |
|
617 ,"<tr>" \ |
|
618 ," <th bgcolor=DDDDDD rowspan=2>Biome Class</th>" \ |
|
619 ," <th bgcolor=DDDDDD colspan=3>"+component(0)+" (m2/m2)</th>" \ |
|
620 ," <th bgcolor=DDDDDD colspan=3>"+component(1)+" (m2/m2)</th>" \ |
|
621 ," <th bgcolor=DDDDDD colspan=3>"+component(2)+" (months)</th>" \ |
|
622 ," <th bgcolor=DDDDDD rowspan=2>Annual Plot</th>" \ |
|
623 ,"</tr>" \ |
|
624 ,"<tr>" \ |
|
625 ," <th bgcolor=DDDDDD >observed</th>" \ |
|
626 ," <th bgcolor=DDDDDD >"+model_name+"</th>" \ |
|
627 ," <th bgcolor=DDDDDD >M_score</th>" \ |
|
628 ," <th bgcolor=DDDDDD >observed</th>" \ |
|
629 ," <th bgcolor=DDDDDD >"+model_name+"</th>" \ |
|
630 ," <th bgcolor=DDDDDD >M_score</th>" \ |
|
631 ," <th bgcolor=DDDDDD >observed</th>" \ |
|
632 ," <th bgcolor=DDDDDD >"+model_name+"</th>" \ |
|
633 ," <th bgcolor=DDDDDD >M_score</th>" \ |
|
634 ,"</tr>" \ |
|
635 /) |
|
636 table_footer = "</table>" |
|
637 row_header = "<tr>" |
|
638 row_footer = "</tr>" |
|
639 |
|
640 lines = new(50000,string) |
|
641 nline = 0 |
|
642 |
|
643 set_line(lines,nline,header) |
|
644 set_line(lines,nline,table_header) |
|
645 ;----------------------------------------------- |
|
646 ;row of table |
|
647 |
|
648 do n = 0,nrow-1 |
|
649 set_line(lines,nline,row_header) |
|
650 |
|
651 txt1 = row_head(n) |
|
652 txt2 = text(n,0) |
|
653 txt3 = text(n,1) |
|
654 txt4 = text(n,2) |
|
655 txt5 = text(n,3) |
|
656 txt6 = text(n,4) |
|
657 txt7 = text(n,5) |
|
658 txt8 = text(n,6) |
|
659 txt9 = text(n,7) |
|
660 txt10 = text(n,8) |
|
661 txt11 = "<a href=./annual_biome_"+n+".png>model_vs_ob</a>" |
|
662 |
|
663 set_line(lines,nline,"<th>"+txt1+"</th>") |
|
664 set_line(lines,nline,"<th>"+txt2+"</th>") |
|
665 set_line(lines,nline,"<th>"+txt3+"</th>") |
|
666 set_line(lines,nline,"<th>"+txt4+"</th>") |
|
667 set_line(lines,nline,"<th>"+txt5+"</th>") |
|
668 set_line(lines,nline,"<th>"+txt6+"</th>") |
|
669 set_line(lines,nline,"<th>"+txt7+"</th>") |
|
670 set_line(lines,nline,"<th>"+txt8+"</th>") |
|
671 set_line(lines,nline,"<th>"+txt9+"</th>") |
|
672 set_line(lines,nline,"<th>"+txt10+"</th>") |
|
673 set_line(lines,nline,"<th>"+txt11+"</th>") |
|
674 |
|
675 set_line(lines,nline,row_footer) |
|
676 end do |
|
677 ;----------------------------------------------- |
|
678 set_line(lines,nline,table_footer) |
|
679 set_line(lines,nline,footer) |
|
680 |
|
681 ; Now write to an HTML file |
|
682 |
|
683 idx = ind(.not.ismissing(lines)) |
|
684 if(.not.any(ismissing(idx))) then |
|
685 asciiwrite(output_html,lines(idx)) |
|
686 else |
|
687 print ("error?") |
|
688 end if |
|
689 |
|
690 delete (yvalues) |
|
691 delete (count) |
|
692 delete (idx) |
|
693 |
|
694 ;************************************************************************ |
|
695 ; go through all ntime |
|
696 ;************************************************************************ |
|
697 |
|
698 ; for 2 data: model and observed |
|
699 data_n = 2 |
|
700 |
|
701 ; using model biome class |
|
702 |
|
703 base = ndtooned(classmod) |
|
704 |
|
705 ; output |
|
706 |
|
707 yvalues = new((/data_n,nx,ntime/),float) |
|
708 |
|
709 ;============================== |
|
710 ; put data into bins |
|
711 ;============================== |
|
712 |
|
713 ; Loop through each range, using base |
|
714 |
|
715 do i=0,nx-1 |
|
716 |
|
717 if (i.ne.(nx-1)) then |
|
718 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1))) |
|
719 else |
|
720 idx = ind(base.ge.range(i)) |
|
721 end if |
|
722 |
|
723 ; loop through each dataset |
|
724 |
|
725 do n = 0,data_n-1 |
|
726 |
|
727 do m = 0,ntime-1 |
|
728 |
|
729 if (n .eq. 0) then |
|
730 data = ndtooned(laiob (m,:,:)) |
|
731 end if |
|
732 |
|
733 if (n .eq. 1) then |
|
734 data = ndtooned(laimod(m,:,:)) |
|
735 end if |
|
736 |
|
737 ; Calculate average |
|
738 |
|
739 if (.not.any(ismissing(idx))) then |
|
740 yvalues(n,i,m) = avg(data(idx)) |
|
741 else |
|
742 yvalues(n,i,m) = yvalues@_FillValue |
|
743 end if |
|
744 |
|
745 ;############################################################# |
|
746 ; using model biome class: |
|
747 ; |
|
748 ; set the following 4 classes to _FillValue: |
|
749 ; (3)Needleleaf Deciduous Boreal Tree, |
|
750 ; (8)Broadleaf Deciduous Boreal Tree, |
|
751 ; (9)Broadleaf Evergreen Shrub, |
|
752 ; (16)Wheat |
|
753 |
|
754 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then |
|
755 yvalues(n,i,m) = yvalues@_FillValue |
|
756 end if |
|
757 ;############################################################# |
|
758 |
|
759 end do ; m-loop |
|
760 |
|
761 delete(data) |
|
762 end do ; n-loop |
|
763 |
|
764 delete(idx) |
|
765 end do ; i-loop |
|
766 |
|
767 good = ind(.not.ismissing(yvalues(0,:,0))) |
|
768 |
|
769 n_biome = dimsizes(good) |
|
770 |
|
771 ;---------------------------------------------------------------- |
|
772 ; data for tseries plot |
|
773 |
|
774 yvalues_g = new((/data_n,n_biome,ntime/),float) |
|
775 |
|
776 yvalues_g@units = "" |
|
777 |
|
778 yvalues_g = yvalues(:,good,:) |
|
779 |
|
780 delete (good) |
|
781 |
|
782 ;******************************************************************* |
|
783 ; res for line plot |
|
784 ;******************************************************************* |
|
785 ; for x-axis in xyplot |
|
786 mon = ispan(1,12,1) |
|
787 mon@long_name = "month" |
|
788 |
|
789 res = True ; plot mods desired |
|
790 res@xyLineThicknesses = (/2.0,2.0,2.0/) ; make 2nd lines thicker |
|
791 res@xyLineColors = (/"blue","red"/) ; change line color |
|
792 |
|
793 ; res@tiMainFontHeightF = 0.025 ; size of title |
|
794 |
|
795 ; Add a boxed legend using the more simple method |
|
796 res@pmLegendDisplayMode = "Always" |
|
797 ; res@pmLegendWidthF = 0.1 |
|
798 res@pmLegendWidthF = 0.08 |
|
799 res@pmLegendHeightF = 0.06 |
|
800 ; res@pmLegendOrthogonalPosF = -1.17 |
|
801 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward) |
|
802 res@pmLegendOrthogonalPosF = -0.30 ;(downward) |
|
803 |
|
804 ; res@pmLegendParallelPosF = 0.18 |
|
805 res@pmLegendParallelPosF = 0.23 ;(rightward) |
|
806 |
|
807 ; res@lgPerimOn = False |
|
808 res@lgLabelFontHeightF = 0.015 |
|
809 res@xyExplicitLegendLabels = (/"observed",model_name/) |
|
810 |
|
811 ;************************************************************ |
|
812 |
|
813 plot_data = new((/2,ntime/),float) |
|
814 plot_data@long_name = "Leaf Area Index" |
|
815 |
|
816 ;---------------------------------------------- |
|
817 ; time series plot : per biome |
|
818 |
|
819 do m = 0, n_biome-1 |
|
820 |
|
821 plot_name = "annual_biome_"+ m |
|
822 |
|
823 wks = gsn_open_wks (plot_type,plot_name) |
|
824 |
|
825 title = "LAI : "+ row_head(m) |
|
826 res@tiMainString = title |
|
827 |
|
828 plot_data(0,:) = yvalues_g(0,m,:) |
|
829 plot_data(1,:) = yvalues_g(1,m,:) |
|
830 |
|
831 plot = gsn_csm_xy(wks,mon,plot_data,res) |
|
832 |
|
833 delete (wks) |
|
834 delete (plot) |
|
835 |
|
836 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
837 "rm "+plot_name+"."+plot_type) |
|
838 end do |
|
839 |
|
840 ;-------------------------------------------- |
|
841 ; time series plot: all biome |
|
842 |
|
843 plot_name = "annual_biome_"+ n_biome |
|
844 |
|
845 wks = gsn_open_wks (plot_type,plot_name) |
|
846 |
|
847 title = "LAI : "+ row_head(n_biome) |
|
848 res@tiMainString = title |
|
849 |
|
850 do k = 0,ntime-1 |
|
851 plot_data(0,k) = avg(yvalues_g(0,:,k)) |
|
852 plot_data(1,k) = avg(yvalues_g(1,:,k)) |
|
853 end do |
|
854 |
|
855 plot = gsn_csm_xy(wks,mon,plot_data,res) |
|
856 |
|
857 delete (wks) |
|
858 delete (plot) |
|
859 |
|
860 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
861 "rm "+plot_name+"."+plot_type) |
|
862 |
|
863 delete (plot_data) |
|
864 |
|
865 ;*************************************************************************** |
|
866 ; write total score to file |
|
867 ;*************************************************************************** |
|
868 |
|
869 asciiwrite("M_save.lai", M_total) |
|
870 |
|
871 ;*************************************************************************** |
|
872 ; output plots |
|
873 ;*************************************************************************** |
|
874 output_dir = model_name+"/lai" |
|
875 |
|
876 system("mv *.png *.html " + output_dir) |
|
877 ;*************************************************************************** |
|
878 end |
|
879 |