Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
1 ;*****************************************************
2 ; combine scatter, histogram, global and zonal plots
3 ; compute all correlation coef and M score
4 ; add 1-to-1 line in scatter plots
5 ;*****************************************************
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
8 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
9 ;*****************************************************
10 procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \
11 ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \
12 ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \
13 ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \
17 ; Calculaee "nice" bins for binning the data in equally spaced ranges.
18 ; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D
19 ; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4
21 nbins = 15 ; Number of bins to use.
23 nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
24 nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
25 range = fspan(nicevals(0),nicevals(1),nvals)
27 ; Use this range information to grab all the values in a
28 ; particular range, and then take an average.
34 ; xvalues = new((/2,nx/),typeof(RAIN1_1D))
35 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
36 dx = xvalues(0,1) - xvalues(0,0) ; range width
37 dx4 = dx/4 ; 1/4 of the range
38 xvalues(1,:) = xvalues(0,:) - dx/5.
39 ; yvalues = new((/2,nx/),typeof(RAIN1_1D))
40 ; mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
41 ; mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
45 ; See if we are doing model or observational data.
55 ; Loop through each range and check for values.
60 ; print("In range ["+range(i)+","+range(i+1)+")")
61 idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
64 ; print("In range ["+range(i)+",)")
65 idx = ind(range(i).le.data)
68 ; Calculate average, and get min and max.
70 if(.not.any(ismissing(idx))) then
71 yvalues(nd,i) = avg(npp_data(idx))
72 mn_yvalues(nd,i) = min(npp_data(idx))
73 mx_yvalues(nd,i) = max(npp_data(idx))
77 yvalues(nd,i) = yvalues@_FillValue
78 mn_yvalues(nd,i) = yvalues@_FillValue
79 mx_yvalues(nd,i) = yvalues@_FillValue
82 ; Print out information.
83 ; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
84 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
86 ; Clean up for next time in loop.
94 ;****************************************************************************
96 ;****************************************************************************
100 plot_type_new = "png"
102 ;************************************************
104 ;************************************************
106 ;film = "i01.06cn_1798-2004_ANN_climo.nc"
110 ;film = "i01.06casa_1798-2004_ANN_climo.nc"
111 ;model_name = "06casa"
114 ;film = "i01.10casa_1948-2004_ANN_climo.nc"
115 ;model_name = "10casa"
118 film = "i01.10cn_1948-2004_ANN_climo.nc"
122 html_name = "table.html." + model_name
123 html_new = html_name +".new"
124 system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new)
125 system("mv -f "+html_new+" "+html_name)
126 ;--------------------------------------------------
127 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
128 fm = addfile (dirm+film,"r")
135 ;************************************************
137 ;************************************************
139 ;(1) data at 81 sites
140 dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
142 f81 = addfile (dir81+fil81,"r")
146 rain81 = tofloat(f81->PREC_ANN)
150 id81@long_name = "SITE_ID"
152 ;change longitude from (-180,180) to (0,360)
153 ;for model data interpolation
155 do i= 0,dimsizes(x81)-1
156 if (x81(i) .lt. 0.) then
157 x81(i) = x81(i)+ 360.
161 ;-------------------------------------
162 ;(2) data at 933 sites
163 dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
164 fil933 = "data.933.nc"
165 f933 = addfile (dir933+fil933,"r")
167 id933 = f933->SITE_ID
168 npp933 = f933->TNPP_C
173 id933@long_name = "SITE_ID"
175 ;change longitude from (-180,180) to (0,360)
176 ;for model data interpolation
178 do i= 0,dimsizes(x933)-1
179 if (x933(i) .lt. 0.) then
180 x933(i) = x933(i)+ 360.
184 ;----------------------------------------
185 ;(3) global data, interpolated into model grid
186 dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
187 filglobe = "Npp_"+model_grid+"_mean.nc"
188 fglobe = addfile (dirglobe+filglobe,"r")
190 nppglobe0 = fglobe->NPP
193 ;***********************************************************************
194 ; interpolate model data into ob sites
195 ;***********************************************************************
197 nppmod = nppmod0(0,:,:)
198 rainmod = rainmod0(0,:,:)
202 nppmod81 =linint2_points(xm,ym,nppmod,True,x81,y81,0)
204 nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)
206 rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)
208 rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)
210 ;**********************************************************
212 ;**********************************************************
213 ; Units for these variables are:
217 ; npp81 : g C/m^2/year
218 ; nppmod81: g C/m^2/s
219 ; nppglobe: g C/m^2/year
221 ; We want to convert these to "m/year" and "g C/m^2/year".
223 nsec_per_year = 60*60*24*365
225 rain81 = rain81 / 1000.
226 rainmod81 = (rainmod81/ 1000.) * nsec_per_year
227 nppmod81 = nppmod81 * nsec_per_year
229 rain933 = rain933 / 1000.
230 rainmod933 = (rainmod933/ 1000.) * nsec_per_year
231 nppmod933 = nppmod933 * nsec_per_year
233 nppmod = nppmod * nsec_per_year
235 npp81@units = "gC/m^2/yr"
236 nppmod81@units = "gC/m^2/yr"
237 npp933@units = "gC/m^2/yr"
238 nppmod933@units = "gC/m^2/yr"
239 nppmod@units = "gC/m^2/yr"
240 nppglobe@units = "gC/m^2/yr"
241 rain81@units = "m/yr"
242 rainmod81@units = "m/yr"
243 rain933@units = "m/yr"
244 rainmod933@units = "m/yr"
246 npp81@long_name = "NPP (gC/m2/year)"
247 npp933@long_name = "NPP (gC/m2/year)"
248 nppmod81@long_name = "NPP (gC/m2/year)"
249 nppmod933@long_name = "NPP (gC/m2/year)"
250 nppmod@long_name = "NPP (gC/m2/year)"
251 nppglobe@long_name = "NPP (gC/m2/year)"
252 rain81@long_name = "PREC (m/year)"
253 rain933@long_name = "PREC (m/year)"
254 rainmod81@long_name = "PREC (m/year)"
255 rainmod933@long_name = "PREC (m/year)"
257 ;*******************************************************************
258 ;(A)-1 table of site81 (1) (the first 41 sites)
259 ;*******************************************************************
264 table_header_name = "Site ID"
266 ; column (not including header column)
267 col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)
268 ncol = dimsizes(col_header)
270 ; row (not including header row)
272 row_header = new ((/nrow/),string )
273 row_header(0:nrow-1) = id81(0:nrow-1)
276 ncr1 = (/1,1/) ; 1 row, 1 column
277 x1 = (/0.005,0.15/) ; Start and end X
278 y1 = (/0.900,0.95/) ; Start and end Y
279 text1 = table_header_name
281 res1@txFontHeightF = 0.015
282 res1@gsFillColor = "CornFlowerBlue"
284 ; Column header (equally space in x2)
285 ncr2 = (/1,ncol/) ; 1 rows, ncol columns
286 x2 = (/x1(1),0.995/) ; start from end of x1
290 res2@txFontHeightF = 0.015
291 res2@gsFillColor = "Gray"
293 ; Row header (equally space in y2)
294 ncr3 = (/nrow,1/) ; nrow rows, 1 columns
296 y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
299 res3@txFontHeightF = 0.010
300 res3@gsFillColor = "Gray"
303 ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
304 x4 = x2 ; Start and end x
305 y4 = y3 ; Start and end Y
306 text4 = new((/nrow,ncol/),string)
308 color_fill4 = new((/nrow,ncol/),string)
309 color_fill4 = "white"
310 ; color_fill4(:,ncol-1) = "grey"
311 ; color_fill4(nrow-1,:) = "green"
313 res4 = True ; Set up resource list
314 ; res4@gsnDebug = True ; Useful to print NDC row,col values used.
315 res4@txFontHeightF = 0.015
316 res4@gsFillColor = color_fill4
319 ;-------------------------------------------------------------------
323 text4(n,0) = sprintf("%5.2f", y81(n))
324 text4(n,1) = sprintf("%5.2f", x81(n))
325 text4(n,2) = sprintf("%5.2f", npp81(n))
326 text4(n,3) = sprintf("%5.2f", rain81(n))
328 ;-------------------------------------------------------------------
330 plot_name = "table_site81_ob1"
332 wks = gsn_open_wks (plot_type,plot_name)
333 ;------------------------------------------
337 gRes@txFontHeightF = 0.02
340 title_text = "Observation at 81 Sites (1)"
342 gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
343 ;------------------------------------------
345 gsn_table(wks,ncr1,x1,y1,text1,res1)
346 gsn_table(wks,ncr2,x2,y2,text2,res2)
347 gsn_table(wks,ncr3,x3,y3,text3,res3)
348 gsn_table(wks,ncr4,x4,y4,text4,res4)
361 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
362 system("rm "+plot_name+"."+plot_type)
363 system("rm "+plot_name+"-1."+plot_type_new)
364 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
366 ;*******************************************************************
367 ;(A)-2 table of site81 (2) (the last 40 sites)
368 ;*******************************************************************
373 table_header_name = "Site ID"
375 ; column (not including header column)
376 col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/)
377 ncol = dimsizes(col_header)
379 ; row (not including header row)
381 row_header = new ((/nrow/),string )
382 row_header(0:nrow-1) = id81(41:41+nrow-1)
385 ncr1 = (/1,1/) ; 1 row, 1 column
386 x1 = (/0.005,0.15/) ; Start and end X
387 y1 = (/0.900,0.95/) ; Start and end Y
388 text1 = table_header_name
390 res1@txFontHeightF = 0.015
391 res1@gsFillColor = "CornFlowerBlue"
393 ; Column header (equally space in x2)
394 ncr2 = (/1,ncol/) ; 1 rows, ncol columns
395 x2 = (/x1(1),0.995/) ; start from end of x1
399 res2@txFontHeightF = 0.015
400 res2@gsFillColor = "Gray"
402 ; Row header (equally space in y2)
403 ncr3 = (/nrow,1/) ; nrow rows, 1 columns
405 y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
408 res3@txFontHeightF = 0.010
409 res3@gsFillColor = "Gray"
412 ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
413 x4 = x2 ; Start and end x
414 y4 = y3 ; Start and end Y
415 text4 = new((/nrow,ncol/),string)
417 color_fill4 = new((/nrow,ncol/),string)
418 color_fill4 = "white"
419 ; color_fill4(:,ncol-1) = "grey"
420 ; color_fill4(nrow-1,:) = "green"
422 res4 = True ; Set up resource list
423 ; res4@gsnDebug = True ; Useful to print NDC row,col values used.
424 res4@txFontHeightF = 0.015
425 res4@gsFillColor = color_fill4
428 ;-------------------------------------------------------------------
432 text4(n,0) = sprintf("%5.2f", y81(n+41))
433 text4(n,1) = sprintf("%5.2f", x81(n+41))
434 text4(n,2) = sprintf("%5.2f", npp81(n+41))
435 text4(n,3) = sprintf("%5.2f", rain81(n+41))
437 ;-------------------------------------------------------------------
439 plot_name = "table_site81_ob2"
441 wks = gsn_open_wks (plot_type,plot_name)
442 ;------------------------------------------
446 gRes@txFontHeightF = 0.02
449 title_text = "Observation at 81 Sites (2)"
451 gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
452 ;------------------------------------------
454 gsn_table(wks,ncr1,x1,y1,text1,res1)
455 gsn_table(wks,ncr2,x2,y2,text2,res2)
456 gsn_table(wks,ncr3,x3,y3,text3,res3)
457 gsn_table(wks,ncr4,x4,y4,text4,res4)
470 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
471 system("rm "+plot_name+"."+plot_type)
472 system("rm "+plot_name+"-1."+plot_type_new)
473 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
475 ;**************************************************************************
476 ;(A)-3 scatter plot, ob 933
477 ;**************************************************************************
479 plot_name = "scatter_ob_933"
480 title = "Observed 933 sites"
482 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
483 res = True ; plot mods desired
484 res@tiMainString = title ; add title
485 res@xyMarkLineModes = "Markers" ; choose which have markers
486 res@xyMarkers = 16 ; choose type of marker
487 res@xyMarkerColor = "red" ; Marker color
488 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
489 res@tmLabelAutoStride = True ; nice tick mark labels
491 plot = gsn_csm_xy (wks,id933,npp933,res) ; create plot
495 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
496 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
497 system("rm "+plot_name+"-*."+plot_type_new)
498 system("rm "+plot_name+"."+plot_type)
500 ;**************************************************************************
501 ;(A)-4 scatter plot, model 933
502 ;**************************************************************************
504 plot_name = "scatter_model_933"
505 title = "Model "+ model_name +" 933 sites"
507 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
508 res = True ; plot mods desired
509 res@tiMainString = title ; add title
510 res@xyMarkLineModes = "Markers" ; choose which have markers
511 res@xyMarkers = 16 ; choose type of marker
512 res@xyMarkerColor = "red" ; Marker color
513 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
514 res@tmLabelAutoStride = True ; nice tick mark labels
516 plot = gsn_csm_xy (wks,id933,nppmod933,res) ; create plot
520 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
521 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
522 system("rm "+plot_name+"-*."+plot_type_new)
523 system("rm "+plot_name+"."+plot_type)
525 ;***************************************************************************
526 ;(A)-5 scatter plot, model vs ob 81
527 ;***************************************************************************
529 plot_name = "scatter_model_vs_ob_81"
530 title = model_name +" vs ob 81 sites"
532 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
533 res = True ; plot mods desired
534 res@tiMainString = title ; add title
535 res@xyMarkLineModes = "Markers" ; choose which have markers
536 res@xyMarkers = 16 ; choose type of marker
537 res@xyMarkerColor = "red" ; Marker color
538 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
539 res@tmLabelAutoStride = True ; nice tick mark labels
542 res@gsnFrame = False ; don't advance frame yet
543 ;-------------------------------
544 ;compute correlation coef. and M
545 ccr81 = esccr(nppmod81,npp81,0)
549 bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81)))
550 M81 = (1. - (bias/dimsizes(y81)))*5.
552 M_npp_S81 = sprintf("%.2f", M81)
553 system("sed s#M_npp_S81#"+M_npp_S81+"# "+html_name+" > "+html_new)
554 system("mv -f "+html_new+" "+html_name)
558 tRes@txFontHeightF = 0.025
560 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
562 gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
563 ;--------------------------------
564 plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot
565 ;-------------------------------
568 dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)
569 ;-------------------------------
575 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
576 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
577 system("rm "+plot_name+"-*."+plot_type_new)
578 system("rm "+plot_name+"."+plot_type)
580 ;***************************************************************************
581 ;(A)-6 scatter plot, model vs ob 933
582 ;***************************************************************************
584 plot_name = "scatter_model_vs_ob_933"
585 title = model_name +" vs ob 933 sites"
587 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
588 res = True ; plot mods desired
589 res@tiMainString = title ; add title
590 res@xyMarkLineModes = "Markers" ; choose which have markers
591 res@xyMarkers = 16 ; choose type of marker
592 res@xyMarkerColor = "red" ; Marker color
593 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
594 res@tmLabelAutoStride = True ; nice tick mark labels
597 res@gsnFrame = False ; don't advance frame yet
598 ;-------------------------------
599 ;compute correlation coef. and M
600 ccr933 = esccr(nppmod933,npp933,0)
604 bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))
605 M933 = (1. - (bias/dimsizes(y933)))*5.
607 M_npp_S933 = sprintf("%.2f", M933)
608 system("sed s#M_npp_S933#"+M_npp_S933+"# "+html_name+" > "+html_new)
609 system("mv -f "+html_new+" "+html_name)
613 tRes@txFontHeightF = 0.025
615 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
617 gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
618 ;--------------------------------
619 plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot
620 ;-------------------------------
623 dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)
624 ;-------------------------------
630 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
631 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
632 system("rm "+plot_name+"-*."+plot_type_new)
633 system("rm "+plot_name+"."+plot_type)
635 ;*******************************************************************
636 ;(A)-7 table of site81 (1) (the first 41 sites), model vs ob
637 ;*******************************************************************
642 table_header_name = "Site ID"
644 ; row (not including header row)
646 row_header = new ((/nrow/),string )
647 row_header(0:nrow-1) = id81(0:nrow-1)
650 ncr1 = (/1,1/) ; 1 row, 1 column
651 x1 = (/0.005,0.15/) ; Start and end X
652 y1 = (/0.85 ,0.95/) ; Start and end Y
653 text1 = table_header_name
655 res1@txFontHeightF = 0.015
656 res1@gsFillColor = "CornFlowerBlue"
658 ; Column header1 (equally space in x2)
664 col_header = (/"Latitude","Longitude"/)
665 ncol = dimsizes(col_header)
667 ncr2 = (/1,ncol/) ; 1 rows, ncol columns
668 x2 = (/x1(1),0.35/) ; start from end of x1
672 res2@txFontHeightF = 0.015
673 res2@gsFillColor = "Gray"
675 ; Column header2 (equally space in x2)
677 col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)
678 col_header2 = (/"ob",model_name \
682 ncol1 = dimsizes(col_header1)
683 ncol2 = dimsizes(col_header2)
685 ncr21 = (/1,ncol1/) ; 1 rows, 4 columns
686 x21 = (/x2(1),0.995/) ; start from end of x2
687 yhalf = (y1(0)+y1(1))*0.5
688 y21 = (/yhalf,y1(1)/) ; top half of y1
691 res21@txFontHeightF = 0.015
692 res21@gsFillColor = "Gray"
694 ncr22 = (/1,ncol2/) ; 1 rows, 12 columns
695 x22 = x21 ; start from end of x1
696 y22 = (/y1(0),yhalf/) ; bottom half of y1
699 res22@txFontHeightF = 0.015
700 res22@gsFillColor = "Gray"
702 ; Row header (equally space in y2)
703 ncr3 = (/nrow,1/) ; nrow rows, 1 columns
705 y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
708 res3@txFontHeightF = 0.010
709 res3@gsFillColor = "Gray"
712 ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
713 x4 = x2 ; Start and end x
714 y4 = y3 ; Start and end Y
715 text4 = new((/nrow,ncol/),string)
717 color_fill4 = new((/nrow,ncol/),string)
718 color_fill4 = "white"
719 ; color_fill4(:,ncol-1) = "grey"
720 ; color_fill4(nrow-1,:) = "green"
722 res4 = True ; Set up resource list
723 ; res4@gsnDebug = True ; Useful to print NDC row,col values used.
724 res4@txFontHeightF = 0.015
725 res4@gsFillColor = color_fill4
730 ncr42 = (/nrow,ncol2/) ; nrow rows, ncol columns
731 x42 = x21 ; Start and end x
732 y42 = y3 ; Start and end Y
733 text42 = new((/nrow,ncol2/),string)
735 color_fill42 = new((/nrow,ncol2/),string)
736 color_fill42 = "white"
737 ; color_fill42(:,ncol-1) = "grey"
738 ; color_fill42(nrow-1,:) = "green"
740 res42 = True ; Set up resource list
741 ; res42@gsnDebug = True ; Useful to print NDC row,col values used.
742 res42@txFontHeightF = 0.015
743 res42@gsFillColor = color_fill42
745 delete (color_fill42)
746 ;-------------------------------------------------------------------
750 text4(n,0) = sprintf("%5.2f", y81(n))
751 text4(n,1) = sprintf("%5.2f", x81(n))
755 text42(n,0) = sprintf("%5.2f", npp81(n))
756 text42(n,1) = sprintf("%5.2f", nppmod81(n))
757 text42(n,2) = sprintf("%5.2f", rain81(n))
758 text42(n,3) = sprintf("%5.2f", rainmod81(n))
760 ;---------------------------------------------------------------------------
762 plot_name = "table_site81_model_vs_ob1"
764 wks = gsn_open_wks (plot_type,plot_name)
765 ;------------------------------------------
769 gRes@txFontHeightF = 0.02
772 title_text = "Model vs Observation at 81 Sites (1)"
774 gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
775 ;------------------------------------------
777 gsn_table(wks,ncr1,x1,y1,text1,res1)
778 gsn_table(wks,ncr2,x2,y2,text2,res2)
779 gsn_table(wks,ncr21,x21,y21,text21,res21)
780 gsn_table(wks,ncr22,x22,y22,text22,res22)
781 gsn_table(wks,ncr3,x3,y3,text3,res3)
782 gsn_table(wks,ncr4,x4,y4,text4,res4)
783 gsn_table(wks,ncr42,x42,y42,text42,res42)
812 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
813 system("rm "+plot_name+"."+plot_type)
814 system("rm "+plot_name+"-1."+plot_type_new)
815 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
817 ;*******************************************************************
818 ;(A)-8 table of site81 (2) (the last 40 sites), model vs ob
819 ;*******************************************************************
824 table_header_name = "Site ID"
826 ; row (not including header row)
828 row_header = new ((/nrow/),string )
829 row_header(0:nrow-1) = id81(41:41+nrow-1)
832 ncr1 = (/1,1/) ; 1 row, 1 column
833 x1 = (/0.005,0.15/) ; Start and end X
834 y1 = (/0.85 ,0.95/) ; Start and end Y
835 text1 = table_header_name
837 res1@txFontHeightF = 0.015
838 res1@gsFillColor = "CornFlowerBlue"
840 ; Column header1 (equally space in x2)
842 col_header = (/"Latitude","Longitude"/)
843 ncol = dimsizes(col_header)
845 ncr2 = (/1,ncol/) ; 1 rows, ncol columns
846 x2 = (/x1(1),0.35/) ; start from end of x1
850 res2@txFontHeightF = 0.015
851 res2@gsFillColor = "Gray"
853 ; Column header2 (equally space in x2)
855 ; col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)
856 col_header2 = (/"ob",model_name \
860 ncol1 = dimsizes(col_header1)
861 ncol2 = dimsizes(col_header2)
863 ncr21 = (/1,ncol1/) ; 1 rows, 4 columns
864 x21 = (/x2(1),0.995/) ; start from end of x2
865 yhalf = (y1(0)+y1(1))*0.5
866 y21 = (/yhalf,y1(1)/) ; top half of y1
869 res21@txFontHeightF = 0.015
870 res21@gsFillColor = "Gray"
872 ncr22 = (/1,ncol2/) ; 1 rows, 12 columns
873 x22 = x21 ; start from end of x1
874 y22 = (/y1(0),yhalf/) ; bottom half of y1
877 res22@txFontHeightF = 0.015
878 res22@gsFillColor = "Gray"
880 ; Row header (equally space in y2)
881 ncr3 = (/nrow,1/) ; nrow rows, 1 columns
883 y3 = (/1.0-table_length,y1(0)/) ; end at start of y1
886 res3@txFontHeightF = 0.010
887 res3@gsFillColor = "Gray"
890 ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns
891 x4 = x2 ; Start and end x
892 y4 = y3 ; Start and end Y
893 text4 = new((/nrow,ncol/),string)
895 color_fill4 = new((/nrow,ncol/),string)
896 color_fill4 = "white"
897 ; color_fill4(:,ncol-1) = "grey"
898 ; color_fill4(nrow-1,:) = "green"
900 res4 = True ; Set up resource list
901 ; res4@gsnDebug = True ; Useful to print NDC row,col values used.
902 res4@txFontHeightF = 0.015
903 res4@gsFillColor = color_fill4
908 ncr42 = (/nrow,ncol2/) ; nrow rows, ncol columns
909 x42 = x21 ; Start and end x
910 y42 = y3 ; Start and end Y
911 text42 = new((/nrow,ncol2/),string)
913 color_fill42 = new((/nrow,ncol2/),string)
914 color_fill42 = "white"
915 ; color_fill42(:,ncol-1) = "grey"
916 ; color_fill42(nrow-1,:) = "green"
918 res42 = True ; Set up resource list
919 ; res42@gsnDebug = True ; Useful to print NDC row,col values used.
920 res42@txFontHeightF = 0.015
921 res42@gsFillColor = color_fill42
923 delete (color_fill42)
924 ;-------------------------------------------------------------------
928 text4(n,0) = sprintf("%5.2f", y81(n+41))
929 text4(n,1) = sprintf("%5.2f", x81(n+41))
933 text42(n,0) = sprintf("%5.2f", npp81(n+41))
934 text42(n,1) = sprintf("%5.2f", nppmod81(n+41))
935 text42(n,2) = sprintf("%5.2f", rain81(n+41))
936 text42(n,3) = sprintf("%5.2f", rainmod81(n+41))
938 ;-------------------------------------------------------------------
940 plot_name = "table_site81_model_vs_ob2"
942 wks = gsn_open_wks (plot_type,plot_name)
943 ;------------------------------------------
947 gRes@txFontHeightF = 0.02
950 title_text = "Model vs Observation at 81 Sites (2)"
952 gsn_text_ndc(wks,title_text,0.50,0.975,gRes)
953 ;------------------------------------------
955 gsn_table(wks,ncr1,x1,y1,text1,res1)
956 gsn_table(wks,ncr2,x2,y2,text2,res2)
957 gsn_table(wks,ncr21,x21,y21,text21,res21)
958 gsn_table(wks,ncr22,x22,y22,text22,res22)
959 gsn_table(wks,ncr3,x3,y3,text3,res3)
960 gsn_table(wks,ncr4,x4,y4,text4,res4)
961 gsn_table(wks,ncr42,x42,y42,text42,res42)
990 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
991 system("rm "+plot_name+"."+plot_type)
992 system("rm "+plot_name+"-1."+plot_type_new)
993 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
995 ;**************************************************************************
997 ;**************************************************************************
999 RAIN1_1D = ndtooned(rain81)
1000 RAIN2_1D = ndtooned(rainmod81)
1001 NPP1_1D = ndtooned(npp81)
1002 NPP2_1D = ndtooned(nppmod81)
1007 xvalues = new((/2,nx/),float)
1008 yvalues = new((/2,nx/),float)
1009 mn_yvalues = new((/2,nx/),float)
1010 mx_yvalues = new((/2,nx/),float)
1011 dx4 = new((/1/),float)
1013 get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
1014 ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
1016 ;----------------------------------------
1017 ;compute correlation coeff and M score
1022 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1026 ccr81h = esccr(uu,vv,0)
1030 bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
1031 M81h = (1.- (bias/dimsizes(uu)))*5.
1033 M_npp_H81 = sprintf("%.2f", M81h)
1034 system("sed s#M_npp_H81#"+M_npp_H81+"# "+html_name+" > "+html_new)
1035 system("mv -f "+html_new+" "+html_name)
1042 ;----------------------------------------------------------------------
1045 resh@gsnMaximize = True
1046 resh@gsnDraw = False
1047 resh@gsnFrame = False
1048 resh@xyMarkLineMode = "Markers"
1049 resh@xyMarkerSizeF = 0.014
1051 resh@xyMarkerColors = (/"Brown","Blue"/)
1052 resh@trYMinF = min(mn_yvalues) - 10.
1053 resh@trYMaxF = max(mx_yvalues) + 10.
1055 resh@tiYAxisString = "NPP (g C/m2/year)"
1056 resh@tiXAxisString = "Precipitation (m/year)"
1058 max_bar = new((/2,nx/),graphic)
1059 min_bar = new((/2,nx/),graphic)
1060 max_cap = new((/2,nx/),graphic)
1061 min_cap = new((/2,nx/),graphic)
1064 line_colors = (/"brown","blue"/)
1066 ;**************************************************************************
1067 ;(B)-1 histogram plot, ob 81 site
1068 ;**************************************************************************
1070 plot_name = "histogram_ob_81"
1071 title = "Observed 81 site"
1072 resh@tiMainString = title
1074 wks = gsn_open_wks (plot_type,plot_name)
1076 xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
1078 ;-------------------------------
1079 ;Attach the vertical bar and the horizontal cap line
1087 lnres@gsLineColor = line_colors(nd)
1090 if(.not.ismissing(mn_yvalues(nd,i)).and. \
1091 .not.ismissing(mx_yvalues(nd,i))) then
1093 ; Attach the vertical bar, both above and below the marker.
1097 y2 = mn_yvalues(nd,i)
1100 min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1102 y2 = mx_yvalues(nd,i)
1105 max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1107 ; Attach the horizontal cap line, both above and below the marker.
1109 x1 = xvalues(nd,i) - dx4
1110 x2 = xvalues(nd,i) + dx4
1111 y1 = mn_yvalues(nd,i)
1114 min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1116 y1 = mx_yvalues(nd,i)
1119 max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1128 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1129 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1130 system("rm "+plot_name+"-*."+plot_type_new)
1131 system("rm "+plot_name+"."+plot_type)
1133 ;****************************************************************************
1134 ;(B)-2 histogram plot, model vs ob 81 site
1135 ;****************************************************************************
1137 plot_name = "histogram_model_vs_ob_81"
1138 title = model_name+ " vs Observed 81 site"
1139 resh@tiMainString = title
1141 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1143 ;-----------------------------
1144 ; Add a boxed legend using the more simple method, which won't have
1145 ; vertical lines going through the markers.
1147 resh@pmLegendDisplayMode = "Always"
1148 ; resh@pmLegendWidthF = 0.1
1149 resh@pmLegendWidthF = 0.08
1150 resh@pmLegendHeightF = 0.05
1151 resh@pmLegendOrthogonalPosF = -1.17
1152 ; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
1153 ; resh@pmLegendParallelPosF = 0.18
1154 resh@pmLegendParallelPosF = 0.88 ;(rightward)
1156 ; resh@lgPerimOn = False
1157 resh@lgLabelFontHeightF = 0.015
1158 resh@xyExplicitLegendLabels = (/"observed",model_name/)
1159 ;-----------------------------
1161 tRes@txFontHeightF = 0.025
1163 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"
1165 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1167 xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
1168 ;-------------------------------
1169 ;Attach the vertical bar and the horizontal cap line
1172 lnres@gsLineColor = line_colors(nd)
1175 if(.not.ismissing(mn_yvalues(nd,i)).and. \
1176 .not.ismissing(mx_yvalues(nd,i))) then
1178 ; Attach the vertical bar, both above and below the marker.
1182 y2 = mn_yvalues(nd,i)
1183 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1185 y2 = mx_yvalues(nd,i)
1186 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1188 ; Attach the horizontal cap line, both above and below the marker.
1190 x1 = xvalues(nd,i) - dx4
1191 x2 = xvalues(nd,i) + dx4
1192 y1 = mn_yvalues(nd,i)
1193 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1195 y1 = mx_yvalues(nd,i)
1196 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1205 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1206 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1207 system("rm "+plot_name+"-*."+plot_type_new)
1208 system("rm "+plot_name+"."+plot_type)
1225 ;**************************************************************************
1227 ;**************************************************************************
1230 RAIN1_1D = ndtooned(rain933)
1231 RAIN2_1D = ndtooned(rainmod933)
1232 NPP1_1D = ndtooned(npp933)
1233 NPP2_1D = ndtooned(nppmod933)
1238 xvalues = new((/2,nx/),float)
1239 yvalues = new((/2,nx/),float)
1240 mn_yvalues = new((/2,nx/),float)
1241 mx_yvalues = new((/2,nx/),float)
1242 dx4 = new((/1/),float)
1244 get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
1245 ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
1247 ;----------------------------------------
1248 ;compute correlation coeff and M score
1253 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1257 ccr933h = esccr(uu,vv,0)
1261 bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
1262 M933h = (1.- (bias/dimsizes(uu)))*5.
1264 M_npp_H933 = sprintf("%.2f", M933h)
1265 system("sed s#M_npp_H933#"+M_npp_H933+"# "+html_name+" > "+html_new)
1266 system("mv -f "+html_new+" "+html_name)
1273 ;----------------------------------------------------------------------
1277 resh@gsnMaximize = True
1278 resh@gsnDraw = False
1279 resh@gsnFrame = False
1280 resh@xyMarkLineMode = "Markers"
1281 resh@xyMarkerSizeF = 0.014
1283 resh@xyMarkerColors = (/"Brown","Blue"/)
1284 resh@trYMinF = min(mn_yvalues) - 10.
1285 resh@trYMaxF = max(mx_yvalues) + 10.
1287 resh@tiYAxisString = "NPP (g C/m2/year)"
1288 resh@tiXAxisString = "Precipitation (m/year)"
1290 max_bar = new((/2,nx/),graphic)
1291 min_bar = new((/2,nx/),graphic)
1292 max_cap = new((/2,nx/),graphic)
1293 min_cap = new((/2,nx/),graphic)
1296 line_colors = (/"brown","blue"/)
1298 ;**************************************************************************
1299 ;(B)-3 histogram plot, ob 933 site
1300 ;**************************************************************************
1302 plot_name = "histogram_ob_933"
1303 title = "Observed 933 site"
1304 resh@tiMainString = title
1306 wks = gsn_open_wks (plot_type,plot_name)
1308 xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
1310 ;-------------------------------
1311 ;Attach the vertical bar and the horizontal cap line
1314 lnres@gsLineColor = line_colors(nd)
1317 if(.not.ismissing(mn_yvalues(nd,i)).and. \
1318 .not.ismissing(mx_yvalues(nd,i))) then
1320 ; Attach the vertical bar, both above and below the marker.
1324 y2 = mn_yvalues(nd,i)
1325 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1327 y2 = mx_yvalues(nd,i)
1328 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1330 ; Attach the horizontal cap line, both above and below the marker.
1332 x1 = xvalues(nd,i) - dx4
1333 x2 = xvalues(nd,i) + dx4
1334 y1 = mn_yvalues(nd,i)
1335 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1337 y1 = mx_yvalues(nd,i)
1338 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1348 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1349 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1350 system("rm "+plot_name+"-*."+plot_type_new)
1351 system("rm "+plot_name+"."+plot_type)
1353 ;****************************************************************************
1354 ;(B)-4 histogram plot, model vs ob 933 site
1355 ;****************************************************************************
1357 plot_name = "histogram_model_vs_ob_933"
1358 title = model_name+ " vs Observed 933 site"
1359 resh@tiMainString = title
1361 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1363 ;-----------------------------
1364 ; Add a boxed legend using the more simple method, which won't have
1365 ; vertical lines going through the markers.
1367 resh@pmLegendDisplayMode = "Always"
1368 ; resh@pmLegendWidthF = 0.1
1369 resh@pmLegendWidthF = 0.08
1370 resh@pmLegendHeightF = 0.05
1371 resh@pmLegendOrthogonalPosF = -1.17
1372 ; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
1373 ; resh@pmLegendParallelPosF = 0.18
1374 resh@pmLegendParallelPosF = 0.88 ;(rightward)
1376 ; resh@lgPerimOn = False
1377 resh@lgLabelFontHeightF = 0.015
1378 resh@xyExplicitLegendLabels = (/"observed",model_name/)
1379 ;-----------------------------
1381 tRes@txFontHeightF = 0.025
1383 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"
1385 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1387 xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
1388 ;-------------------------------
1389 ;Attach the vertical bar and the horizontal cap line
1392 lnres@gsLineColor = line_colors(nd)
1395 if(.not.ismissing(mn_yvalues(nd,i)).and. \
1396 .not.ismissing(mx_yvalues(nd,i))) then
1398 ; Attach the vertical bar, both above and below the marker.
1402 y2 = mn_yvalues(nd,i)
1403 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1405 y2 = mx_yvalues(nd,i)
1406 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1408 ; Attach the horizontal cap line, both above and below the marker.
1410 x1 = xvalues(nd,i) - dx4
1411 x2 = xvalues(nd,i) + dx4
1412 y1 = mn_yvalues(nd,i)
1413 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1415 y1 = mx_yvalues(nd,i)
1416 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1426 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1427 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1428 system("rm "+plot_name+"-*."+plot_type_new)
1429 system("rm "+plot_name+"."+plot_type)
1431 ;***************************************************************************
1433 ;***************************************************************************
1436 resg = True ; Use plot options
1437 resg@cnFillOn = True ; Turn on color fill
1438 resg@gsnSpreadColors = True ; use full colormap
1439 ; resg@cnFillMode = "RasterFill" ; Turn on raster color
1440 ; resg@lbLabelAutoStride = True
1441 resg@cnLinesOn = False ; Turn off contourn lines
1442 resg@mpFillOn = False ; Turn off map fill
1444 resg@gsnSpreadColors = True ; use full colormap
1445 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1446 resg@cnMinLevelValF = 0. ; Min level
1447 resg@cnMaxLevelValF = 2200. ; Max level
1448 resg@cnLevelSpacingF = 200. ; interval
1450 ;****************************************************************************
1451 ;(C)-1 global contour plot, ob
1452 ;****************************************************************************
1454 delta = 0.00000000001
1455 nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)
1457 plot_name = "global_ob"
1458 title = "Observed MODIS MOD 17"
1459 resg@tiMainString = title
1461 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1462 gsn_define_colormap(wks,"gui_default") ; choose colormap
1464 plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1469 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1470 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1471 system("rm "+plot_name+"-*."+plot_type_new)
1472 system("rm "+plot_name+"."+plot_type)
1474 ;****************************************************************************
1475 ;(C)-2 global contour plot, model
1476 ;****************************************************************************
1478 plot_name = "global_model"
1479 title = "Model "+ model_name
1480 resg@tiMainString = title
1482 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1483 gsn_define_colormap(wks,"gui_default") ; choose colormap
1485 plot = gsn_csm_contour_map_ce(wks,nppmod,resg)
1490 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1491 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1492 system("rm "+plot_name+"-*."+plot_type_new)
1493 system("rm "+plot_name+"."+plot_type)
1495 ;****************************************************************************
1496 ;(C)-3 global contour plot, model vs ob
1497 ;****************************************************************************
1499 plot_name = "global_model_vs_ob"
1501 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1502 gsn_define_colormap(wks,"gui_default") ; choose colormap
1505 plot=new(3,graphic) ; create graphic array
1507 resg@gsnFrame = False ; Do not draw plot
1508 resg@gsnDraw = False ; Do not advance frame
1510 ; compute correlation coef and M score
1512 uu1 = ndtooned(nppmod)
1513 vv1 = ndtooned(nppglobe)
1516 good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))
1521 ccrG = esccr(ug,vg,0)
1524 MG = (ccrG*ccrG)* 5.0
1526 M_npp_G = sprintf("%.2f", MG)
1527 system("sed s#M_npp_G#"+M_npp_G+"# "+html_name+" > "+html_new)
1528 system("mv -f "+html_new+" "+html_name)
1531 ; plot correlation coef
1534 gRes@txFontHeightF = 0.02
1537 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
1539 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1540 ;--------------------------------------------------------------------
1543 title = "Observed MODIS MOD 17"
1544 resg@tiMainString = title
1546 plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1550 title = "Model "+ model_name
1551 resg@tiMainString = title
1553 plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)
1558 zz = nppmod - nppglobe
1559 title = "Model_"+model_name+" - Observed"
1561 resg@cnMinLevelValF = -500 ; Min level
1562 resg@cnMaxLevelValF = 500. ; Max level
1563 resg@cnLevelSpacingF = 50. ; interval
1564 resg@tiMainString = title
1566 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1568 pres = True ; panel plot mods desired
1569 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1570 ; indiv. plots in panel
1571 pres@gsnMaximize = True ; fill the page
1573 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1579 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1580 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1581 system("rm "+plot_name+"-*."+plot_type_new)
1582 system("rm "+plot_name+"."+plot_type)
1584 ;***************************************************************************
1585 ;(D)-1 zonal line plot, ob
1586 ;***************************************************************************
1588 vv = zonalAve(nppglobe)
1589 vv@long_name = nppglobe@long_name
1591 plot_name = "zonal_ob"
1592 title = "Observed MODIS MOD 17"
1595 resz@tiMainString = title
1597 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1599 plot = gsn_csm_xy (wks,ym,vv,resz)
1604 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1605 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1606 system("rm "+plot_name+"-*."+plot_type_new)
1607 system("rm "+plot_name+"."+plot_type)
1609 ;****************************************************************************
1610 ;(D)-2 zonal line plot, model vs ob
1611 ;****************************************************************************
1613 s = new ((/2,dimsizes(ym)/), float)
1615 s(0,:) = zonalAve(nppglobe)
1616 s(1,:) = zonalAve(nppmod)
1618 s@long_name = nppglobe@long_name
1619 ;-------------------------------------------
1620 ; compute correlation coef and M score
1622 ccrZ = esccr(s(1,:), s(0,:),0)
1625 MZ = (ccrZ*ccrZ)* 5.0
1627 M_npp_Z = sprintf("%.2f", MZ)
1628 system("sed s#M_npp_Z#"+M_npp_Z+"# "+html_name+" > "+html_new)
1629 system("mv -f "+html_new+" "+html_name)
1631 ;-------------------------------------------
1632 plot_name = "zonal_model_vs_ob"
1633 title = "Zonal Average"
1634 resz@tiMainString = title
1636 wks = gsn_open_wks (plot_type,plot_name)
1638 ; resz@vpHeightF = 0.4 ; change aspect ratio of plot
1639 ; resz@vpWidthF = 0.7
1641 resz@xyMonoLineColor = "False" ; want colored lines
1642 resz@xyLineColors = (/"black","red"/) ; colors chosen
1643 ; resz@xyLineThicknesses = (/3.,3./) ; line thicknesses
1644 resz@xyLineThicknesses = (/2.,2./) ; line thicknesses
1645 resz@xyDashPatterns = (/0.,0./) ; make all lines solid
1647 resz@tiYAxisString = s@long_name ; add a axis title
1648 resz@txFontHeightF = 0.0195 ; change title font heights
1651 resz@pmLegendDisplayMode = "Always" ; turn on legend
1652 resz@pmLegendSide = "Top" ; Change location of
1653 ; resz@pmLegendParallelPosF = .45 ; move units right
1654 resz@pmLegendParallelPosF = .82 ; move units right
1655 resz@pmLegendOrthogonalPosF = -0.4 ; move units down
1657 resz@pmLegendWidthF = 0.10 ; Change width and
1658 resz@pmLegendHeightF = 0.10 ; height of legend.
1659 resz@lgLabelFontHeightF = .02 ; change font height
1660 ; resz@lgTitleOn = True ; turn on legend title
1661 ; resz@lgTitleString = "Example" ; create legend title
1662 ; resz@lgTitleFontHeightF = .025 ; font of legend title
1663 resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels
1664 ;--------------------------------------------------------------------
1666 zRes@txFontHeightF = 0.025
1668 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"
1670 gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)
1671 ;--------------------------------------------------------------------
1673 plot = gsn_csm_xy (wks,ym,s,resz)
1678 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1679 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1680 system("rm "+plot_name+"-*."+plot_type_new)
1681 system("rm "+plot_name+"."+plot_type)
1683 ;***************************************************************************
1684 ; tar up output plots
1685 ;***************************************************************************
1687 temp_name = "npp." + model_name
1688 system("mkdir -p " + temp_name)
1689 system("cp "+ html_name + " " +temp_name+"/table.html")
1690 system("mv *.png " + temp_name)
1691 system("tar cf "+ temp_name +".tar " + temp_name)
1693 ;***************************************************************************