diff -r 000000000000 -r 0c6405ab2ff4 npp/99.all.ncl.0 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/npp/99.all.ncl.0 Mon Jan 26 22:08:20 2009 -0500 @@ -0,0 +1,1683 @@ +; **************************************************** +; combine scatter, histogram, global and zonal plots +; compute all correlation coef and M score +; add 1-to-1 line in scatter plots +; **************************************************** +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" +;************************************************ +procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \ + ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \ + ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \ + ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \ + ,dx4[1]:numeric \ + ) +begin +; Calculaee "nice" bins for binning the data in equally spaced ranges. +; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D +; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4 + + nbins = 15 ; Number of bins to use. + + nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True) + nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1) + range = fspan(nicevals(0),nicevals(1),nvals) + +; Use this range information to grab all the values in a +; particular range, and then take an average. + + nr = dimsizes(range) + nx = nr-1 +; print (nx) + +; xvalues = new((/2,nx/),typeof(RAIN1_1D)) + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2. + dx = xvalues(0,1) - xvalues(0,0) ; range width + dx4 = dx/4 ; 1/4 of the range + xvalues(1,:) = xvalues(0,:) - dx/5. +; yvalues = new((/2,nx/),typeof(RAIN1_1D)) +; mn_yvalues = new((/2,nx/),typeof(RAIN1_1D)) +; mx_yvalues = new((/2,nx/),typeof(RAIN1_1D)) + + do nd=0,1 + +; See if we are doing model or observational data. + + if(nd.eq.0) then + data = RAIN1_1D + npp_data = NPP1_1D + else + data = RAIN2_1D + npp_data = NPP2_1D + end if + +; Loop through each range and check for values. + + do i=0,nr-2 + if (i.ne.(nr-2)) then +; print("") +; print("In range ["+range(i)+","+range(i+1)+")") + idx = ind((range(i).le.data).and.(data.lt.range(i+1))) + else +; print("") +; print("In range ["+range(i)+",)") + idx = ind(range(i).le.data) + end if + +; Calculate average, and get min and max. + + if(.not.any(ismissing(idx))) then + yvalues(nd,i) = avg(npp_data(idx)) + mn_yvalues(nd,i) = min(npp_data(idx)) + mx_yvalues(nd,i) = max(npp_data(idx)) + count = dimsizes(idx) + else + count = 0 + yvalues(nd,i) = yvalues@_FillValue + mn_yvalues(nd,i) = yvalues@_FillValue + mx_yvalues(nd,i) = yvalues@_FillValue + end if + +; Print out information. + +; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i)) +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) + + +; Clean up for next time in loop. + + delete(idx) + end do + delete(data) + delete(npp_data) + end do +end + +; Main code. +begin + + plot_type = "ps" + plot_type_new = "png" +;************************************************ +; read in model data +;************************************************ +;film = "i01.06cn_1798-2004_ANN_climo.nc" +;model_name = "06cn" +;model_grid = "T42" + +;film = "i01.06casa_1798-2004_ANN_climo.nc" +;model_name = "06casa" +;model_grid = "T42" + + film = "i01.10casa_1948-2004_ANN_climo.nc" + model_name = "10casa" + model_grid = "T42" + +;film = "i01.10cn_1948-2004_ANN_climo.nc" +;model_name = "10cn" +;model_grid = "T42" +;-------------------------------------------------- + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" + fm = addfile (dirm+film,"r") + + nppmod0 = fm->NPP + rainmod0 = fm->RAIN + xm = fm->lon + ym = fm->lat + +;************************************************ +; read in ob data +;************************************************ +;(1) + dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/" + fil81 = "data.81.nc" + f81 = addfile (dir81+fil81,"r") + + id81 = f81->SITE_ID + npp81 = f81->TNPP_C + rain81 = tofloat(f81->PREC_ANN) + x81 = f81->LONG_DD + y81 = f81->LAT_DD + + id81@long_name = "SITE_ID" + +;change longitude from (-180,180) to (0,360) +;for model data interpolation + + do i= 0,dimsizes(x81)-1 + if (x81(i) .lt. 0.) then + x81(i) = x81(i)+ 360. + end if + end do +;print (x81) +;------------------------------------- +;(2) + dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/" + fil933 = "data.933.nc" + f933 = addfile (dir933+fil933,"r") + + id933 = f933->SITE_ID + npp933 = f933->TNPP_C + rain933 = f933->PREC + x933 = f933->LONG_DD + y933 = f933->LAT_DD + + id933@long_name = "SITE_ID" + +;change longitude from (-180,180) to (0,360) +;for model data interpolation + + do i= 0,dimsizes(x933)-1 + if (x933(i) .lt. 0.) then + x933(i) = x933(i)+ 360. + end if + end do +;print (x933) +;---------------------------------------- +;(3) + dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/" + filglobe = "Npp_"+model_grid+"_mean.nc" + fglobe = addfile (dirglobe+filglobe,"r") + + nppglobe0 = fglobe->NPP + nppglobe = nppglobe0 +;*********************************************************************** +; interpolate model data +;*********************************************************************** + nppmod = nppmod0(0,:,:) + rainmod = rainmod0(0,:,:) + delete (nppmod0) + delete (rainmod0) + + nppmod81 =linint2_points(xm,ym,nppmod,True,x81,y81,0) + + nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0) + + rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0) + + rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0) + +;************************************************ +; Units for these variables are: +; +; rain81 : mm/year +; rainmod : mm/s +; npp81 : g C/m^2/year +; nppmod81: g C/m^2/s +; nppglobe: g C/m^2/year +; +; We want to convert these to "m/year" and "g C/m^2/year". +; + nsec_per_year = 60*60*24*365 + + rain81 = rain81 / 1000. + rainmod81 = (rainmod81/ 1000.) * nsec_per_year + nppmod81 = nppmod81 * nsec_per_year + + rain933 = rain933 / 1000. + rainmod933 = (rainmod933/ 1000.) * nsec_per_year + nppmod933 = nppmod933 * nsec_per_year + + nppmod = nppmod * nsec_per_year + + npp81@units = "gC/m^2/yr" + nppmod81@units = "gC/m^2/yr" + npp933@units = "gC/m^2/yr" + nppmod933@units = "gC/m^2/yr" + nppmod@units = "gC/m^2/yr" + nppglobe@units = "gC/m^2/yr" + rain81@units = "m/yr" + rainmod81@units = "m/yr" + rain933@units = "m/yr" + rainmod933@units = "m/yr" + + npp81@long_name = "NPP (gC/m2/year)" + npp933@long_name = "NPP (gC/m2/year)" + nppmod81@long_name = "NPP (gC/m2/year)" + nppmod933@long_name = "NPP (gC/m2/year)" + nppmod@long_name = "NPP (gC/m2/year)" + nppglobe@long_name = "NPP (gC/m2/year)" + rain81@long_name = "PREC (m/year)" + rain933@long_name = "PREC (m/year)" + rainmod81@long_name = "PREC (m/year)" + rainmod933@long_name = "PREC (m/year)" + +;(A)-1 plot scatter ob 81 + +;plot_name = "scatter_ob_81" +;title = "Observed 81 sites" + +;wks = gsn_open_wks (plot_type,plot_name) ; open workstation +;res = True ; plot mods desired +;res@tiMainString = title ; add title +;res@xyMarkLineModes = "Markers" ; choose which have markers +;res@xyMarkers = 16 ; choose type of marker +;res@xyMarkerColor = "red" ; Marker color +;res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) +;res@tmLabelAutoStride = True ; nice tick mark labels + +;plot = gsn_csm_xy (wks,id81,npp81,res) ; create plot +;frame(wks) + +;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +;system("rm "+plot_name+"."+plot_type) +;system("rm "+plot_name+"-1."+plot_type_new) +;system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + +;clear (wks) +;******************************************************************* +;(A)-2 for table of site(1) (the first 41 sites) +;******************************************************************* + + table_length = 0.95 + +; table header name + table_header_name = "Site ID" + +; column (not including header column) + col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/) + ncol = dimsizes(col_header) + +; row (not including header row) + nrow = 41 + row_header = new ((/nrow/),string ) + row_header(0:nrow-1) = id81(0:nrow-1) + +; Table header + ncr1 = (/1,1/) ; 1 row, 1 column + x1 = (/0.005,0.15/) ; Start and end X + y1 = (/0.900,0.95/) ; Start and end Y + text1 = table_header_name + res1 = True + res1@txFontHeightF = 0.015 + res1@gsFillColor = "CornFlowerBlue" + +; Column header (equally space in x2) + ncr2 = (/1,ncol/) ; 1 rows, ncol columns + x2 = (/x1(1),0.995/) ; start from end of x1 + y2 = y1 ; same as y1 + text2 = col_header + res2 = True + res2@txFontHeightF = 0.015 + res2@gsFillColor = "Gray" + +; Row header (equally space in y2) + ncr3 = (/nrow,1/) ; nrow rows, 1 columns + x3 = x1 ; same as x1 + y3 = (/1.0-table_length,y1(0)/) ; end at start of y1 + text3 = row_header + res3 = True + res3@txFontHeightF = 0.010 + res3@gsFillColor = "Gray" + +; Main table body + ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns + x4 = x2 ; Start and end x + y4 = y3 ; Start and end Y + text4 = new((/nrow,ncol/),string) + + color_fill4 = new((/nrow,ncol/),string) + color_fill4 = "white" +; color_fill4(:,ncol-1) = "grey" +; color_fill4(nrow-1,:) = "green" + + res4 = True ; Set up resource list +; res4@gsnDebug = True ; Useful to print NDC row,col values used. + res4@txFontHeightF = 0.015 + res4@gsFillColor = color_fill4 + + delete (color_fill4) +;------------------------------------------------------------------- +; for table value + + do n = 0,nrow-1 + text4(n,0) = sprintf("%5.2f", y81(n)) + text4(n,1) = sprintf("%5.2f", x81(n)) + text4(n,2) = sprintf("%5.2f", npp81(n)) + text4(n,3) = sprintf("%5.2f", rain81(n)) + end do +;--------------------------------------------------------------------------- + + plot_name = "table_site81_ob1" + + wks = gsn_open_wks (plot_type,plot_name) +;------------------------------------------ +; for table title + + gRes = True + gRes@txFontHeightF = 0.02 +; gRes@txAngleF = 90 + + title_text = "Observation at 81 Sites (1)" + + gsn_text_ndc(wks,title_text,0.50,0.975,gRes) +;------------------------------------------ + + gsn_table(wks,ncr1,x1,y1,text1,res1) + gsn_table(wks,ncr2,x2,y2,text2,res2) + gsn_table(wks,ncr3,x3,y3,text3,res3) + gsn_table(wks,ncr4,x4,y4,text4,res4) + + frame(wks) + clear (wks) + + delete (row_header) + delete (res3) + delete (ncr3) + delete (text3) + delete (res4) + delete (ncr4) + delete (text4) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + +;******************************************************************* +;(A)-3 for table of site(2) (the second 40 sites) +;******************************************************************* + + table_length = 0.95 + +; table header name + table_header_name = "Site ID" + +; column (not including header column) + col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/) + ncol = dimsizes(col_header) + +; row (not including header row) + nrow = 40 + row_header = new ((/nrow/),string ) + row_header(0:nrow-1) = id81(41:41+nrow-1) + +; Table header + ncr1 = (/1,1/) ; 1 row, 1 column + x1 = (/0.005,0.15/) ; Start and end X + y1 = (/0.900,0.95/) ; Start and end Y + text1 = table_header_name + res1 = True + res1@txFontHeightF = 0.015 + res1@gsFillColor = "CornFlowerBlue" + +; Column header (equally space in x2) + ncr2 = (/1,ncol/) ; 1 rows, ncol columns + x2 = (/x1(1),0.995/) ; start from end of x1 + y2 = y1 ; same as y1 + text2 = col_header + res2 = True + res2@txFontHeightF = 0.015 + res2@gsFillColor = "Gray" + +; Row header (equally space in y2) + ncr3 = (/nrow,1/) ; nrow rows, 1 columns + x3 = x1 ; same as x1 + y3 = (/1.0-table_length,y1(0)/) ; end at start of y1 + text3 = row_header + res3 = True + res3@txFontHeightF = 0.010 + res3@gsFillColor = "Gray" + +; Main table body + ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns + x4 = x2 ; Start and end x + y4 = y3 ; Start and end Y + text4 = new((/nrow,ncol/),string) + + color_fill4 = new((/nrow,ncol/),string) + color_fill4 = "white" +; color_fill4(:,ncol-1) = "grey" +; color_fill4(nrow-1,:) = "green" + + res4 = True ; Set up resource list +; res4@gsnDebug = True ; Useful to print NDC row,col values used. + res4@txFontHeightF = 0.015 + res4@gsFillColor = color_fill4 + + delete (color_fill4) +;------------------------------------------------------------------- +; for table value + + do n = 0,nrow-1 + text4(n,0) = sprintf("%5.2f", y81(n+41)) + text4(n,1) = sprintf("%5.2f", x81(n+41)) + text4(n,2) = sprintf("%5.2f", npp81(n+41)) + text4(n,3) = sprintf("%5.2f", rain81(n+41)) + end do +;--------------------------------------------------------------------------- + + plot_name = "table_site81_ob2" + + wks = gsn_open_wks (plot_type,plot_name) +;------------------------------------------ +; for table title + + gRes = True + gRes@txFontHeightF = 0.02 +; gRes@txAngleF = 90 + + title_text = "Observation at 81 Sites (2)" + + gsn_text_ndc(wks,title_text,0.50,0.975,gRes) +;------------------------------------------ + + gsn_table(wks,ncr1,x1,y1,text1,res1) + gsn_table(wks,ncr2,x2,y2,text2,res2) + gsn_table(wks,ncr3,x3,y3,text3,res3) + gsn_table(wks,ncr4,x4,y4,text4,res4) + + frame(wks) + clear (wks) + + delete (row_header) + delete (res3) + delete (ncr3) + delete (text3) + delete (res4) + delete (ncr4) + delete (text4) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) +;************************************************************************** +;(A)-4 plot scatter ob 933 + + plot_name = "scatter_ob_933" + title = "Observed 933 sites" + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + res = True ; plot mods desired + res@tiMainString = title ; add title + res@xyMarkLineModes = "Markers" ; choose which have markers + res@xyMarkers = 16 ; choose type of marker + res@xyMarkerColor = "red" ; Marker color + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) + res@tmLabelAutoStride = True ; nice tick mark labels + + plot = gsn_csm_xy (wks,id933,npp933,res) ; create plot + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) + +;(A)-5 plot scatter model 81 + +;plot_name = "scatter_model_81" +;title = "Model "+ model_name +" 81 sites" + +;wks = gsn_open_wks (plot_type,plot_name) ; open workstation +;res = True ; plot mods desired +;res@tiMainString = title ; add title +;res@xyMarkLineModes = "Markers" ; choose which have markers +;res@xyMarkers = 16 ; choose type of marker +;res@xyMarkerColor = "red" ; Marker color +;res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) +;res@tmLabelAutoStride = True ; nice tick mark labels + +;plot = gsn_csm_xy (wks,id81,nppmod81,res) ; create plot +;frame(wks) + +;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +;system("rm "+plot_name+"."+plot_type) +;system("rm "+plot_name+"-1."+plot_type_new) +;system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + +;clear (wks) + +;(A)-6 plot scatter model 933 + + plot_name = "scatter_model_933" + title = "Model "+ model_name +" 933 sites" + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + res = True ; plot mods desired + res@tiMainString = title ; add title + res@xyMarkLineModes = "Markers" ; choose which have markers + res@xyMarkers = 16 ; choose type of marker + res@xyMarkerColor = "red" ; Marker color + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) + res@tmLabelAutoStride = True ; nice tick mark labels + + plot = gsn_csm_xy (wks,id933,nppmod933,res) ; create plot + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) + +;(A)-7 scatter model vs ob 81 + + plot_name = "scatter_model_vs_ob_81" + title = model_name +" vs ob 81 sites" + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + res = True ; plot mods desired + res@tiMainString = title ; add title + res@xyMarkLineModes = "Markers" ; choose which have markers + res@xyMarkers = 16 ; choose type of marker + res@xyMarkerColor = "red" ; Marker color + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) + res@tmLabelAutoStride = True ; nice tick mark labels + + res@gsnDraw = False + res@gsnFrame = False ; don't advance frame yet +;------------------------------- +;compute correlation coef. and M + ccr81 = esccr(nppmod81,npp81,0) +;print (ccr81) + +;new eq + bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81))) + M81 = (1. - (bias/dimsizes(y81)))*5. + print (M81) + + tRes = True + tRes@txFontHeightF = 0.025 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")" + + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes) +;-------------------------------- + plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot +;------------------------------- +; add polyline + + dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False) +;------------------------------- + draw (plot) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +;system("rm "+plot_name+"-1."+plot_type_new) +;system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) + delete (dum) + +;(A)-8 scatter model vs ob 933 + + plot_name = "scatter_model_vs_ob_933" + title = model_name +" vs ob 933 sites" + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + res = True ; plot mods desired + res@tiMainString = title ; add title + res@xyMarkLineModes = "Markers" ; choose which have markers + res@xyMarkers = 16 ; choose type of marker + res@xyMarkerColor = "red" ; Marker color + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) + res@tmLabelAutoStride = True ; nice tick mark labels + + res@gsnDraw = False + res@gsnFrame = False ; don't advance frame yet +;------------------------------- +;compute correlation coef. and M + ccr933 = esccr(nppmod933,npp933,0) +;print (ccr933) + +;new eq + bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933))) + M933 = (1. - (bias/dimsizes(y933)))*5. + print (M933) + + tRes = True + tRes@txFontHeightF = 0.025 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")" + + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes) +;-------------------------------- + plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot +;------------------------------- +; add polyline + + dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False) +;------------------------------- + draw (plot) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +;system("rm "+plot_name+"-1."+plot_type_new) +;system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) + delete (dum) + +;******************************************************************* +;(A)-9 for table of site(1) (the first 41 sites), model vs ob +;******************************************************************* + + table_length = 0.95 + +; table header name + table_header_name = "Site ID" + +; row (not including header row) + nrow = 41 + row_header = new ((/nrow/),string ) + row_header(0:nrow-1) = id81(0:nrow-1) + +; Table header + ncr1 = (/1,1/) ; 1 row, 1 column + x1 = (/0.005,0.15/) ; Start and end X + y1 = (/0.85 ,0.95/) ; Start and end Y + text1 = table_header_name + res1 = True + res1@txFontHeightF = 0.015 + res1@gsFillColor = "CornFlowerBlue" + +; Column header1 (equally space in x2) + + delete (col_header) + delete (ncr2) + delete (text2) + + col_header = (/"Latitude","Longitude"/) + ncol = dimsizes(col_header) + + ncr2 = (/1,ncol/) ; 1 rows, ncol columns + x2 = (/x1(1),0.35/) ; start from end of x1 + y2 = y1 ; same as y1 + text2 = col_header + res2 = True + res2@txFontHeightF = 0.015 + res2@gsFillColor = "Gray" + +; Column header2 (equally space in x2) + + col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/) + col_header2 = (/"ob",model_name \ + ,"ob",model_name \ + /) + + ncol1 = dimsizes(col_header1) + ncol2 = dimsizes(col_header2) + + ncr21 = (/1,ncol1/) ; 1 rows, 4 columns + x21 = (/x2(1),0.995/) ; start from end of x2 + yhalf = (y1(0)+y1(1))*0.5 + y21 = (/yhalf,y1(1)/) ; top half of y1 + text21 = col_header1 + res21 = True + res21@txFontHeightF = 0.015 + res21@gsFillColor = "Gray" + + ncr22 = (/1,ncol2/) ; 1 rows, 12 columns + x22 = x21 ; start from end of x1 + y22 = (/y1(0),yhalf/) ; bottom half of y1 + text22 = col_header2 + res22 = True + res22@txFontHeightF = 0.015 + res22@gsFillColor = "Gray" + +; Row header (equally space in y2) + ncr3 = (/nrow,1/) ; nrow rows, 1 columns + x3 = x1 ; same as x1 + y3 = (/1.0-table_length,y1(0)/) ; end at start of y1 + text3 = row_header + res3 = True + res3@txFontHeightF = 0.010 + res3@gsFillColor = "Gray" + +; Main table body-1 + ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns + x4 = x2 ; Start and end x + y4 = y3 ; Start and end Y + text4 = new((/nrow,ncol/),string) + + color_fill4 = new((/nrow,ncol/),string) + color_fill4 = "white" +; color_fill4(:,ncol-1) = "grey" +; color_fill4(nrow-1,:) = "green" + + res4 = True ; Set up resource list +; res4@gsnDebug = True ; Useful to print NDC row,col values used. + res4@txFontHeightF = 0.015 + res4@gsFillColor = color_fill4 + + delete (color_fill4) + +; Main table body-2 + ncr42 = (/nrow,ncol2/) ; nrow rows, ncol columns + x42 = x21 ; Start and end x + y42 = y3 ; Start and end Y + text42 = new((/nrow,ncol2/),string) + + color_fill42 = new((/nrow,ncol2/),string) + color_fill42 = "white" +; color_fill42(:,ncol-1) = "grey" +; color_fill42(nrow-1,:) = "green" + + res42 = True ; Set up resource list +; res42@gsnDebug = True ; Useful to print NDC row,col values used. + res42@txFontHeightF = 0.015 + res42@gsFillColor = color_fill42 + + delete (color_fill42) +;------------------------------------------------------------------- +; for table value + + do n = 0,nrow-1 + text4(n,0) = sprintf("%5.2f", y81(n)) + text4(n,1) = sprintf("%5.2f", x81(n)) + end do + + do n = 0,nrow-1 + text42(n,0) = sprintf("%5.2f", npp81(n)) + text42(n,1) = sprintf("%5.2f", nppmod81(n)) + text42(n,2) = sprintf("%5.2f", rain81(n)) + text42(n,3) = sprintf("%5.2f", rainmod81(n)) + end do +;--------------------------------------------------------------------------- + + plot_name = "table_site81_model_vs_ob1" + + wks = gsn_open_wks (plot_type,plot_name) +;------------------------------------------ +; for table title + + gRes = True + gRes@txFontHeightF = 0.02 +; gRes@txAngleF = 90 + + title_text = "Model vs Observation at 81 Sites (1)" + + gsn_text_ndc(wks,title_text,0.50,0.975,gRes) +;------------------------------------------ + + gsn_table(wks,ncr1,x1,y1,text1,res1) + gsn_table(wks,ncr2,x2,y2,text2,res2) + gsn_table(wks,ncr21,x21,y21,text21,res21) + gsn_table(wks,ncr22,x22,y22,text22,res22) + gsn_table(wks,ncr3,x3,y3,text3,res3) + gsn_table(wks,ncr4,x4,y4,text4,res4) + gsn_table(wks,ncr42,x42,y42,text42,res42) + + frame(wks) + clear (wks) + + delete (col_header) + delete (row_header) + delete (res1) + delete (ncr1) + delete (text1) + delete (res2) + delete (ncr2) + delete (text2) + delete (res21) + delete (ncr21) + delete (text21) + delete (res22) + delete (ncr22) + delete (text22) + delete (res3) + delete (ncr3) + delete (text3) + delete (res4) + delete (ncr4) + delete (text4) + delete (res42) + delete (ncr42) + delete (text42) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + +;******************************************************************* +;(A)-10 for table of site(2) (the last 40 sites), model vs ob +;******************************************************************* + + table_length = 0.95 + +; table header name + table_header_name = "Site ID" + +; row (not including header row) + nrow = 40 + row_header = new ((/nrow/),string ) + row_header(0:nrow-1) = id81(41:41+nrow-1) + +; Table header + ncr1 = (/1,1/) ; 1 row, 1 column + x1 = (/0.005,0.15/) ; Start and end X + y1 = (/0.85 ,0.95/) ; Start and end Y + text1 = table_header_name + res1 = True + res1@txFontHeightF = 0.015 + res1@gsFillColor = "CornFlowerBlue" + +; Column header1 (equally space in x2) + + col_header = (/"Latitude","Longitude"/) + ncol = dimsizes(col_header) + + ncr2 = (/1,ncol/) ; 1 rows, ncol columns + x2 = (/x1(1),0.35/) ; start from end of x1 + y2 = y1 ; same as y1 + text2 = col_header + res2 = True + res2@txFontHeightF = 0.015 + res2@gsFillColor = "Gray" + +; Column header2 (equally space in x2) + +; col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/) + col_header2 = (/"ob",model_name \ + ,"ob",model_name \ + /) + + ncol1 = dimsizes(col_header1) + ncol2 = dimsizes(col_header2) + + ncr21 = (/1,ncol1/) ; 1 rows, 4 columns + x21 = (/x2(1),0.995/) ; start from end of x2 + yhalf = (y1(0)+y1(1))*0.5 + y21 = (/yhalf,y1(1)/) ; top half of y1 + text21 = col_header1 + res21 = True + res21@txFontHeightF = 0.015 + res21@gsFillColor = "Gray" + + ncr22 = (/1,ncol2/) ; 1 rows, 12 columns + x22 = x21 ; start from end of x1 + y22 = (/y1(0),yhalf/) ; bottom half of y1 + text22 = col_header2 + res22 = True + res22@txFontHeightF = 0.015 + res22@gsFillColor = "Gray" + +; Row header (equally space in y2) + ncr3 = (/nrow,1/) ; nrow rows, 1 columns + x3 = x1 ; same as x1 + y3 = (/1.0-table_length,y1(0)/) ; end at start of y1 + text3 = row_header + res3 = True + res3@txFontHeightF = 0.010 + res3@gsFillColor = "Gray" + +; Main table body-1 + ncr4 = (/nrow,ncol/) ; nrow rows, ncol columns + x4 = x2 ; Start and end x + y4 = y3 ; Start and end Y + text4 = new((/nrow,ncol/),string) + + color_fill4 = new((/nrow,ncol/),string) + color_fill4 = "white" +; color_fill4(:,ncol-1) = "grey" +; color_fill4(nrow-1,:) = "green" + + res4 = True ; Set up resource list +; res4@gsnDebug = True ; Useful to print NDC row,col values used. + res4@txFontHeightF = 0.015 + res4@gsFillColor = color_fill4 + + delete (color_fill4) + +; Main table body-2 + ncr42 = (/nrow,ncol2/) ; nrow rows, ncol columns + x42 = x21 ; Start and end x + y42 = y3 ; Start and end Y + text42 = new((/nrow,ncol2/),string) + + color_fill42 = new((/nrow,ncol2/),string) + color_fill42 = "white" +; color_fill42(:,ncol-1) = "grey" +; color_fill42(nrow-1,:) = "green" + + res42 = True ; Set up resource list +; res42@gsnDebug = True ; Useful to print NDC row,col values used. + res42@txFontHeightF = 0.015 + res42@gsFillColor = color_fill42 + + delete (color_fill42) +;------------------------------------------------------------------- +; for table value + + do n = 0,nrow-1 + text4(n,0) = sprintf("%5.2f", y81(n+41)) + text4(n,1) = sprintf("%5.2f", x81(n+41)) + end do + + do n = 0,nrow-1 + text42(n,0) = sprintf("%5.2f", npp81(n+41)) + text42(n,1) = sprintf("%5.2f", nppmod81(n+41)) + text42(n,2) = sprintf("%5.2f", rain81(n+41)) + text42(n,3) = sprintf("%5.2f", rainmod81(n+41)) + end do +;--------------------------------------------------------------------------- + + plot_name = "table_site81_model_vs_ob2" + + wks = gsn_open_wks (plot_type,plot_name) +;------------------------------------------ +; for table title + + gRes = True + gRes@txFontHeightF = 0.02 +; gRes@txAngleF = 90 + + title_text = "Model vs Observation at 81 Sites (2)" + + gsn_text_ndc(wks,title_text,0.50,0.975,gRes) +;------------------------------------------ + + gsn_table(wks,ncr1,x1,y1,text1,res1) + gsn_table(wks,ncr2,x2,y2,text2,res2) + gsn_table(wks,ncr21,x21,y21,text21,res21) + gsn_table(wks,ncr22,x22,y22,text22,res22) + gsn_table(wks,ncr3,x3,y3,text3,res3) + gsn_table(wks,ncr4,x4,y4,text4,res4) + gsn_table(wks,ncr42,x42,y42,text42,res42) + + frame(wks) + clear (wks) + + delete (col_header) + delete (row_header) + delete (res1) + delete (ncr1) + delete (text1) + delete (res2) + delete (ncr2) + delete (text2) + delete (res21) + delete (ncr21) + delete (text21) + delete (res22) + delete (ncr22) + delete (text22) + delete (res3) + delete (ncr3) + delete (text3) + delete (res4) + delete (ncr4) + delete (text4) + delete (res42) + delete (ncr42) + delete (text42) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + +;************************************************************************** +;(B) histogram 81 +;***********************************************************************- +;get data + + RAIN1_1D = ndtooned(rain81) + RAIN2_1D = ndtooned(rainmod81) + NPP1_1D = ndtooned(npp81) + NPP2_1D = ndtooned(nppmod81) + + nx = 8 + + xvalues = new((/2,nx/),float) + yvalues = new((/2,nx/),float) + mn_yvalues = new((/2,nx/),float) + mx_yvalues = new((/2,nx/),float) + dx4 = new((/1/),float) + + get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \ + ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4) + +;---------------------------------------- +;compute correlation coeff and M score + + u = yvalues(0,:) + v = yvalues(1,:) + + good = ind(.not.ismissing(u) .and. .not.ismissing(v)) + uu = u(good) + vv = v(good) + + ccr81h = esccr(uu,vv,0) +;print (ccr81h) + +;new eq + bias = sum(abs(vv-uu)/(abs(vv)+abs(uu))) + M81h = (1.- (bias/dimsizes(uu)))*5. + print (M81h) + + delete (u) + delete (v) + delete (uu) + delete (vv) +;---------------------------------------------------------------------- +; histogram res + + resh = True + resh@gsnMaximize = True + resh@gsnDraw = False + resh@gsnFrame = False + resh@xyMarkLineMode = "Markers" + resh@xyMarkerSizeF = 0.014 + resh@xyMarker = 16 + resh@xyMarkerColors = (/"Brown","Blue"/) + resh@trYMinF = min(mn_yvalues) - 10. + resh@trYMaxF = max(mx_yvalues) + 10. + + resh@tiYAxisString = "NPP (g C/m2/year)" + resh@tiXAxisString = "Precipitation (m/year)" + + max_bar = new((/2,nx/),graphic) + min_bar = new((/2,nx/),graphic) + max_cap = new((/2,nx/),graphic) + min_cap = new((/2,nx/),graphic) + + lnres = True + line_colors = (/"brown","blue"/) +;================================================================= +;(B)-1 histogram ob 81 site only +; + plot_name = "histogram_ob_81" + title = "Observed 81 site" + resh@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) + + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh) + +;------------------------------- +;Attach the vertical bar and the horizontal cap line + + delete (x1) + delete (x2) + delete (y1) + delete (y2) + + do nd=0,0 + lnres@gsLineColor = line_colors(nd) + do i=0,nx-1 + + if(.not.ismissing(mn_yvalues(nd,i)).and. \ + .not.ismissing(mx_yvalues(nd,i))) then +; +; Attach the vertical bar, both above and below the marker. +; + x1 = xvalues(nd,i) + y1 = yvalues(nd,i) + y2 = mn_yvalues(nd,i) + plx = (/x1,x1/) + ply = (/y1,y2/) + min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres) + + y2 = mx_yvalues(nd,i) + plx = (/x1,x1/) + ply = (/y1,y2/) + max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres) +; +; Attach the horizontal cap line, both above and below the marker. +; + x1 = xvalues(nd,i) - dx4 + x2 = xvalues(nd,i) + dx4 + y1 = mn_yvalues(nd,i) + plx = (/x1,x2/) + ply = (/y1,y1/) + min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres) + + y1 = mx_yvalues(nd,i) + plx = (/x1,x2/) + ply = (/y1,y1/) + max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres) + end if + end do + end do + + draw(xy) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;=========================================================================== +;(B)-2 histogram model vs ob 81 site + + plot_name = "histogram_mod-ob_81" + title = model_name+ " vs Observed 81 site" + resh@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + +;----------------------------- +; Add a boxed legend using the more simple method, which won't have +; vertical lines going through the markers. + + resh@pmLegendDisplayMode = "Always" +; resh@pmLegendWidthF = 0.1 + resh@pmLegendWidthF = 0.08 + resh@pmLegendHeightF = 0.05 + resh@pmLegendOrthogonalPosF = -1.17 +; resh@pmLegendOrthogonalPosF = -1.00 ;(downward) +; resh@pmLegendParallelPosF = 0.18 + resh@pmLegendParallelPosF = 0.88 ;(rightward) + +; resh@lgPerimOn = False + resh@lgLabelFontHeightF = 0.015 + resh@xyExplicitLegendLabels = (/"observed",model_name/) +;----------------------------- + tRes = True + tRes@txFontHeightF = 0.025 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")" + + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) + + xy = gsn_csm_xy(wks,xvalues,yvalues,resh) +;------------------------------- +;Attach the vertical bar and the horizontal cap line + + do nd=0,1 + lnres@gsLineColor = line_colors(nd) + do i=0,nx-1 + + if(.not.ismissing(mn_yvalues(nd,i)).and. \ + .not.ismissing(mx_yvalues(nd,i))) then +; +; Attach the vertical bar, both above and below the marker. +; + x1 = xvalues(nd,i) + y1 = yvalues(nd,i) + y2 = mn_yvalues(nd,i) + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + + y2 = mx_yvalues(nd,i) + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) +; +; Attach the horizontal cap line, both above and below the marker. +; + x1 = xvalues(nd,i) - dx4 + x2 = xvalues(nd,i) + dx4 + y1 = mn_yvalues(nd,i) + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + + y1 = mx_yvalues(nd,i) + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + end if + end do + end do + draw(xy) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) + + delete (RAIN1_1D) + delete (RAIN2_1D) + delete (NPP1_1D) + delete (NPP2_1D) +;delete (range) + delete (xvalues) + delete (yvalues) + delete (mn_yvalues) + delete (mx_yvalues) + delete (good) + delete (max_bar) + delete (min_bar) + delete (max_cap) + delete (min_cap) +;************************************************************************** +;(B) histogram 933 + +;-------------------------------------------------------------------- +;get data + + RAIN1_1D = ndtooned(rain933) + RAIN2_1D = ndtooned(rainmod933) + NPP1_1D = ndtooned(npp933) + NPP2_1D = ndtooned(nppmod933) + + nx = 10 + xvalues = new((/2,nx/),float) + yvalues = new((/2,nx/),float) + mn_yvalues = new((/2,nx/),float) + mx_yvalues = new((/2,nx/),float) + dx4 = new((/1/),float) + + get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \ + ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4) + +;---------------------------------------- +;compute correlation coeff and M score + + u = yvalues(0,:) + v = yvalues(1,:) + + good = ind(.not.ismissing(u) .and. .not.ismissing(v)) + uu = u(good) + vv = v(good) + + ccr933h = esccr(uu,vv,0) +;print (ccr933h) + +;new eq + bias = sum(abs(vv-uu)/(abs(vv)+abs(uu))) + M933h = (1.- (bias/dimsizes(uu)))*5. + print (M933h) + + delete (u) + delete (v) + delete (uu) + delete (vv) +;---------------------------------------------------------------------- +; histogram res + + delete (resh) + resh = True + resh@gsnMaximize = True + resh@gsnDraw = False + resh@gsnFrame = False + resh@xyMarkLineMode = "Markers" + resh@xyMarkerSizeF = 0.014 + resh@xyMarker = 16 + resh@xyMarkerColors = (/"Brown","Blue"/) + resh@trYMinF = min(mn_yvalues) - 10. + resh@trYMaxF = max(mx_yvalues) + 10. + + resh@tiYAxisString = "NPP (g C/m2/year)" + resh@tiXAxisString = "Precipitation (m/year)" + + max_bar = new((/2,nx/),graphic) + min_bar = new((/2,nx/),graphic) + max_cap = new((/2,nx/),graphic) + min_cap = new((/2,nx/),graphic) + + lnres = True + line_colors = (/"brown","blue"/) +;================================================================= +;(B)-3 histogram ob 933 site only + + plot_name = "histogram_ob_933" + title = "Observed 933 site" + resh@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) + + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh) + +;------------------------------- +;Attach the vertical bar and the horizontal cap line + + do nd=0,0 + lnres@gsLineColor = line_colors(nd) + do i=0,nx-1 + + if(.not.ismissing(mn_yvalues(nd,i)).and. \ + .not.ismissing(mx_yvalues(nd,i))) then + +; Attach the vertical bar, both above and below the marker. + + x1 = xvalues(nd,i) + y1 = yvalues(nd,i) + y2 = mn_yvalues(nd,i) + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + + y2 = mx_yvalues(nd,i) + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + +; Attach the horizontal cap line, both above and below the marker. + + x1 = xvalues(nd,i) - dx4 + x2 = xvalues(nd,i) + dx4 + y1 = mn_yvalues(nd,i) + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + + y1 = mx_yvalues(nd,i) + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + end if + end do + end do + + draw(xy) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + delete (xy) + clear (wks) + +;=========================================================================== +;(B)-4 histogram model vs ob 933 site + + plot_name = "histogram_mod-ob_933" + title = model_name+ " vs Observed 933 site" + resh@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + +;----------------------------- +; Add a boxed legend using the more simple method, which won't have +; vertical lines going through the markers. + + resh@pmLegendDisplayMode = "Always" +; resh@pmLegendWidthF = 0.1 + resh@pmLegendWidthF = 0.08 + resh@pmLegendHeightF = 0.05 + resh@pmLegendOrthogonalPosF = -1.17 +; resh@pmLegendOrthogonalPosF = -1.00 ;(downward) +; resh@pmLegendParallelPosF = 0.18 + resh@pmLegendParallelPosF = 0.88 ;(rightward) + +; resh@lgPerimOn = False + resh@lgLabelFontHeightF = 0.015 + resh@xyExplicitLegendLabels = (/"observed",model_name/) +;----------------------------- + tRes = True + tRes@txFontHeightF = 0.025 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")" + + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) + + xy = gsn_csm_xy(wks,xvalues,yvalues,resh) +;------------------------------- +;Attach the vertical bar and the horizontal cap line + + do nd=0,1 + lnres@gsLineColor = line_colors(nd) + do i=0,nx-1 + + if(.not.ismissing(mn_yvalues(nd,i)).and. \ + .not.ismissing(mx_yvalues(nd,i))) then +; +; Attach the vertical bar, both above and below the marker. +; + x1 = xvalues(nd,i) + y1 = yvalues(nd,i) + y2 = mn_yvalues(nd,i) + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + + y2 = mx_yvalues(nd,i) + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) +; +; Attach the horizontal cap line, both above and below the marker. +; + x1 = xvalues(nd,i) - dx4 + x2 = xvalues(nd,i) + dx4 + y1 = mn_yvalues(nd,i) + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + + y1 = mx_yvalues(nd,i) + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + end if + end do + end do + draw(xy) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + delete (xy) + clear (wks) +;***********************************************************************- +; global contour + +;res + resg = True ; Use plot options + resg@cnFillOn = True ; Turn on color fill + resg@gsnSpreadColors = True ; use full colormap +; resg@cnFillMode = "RasterFill" ; Turn on raster color +; resg@lbLabelAutoStride = True + resg@cnLinesOn = False ; Turn off contourn lines + resg@mpFillOn = False ; Turn off map fill + + resg@gsnSpreadColors = True ; use full colormap + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 0. ; Min level + resg@cnMaxLevelValF = 2200. ; Max level + resg@cnLevelSpacingF = 200. ; interval +;------------------------------------------------------------------------ +;(C)-1 global contour ob + + delta = 0.00000000001 + nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe) + + plot_name = "global_ob" + title = "Observed MODIS MOD 17" + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot = gsn_csm_contour_map_ce(wks,nppglobe,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;(C)-2 global contour model + + plot_name = "global_model" + title = "Model "+ model_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot = gsn_csm_contour_map_ce(wks,nppmod,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;(C)-3 global contour model vs ob + + plot_name = "global_model_vs_ob" + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + delete (plot) + plot=new(3,graphic) ; create graphic array + + resg@gsnFrame = False ; Do not draw plot + resg@gsnDraw = False ; Do not advance frame + +;(d) compute correlation coef and M score + + uu1 = ndtooned(nppmod) + vv1 = ndtooned(nppglobe) + + delete (good) + good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1)) + + ug = uu1(good) + vg = vv1(good) + + ccrG = esccr(ug,vg,0) +; print (ccrG) + + MG = (ccrG*ccrG)* 5.0 + print (MG) + +; plot correlation coef + + gRes = True + gRes@txFontHeightF = 0.02 + gRes@txAngleF = 90 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")" + + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) +;-------------------------------------------------------------------- + +;(a) ob + + title = "Observed MODIS MOD 17" + resg@tiMainString = title + + plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg) + +;(b) model + + title = "Model "+ model_name + resg@tiMainString = title + + plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg) + +;(c) model-ob + + zz = nppmod + zz = nppmod - nppglobe + title = "Model_"+model_name+" - Observed" + + resg@cnMinLevelValF = -500 ; Min level + resg@cnMaxLevelValF = 500. ; Max level + resg@cnLevelSpacingF = 50. ; interval + resg@tiMainString = title + + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) + + pres = True ; panel plot mods desired + pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around + ; indiv. plots in panel + pres@gsnMaximize = True ; fill the page + + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + frame (wks) + clear (wks) + + delete (plot) +;********************************************************************** +;(D)-1 zonal line plot ob + + resz = True + + vv = zonalAve(nppglobe) + vv@long_name = nppglobe@long_name + + plot_name = "zonal_ob" + title = "Observed MODIS MOD 17" + resz@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + + plot = gsn_csm_xy (wks,ym,vv,resz) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;--------------------------------------------------------------------- +;(D)-2 zonal line plot model vs ob + + s = new ((/2,dimsizes(ym)/), float) + + s(0,:) = zonalAve(nppglobe) + s(1,:) = zonalAve(nppmod) + + s@long_name = nppglobe@long_name +;------------------------------------------- +;(d) compute correlation coef and M score + + ccrZ = esccr(s(1,:), s(0,:),0) +; print (ccrZ) + + MZ = (ccrZ*ccrZ)* 5.0 + print (MZ) +;------------------------------------------- + plot_name = "zonal_model_vs_ob" + title = "Zonal Average" + resz@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + +; resz@vpHeightF = 0.4 ; change aspect ratio of plot +; resz@vpWidthF = 0.7 + + resz@xyMonoLineColor = "False" ; want colored lines + resz@xyLineColors = (/"black","red"/) ; colors chosen +; resz@xyLineThicknesses = (/3.,3./) ; line thicknesses + resz@xyLineThicknesses = (/2.,2./) ; line thicknesses + resz@xyDashPatterns = (/0.,0./) ; make all lines solid + + resz@tiYAxisString = s@long_name ; add a axis title + resz@txFontHeightF = 0.0195 ; change title font heights + +; Legent + resz@pmLegendDisplayMode = "Always" ; turn on legend + resz@pmLegendSide = "Top" ; Change location of +; resz@pmLegendParallelPosF = .45 ; move units right + resz@pmLegendParallelPosF = .82 ; move units right + resz@pmLegendOrthogonalPosF = -0.4 ; move units down + + resz@pmLegendWidthF = 0.10 ; Change width and + resz@pmLegendHeightF = 0.10 ; height of legend. + resz@lgLabelFontHeightF = .02 ; change font height +; resz@lgTitleOn = True ; turn on legend title +; resz@lgTitleString = "Example" ; create legend title +; resz@lgTitleFontHeightF = .025 ; font of legend title + resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels +;-------------------------------------------------------------------- + zRes = True + zRes@txFontHeightF = 0.025 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")" + + gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes) +;-------------------------------------------------------------------- + + plot = gsn_csm_xy (wks,ym,s,resz) ; create plot + + frame(wks) ; advance frame + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + system("rm "+plot_name+"."+plot_type) + system("rm "+plot_name+"-1."+plot_type_new) + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ + temp_name = "temp." + model_name + system("mkdir -p " + temp_name) + system("mv *.png " + temp_name) + system("tar cf "+ temp_name +".tar " + temp_name) +;------------------------------------------------------------------------ +end +