Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
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)
8 ; add line to ascci/html file
10 nnewlines = dimsizes(newlines)
11 if(nline+nnewlines-1.ge.dimsizes(lines))
12 print("set_line: bad index, not setting anything.")
15 lines(nline:nline+nnewlines-1) = newlines
16 ; print ("lines = " + lines(nline:nline+nnewlines-1))
17 nline = nline + nnewlines
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 \
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
32 nbins = 15 ; Number of bins to use.
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)
38 ; Use this range information to grab all the values in a
39 ; particular range, and then take an average.
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))
56 ; See if we are doing model or observational data.
66 ; Loop through each range and check for values.
71 ; print("In range ["+range(i)+","+range(i+1)+")")
72 idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
75 ; print("In range ["+range(i)+",)")
76 idx = ind(range(i).le.data)
79 ; Calculate average, and get min and max.
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))
88 yvalues(nd,i) = yvalues@_FillValue
89 mn_yvalues(nd,i) = yvalues@_FillValue
90 mx_yvalues(nd,i) = yvalues@_FillValue
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))
97 ; Clean up for next time in loop.
105 ;****************************************************************************
107 ;****************************************************************************
111 plot_type_new = "png"
113 ;-----------------------------------------------------
114 ; edit table.html of current model for movel1_vs_model2
116 if (isvar("compare")) then
117 html_name2 = compare+"/table.html"
118 html_new2 = html_name2 +".new"
120 ;------------------------------------------------------
121 ; edit table.html for current model
123 html_name = model_name+"/table.html"
124 html_new = html_name +".new"
126 ;------------------------------------------------------
129 fm = addfile (dirm+film1,"r")
138 ;----------------------------------------------------
141 ;(1) data at 81 sites
143 dir81 = diro + "npp/"
145 f81 = addfile (dir81+fil81,"r")
149 rain81 = tofloat(f81->PREC_ANN)
155 id81@long_name = "SITE_ID"
157 ;change longitude from (-180,180) to (0,360)
158 ;for model data interpolation
160 do i= 0,dimsizes(x81)-1
161 if (x81(i) .lt. 0.) then
162 x81(i) = x81(i)+ 360.
166 ;-------------------------------------
167 ;(2) data at 933 sites
169 dir933 = diro + "npp/"
170 fil933 = "data.933.nc"
171 f933 = addfile (dir933+fil933,"r")
173 id933 = f933->SITE_ID
174 npp933 = f933->TNPP_C
181 id933@long_name = "SITE_ID"
183 ;change longitude from (-180,180) to (0,360)
184 ;for model data interpolation
186 do i= 0,dimsizes(x933)-1
187 if (x933(i) .lt. 0.) then
188 x933(i) = x933(i)+ 360.
192 ;----------------------------------------
193 ;(3) global data, interpolated into model grid
196 ;filg = "Npp_"+model_grid+"_mean.nc"
197 filg = "npp_"+model_grid+"_mean_2000-2004.nc"
198 fglobe = addfile (dirg+filg,"r")
200 nppglobe0 = fglobe->NPP
205 ;***********************************************************************
206 ; interpolate model data into ob sites
207 ;***********************************************************************
209 nppmod = nppmod0(0,:,:)
210 rainmod = rainmod0(0,:,:)
214 nppmod81 =linint2_points(xm,ym,nppmod,True,x81,y81,0)
216 nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)
218 rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)
220 rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)
222 ;**********************************************************
224 ;**********************************************************
225 ; Units for these variables are:
229 ; npp81 : gC/m^2/year
231 ; nppglobe: gC/m^2/year
233 ; We want to convert these to "m/year" and "gC/m^2/year".
235 nsec_per_year = 60*60*24*365
237 rain81 = rain81 / 1000.
238 rainmod81 = (rainmod81/ 1000.) * nsec_per_year
239 nppmod81 = nppmod81 * nsec_per_year
241 rain933 = rain933 / 1000.
242 rainmod933 = (rainmod933/ 1000.) * nsec_per_year
243 nppmod933 = nppmod933 * nsec_per_year
245 nppmod = nppmod * nsec_per_year
247 npp81@units = "gC m~S~-2~N~ y~S~-1~N~"
248 nppmod81@units = "gC m~S~-2~N~ y~S~-1~N~"
249 npp933@units = "gC m~S~-2~N~ y~S~-1~N~"
250 nppmod933@units = "gC m~S~-2~N~ y~S~-1~N~"
251 nppmod@units = "gC m~S~-2~N~ y~S~-1~N~"
252 nppglobe@units = "gC m~S~-2~N~ y~S~-1~N~"
253 rain81@units = "m y~S~-1~N~"
254 rainmod81@units = "m y~S~-1~N~"
255 rain933@units = "m y~S~-1~N~"
256 rainmod933@units = "m y~S~-1~N~"
258 npp81@long_name = "Obs. NPP (gC m~S~-2~N~ y~S~-1~N~)"
259 npp933@long_name = "Obs. NPP (gC m~S~-2~N~ y~S~-1~N~)"
260 nppmod81@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"
261 nppmod933@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"
262 nppmod@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"
263 nppglobe@long_name = "NPP (gC m~S~-2~N~ y~S~-1~N~)"
264 rain81@long_name = "PREC (m y~S~-1~N~)"
265 rain933@long_name = "PREC (m y~S~-1~N~)"
266 rainmod81@long_name = "PREC (m y~S~-1~N~)"
267 rainmod933@long_name = "PREC (m y~S~-1~N~)"
269 ; change longitude from 0-360 to -180-180
270 x81 = where(x81 .gt. 180., x81 -360., x81 )
271 x933 = where(x933 .gt. 180., x933-360., x933)
273 ;*******************************************************************
274 ;(A)-1 html table of site81 -- observed
275 ;*******************************************************************
277 output_html = "table_site81_ob.html"
279 header = (/"<HTML>" \
281 ,"<TITLE>CLAMP metrics</TITLE>" \
283 ,"<H1>EMDI Observations Class A (81 sites)</H1>" \
288 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
290 ," <th bgcolor=DDDDDD >Site ID</th>" \
291 ," <th bgcolor=DDDDDD >Latitude</th>" \
292 ," <th bgcolor=DDDDDD >Longitude</th>" \
293 ," <th bgcolor=DDDDDD >NPP (gC m<small><sup>-2</sup></small> y<small><sup>-1</sup></small>)</th>" \
294 ," <th bgcolor=DDDDDD >PPT (m y<small><sup>-1</sup></small>)</th>" \
297 table_footer = "</table>"
301 lines = new(50000,string)
304 set_line(lines,nline,header)
305 set_line(lines,nline,table_header)
306 ;-----------------------------------------------
309 nrow = dimsizes(id81)
311 set_line(lines,nline,row_header)
313 txt1 = sprintf("%5.0f", id81(n))
314 txt2 = sprintf("%5.2f", y81(n))
315 txt3 = sprintf("%5.2f", x81(n))
316 txt4 = sprintf("%5.2f", npp81(n))
317 txt5 = sprintf("%5.2f", rain81(n))
319 set_line(lines,nline,"<th>"+txt1+"</th>")
320 set_line(lines,nline,"<th>"+txt2+"</th>")
321 set_line(lines,nline,"<th>"+txt3+"</th>")
322 set_line(lines,nline,"<th>"+txt4+"</th>")
323 set_line(lines,nline,"<th>"+txt5+"</th>")
325 set_line(lines,nline,row_footer)
327 ;-----------------------------------------------
328 set_line(lines,nline,table_footer)
329 set_line(lines,nline,footer)
331 ; Now write to an HTML file.
332 idx = ind(.not.ismissing(lines))
333 if(.not.any(ismissing(idx))) then
334 asciiwrite(output_html,lines(idx))
339 ;*******************************************************************
340 ;(A)-2 html table of site933 -- observed
341 ;*******************************************************************
342 output_html = "table_site933_ob.html"
344 header = (/"<HTML>" \
346 ,"<TITLE>CLAMP metrics</TITLE>" \
348 ,"<H1>EMDI Observations Class B (933 sites)</H1>" \
353 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
355 ," <th bgcolor=DDDDDD >Site ID</th>" \
356 ," <th bgcolor=DDDDDD >Latitude</th>" \
357 ," <th bgcolor=DDDDDD >Longitude</th>" \
358 ," <th bgcolor=DDDDDD >NPP (gC m<small><sup>-2</sup></small> y<small><sup>-1</sup></small>)</th>" \
359 ," <th bgcolor=DDDDDD >PPT (m y<small><sup>-1</sup></small>)</th>" \
362 table_footer = "</table>"
367 lines = new(50000,string)
370 set_line(lines,nline,header)
371 set_line(lines,nline,table_header)
372 ;-----------------------------------------------
375 nrow = dimsizes(id933)
377 set_line(lines,nline,row_header)
379 txt1 = sprintf("%5.0f", id933(n))
380 txt2 = sprintf("%5.2f", y933(n))
381 txt3 = sprintf("%5.2f", x933(n))
382 txt4 = sprintf("%5.2f", npp933(n))
383 txt5 = sprintf("%5.2f", rain933(n))
385 set_line(lines,nline,"<th>"+txt1+"</th>")
386 set_line(lines,nline,"<th>"+txt2+"</th>")
387 set_line(lines,nline,"<th>"+txt3+"</th>")
388 set_line(lines,nline,"<th>"+txt4+"</th>")
389 set_line(lines,nline,"<th>"+txt5+"</th>")
391 set_line(lines,nline,row_footer)
393 ;-----------------------------------------------
394 set_line(lines,nline,table_footer)
395 set_line(lines,nline,footer)
397 ; Now write to an HTML file.
399 idx = ind(.not.ismissing(lines))
400 if(.not.any(ismissing(idx))) then
401 asciiwrite(output_html,lines(idx))
406 ;*******************************************************************
407 ;(A)-3 html table of site81 -- model vs ob
408 ;*******************************************************************
409 output_html = "table_site81_model_vs_ob.html"
411 header_text = "<H1>Model "+model_name+" vs Class A Observations (81 sites)</H1>"
413 header = (/"<HTML>" \
415 ,"<TITLE>CLAMP metrics</TITLE>" \
421 delete (table_header)
422 table_header_text = " <th bgcolor=DDDDDD >"+model_name+"</th>"
424 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
426 ," <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \
427 ," <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \
428 ," <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \
429 ," <th bgcolor=DDDDDD colspan=2>NPP (gC m<small><sup>-2</sup></small> y<small><sup>-1</sup></small>)</th>" \
430 ," <th bgcolor=DDDDDD colspan=2>RAIN (m y<small><sup>-1</sup></small>)</th>" \
433 ," <th bgcolor=DDDDDD >observed</th>" \
434 , table_header_text \
435 ," <th bgcolor=DDDDDD >observed</th>" \
436 , table_header_text \
439 table_footer = "</table>"
444 lines = new(50000,string)
447 set_line(lines,nline,header)
448 set_line(lines,nline,table_header)
449 ;-----------------------------------------------
452 nrow = dimsizes(id81)
454 set_line(lines,nline,row_header)
456 txt1 = sprintf("%5.0f", id81(n))
457 txt2 = sprintf("%5.2f", y81(n))
458 txt3 = sprintf("%5.2f", x81(n))
459 txt4 = sprintf("%5.2f", npp81(n))
460 txt5 = sprintf("%5.2f", nppmod81(n))
461 txt6 = sprintf("%5.2f", rain81(n))
462 txt7 = sprintf("%5.2f", rainmod81(n))
464 set_line(lines,nline,"<th>"+txt1+"</th>")
465 set_line(lines,nline,"<th>"+txt2+"</th>")
466 set_line(lines,nline,"<th>"+txt3+"</th>")
467 set_line(lines,nline,"<th>"+txt4+"</th>")
468 set_line(lines,nline,"<th>"+txt5+"</th>")
469 set_line(lines,nline,"<th>"+txt6+"</th>")
470 set_line(lines,nline,"<th>"+txt7+"</th>")
472 set_line(lines,nline,row_footer)
474 ;-----------------------------------------------
475 set_line(lines,nline,table_footer)
476 set_line(lines,nline,footer)
478 ; Now write to an HTML file.
480 idx = ind(.not.ismissing(lines))
481 if(.not.any(ismissing(idx))) then
482 asciiwrite(output_html,lines(idx))
487 ;*******************************************************************
488 ;(A)-4 html table of site933 -- model vs ob
489 ;*******************************************************************
490 output_html = "table_site933_model_vs_ob.html"
492 header_text = "<H1>Model "+model_name+" vs Class B Observations (933 sites)</H1>"
494 header = (/"<HTML>" \
496 ,"<TITLE>CLAMP metrics</TITLE>" \
502 delete (table_header)
503 table_header_text = " <th bgcolor=DDDDDD >"+model_name+"</th>"
505 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
507 ," <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \
508 ," <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \
509 ," <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \
510 ," <th bgcolor=DDDDDD colspan=2>NPP (gC m<small><sup>-2</sup></small> y<small><sup>-1</sup></small>)</th>" \
511 ," <th bgcolor=DDDDDD colspan=2>RAIN (m y<small><sup>-1</sup></small>)</th>" \
514 ," <th bgcolor=DDDDDD >observed</th>" \
515 , table_header_text \
516 ," <th bgcolor=DDDDDD >observed</th>" \
517 , table_header_text \
520 table_footer = "</table>"
525 lines = new(50000,string)
528 set_line(lines,nline,header)
529 set_line(lines,nline,table_header)
530 ;-----------------------------------------------
533 nrow = dimsizes(id933)
535 set_line(lines,nline,row_header)
537 txt1 = sprintf("%5.0f", id933(n))
538 txt2 = sprintf("%5.2f", y933(n))
539 txt3 = sprintf("%5.2f", x933(n))
540 txt4 = sprintf("%5.2f", npp933(n))
541 txt5 = sprintf("%5.2f", nppmod933(n))
542 txt6 = sprintf("%5.2f", rain933(n))
543 txt7 = sprintf("%5.2f", rainmod933(n))
545 set_line(lines,nline,"<th>"+txt1+"</th>")
546 set_line(lines,nline,"<th>"+txt2+"</th>")
547 set_line(lines,nline,"<th>"+txt3+"</th>")
548 set_line(lines,nline,"<th>"+txt4+"</th>")
549 set_line(lines,nline,"<th>"+txt5+"</th>")
550 set_line(lines,nline,"<th>"+txt6+"</th>")
551 set_line(lines,nline,"<th>"+txt7+"</th>")
553 set_line(lines,nline,row_footer)
555 ;-----------------------------------------------
556 set_line(lines,nline,table_footer)
557 set_line(lines,nline,footer)
559 ; Now write to an HTML file.
561 idx = ind(.not.ismissing(lines))
562 if(.not.any(ismissing(idx))) then
563 asciiwrite(output_html,lines(idx))
568 ;***************************************************************************
569 ;(A)-5 scatter plot, model vs ob 81
570 ;***************************************************************************
572 plot_name = "scatter_model_vs_ob_81"
573 ;title = model_name +" vs Class A observations (81 sites)"
574 title = model_name +" (1975-2000) vs Class A observations (81 sites)"
576 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
577 res = True ; plot mods desired
578 res@tiMainString = title ; add title
579 res@xyMarkLineModes = "Markers" ; choose which have markers
580 res@xyMarkers = 16 ; choose type of marker
581 res@xyMarkerColor = "red" ; Marker color
582 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
583 res@tmLabelAutoStride = True ; nice tick mark labels
586 res@gsnFrame = False ; don't advance frame yet
587 ;-------------------------------
588 ;compute correlation coef. and M
589 ccr81 = esccr(nppmod81,npp81,0)
594 bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81)))
595 M81s = (1. - (bias/dimsizes(y81)))*score_max
596 M_npp_S81 = sprintf("%.2f", M81s)
598 if (isvar("compare")) then
599 system("sed -e '1,/M_npp_S81/s/M_npp_S81/"+M_npp_S81+"/' "+html_name2+" > "+html_new2+";"+ \
600 "mv -f "+html_new2+" "+html_name2)
603 system("sed s#M_npp_S81#"+M_npp_S81+"# "+html_name+" > "+html_new+";"+ \
604 "mv -f "+html_new+" "+html_name)
607 tRes@txFontHeightF = 0.025
609 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
611 gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
612 ;--------------------------------
613 plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot
614 ;-------------------------------
617 dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)
618 ;-------------------------------
623 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
624 "rm "+plot_name+"."+plot_type)
625 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
627 ;***************************************************************************
628 ;(A)-6 scatter plot, model vs ob 933
629 ;***************************************************************************
631 plot_name = "scatter_model_vs_ob_933"
632 title = model_name +" (1975-2000) vs Class B Observations (933 sites)"
634 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
635 res = True ; plot mods desired
636 res@tiMainString = title ; add title
637 res@xyMarkLineModes = "Markers" ; choose which have markers
638 res@xyMarkers = 16 ; choose type of marker
639 res@xyMarkerColor = "red" ; Marker color
640 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
641 res@tmLabelAutoStride = True ; nice tick mark labels
644 res@gsnFrame = False ; don't advance frame yet
645 ;-------------------------------
646 ;compute correlation coef. and M
648 ccr933 = esccr(nppmod933,npp933,0)
652 bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))
653 M933s = (1. - (bias/dimsizes(y933)))*score_max
654 M_npp_S933 = sprintf("%.2f", M933s)
656 if (isvar("compare")) then
657 system("sed -e '1,/M_npp_S933/s/M_npp_S933/"+M_npp_S933+"/' "+html_name2+" > "+html_new2+";"+ \
658 "mv -f "+html_new2+" "+html_name2)
661 system("sed s#M_npp_S933#"+M_npp_S933+"# "+html_name+" > "+html_new+";"+ \
662 "mv -f "+html_new+" "+html_name)
665 tRes@txFontHeightF = 0.025
667 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
669 gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
670 ;--------------------------------
671 plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot
672 ;-------------------------------
675 dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)
676 ;-------------------------------
681 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
682 "rm "+plot_name+"."+plot_type)
683 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
685 ;**************************************************************************
687 ;**************************************************************************
689 RAIN1_1D = ndtooned(rain81)
690 RAIN2_1D = ndtooned(rainmod81)
691 NPP1_1D = ndtooned(npp81)
692 NPP2_1D = ndtooned(nppmod81)
697 xvalues = new((/2,nx/),float)
698 yvalues = new((/2,nx/),float)
699 mn_yvalues = new((/2,nx/),float)
700 mx_yvalues = new((/2,nx/),float)
701 dx4 = new((/1/),float)
703 get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
704 ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
706 ;----------------------------------------
707 ;compute correlation coeff and M score
712 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
716 ccr81h = esccr(uu,vv,0)
720 bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
721 M81h = (1.- (bias/dimsizes(uu)))*score_max
722 M_npp_H81 = sprintf("%.2f", M81h)
724 if (isvar("compare")) then
725 system("sed -e '1,/M_npp_H81/s/M_npp_H81/"+M_npp_H81+"/' "+html_name2+" > "+html_new2+";"+ \
726 "mv -f "+html_new2+" "+html_name2)
729 system("sed s#M_npp_H81#"+M_npp_H81+"# "+html_name+" > "+html_new+";"+ \
730 "mv -f "+html_new+" "+html_name)
736 ;----------------------------------------------------------------------
739 resh@gsnMaximize = True
741 resh@gsnFrame = False
742 resh@xyMarkLineMode = "Markers"
743 resh@xyMarkerSizeF = 0.014
745 resh@xyMarkerColors = (/"Brown","Blue"/)
746 resh@trYMinF = min(mn_yvalues) - 10.
747 resh@trYMaxF = max(mx_yvalues) + 10.
749 resh@tiYAxisString = "NPP (gC m~S~-2~N~ y~S~-1~N~)"
750 resh@tiXAxisString = "Precipitation (m y~S~-1~N~)"
752 max_bar = new((/2,nx/),graphic)
753 min_bar = new((/2,nx/),graphic)
754 max_cap = new((/2,nx/),graphic)
755 min_cap = new((/2,nx/),graphic)
758 line_colors = (/"brown","blue"/)
760 ;**************************************************************************
761 ;(B)-1 histogram plot, ob 81 site
762 ;**************************************************************************
764 plot_name = "histogram_ob_81"
765 title = "Class A Observations (81 sites)"
766 resh@tiMainString = title
768 wks = gsn_open_wks (plot_type,plot_name)
770 xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
772 ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009
773 print("Observations: " + xvalues(0,:) + " " + yvalues(0,:))
775 ;-------------------------------
776 ;Attach the vertical bar and the horizontal cap line
779 lnres@gsLineColor = line_colors(nd)
782 if(.not.ismissing(mn_yvalues(nd,i)).and. \
783 .not.ismissing(mx_yvalues(nd,i))) then
785 ; Attach the vertical bar, both above and below the marker.
789 y2 = mn_yvalues(nd,i)
792 min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
794 y2 = mx_yvalues(nd,i)
797 max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
799 ; Attach the horizontal cap line, both above and below the marker.
801 x1 = xvalues(nd,i) - dx4
802 x2 = xvalues(nd,i) + dx4
803 y1 = mn_yvalues(nd,i)
806 min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
808 y1 = mx_yvalues(nd,i)
811 max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
819 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
820 "rm "+plot_name+"."+plot_type)
821 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
823 ;****************************************************************************
824 ;(B)-2 histogram plot, model vs ob 81 site
825 ;****************************************************************************
827 plot_name = "histogram_model_vs_ob_81"
828 title = model_name+ " (1975-2000) vs Class A Observations (81 sites)"
829 resh@tiMainString = title
831 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
833 ;-----------------------------
834 ; Add a boxed legend using the more simple method, which won't have
835 ; vertical lines going through the markers.
837 resh@pmLegendDisplayMode = "Always"
838 ; resh@pmLegendWidthF = 0.1
839 resh@pmLegendWidthF = 0.08
840 resh@pmLegendHeightF = 0.05
841 resh@pmLegendOrthogonalPosF = -1.17
842 ; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
843 ; resh@pmLegendParallelPosF = 0.18
844 resh@pmLegendParallelPosF = 0.88 ;(rightward)
846 ; resh@lgPerimOn = False
847 resh@lgLabelFontHeightF = 0.015
848 resh@xyExplicitLegendLabels = (/"observed",model_name/)
849 ;-----------------------------
851 tRes@txFontHeightF = 0.025
853 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"
855 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
857 xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
859 ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009
860 print(model_name + ": " + xvalues(1,:) + " " + yvalues(1,:))
861 ;print("All: " + xvalues(:,:) + " " + yvalues(:,:))
862 ;-------------------------------
863 ;Attach the vertical bar and the horizontal cap line
866 lnres@gsLineColor = line_colors(nd)
869 if(.not.ismissing(mn_yvalues(nd,i)).and. \
870 .not.ismissing(mx_yvalues(nd,i))) then
872 ; Attach the vertical bar, both above and below the marker.
876 y2 = mn_yvalues(nd,i)
877 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
879 y2 = mx_yvalues(nd,i)
880 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
882 ; Attach the horizontal cap line, both above and below the marker.
884 x1 = xvalues(nd,i) - dx4
885 x2 = xvalues(nd,i) + dx4
886 y1 = mn_yvalues(nd,i)
887 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
889 y1 = mx_yvalues(nd,i)
890 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
898 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
899 "rm "+plot_name+"."+plot_type)
900 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
917 ;**************************************************************************
919 ;**************************************************************************
922 RAIN1_1D = ndtooned(rain933)
923 RAIN2_1D = ndtooned(rainmod933)
924 NPP1_1D = ndtooned(npp933)
925 NPP2_1D = ndtooned(nppmod933)
930 xvalues = new((/2,nx/),float)
931 yvalues = new((/2,nx/),float)
932 mn_yvalues = new((/2,nx/),float)
933 mx_yvalues = new((/2,nx/),float)
934 dx4 = new((/1/),float)
936 get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
937 ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
939 ;----------------------------------------
940 ;compute correlation coeff and M score
945 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
949 ccr933h = esccr(uu,vv,0)
953 bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
954 M933h = (1.- (bias/dimsizes(uu)))*score_max
955 M_npp_H933 = sprintf("%.2f", M933h)
957 if (isvar("compare")) then
958 system("sed -e '1,/M_npp_H933/s/M_npp_H933/"+M_npp_H933+"/' "+html_name2+" > "+html_new2+";"+ \
959 "mv -f "+html_new2+" "+html_name2)
962 system("sed s#M_npp_H933#"+M_npp_H933+"# "+html_name+" > "+html_new+";"+ \
963 "mv -f "+html_new+" "+html_name)
969 ;----------------------------------------------------------------------
973 resh@gsnMaximize = True
975 resh@gsnFrame = False
976 resh@xyMarkLineMode = "Markers"
977 resh@xyMarkerSizeF = 0.014
979 resh@xyMarkerColors = (/"Brown","Blue"/)
980 resh@trYMinF = min(mn_yvalues) - 10.
981 resh@trYMaxF = max(mx_yvalues) + 10.
983 resh@tiYAxisString = "NPP (gC m~S~-2~N~ y~S~-1~N~)"
984 resh@tiXAxisString = "Precipitation (m y~S~-1~N~)"
986 max_bar = new((/2,nx/),graphic)
987 min_bar = new((/2,nx/),graphic)
988 max_cap = new((/2,nx/),graphic)
989 min_cap = new((/2,nx/),graphic)
992 line_colors = (/"brown","blue"/)
994 ;**************************************************************************
995 ;(B)-3 histogram plot, ob 933 site
996 ;**************************************************************************
998 plot_name = "histogram_ob_933"
999 title = "Class B Observations (933 sites)"
1000 resh@tiMainString = title
1002 wks = gsn_open_wks (plot_type,plot_name)
1004 xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
1006 ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009
1007 print("Observations: " + xvalues(0,:) + " " + yvalues(0,:))
1009 ;-------------------------------
1010 ;Attach the vertical bar and the horizontal cap line
1013 lnres@gsLineColor = line_colors(nd)
1016 if(.not.ismissing(mn_yvalues(nd,i)).and. \
1017 .not.ismissing(mx_yvalues(nd,i))) then
1019 ; Attach the vertical bar, both above and below the marker.
1023 y2 = mn_yvalues(nd,i)
1024 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1026 y2 = mx_yvalues(nd,i)
1027 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1029 ; Attach the horizontal cap line, both above and below the marker.
1031 x1 = xvalues(nd,i) - dx4
1032 x2 = xvalues(nd,i) + dx4
1033 y1 = mn_yvalues(nd,i)
1034 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1036 y1 = mx_yvalues(nd,i)
1037 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1046 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1047 "rm "+plot_name+"."+plot_type)
1048 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1050 ;****************************************************************************
1051 ;(B)-4 histogram plot, model vs ob 933 site
1052 ;****************************************************************************
1054 plot_name = "histogram_model_vs_ob_933"
1055 title = model_name+ " (1975-2000) vs Class B Observations (933 sites)"
1056 resh@tiMainString = title
1058 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1060 ;-----------------------------
1061 ; Add a boxed legend using the more simple method, which won't have
1062 ; vertical lines going through the markers.
1064 resh@pmLegendDisplayMode = "Always"
1065 ; resh@pmLegendWidthF = 0.1
1066 resh@pmLegendWidthF = 0.08
1067 resh@pmLegendHeightF = 0.05
1068 resh@pmLegendOrthogonalPosF = -1.17
1069 ; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
1070 ; resh@pmLegendParallelPosF = 0.18
1071 resh@pmLegendParallelPosF = 0.88 ;(rightward)
1073 ; resh@lgPerimOn = False
1074 resh@lgLabelFontHeightF = 0.015
1075 resh@xyExplicitLegendLabels = (/"observed",model_name/)
1076 ;-----------------------------
1078 tRes@txFontHeightF = 0.025
1080 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"
1082 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1084 xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
1086 ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009
1087 print("Observations: " + xvalues(1,:) + " " + yvalues(1,:))
1088 ;-------------------------------
1089 ;Attach the vertical bar and the horizontal cap line
1092 lnres@gsLineColor = line_colors(nd)
1095 if(.not.ismissing(mn_yvalues(nd,i)).and. \
1096 .not.ismissing(mx_yvalues(nd,i))) then
1098 ; Attach the vertical bar, both above and below the marker.
1102 y2 = mn_yvalues(nd,i)
1103 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1105 y2 = mx_yvalues(nd,i)
1106 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1108 ; Attach the horizontal cap line, both above and below the marker.
1110 x1 = xvalues(nd,i) - dx4
1111 x2 = xvalues(nd,i) + dx4
1112 y1 = mn_yvalues(nd,i)
1113 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1115 y1 = mx_yvalues(nd,i)
1116 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1125 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1126 "rm "+plot_name+"."+plot_type)
1127 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1129 ;***************************************************************************
1130 ; Read 2000-2004 dataset for MODIS comparison
1131 ;***************************************************************************
1133 ;------------------------------------------------------
1136 fm = addfile (dirm+film9,"r")
1144 nppmod = nppmod0(0,:,:)
1146 nppmod = nppmod * nsec_per_year
1147 nppmod@units = "gC m~S~-2~N~ y~S~-1~N~"
1148 nppmod@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)"
1150 ;***************************************************************************
1152 ;***************************************************************************
1155 resg = True ; Use plot options
1156 resg@cnFillOn = True ; Turn on color fill
1157 resg@gsnSpreadColors = True ; use full colormap
1158 ; resg@cnFillMode = "RasterFill" ; Turn on raster color
1159 ; resg@lbLabelAutoStride = True
1160 resg@cnLinesOn = False ; Turn off contourn lines
1161 resg@mpFillOn = False ; Turn off map fill
1163 resg@gsnSpreadColors = True ; use full colormap
1164 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1165 resg@cnMinLevelValF = 0. ; Min level
1166 resg@cnMaxLevelValF = 2200. ; Max level
1167 resg@cnLevelSpacingF = 200. ; interval
1169 ;****************************************************************************
1170 ;(C)-1 global contour plot, ob
1171 ;****************************************************************************
1173 delta = 0.00000000001
1174 nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)
1176 plot_name = "global_ob"
1177 title = "MODIS MOD17A3 (2000-2004)"
1178 resg@tiMainString = title
1180 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1181 gsn_define_colormap(wks,"gui_default") ; choose colormap
1183 plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1187 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1188 "rm "+plot_name+"."+plot_type)
1189 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1190 ;****************************************************************************
1191 ;(C)-2 global contour plot, model
1192 ;****************************************************************************
1194 plot_name = "global_model"
1195 title = "Model "+ model_name + " (2000-2004)"
1196 resg@tiMainString = title
1198 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1199 gsn_define_colormap(wks,"gui_default") ; choose colormap
1201 plot = gsn_csm_contour_map_ce(wks,nppmod,resg)
1205 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1206 "rm "+plot_name+"."+plot_type)
1207 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1209 ;****************************************************************************
1210 ;(C)-3 global contour plot, model vs ob
1211 ;****************************************************************************
1213 plot_name = "global_model_vs_ob"
1215 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1216 gsn_define_colormap(wks,"gui_default") ; choose colormap
1219 plot=new(3,graphic) ; create graphic array
1221 resg@gsnFrame = False ; Do not draw plot
1222 resg@gsnDraw = False ; Do not advance frame
1224 ; compute correlation coef and M score
1226 uu1 = ndtooned(nppmod)
1227 vv1 = ndtooned(nppglobe)
1230 good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))
1235 ccrG = esccr(ug,vg,0)
1239 MG = (ccrG*ccrG)* score_max
1240 M_npp_G = sprintf("%.2f", MG)
1242 if (isvar("compare")) then
1243 system("sed -e '1,/M_npp_G/s/M_npp_G/"+M_npp_G+"/' "+html_name2+" > "+html_new2+";"+ \
1244 "mv -f "+html_new2+" "+html_name2)
1247 system("sed s#M_npp_G#"+M_npp_G+"# "+html_name+" > "+html_new+";"+ \
1248 "mv -f "+html_new+" "+html_name)
1250 ; plot correlation coef
1253 gRes@txFontHeightF = 0.02
1256 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
1258 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1259 ;--------------------------------------------------------------------
1262 title = "MODIS MOD17A3 (2000-2004)"
1263 resg@tiMainString = title
1265 plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1269 title = "Model "+ model_name + " (2000-2004)"
1270 resg@tiMainString = title
1272 plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)
1277 zz = nppmod - nppglobe
1278 title = "Model "+model_name+" - MODIS MOD17A3"
1280 resg@cnMinLevelValF = -500 ; Min level
1281 resg@cnMaxLevelValF = 500. ; Max level
1282 resg@cnLevelSpacingF = 50. ; interval
1283 resg@tiMainString = title
1285 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1287 pres = True ; panel plot mods desired
1288 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1289 ; indiv. plots in panel
1290 pres@gsnMaximize = True ; fill the page
1292 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1297 system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1298 "rm "+plot_name+"."+plot_type)
1299 ;system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1301 ;***************************************************************************
1302 ;(D)-1 zonal line plot, ob
1303 ;***************************************************************************
1305 vv = zonalAve(nppglobe)
1306 vv@long_name = nppglobe@long_name
1308 plot_name = "zonal_ob"
1309 title = "MODIS MOD17A3 (2000-2004)"
1312 resz@tiMainString = title
1314 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1316 plot = gsn_csm_xy (wks,ym,vv,resz)
1320 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1321 "rm "+plot_name+"."+plot_type)
1322 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1324 ;****************************************************************************
1325 ;(D)-2 zonal line plot, model vs ob
1326 ;****************************************************************************
1328 s = new ((/2,dimsizes(ym)/), float)
1330 s(0,:) = zonalAve(nppglobe)
1331 s(1,:) = zonalAve(nppmod)
1333 s@long_name = nppglobe@long_name
1334 ;-------------------------------------------
1335 ; compute correlation coef and M score
1339 ccrZ = esccr(s(1,:), s(0,:),0)
1340 MZ = (ccrZ*ccrZ)* score_max
1341 M_npp_Z = sprintf("%.2f", MZ)
1343 if (isvar("compare")) then
1344 system("sed -e '1,/M_npp_Z/s/M_npp_Z/"+M_npp_Z+"/' "+html_name2+" > "+html_new2+";"+ \
1345 "mv -f "+html_new2+" "+html_name2)
1348 system("sed s#M_npp_Z#"+M_npp_Z+"# "+html_name+" > "+html_new+";"+ \
1349 "mv -f "+html_new+" "+html_name)
1350 ;-------------------------------------------
1351 plot_name = "zonal_model_vs_ob"
1352 title = "Zonal Average (2000-2004)"
1353 resz@tiMainString = title
1355 wks = gsn_open_wks (plot_type,plot_name)
1357 ; resz@vpHeightF = 0.4 ; change aspect ratio of plot
1358 ; resz@vpWidthF = 0.7
1360 resz@xyMonoLineColor = "False" ; want colored lines
1361 resz@xyLineColors = (/"black","red"/) ; colors chosen
1362 ; resz@xyLineThicknesses = (/3.,3./) ; line thicknesses
1363 resz@xyLineThicknesses = (/2.,2./) ; line thicknesses
1364 resz@xyDashPatterns = (/0.,0./) ; make all lines solid
1366 resz@tiYAxisString = s@long_name ; add a axis title
1367 resz@txFontHeightF = 0.0195 ; change title font heights
1370 resz@pmLegendDisplayMode = "Always" ; turn on legend
1371 resz@pmLegendSide = "Top" ; Change location of
1372 ; resz@pmLegendParallelPosF = .45 ; move units right
1373 resz@pmLegendParallelPosF = .75 ; move units right
1374 resz@pmLegendOrthogonalPosF = -0.4 ; move units down
1376 resz@pmLegendWidthF = 0.10 ; Change width and
1377 resz@pmLegendHeightF = 0.10 ; height of legend.
1378 resz@lgLabelFontHeightF = .015 ; change font height
1379 ; resz@lgTitleOn = True ; turn on legend title
1380 ; resz@lgTitleString = "Example" ; create legend title
1381 ; resz@lgTitleFontHeightF = .025 ; font of legend title
1382 resz@xyExplicitLegendLabels = (/"MODIS MOD17A3",model_name/) ; explicit labels
1383 ;--------------------------------------------------------------------
1385 zRes@txFontHeightF = 0.025
1387 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"
1389 gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)
1390 ;--------------------------------------------------------------------
1392 plot = gsn_csm_xy (wks,ym,s,resz)
1396 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1397 "rm "+plot_name+"."+plot_type)
1398 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1400 ;***************************************************************************
1401 ; add total score and write to file
1402 ;***************************************************************************
1403 M_total = M81s + M81h + M933s + M933h + MG + MZ
1405 asciiwrite("M_save.npp", M_total)
1407 ;***************************************************************************
1409 ;***************************************************************************
1410 output_dir = model_name+"/npp"
1412 system("mv table_*.html " + output_dir +";"+ \
1413 "mv *.png " + output_dir)
1415 ;***************************************************************************