diff -r 000000000000 -r 0c6405ab2ff4 npp/18.scatter_bias.ncl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/npp/18.scatter_bias.ncl Mon Jan 26 22:08:20 2009 -0500 @@ -0,0 +1,1096 @@ +; **************************************************** +; combine scatter, histogram, global and zonal plots +; **************************************************** +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" +;************************************************ +procedure pminmax(data:numeric,name:string) +begin + print ("min/max " + name + " = " + min(data) + "/" + max(data)) + if(isatt(data,"units")) then + print (name + " units = " + data@units) + end if +end + +; Main code. +begin + + plot_type = "ps" + plot_type_new = "png" + +;************************************************ +; 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_T31_mean.nc" + fglobe = addfile (dirglobe+filglobe,"r") + + nppglobe0 = fglobe->NPP + nppglobe = nppglobe0 +;************************************************ +; read in model data +;************************************************ + model_name = "b30.061n" + + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" + film = "b30.061n_1995-2004_ANN_climo_lnd.nc" + fm = addfile (dirm+film,"r") + + nppmod0 = fm->NPP + rainmod0 = fm->RAIN + xm = fm->lon + ym = fm->lat + + 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) 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) + +;(B) 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) + +;(C) 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) + +;(D) 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) + +;(E) 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 +;------------------------------- +;compute correlation coef. and M + ccr81 = esccr(nppmod81,npp81,0) + print (ccr81) + +;new eq + bias = sum(abs(nppmod81-npp81)/(nppmod81+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 + 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) + +;(F) 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 +;------------------------------- +;compute correlation coef. and M + ccr933 = esccr(nppmod933,npp933,0) + print (ccr933) + +;new eq + bias = sum(abs(nppmod933-npp933)/(nppmod933+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 + 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) +;************************************************************************** +;(G) histogram 81 + +;-------------------------------------------------------------------- +;get data + + RAIN1_1D = ndtooned(rain81) + RAIN2_1D = ndtooned(rainmod81) + NPP1_1D = ndtooned(npp81) + NPP2_1D = ndtooned(nppmod81) + +; Calculaee "nice" bins for binning the data in equally spaced ranges. + + 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 + 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 +;---------------------------------------- +;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)/(vv+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"/) +;================================================================= +; 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 + + 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) + + clear (wks) +;=========================================================================== +; 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) +;************************************************************************** +;(H) histogram 933 + +;-------------------------------------------------------------------- +;get data + + RAIN1_1D = ndtooned(rain933) + RAIN2_1D = ndtooned(rainmod933) + NPP1_1D = ndtooned(npp933) + NPP2_1D = ndtooned(nppmod933) + +; Calculate "nice" bins for binning the data in equally spaced ranges. + + 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 + 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 +;---------------------------------------- +;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)/(vv+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"/) +;================================================================= +; 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) + +;=========================================================================== +; 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 +;------------------------------------------------------------------------ +;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) +;------------------------------------------------------------------------ +;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) +;------------------------------------------------------------------------ +;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) +;--------------------------------------------------------------------- +; 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) +;--------------------------------------------------------------------- +; 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) +end