|
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 procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \ |
|
22 ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \ |
|
23 ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \ |
|
24 ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \ |
|
25 ,dx4[1]:numeric \ |
|
26 ) |
|
27 begin |
|
28 ; Calculaee "nice" bins for binning the data in equally spaced ranges. |
|
29 ; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D |
|
30 ; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4 |
|
31 |
|
32 nbins = 15 ; Number of bins to use. |
|
33 |
|
34 nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True) |
|
35 nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1) |
|
36 range = fspan(nicevals(0),nicevals(1),nvals) |
|
37 |
|
38 ; Use this range information to grab all the values in a |
|
39 ; particular range, and then take an average. |
|
40 |
|
41 nr = dimsizes(range) |
|
42 nx = nr-1 |
|
43 ; print (nx) |
|
44 |
|
45 ; xvalues = new((/2,nx/),typeof(RAIN1_1D)) |
|
46 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2. |
|
47 dx = xvalues(0,1) - xvalues(0,0) ; range width |
|
48 dx4 = dx/4 ; 1/4 of the range |
|
49 xvalues(1,:) = xvalues(0,:) - dx/5. |
|
50 ; yvalues = new((/2,nx/),typeof(RAIN1_1D)) |
|
51 ; mn_yvalues = new((/2,nx/),typeof(RAIN1_1D)) |
|
52 ; mx_yvalues = new((/2,nx/),typeof(RAIN1_1D)) |
|
53 |
|
54 do nd=0,1 |
|
55 |
|
56 ; See if we are doing model or observational data. |
|
57 |
|
58 if(nd.eq.0) then |
|
59 data = RAIN1_1D |
|
60 npp_data = NPP1_1D |
|
61 else |
|
62 data = RAIN2_1D |
|
63 npp_data = NPP2_1D |
|
64 end if |
|
65 |
|
66 ; Loop through each range and check for values. |
|
67 |
|
68 do i=0,nr-2 |
|
69 if (i.ne.(nr-2)) then |
|
70 ; print("") |
|
71 ; print("In range ["+range(i)+","+range(i+1)+")") |
|
72 idx = ind((range(i).le.data).and.(data.lt.range(i+1))) |
|
73 else |
|
74 ; print("") |
|
75 ; print("In range ["+range(i)+",)") |
|
76 idx = ind(range(i).le.data) |
|
77 end if |
|
78 |
|
79 ; Calculate average, and get min and max. |
|
80 |
|
81 if(.not.any(ismissing(idx))) then |
|
82 yvalues(nd,i) = avg(npp_data(idx)) |
|
83 mn_yvalues(nd,i) = min(npp_data(idx)) |
|
84 mx_yvalues(nd,i) = max(npp_data(idx)) |
|
85 count = dimsizes(idx) |
|
86 else |
|
87 count = 0 |
|
88 yvalues(nd,i) = yvalues@_FillValue |
|
89 mn_yvalues(nd,i) = yvalues@_FillValue |
|
90 mx_yvalues(nd,i) = yvalues@_FillValue |
|
91 end if |
|
92 |
|
93 ; Print out information. |
|
94 ; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i)) |
|
95 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) |
|
96 |
|
97 ; Clean up for next time in loop. |
|
98 |
|
99 delete(idx) |
|
100 end do |
|
101 delete(data) |
|
102 delete(npp_data) |
|
103 end do |
|
104 end |
|
105 ;**************************************************************************** |
|
106 ; Main code. |
|
107 ;**************************************************************************** |
|
108 begin |
|
109 |
|
110 plot_type = "ps" |
|
111 plot_type_new = "png" |
|
112 |
|
113 ;----------------------------------------------------- |
|
114 ; edit table.html of current model for movel1_vs_model2 |
|
115 |
|
116 if (isvar("compare")) then |
|
117 html_name2 = compare+"/table.html" |
|
118 html_new2 = html_name2 +".new" |
|
119 end if |
|
120 ;------------------------------------------------------ |
|
121 ; edit table.html for current model |
|
122 |
|
123 html_name = model_name+"/table.html" |
|
124 html_new = html_name +".new" |
|
125 |
|
126 ;------------------------------------------------------ |
|
127 ; read model data |
|
128 |
|
129 fm = addfile (dirm+film1,"r") |
|
130 |
|
131 nppmod0 = fm->NPP |
|
132 rainmod0 = fm->RAIN |
|
133 xm = fm->lon |
|
134 ym = fm->lat |
|
135 |
|
136 delete (fm) |
|
137 |
|
138 ;---------------------------------------------------- |
|
139 ; read ob data |
|
140 |
|
141 ;(1) data at 81 sites |
|
142 |
|
143 dir81 = diro + "npp/" |
|
144 fil81 = "data.81.nc" |
|
145 f81 = addfile (dir81+fil81,"r") |
|
146 |
|
147 id81 = f81->SITE_ID |
|
148 npp81 = f81->TNPP_C |
|
149 rain81 = tofloat(f81->PREC_ANN) |
|
150 x81 = f81->LONG_DD |
|
151 y81 = f81->LAT_DD |
|
152 |
|
153 delete (f81) |
|
154 |
|
155 id81@long_name = "SITE_ID" |
|
156 |
|
157 ;change longitude from (-180,180) to (0,360) |
|
158 ;for model data interpolation |
|
159 |
|
160 do i= 0,dimsizes(x81)-1 |
|
161 if (x81(i) .lt. 0.) then |
|
162 x81(i) = x81(i)+ 360. |
|
163 end if |
|
164 end do |
|
165 |
|
166 ;------------------------------------- |
|
167 ;(2) data at 933 sites |
|
168 |
|
169 dir933 = diro + "npp/" |
|
170 fil933 = "data.933.nc" |
|
171 f933 = addfile (dir933+fil933,"r") |
|
172 |
|
173 id933 = f933->SITE_ID |
|
174 npp933 = f933->TNPP_C |
|
175 rain933 = f933->PREC |
|
176 x933 = f933->LONG_DD |
|
177 y933 = f933->LAT_DD |
|
178 |
|
179 delete (f933) |
|
180 |
|
181 id933@long_name = "SITE_ID" |
|
182 |
|
183 ;change longitude from (-180,180) to (0,360) |
|
184 ;for model data interpolation |
|
185 |
|
186 do i= 0,dimsizes(x933)-1 |
|
187 if (x933(i) .lt. 0.) then |
|
188 x933(i) = x933(i)+ 360. |
|
189 end if |
|
190 end do |
|
191 |
|
192 ;---------------------------------------- |
|
193 ;(3) global data, interpolated into model grid |
|
194 |
|
195 dirg = diro + "npp/" |
|
196 filg = "Npp_"+model_grid+"_mean.nc" |
|
197 fglobe = addfile (dirg+filg,"r") |
|
198 |
|
199 nppglobe0 = fglobe->NPP |
|
200 nppglobe = nppglobe0 |
|
201 |
|
202 delete (fglobe) |
|
203 |
|
204 ;*********************************************************************** |
|
205 ; interpolate model data into ob sites |
|
206 ;*********************************************************************** |
|
207 |
|
208 nppmod = nppmod0(0,:,:) |
|
209 rainmod = rainmod0(0,:,:) |
|
210 delete (nppmod0) |
|
211 delete (rainmod0) |
|
212 |
|
213 nppmod81 =linint2_points(xm,ym,nppmod,True,x81,y81,0) |
|
214 |
|
215 nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0) |
|
216 |
|
217 rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0) |
|
218 |
|
219 rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0) |
|
220 |
|
221 ;********************************************************** |
|
222 ; unified units |
|
223 ;********************************************************** |
|
224 ; Units for these variables are: |
|
225 ; |
|
226 ; rain81 : mm/year |
|
227 ; rainmod : mm/s |
|
228 ; npp81 : g C/m^2/year |
|
229 ; nppmod81: g C/m^2/s |
|
230 ; nppglobe: g C/m^2/year |
|
231 ; |
|
232 ; We want to convert these to "m/year" and "g C/m^2/year". |
|
233 |
|
234 nsec_per_year = 60*60*24*365 |
|
235 |
|
236 rain81 = rain81 / 1000. |
|
237 rainmod81 = (rainmod81/ 1000.) * nsec_per_year |
|
238 nppmod81 = nppmod81 * nsec_per_year |
|
239 |
|
240 rain933 = rain933 / 1000. |
|
241 rainmod933 = (rainmod933/ 1000.) * nsec_per_year |
|
242 nppmod933 = nppmod933 * nsec_per_year |
|
243 |
|
244 nppmod = nppmod * nsec_per_year |
|
245 |
|
246 npp81@units = "gC/m^2/yr" |
|
247 nppmod81@units = "gC/m^2/yr" |
|
248 npp933@units = "gC/m^2/yr" |
|
249 nppmod933@units = "gC/m^2/yr" |
|
250 nppmod@units = "gC/m^2/yr" |
|
251 nppglobe@units = "gC/m^2/yr" |
|
252 rain81@units = "m/yr" |
|
253 rainmod81@units = "m/yr" |
|
254 rain933@units = "m/yr" |
|
255 rainmod933@units = "m/yr" |
|
256 |
|
257 npp81@long_name = "Obs. NPP (gC/m2/year)" |
|
258 npp933@long_name = "Obs. NPP (gC/m2/year)" |
|
259 nppmod81@long_name = "Model NPP (gC/m2/year)" |
|
260 nppmod933@long_name = "Model NPP (gC/m2/year)" |
|
261 nppmod@long_name = "Model NPP (gC/m2/year)" |
|
262 nppglobe@long_name = "NPP (gC/m2/year)" |
|
263 rain81@long_name = "PREC (m/year)" |
|
264 rain933@long_name = "PREC (m/year)" |
|
265 rainmod81@long_name = "PREC (m/year)" |
|
266 rainmod933@long_name = "PREC (m/year)" |
|
267 |
|
268 ; change longitude from 0-360 to -180-180 |
|
269 x81 = where(x81 .gt. 180., x81 -360., x81 ) |
|
270 x933 = where(x933 .gt. 180., x933-360., x933) |
|
271 |
|
272 ;******************************************************************* |
|
273 ;(A)-1 html table of site81 -- observed |
|
274 ;******************************************************************* |
|
275 |
|
276 output_html = "table_site81_ob.html" |
|
277 |
|
278 header = (/"<HTML>" \ |
|
279 ,"<HEAD>" \ |
|
280 ,"<TITLE>CLAMP metrics</TITLE>" \ |
|
281 ,"</HEAD>" \ |
|
282 ,"<H1>EMDI Observations Class A (81 sites)</H1>" \ |
|
283 /) |
|
284 footer = "</HTML>" |
|
285 |
|
286 table_header = (/ \ |
|
287 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \ |
|
288 ,"<tr>" \ |
|
289 ," <th bgcolor=DDDDDD >Site ID</th>" \ |
|
290 ," <th bgcolor=DDDDDD >Latitude</th>" \ |
|
291 ," <th bgcolor=DDDDDD >Longitude</th>" \ |
|
292 ," <th bgcolor=DDDDDD >NPP(gC/m2/year)</th>" \ |
|
293 ," <th bgcolor=DDDDDD >PPT(m/year)</th>" \ |
|
294 ,"</tr>" \ |
|
295 /) |
|
296 table_footer = "</table>" |
|
297 row_header = "<tr>" |
|
298 row_footer = "</tr>" |
|
299 |
|
300 lines = new(50000,string) |
|
301 nline = 0 |
|
302 |
|
303 set_line(lines,nline,header) |
|
304 set_line(lines,nline,table_header) |
|
305 ;----------------------------------------------- |
|
306 ;row of table |
|
307 |
|
308 nrow = dimsizes(id81) |
|
309 do n = 0,nrow-1 |
|
310 set_line(lines,nline,row_header) |
|
311 |
|
312 txt1 = sprintf("%5.0f", id81(n)) |
|
313 txt2 = sprintf("%5.2f", y81(n)) |
|
314 txt3 = sprintf("%5.2f", x81(n)) |
|
315 txt4 = sprintf("%5.2f", npp81(n)) |
|
316 txt5 = sprintf("%5.2f", rain81(n)) |
|
317 |
|
318 set_line(lines,nline,"<th>"+txt1+"</th>") |
|
319 set_line(lines,nline,"<th>"+txt2+"</th>") |
|
320 set_line(lines,nline,"<th>"+txt3+"</th>") |
|
321 set_line(lines,nline,"<th>"+txt4+"</th>") |
|
322 set_line(lines,nline,"<th>"+txt5+"</th>") |
|
323 |
|
324 set_line(lines,nline,row_footer) |
|
325 end do |
|
326 ;----------------------------------------------- |
|
327 set_line(lines,nline,table_footer) |
|
328 set_line(lines,nline,footer) |
|
329 |
|
330 ; Now write to an HTML file. |
|
331 idx = ind(.not.ismissing(lines)) |
|
332 if(.not.any(ismissing(idx))) then |
|
333 asciiwrite(output_html,lines(idx)) |
|
334 else |
|
335 print ("error?") |
|
336 end if |
|
337 |
|
338 ;******************************************************************* |
|
339 ;(A)-2 html table of site933 -- observed |
|
340 ;******************************************************************* |
|
341 output_html = "table_site933_ob.html" |
|
342 |
|
343 header = (/"<HTML>" \ |
|
344 ,"<HEAD>" \ |
|
345 ,"<TITLE>CLAMP metrics</TITLE>" \ |
|
346 ,"</HEAD>" \ |
|
347 ,"<H1>EMDI Observations Class B (933 sites)</H1>" \ |
|
348 /) |
|
349 footer = "</HTML>" |
|
350 |
|
351 table_header = (/ \ |
|
352 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \ |
|
353 ,"<tr>" \ |
|
354 ," <th bgcolor=DDDDDD >Site ID</th>" \ |
|
355 ," <th bgcolor=DDDDDD >Latitude</th>" \ |
|
356 ," <th bgcolor=DDDDDD >Longitude</th>" \ |
|
357 ," <th bgcolor=DDDDDD >NPP(gC/m2/year)</th>" \ |
|
358 ," <th bgcolor=DDDDDD >PPT(m/year)</th>" \ |
|
359 ,"</tr>" \ |
|
360 /) |
|
361 table_footer = "</table>" |
|
362 row_header = "<tr>" |
|
363 row_footer = "</tr>" |
|
364 |
|
365 delete (lines) |
|
366 lines = new(50000,string) |
|
367 nline = 0 |
|
368 |
|
369 set_line(lines,nline,header) |
|
370 set_line(lines,nline,table_header) |
|
371 ;----------------------------------------------- |
|
372 ;row of table |
|
373 |
|
374 nrow = dimsizes(id933) |
|
375 do n = 0,nrow-1 |
|
376 set_line(lines,nline,row_header) |
|
377 |
|
378 txt1 = sprintf("%5.0f", id933(n)) |
|
379 txt2 = sprintf("%5.2f", y933(n)) |
|
380 txt3 = sprintf("%5.2f", x933(n)) |
|
381 txt4 = sprintf("%5.2f", npp933(n)) |
|
382 txt5 = sprintf("%5.2f", rain933(n)) |
|
383 |
|
384 set_line(lines,nline,"<th>"+txt1+"</th>") |
|
385 set_line(lines,nline,"<th>"+txt2+"</th>") |
|
386 set_line(lines,nline,"<th>"+txt3+"</th>") |
|
387 set_line(lines,nline,"<th>"+txt4+"</th>") |
|
388 set_line(lines,nline,"<th>"+txt5+"</th>") |
|
389 |
|
390 set_line(lines,nline,row_footer) |
|
391 end do |
|
392 ;----------------------------------------------- |
|
393 set_line(lines,nline,table_footer) |
|
394 set_line(lines,nline,footer) |
|
395 |
|
396 ; Now write to an HTML file. |
|
397 delete (idx) |
|
398 idx = ind(.not.ismissing(lines)) |
|
399 if(.not.any(ismissing(idx))) then |
|
400 asciiwrite(output_html,lines(idx)) |
|
401 else |
|
402 print ("error?") |
|
403 end if |
|
404 |
|
405 ;******************************************************************* |
|
406 ;(A)-3 html table of site81 -- model vs ob |
|
407 ;******************************************************************* |
|
408 output_html = "table_site81_model_vs_ob.html" |
|
409 |
|
410 header_text = "<H1>Model "+model_name+" vs Class A Observations (81 sites)</H1>" |
|
411 |
|
412 header = (/"<HTML>" \ |
|
413 ,"<HEAD>" \ |
|
414 ,"<TITLE>CLAMP metrics</TITLE>" \ |
|
415 ,"</HEAD>" \ |
|
416 ,header_text \ |
|
417 /) |
|
418 footer = "</HTML>" |
|
419 |
|
420 delete (table_header) |
|
421 table_header_text = " <th bgcolor=DDDDDD >"+model_name+"</th>" |
|
422 table_header = (/ \ |
|
423 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \ |
|
424 ,"<tr>" \ |
|
425 ," <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \ |
|
426 ," <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \ |
|
427 ," <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \ |
|
428 ," <th bgcolor=DDDDDD colspan=2>NPP(gC/m2.year)</th>" \ |
|
429 ," <th bgcolor=DDDDDD colspan=2>RAIN(m/year)</th>" \ |
|
430 ,"</tr>" \ |
|
431 ,"<tr>" \ |
|
432 ," <th bgcolor=DDDDDD >observed</th>" \ |
|
433 , table_header_text \ |
|
434 ," <th bgcolor=DDDDDD >observed</th>" \ |
|
435 , table_header_text \ |
|
436 ,"</tr>" \ |
|
437 /) |
|
438 table_footer = "</table>" |
|
439 row_header = "<tr>" |
|
440 row_footer = "</tr>" |
|
441 |
|
442 delete (lines) |
|
443 lines = new(50000,string) |
|
444 nline = 0 |
|
445 |
|
446 set_line(lines,nline,header) |
|
447 set_line(lines,nline,table_header) |
|
448 ;----------------------------------------------- |
|
449 ;row of table |
|
450 |
|
451 nrow = dimsizes(id81) |
|
452 do n = 0,nrow-1 |
|
453 set_line(lines,nline,row_header) |
|
454 |
|
455 txt1 = sprintf("%5.0f", id81(n)) |
|
456 txt2 = sprintf("%5.2f", y81(n)) |
|
457 txt3 = sprintf("%5.2f", x81(n)) |
|
458 txt4 = sprintf("%5.2f", npp81(n)) |
|
459 txt5 = sprintf("%5.2f", nppmod81(n)) |
|
460 txt6 = sprintf("%5.2f", rain81(n)) |
|
461 txt7 = sprintf("%5.2f", rainmod81(n)) |
|
462 |
|
463 set_line(lines,nline,"<th>"+txt1+"</th>") |
|
464 set_line(lines,nline,"<th>"+txt2+"</th>") |
|
465 set_line(lines,nline,"<th>"+txt3+"</th>") |
|
466 set_line(lines,nline,"<th>"+txt4+"</th>") |
|
467 set_line(lines,nline,"<th>"+txt5+"</th>") |
|
468 set_line(lines,nline,"<th>"+txt6+"</th>") |
|
469 set_line(lines,nline,"<th>"+txt7+"</th>") |
|
470 |
|
471 set_line(lines,nline,row_footer) |
|
472 end do |
|
473 ;----------------------------------------------- |
|
474 set_line(lines,nline,table_footer) |
|
475 set_line(lines,nline,footer) |
|
476 |
|
477 ; Now write to an HTML file. |
|
478 delete (idx) |
|
479 idx = ind(.not.ismissing(lines)) |
|
480 if(.not.any(ismissing(idx))) then |
|
481 asciiwrite(output_html,lines(idx)) |
|
482 else |
|
483 print ("error?") |
|
484 end if |
|
485 |
|
486 ;******************************************************************* |
|
487 ;(A)-4 html table of site933 -- model vs ob |
|
488 ;******************************************************************* |
|
489 output_html = "table_site933_model_vs_ob.html" |
|
490 |
|
491 header_text = "<H1>Model "+model_name+" vs Class B Observations (933 sites)</H1>" |
|
492 |
|
493 header = (/"<HTML>" \ |
|
494 ,"<HEAD>" \ |
|
495 ,"<TITLE>CLAMP metrics</TITLE>" \ |
|
496 ,"</HEAD>" \ |
|
497 ,header_text \ |
|
498 /) |
|
499 footer = "</HTML>" |
|
500 |
|
501 delete (table_header) |
|
502 table_header_text = " <th bgcolor=DDDDDD >"+model_name+"</th>" |
|
503 table_header = (/ \ |
|
504 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \ |
|
505 ,"<tr>" \ |
|
506 ," <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \ |
|
507 ," <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \ |
|
508 ," <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \ |
|
509 ," <th bgcolor=DDDDDD colspan=2>NPP(gC/m2.year)</th>" \ |
|
510 ," <th bgcolor=DDDDDD colspan=2>RAIN(m/year)</th>" \ |
|
511 ,"</tr>" \ |
|
512 ,"<tr>" \ |
|
513 ," <th bgcolor=DDDDDD >observed</th>" \ |
|
514 , table_header_text \ |
|
515 ," <th bgcolor=DDDDDD >observed</th>" \ |
|
516 , table_header_text \ |
|
517 ,"</tr>" \ |
|
518 /) |
|
519 table_footer = "</table>" |
|
520 row_header = "<tr>" |
|
521 row_footer = "</tr>" |
|
522 |
|
523 delete (lines) |
|
524 lines = new(50000,string) |
|
525 nline = 0 |
|
526 |
|
527 set_line(lines,nline,header) |
|
528 set_line(lines,nline,table_header) |
|
529 ;----------------------------------------------- |
|
530 ;row of table |
|
531 |
|
532 nrow = dimsizes(id933) |
|
533 do n = 0,nrow-1 |
|
534 set_line(lines,nline,row_header) |
|
535 |
|
536 txt1 = sprintf("%5.0f", id933(n)) |
|
537 txt2 = sprintf("%5.2f", y933(n)) |
|
538 txt3 = sprintf("%5.2f", x933(n)) |
|
539 txt4 = sprintf("%5.2f", npp933(n)) |
|
540 txt5 = sprintf("%5.2f", nppmod933(n)) |
|
541 txt6 = sprintf("%5.2f", rain933(n)) |
|
542 txt7 = sprintf("%5.2f", rainmod933(n)) |
|
543 |
|
544 set_line(lines,nline,"<th>"+txt1+"</th>") |
|
545 set_line(lines,nline,"<th>"+txt2+"</th>") |
|
546 set_line(lines,nline,"<th>"+txt3+"</th>") |
|
547 set_line(lines,nline,"<th>"+txt4+"</th>") |
|
548 set_line(lines,nline,"<th>"+txt5+"</th>") |
|
549 set_line(lines,nline,"<th>"+txt6+"</th>") |
|
550 set_line(lines,nline,"<th>"+txt7+"</th>") |
|
551 |
|
552 set_line(lines,nline,row_footer) |
|
553 end do |
|
554 ;----------------------------------------------- |
|
555 set_line(lines,nline,table_footer) |
|
556 set_line(lines,nline,footer) |
|
557 |
|
558 ; Now write to an HTML file. |
|
559 delete (idx) |
|
560 idx = ind(.not.ismissing(lines)) |
|
561 if(.not.any(ismissing(idx))) then |
|
562 asciiwrite(output_html,lines(idx)) |
|
563 else |
|
564 print ("error?") |
|
565 end if |
|
566 |
|
567 ;*************************************************************************** |
|
568 ;(A)-5 scatter plot, model vs ob 81 |
|
569 ;*************************************************************************** |
|
570 |
|
571 plot_name = "scatter_model_vs_ob_81" |
|
572 title = model_name +" vs Class A observations (81 sites)" |
|
573 |
|
574 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
575 res = True ; plot mods desired |
|
576 res@tiMainString = title ; add title |
|
577 res@xyMarkLineModes = "Markers" ; choose which have markers |
|
578 res@xyMarkers = 16 ; choose type of marker |
|
579 res@xyMarkerColor = "red" ; Marker color |
|
580 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) |
|
581 res@tmLabelAutoStride = True ; nice tick mark labels |
|
582 |
|
583 res@gsnDraw = False |
|
584 res@gsnFrame = False ; don't advance frame yet |
|
585 ;------------------------------- |
|
586 ;compute correlation coef. and M |
|
587 ccr81 = esccr(nppmod81,npp81,0) |
|
588 ;print (ccr81) |
|
589 |
|
590 score_max = 2.5 |
|
591 |
|
592 bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81))) |
|
593 M81s = (1. - (bias/dimsizes(y81)))*score_max |
|
594 M_npp_S81 = sprintf("%.2f", M81s) |
|
595 |
|
596 if (isvar("compare")) then |
|
597 system("sed -e '1,/M_npp_S81/s/M_npp_S81/"+M_npp_S81+"/' "+html_name2+" > "+html_new2+";"+ \ |
|
598 "mv -f "+html_new2+" "+html_name2) |
|
599 end if |
|
600 |
|
601 system("sed s#M_npp_S81#"+M_npp_S81+"# "+html_name+" > "+html_new+";"+ \ |
|
602 "mv -f "+html_new+" "+html_name) |
|
603 |
|
604 tRes = True |
|
605 tRes@txFontHeightF = 0.025 |
|
606 |
|
607 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")" |
|
608 |
|
609 gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes) |
|
610 ;-------------------------------- |
|
611 plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot |
|
612 ;------------------------------- |
|
613 ; add polyline |
|
614 |
|
615 dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False) |
|
616 ;------------------------------- |
|
617 draw (plot) |
|
618 delete (wks) |
|
619 delete (dum) |
|
620 |
|
621 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
622 "rm "+plot_name+"."+plot_type) |
|
623 |
|
624 ;*************************************************************************** |
|
625 ;(A)-6 scatter plot, model vs ob 933 |
|
626 ;*************************************************************************** |
|
627 |
|
628 plot_name = "scatter_model_vs_ob_933" |
|
629 title = model_name +" vs Class B Observations (933 sites)" |
|
630 |
|
631 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
632 res = True ; plot mods desired |
|
633 res@tiMainString = title ; add title |
|
634 res@xyMarkLineModes = "Markers" ; choose which have markers |
|
635 res@xyMarkers = 16 ; choose type of marker |
|
636 res@xyMarkerColor = "red" ; Marker color |
|
637 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) |
|
638 res@tmLabelAutoStride = True ; nice tick mark labels |
|
639 |
|
640 res@gsnDraw = False |
|
641 res@gsnFrame = False ; don't advance frame yet |
|
642 ;------------------------------- |
|
643 ;compute correlation coef. and M |
|
644 |
|
645 ccr933 = esccr(nppmod933,npp933,0) |
|
646 |
|
647 score_max = 2.5 |
|
648 |
|
649 bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933))) |
|
650 M933s = (1. - (bias/dimsizes(y933)))*score_max |
|
651 M_npp_S933 = sprintf("%.2f", M933s) |
|
652 |
|
653 if (isvar("compare")) then |
|
654 system("sed -e '1,/M_npp_S933/s/M_npp_S933/"+M_npp_S933+"/' "+html_name2+" > "+html_new2+";"+ \ |
|
655 "mv -f "+html_new2+" "+html_name2) |
|
656 end if |
|
657 |
|
658 system("sed s#M_npp_S933#"+M_npp_S933+"# "+html_name+" > "+html_new+";"+ \ |
|
659 "mv -f "+html_new+" "+html_name) |
|
660 |
|
661 tRes = True |
|
662 tRes@txFontHeightF = 0.025 |
|
663 |
|
664 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")" |
|
665 |
|
666 gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes) |
|
667 ;-------------------------------- |
|
668 plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot |
|
669 ;------------------------------- |
|
670 ; add polyline |
|
671 |
|
672 dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False) |
|
673 ;------------------------------- |
|
674 draw (plot) |
|
675 delete (wks) |
|
676 delete (dum) |
|
677 |
|
678 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
679 "rm "+plot_name+"."+plot_type) |
|
680 |
|
681 ;************************************************************************** |
|
682 ;(B) histogram 81 |
|
683 ;************************************************************************** |
|
684 ;get data |
|
685 RAIN1_1D = ndtooned(rain81) |
|
686 RAIN2_1D = ndtooned(rainmod81) |
|
687 NPP1_1D = ndtooned(npp81) |
|
688 NPP2_1D = ndtooned(nppmod81) |
|
689 |
|
690 ; number of bin |
|
691 nx = 8 |
|
692 |
|
693 xvalues = new((/2,nx/),float) |
|
694 yvalues = new((/2,nx/),float) |
|
695 mn_yvalues = new((/2,nx/),float) |
|
696 mx_yvalues = new((/2,nx/),float) |
|
697 dx4 = new((/1/),float) |
|
698 |
|
699 get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \ |
|
700 ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4) |
|
701 |
|
702 ;---------------------------------------- |
|
703 ;compute correlation coeff and M score |
|
704 |
|
705 u = yvalues(0,:) |
|
706 v = yvalues(1,:) |
|
707 |
|
708 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
709 uu = u(good) |
|
710 vv = v(good) |
|
711 |
|
712 ccr81h = esccr(uu,vv,0) |
|
713 |
|
714 score_max = 2.5 |
|
715 |
|
716 bias = sum(abs(vv-uu)/(abs(vv)+abs(uu))) |
|
717 M81h = (1.- (bias/dimsizes(uu)))*score_max |
|
718 M_npp_H81 = sprintf("%.2f", M81h) |
|
719 |
|
720 if (isvar("compare")) then |
|
721 system("sed -e '1,/M_npp_H81/s/M_npp_H81/"+M_npp_H81+"/' "+html_name2+" > "+html_new2+";"+ \ |
|
722 "mv -f "+html_new2+" "+html_name2) |
|
723 end if |
|
724 |
|
725 system("sed s#M_npp_H81#"+M_npp_H81+"# "+html_name+" > "+html_new+";"+ \ |
|
726 "mv -f "+html_new+" "+html_name) |
|
727 |
|
728 delete (u) |
|
729 delete (v) |
|
730 delete (uu) |
|
731 delete (vv) |
|
732 ;---------------------------------------------------------------------- |
|
733 ; histogram res |
|
734 resh = True |
|
735 resh@gsnMaximize = True |
|
736 resh@gsnDraw = False |
|
737 resh@gsnFrame = False |
|
738 resh@xyMarkLineMode = "Markers" |
|
739 resh@xyMarkerSizeF = 0.014 |
|
740 resh@xyMarker = 16 |
|
741 resh@xyMarkerColors = (/"Brown","Blue"/) |
|
742 resh@trYMinF = min(mn_yvalues) - 10. |
|
743 resh@trYMaxF = max(mx_yvalues) + 10. |
|
744 |
|
745 resh@tiYAxisString = "NPP (g C/m2/year)" |
|
746 resh@tiXAxisString = "Precipitation (m/year)" |
|
747 |
|
748 max_bar = new((/2,nx/),graphic) |
|
749 min_bar = new((/2,nx/),graphic) |
|
750 max_cap = new((/2,nx/),graphic) |
|
751 min_cap = new((/2,nx/),graphic) |
|
752 |
|
753 lnres = True |
|
754 line_colors = (/"brown","blue"/) |
|
755 |
|
756 ;************************************************************************** |
|
757 ;(B)-1 histogram plot, ob 81 site |
|
758 ;************************************************************************** |
|
759 |
|
760 plot_name = "histogram_ob_81" |
|
761 title = "Class A Observations (81 sites)" |
|
762 resh@tiMainString = title |
|
763 |
|
764 wks = gsn_open_wks (plot_type,plot_name) |
|
765 |
|
766 xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh) |
|
767 |
|
768 ;------------------------------- |
|
769 ;Attach the vertical bar and the horizontal cap line |
|
770 |
|
771 do nd=0,0 |
|
772 lnres@gsLineColor = line_colors(nd) |
|
773 do i=0,nx-1 |
|
774 |
|
775 if(.not.ismissing(mn_yvalues(nd,i)).and. \ |
|
776 .not.ismissing(mx_yvalues(nd,i))) then |
|
777 ; |
|
778 ; Attach the vertical bar, both above and below the marker. |
|
779 ; |
|
780 x1 = xvalues(nd,i) |
|
781 y1 = yvalues(nd,i) |
|
782 y2 = mn_yvalues(nd,i) |
|
783 plx = (/x1,x1/) |
|
784 ply = (/y1,y2/) |
|
785 min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres) |
|
786 |
|
787 y2 = mx_yvalues(nd,i) |
|
788 plx = (/x1,x1/) |
|
789 ply = (/y1,y2/) |
|
790 max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres) |
|
791 ; |
|
792 ; Attach the horizontal cap line, both above and below the marker. |
|
793 ; |
|
794 x1 = xvalues(nd,i) - dx4 |
|
795 x2 = xvalues(nd,i) + dx4 |
|
796 y1 = mn_yvalues(nd,i) |
|
797 plx = (/x1,x2/) |
|
798 ply = (/y1,y1/) |
|
799 min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres) |
|
800 |
|
801 y1 = mx_yvalues(nd,i) |
|
802 plx = (/x1,x2/) |
|
803 ply = (/y1,y1/) |
|
804 max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres) |
|
805 end if |
|
806 end do |
|
807 end do |
|
808 |
|
809 draw(xy) |
|
810 delete (wks) |
|
811 |
|
812 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
813 "rm "+plot_name+"."+plot_type) |
|
814 |
|
815 ;**************************************************************************** |
|
816 ;(B)-2 histogram plot, model vs ob 81 site |
|
817 ;**************************************************************************** |
|
818 |
|
819 plot_name = "histogram_model_vs_ob_81" |
|
820 title = model_name+ " vs Class A Observations (81 sites)" |
|
821 resh@tiMainString = title |
|
822 |
|
823 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
824 |
|
825 ;----------------------------- |
|
826 ; Add a boxed legend using the more simple method, which won't have |
|
827 ; vertical lines going through the markers. |
|
828 |
|
829 resh@pmLegendDisplayMode = "Always" |
|
830 ; resh@pmLegendWidthF = 0.1 |
|
831 resh@pmLegendWidthF = 0.08 |
|
832 resh@pmLegendHeightF = 0.05 |
|
833 resh@pmLegendOrthogonalPosF = -1.17 |
|
834 ; resh@pmLegendOrthogonalPosF = -1.00 ;(downward) |
|
835 ; resh@pmLegendParallelPosF = 0.18 |
|
836 resh@pmLegendParallelPosF = 0.88 ;(rightward) |
|
837 |
|
838 ; resh@lgPerimOn = False |
|
839 resh@lgLabelFontHeightF = 0.015 |
|
840 resh@xyExplicitLegendLabels = (/"observed",model_name/) |
|
841 ;----------------------------- |
|
842 tRes = True |
|
843 tRes@txFontHeightF = 0.025 |
|
844 |
|
845 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")" |
|
846 |
|
847 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) |
|
848 |
|
849 xy = gsn_csm_xy(wks,xvalues,yvalues,resh) |
|
850 ;------------------------------- |
|
851 ;Attach the vertical bar and the horizontal cap line |
|
852 |
|
853 do nd=0,1 |
|
854 lnres@gsLineColor = line_colors(nd) |
|
855 do i=0,nx-1 |
|
856 |
|
857 if(.not.ismissing(mn_yvalues(nd,i)).and. \ |
|
858 .not.ismissing(mx_yvalues(nd,i))) then |
|
859 ; |
|
860 ; Attach the vertical bar, both above and below the marker. |
|
861 ; |
|
862 x1 = xvalues(nd,i) |
|
863 y1 = yvalues(nd,i) |
|
864 y2 = mn_yvalues(nd,i) |
|
865 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
866 |
|
867 y2 = mx_yvalues(nd,i) |
|
868 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
869 ; |
|
870 ; Attach the horizontal cap line, both above and below the marker. |
|
871 ; |
|
872 x1 = xvalues(nd,i) - dx4 |
|
873 x2 = xvalues(nd,i) + dx4 |
|
874 y1 = mn_yvalues(nd,i) |
|
875 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
876 |
|
877 y1 = mx_yvalues(nd,i) |
|
878 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
879 end if |
|
880 end do |
|
881 end do |
|
882 |
|
883 draw(xy) |
|
884 delete(wks) |
|
885 |
|
886 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
887 "rm "+plot_name+"."+plot_type) |
|
888 |
|
889 delete (RAIN1_1D) |
|
890 delete (RAIN2_1D) |
|
891 delete (NPP1_1D) |
|
892 delete (NPP2_1D) |
|
893 ;delete (range) |
|
894 delete (xvalues) |
|
895 delete (yvalues) |
|
896 delete (mn_yvalues) |
|
897 delete (mx_yvalues) |
|
898 delete (good) |
|
899 delete (max_bar) |
|
900 delete (min_bar) |
|
901 delete (max_cap) |
|
902 delete (min_cap) |
|
903 |
|
904 ;************************************************************************** |
|
905 ;(B) histogram 933 |
|
906 ;************************************************************************** |
|
907 |
|
908 ; get data |
|
909 RAIN1_1D = ndtooned(rain933) |
|
910 RAIN2_1D = ndtooned(rainmod933) |
|
911 NPP1_1D = ndtooned(npp933) |
|
912 NPP2_1D = ndtooned(nppmod933) |
|
913 |
|
914 ; number of bin |
|
915 nx = 10 |
|
916 |
|
917 xvalues = new((/2,nx/),float) |
|
918 yvalues = new((/2,nx/),float) |
|
919 mn_yvalues = new((/2,nx/),float) |
|
920 mx_yvalues = new((/2,nx/),float) |
|
921 dx4 = new((/1/),float) |
|
922 |
|
923 get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \ |
|
924 ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4) |
|
925 |
|
926 ;---------------------------------------- |
|
927 ;compute correlation coeff and M score |
|
928 |
|
929 u = yvalues(0,:) |
|
930 v = yvalues(1,:) |
|
931 |
|
932 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
933 uu = u(good) |
|
934 vv = v(good) |
|
935 |
|
936 ccr933h = esccr(uu,vv,0) |
|
937 |
|
938 score_max = 2.5 |
|
939 |
|
940 bias = sum(abs(vv-uu)/(abs(vv)+abs(uu))) |
|
941 M933h = (1.- (bias/dimsizes(uu)))*score_max |
|
942 M_npp_H933 = sprintf("%.2f", M933h) |
|
943 |
|
944 if (isvar("compare")) then |
|
945 system("sed -e '1,/M_npp_H933/s/M_npp_H933/"+M_npp_H933+"/' "+html_name2+" > "+html_new2+";"+ \ |
|
946 "mv -f "+html_new2+" "+html_name2) |
|
947 end if |
|
948 |
|
949 system("sed s#M_npp_H933#"+M_npp_H933+"# "+html_name+" > "+html_new+";"+ \ |
|
950 "mv -f "+html_new+" "+html_name) |
|
951 |
|
952 delete (u) |
|
953 delete (v) |
|
954 delete (uu) |
|
955 delete (vv) |
|
956 ;---------------------------------------------------------------------- |
|
957 ; histogram res |
|
958 delete (resh) |
|
959 resh = True |
|
960 resh@gsnMaximize = True |
|
961 resh@gsnDraw = False |
|
962 resh@gsnFrame = False |
|
963 resh@xyMarkLineMode = "Markers" |
|
964 resh@xyMarkerSizeF = 0.014 |
|
965 resh@xyMarker = 16 |
|
966 resh@xyMarkerColors = (/"Brown","Blue"/) |
|
967 resh@trYMinF = min(mn_yvalues) - 10. |
|
968 resh@trYMaxF = max(mx_yvalues) + 10. |
|
969 |
|
970 resh@tiYAxisString = "NPP (g C/m2/year)" |
|
971 resh@tiXAxisString = "Precipitation (m/year)" |
|
972 |
|
973 max_bar = new((/2,nx/),graphic) |
|
974 min_bar = new((/2,nx/),graphic) |
|
975 max_cap = new((/2,nx/),graphic) |
|
976 min_cap = new((/2,nx/),graphic) |
|
977 |
|
978 lnres = True |
|
979 line_colors = (/"brown","blue"/) |
|
980 |
|
981 ;************************************************************************** |
|
982 ;(B)-3 histogram plot, ob 933 site |
|
983 ;************************************************************************** |
|
984 |
|
985 plot_name = "histogram_ob_933" |
|
986 title = "Class B Observations (933 sites)" |
|
987 resh@tiMainString = title |
|
988 |
|
989 wks = gsn_open_wks (plot_type,plot_name) |
|
990 |
|
991 xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh) |
|
992 |
|
993 ;------------------------------- |
|
994 ;Attach the vertical bar and the horizontal cap line |
|
995 |
|
996 do nd=0,0 |
|
997 lnres@gsLineColor = line_colors(nd) |
|
998 do i=0,nx-1 |
|
999 |
|
1000 if(.not.ismissing(mn_yvalues(nd,i)).and. \ |
|
1001 .not.ismissing(mx_yvalues(nd,i))) then |
|
1002 |
|
1003 ; Attach the vertical bar, both above and below the marker. |
|
1004 |
|
1005 x1 = xvalues(nd,i) |
|
1006 y1 = yvalues(nd,i) |
|
1007 y2 = mn_yvalues(nd,i) |
|
1008 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
1009 |
|
1010 y2 = mx_yvalues(nd,i) |
|
1011 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
1012 |
|
1013 ; Attach the horizontal cap line, both above and below the marker. |
|
1014 |
|
1015 x1 = xvalues(nd,i) - dx4 |
|
1016 x2 = xvalues(nd,i) + dx4 |
|
1017 y1 = mn_yvalues(nd,i) |
|
1018 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
1019 |
|
1020 y1 = mx_yvalues(nd,i) |
|
1021 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
1022 end if |
|
1023 end do |
|
1024 end do |
|
1025 |
|
1026 draw(xy) |
|
1027 delete (xy) |
|
1028 delete (wks) |
|
1029 |
|
1030 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1031 "rm "+plot_name+"."+plot_type) |
|
1032 |
|
1033 ;**************************************************************************** |
|
1034 ;(B)-4 histogram plot, model vs ob 933 site |
|
1035 ;**************************************************************************** |
|
1036 |
|
1037 plot_name = "histogram_model_vs_ob_933" |
|
1038 title = model_name+ " vs Class B Observations (933 sites)" |
|
1039 resh@tiMainString = title |
|
1040 |
|
1041 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1042 |
|
1043 ;----------------------------- |
|
1044 ; Add a boxed legend using the more simple method, which won't have |
|
1045 ; vertical lines going through the markers. |
|
1046 |
|
1047 resh@pmLegendDisplayMode = "Always" |
|
1048 ; resh@pmLegendWidthF = 0.1 |
|
1049 resh@pmLegendWidthF = 0.08 |
|
1050 resh@pmLegendHeightF = 0.05 |
|
1051 resh@pmLegendOrthogonalPosF = -1.17 |
|
1052 ; resh@pmLegendOrthogonalPosF = -1.00 ;(downward) |
|
1053 ; resh@pmLegendParallelPosF = 0.18 |
|
1054 resh@pmLegendParallelPosF = 0.88 ;(rightward) |
|
1055 |
|
1056 ; resh@lgPerimOn = False |
|
1057 resh@lgLabelFontHeightF = 0.015 |
|
1058 resh@xyExplicitLegendLabels = (/"observed",model_name/) |
|
1059 ;----------------------------- |
|
1060 tRes = True |
|
1061 tRes@txFontHeightF = 0.025 |
|
1062 |
|
1063 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")" |
|
1064 |
|
1065 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) |
|
1066 |
|
1067 xy = gsn_csm_xy(wks,xvalues,yvalues,resh) |
|
1068 ;------------------------------- |
|
1069 ;Attach the vertical bar and the horizontal cap line |
|
1070 |
|
1071 do nd=0,1 |
|
1072 lnres@gsLineColor = line_colors(nd) |
|
1073 do i=0,nx-1 |
|
1074 |
|
1075 if(.not.ismissing(mn_yvalues(nd,i)).and. \ |
|
1076 .not.ismissing(mx_yvalues(nd,i))) then |
|
1077 ; |
|
1078 ; Attach the vertical bar, both above and below the marker. |
|
1079 ; |
|
1080 x1 = xvalues(nd,i) |
|
1081 y1 = yvalues(nd,i) |
|
1082 y2 = mn_yvalues(nd,i) |
|
1083 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
1084 |
|
1085 y2 = mx_yvalues(nd,i) |
|
1086 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
1087 ; |
|
1088 ; Attach the horizontal cap line, both above and below the marker. |
|
1089 ; |
|
1090 x1 = xvalues(nd,i) - dx4 |
|
1091 x2 = xvalues(nd,i) + dx4 |
|
1092 y1 = mn_yvalues(nd,i) |
|
1093 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
1094 |
|
1095 y1 = mx_yvalues(nd,i) |
|
1096 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
1097 end if |
|
1098 end do |
|
1099 end do |
|
1100 |
|
1101 draw(xy) |
|
1102 delete(xy) |
|
1103 delete(wks) |
|
1104 |
|
1105 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1106 "rm "+plot_name+"."+plot_type) |
|
1107 |
|
1108 ;*************************************************************************** |
|
1109 ;(C) global contour |
|
1110 ;*************************************************************************** |
|
1111 |
|
1112 ;res |
|
1113 resg = True ; Use plot options |
|
1114 resg@cnFillOn = True ; Turn on color fill |
|
1115 resg@gsnSpreadColors = True ; use full colormap |
|
1116 ; resg@cnFillMode = "RasterFill" ; Turn on raster color |
|
1117 ; resg@lbLabelAutoStride = True |
|
1118 resg@cnLinesOn = False ; Turn off contourn lines |
|
1119 resg@mpFillOn = False ; Turn off map fill |
|
1120 |
|
1121 resg@gsnSpreadColors = True ; use full colormap |
|
1122 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals |
|
1123 resg@cnMinLevelValF = 0. ; Min level |
|
1124 resg@cnMaxLevelValF = 2200. ; Max level |
|
1125 resg@cnLevelSpacingF = 200. ; interval |
|
1126 |
|
1127 ;**************************************************************************** |
|
1128 ;(C)-1 global contour plot, ob |
|
1129 ;**************************************************************************** |
|
1130 |
|
1131 delta = 0.00000000001 |
|
1132 nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe) |
|
1133 |
|
1134 plot_name = "global_ob" |
|
1135 title = "Observed MODIS MOD 17" |
|
1136 resg@tiMainString = title |
|
1137 |
|
1138 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1139 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1140 |
|
1141 plot = gsn_csm_contour_map_ce(wks,nppglobe,resg) |
|
1142 |
|
1143 delete (wks) |
|
1144 |
|
1145 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1146 "rm "+plot_name+"."+plot_type) |
|
1147 ;**************************************************************************** |
|
1148 ;(C)-2 global contour plot, model |
|
1149 ;**************************************************************************** |
|
1150 |
|
1151 plot_name = "global_model" |
|
1152 title = "Model "+ model_name |
|
1153 resg@tiMainString = title |
|
1154 |
|
1155 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1156 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1157 |
|
1158 plot = gsn_csm_contour_map_ce(wks,nppmod,resg) |
|
1159 |
|
1160 delete (wks) |
|
1161 |
|
1162 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1163 "rm "+plot_name+"."+plot_type) |
|
1164 |
|
1165 ;**************************************************************************** |
|
1166 ;(C)-3 global contour plot, model vs ob |
|
1167 ;**************************************************************************** |
|
1168 |
|
1169 plot_name = "global_model_vs_ob" |
|
1170 |
|
1171 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1172 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
1173 |
|
1174 delete (plot) |
|
1175 plot=new(3,graphic) ; create graphic array |
|
1176 |
|
1177 resg@gsnFrame = False ; Do not draw plot |
|
1178 resg@gsnDraw = False ; Do not advance frame |
|
1179 |
|
1180 ; compute correlation coef and M score |
|
1181 |
|
1182 uu1 = ndtooned(nppmod) |
|
1183 vv1 = ndtooned(nppglobe) |
|
1184 |
|
1185 delete (good) |
|
1186 good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1)) |
|
1187 |
|
1188 ug = uu1(good) |
|
1189 vg = vv1(good) |
|
1190 |
|
1191 ccrG = esccr(ug,vg,0) |
|
1192 |
|
1193 score_max = 5.0 |
|
1194 |
|
1195 MG = (ccrG*ccrG)* score_max |
|
1196 M_npp_G = sprintf("%.2f", MG) |
|
1197 |
|
1198 if (isvar("compare")) then |
|
1199 system("sed -e '1,/M_npp_G/s/M_npp_G/"+M_npp_G+"/' "+html_name2+" > "+html_new2+";"+ \ |
|
1200 "mv -f "+html_new2+" "+html_name2) |
|
1201 end if |
|
1202 |
|
1203 system("sed s#M_npp_G#"+M_npp_G+"# "+html_name+" > "+html_new+";"+ \ |
|
1204 "mv -f "+html_new+" "+html_name) |
|
1205 |
|
1206 ; plot correlation coef |
|
1207 |
|
1208 gRes = True |
|
1209 gRes@txFontHeightF = 0.02 |
|
1210 gRes@txAngleF = 90 |
|
1211 |
|
1212 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")" |
|
1213 |
|
1214 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
|
1215 ;-------------------------------------------------------------------- |
|
1216 ;(a) ob |
|
1217 |
|
1218 title = "Observed MODIS MOD 17" |
|
1219 resg@tiMainString = title |
|
1220 |
|
1221 plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg) |
|
1222 |
|
1223 ;(b) model |
|
1224 |
|
1225 title = "Model "+ model_name |
|
1226 resg@tiMainString = title |
|
1227 |
|
1228 plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg) |
|
1229 |
|
1230 ;(c) model-ob |
|
1231 |
|
1232 zz = nppmod |
|
1233 zz = nppmod - nppglobe |
|
1234 title = "Model_"+model_name+" - Observed" |
|
1235 |
|
1236 resg@cnMinLevelValF = -500 ; Min level |
|
1237 resg@cnMaxLevelValF = 500. ; Max level |
|
1238 resg@cnLevelSpacingF = 50. ; interval |
|
1239 resg@tiMainString = title |
|
1240 |
|
1241 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) |
|
1242 |
|
1243 pres = True ; panel plot mods desired |
|
1244 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around |
|
1245 ; indiv. plots in panel |
|
1246 pres@gsnMaximize = True ; fill the page |
|
1247 |
|
1248 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
|
1249 |
|
1250 delete (wks) |
|
1251 delete (plot) |
|
1252 |
|
1253 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1254 "rm "+plot_name+"."+plot_type) |
|
1255 |
|
1256 ;*************************************************************************** |
|
1257 ;(D)-1 zonal line plot, ob |
|
1258 ;*************************************************************************** |
|
1259 |
|
1260 vv = zonalAve(nppglobe) |
|
1261 vv@long_name = nppglobe@long_name |
|
1262 |
|
1263 plot_name = "zonal_ob" |
|
1264 title = "Observed MODIS MOD 17" |
|
1265 |
|
1266 resz = True |
|
1267 resz@tiMainString = title |
|
1268 |
|
1269 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
1270 |
|
1271 plot = gsn_csm_xy (wks,ym,vv,resz) |
|
1272 |
|
1273 delete (wks) |
|
1274 |
|
1275 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1276 "rm "+plot_name+"."+plot_type) |
|
1277 |
|
1278 ;**************************************************************************** |
|
1279 ;(D)-2 zonal line plot, model vs ob |
|
1280 ;**************************************************************************** |
|
1281 |
|
1282 s = new ((/2,dimsizes(ym)/), float) |
|
1283 |
|
1284 s(0,:) = zonalAve(nppglobe) |
|
1285 s(1,:) = zonalAve(nppmod) |
|
1286 |
|
1287 s@long_name = nppglobe@long_name |
|
1288 ;------------------------------------------- |
|
1289 ; compute correlation coef and M score |
|
1290 |
|
1291 score_max = 5.0 |
|
1292 |
|
1293 ccrZ = esccr(s(1,:), s(0,:),0) |
|
1294 MZ = (ccrZ*ccrZ)* score_max |
|
1295 M_npp_Z = sprintf("%.2f", MZ) |
|
1296 |
|
1297 if (isvar("compare")) then |
|
1298 system("sed -e '1,/M_npp_Z/s/M_npp_Z/"+M_npp_Z+"/' "+html_name2+" > "+html_new2+";"+ \ |
|
1299 "mv -f "+html_new2+" "+html_name2) |
|
1300 end if |
|
1301 |
|
1302 system("sed s#M_npp_Z#"+M_npp_Z+"# "+html_name+" > "+html_new+";"+ \ |
|
1303 "mv -f "+html_new+" "+html_name) |
|
1304 ;------------------------------------------- |
|
1305 plot_name = "zonal_model_vs_ob" |
|
1306 title = "Zonal Average" |
|
1307 resz@tiMainString = title |
|
1308 |
|
1309 wks = gsn_open_wks (plot_type,plot_name) |
|
1310 |
|
1311 ; resz@vpHeightF = 0.4 ; change aspect ratio of plot |
|
1312 ; resz@vpWidthF = 0.7 |
|
1313 |
|
1314 resz@xyMonoLineColor = "False" ; want colored lines |
|
1315 resz@xyLineColors = (/"black","red"/) ; colors chosen |
|
1316 ; resz@xyLineThicknesses = (/3.,3./) ; line thicknesses |
|
1317 resz@xyLineThicknesses = (/2.,2./) ; line thicknesses |
|
1318 resz@xyDashPatterns = (/0.,0./) ; make all lines solid |
|
1319 |
|
1320 resz@tiYAxisString = s@long_name ; add a axis title |
|
1321 resz@txFontHeightF = 0.0195 ; change title font heights |
|
1322 |
|
1323 ; Legent |
|
1324 resz@pmLegendDisplayMode = "Always" ; turn on legend |
|
1325 resz@pmLegendSide = "Top" ; Change location of |
|
1326 ; resz@pmLegendParallelPosF = .45 ; move units right |
|
1327 resz@pmLegendParallelPosF = .82 ; move units right |
|
1328 resz@pmLegendOrthogonalPosF = -0.4 ; move units down |
|
1329 |
|
1330 resz@pmLegendWidthF = 0.10 ; Change width and |
|
1331 resz@pmLegendHeightF = 0.10 ; height of legend. |
|
1332 resz@lgLabelFontHeightF = .02 ; change font height |
|
1333 ; resz@lgTitleOn = True ; turn on legend title |
|
1334 ; resz@lgTitleString = "Example" ; create legend title |
|
1335 ; resz@lgTitleFontHeightF = .025 ; font of legend title |
|
1336 resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels |
|
1337 ;-------------------------------------------------------------------- |
|
1338 zRes = True |
|
1339 zRes@txFontHeightF = 0.025 |
|
1340 |
|
1341 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")" |
|
1342 |
|
1343 gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes) |
|
1344 ;-------------------------------------------------------------------- |
|
1345 |
|
1346 plot = gsn_csm_xy (wks,ym,s,resz) |
|
1347 |
|
1348 delete (wks) |
|
1349 |
|
1350 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
|
1351 "rm "+plot_name+"."+plot_type) |
|
1352 |
|
1353 ;*************************************************************************** |
|
1354 ; add total score and write to file |
|
1355 ;*************************************************************************** |
|
1356 M_total = M81s + M81h + M933s + M933h + MG + MZ |
|
1357 |
|
1358 asciiwrite("M_save.npp", M_total) |
|
1359 |
|
1360 ;*************************************************************************** |
|
1361 ; output plots |
|
1362 ;*************************************************************************** |
|
1363 output_dir = model_name+"/npp" |
|
1364 |
|
1365 system("mv table_*.html " + output_dir +";"+ \ |
|
1366 "mv *.png " + output_dir) |
|
1367 |
|
1368 ;*************************************************************************** |
|
1369 end |
|
1370 |