forrest@0: ; *********************************************** forrest@0: ; combine all scatter forrest@0: ; *********************************************** forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" forrest@0: ;************************************************ forrest@0: procedure pminmax(data:numeric,name:string) forrest@0: begin forrest@0: print ("min/max " + name + " = " + min(data) + "/" + max(data)) forrest@0: if(isatt(data,"units")) then forrest@0: print (name + " units = " + data@units) forrest@0: end if forrest@0: end forrest@0: forrest@0: ; Main code. forrest@0: begin forrest@0: forrest@0: plot_type = "ps" forrest@0: plot_type_new = "png" forrest@0: forrest@0: ;************************************************ forrest@0: ; read in ob data forrest@0: ;************************************************ forrest@0: dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/" forrest@0: fil81 = "data.81.nc" forrest@0: f81 = addfile (dir81+fil81,"r") forrest@0: forrest@0: id81 = f81->SITE_ID forrest@0: npp81 = f81->TNPP_C forrest@0: rain81 = tofloat(f81->PREC_ANN) forrest@0: x81 = f81->LONG_DD forrest@0: y81 = f81->LAT_DD forrest@0: forrest@0: id81@long_name = "SITE_ID" forrest@0: npp81@long_name = "NPP (gC/m2/year)" forrest@0: forrest@0: ;change longitude from (-180,180) to (0,360) forrest@0: ;for model data interpolation forrest@0: forrest@0: do i= 0,dimsizes(x81)-1 forrest@0: if (x81(i) .lt. 0.) then forrest@0: x81(i) = x81(i)+ 360. forrest@0: end if forrest@0: end do forrest@0: ;print (x81) forrest@0: ;------------------------------------- forrest@0: dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/" forrest@0: fil933 = "data.933.nc" forrest@0: f933 = addfile (dir933+fil933,"r") forrest@0: forrest@0: id933 = f933->SITE_ID forrest@0: npp933 = f933->TNPP_C forrest@0: rain933 = f933->PREC forrest@0: x933 = f933->LONG_DD forrest@0: y933 = f933->LAT_DD forrest@0: forrest@0: id933@long_name = "SITE_ID" forrest@0: npp933@long_name = "NPP (gC/m2/year)" forrest@0: forrest@0: ;change longitude from (-180,180) to (0,360) forrest@0: ;for model data interpolation forrest@0: forrest@0: do i= 0,dimsizes(x933)-1 forrest@0: if (x933(i) .lt. 0.) then forrest@0: x933(i) = x933(i)+ 360. forrest@0: end if forrest@0: end do forrest@0: ;print (x933) forrest@0: ;************************************************ forrest@0: ; read in model data forrest@0: ;************************************************ forrest@0: model_name = "b30.061n" forrest@0: forrest@0: dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" forrest@0: film = "b30.061n_1995-2004_ANN_climo_lnd.nc" forrest@0: fm = addfile (dirm+film,"r") forrest@0: forrest@0: nppmod = fm->NPP forrest@0: rainmod = fm->RAIN forrest@0: xm = fm->lon forrest@0: ym = fm->lat forrest@0: forrest@0: nppmod81 =linint2_points(xm,ym,nppmod(0,:,:),True,x81,y81,0) forrest@0: forrest@0: nppmod933 =linint2_points(xm,ym,nppmod(0,:,:),True,x933,y933,0) forrest@0: forrest@0: rainmod81 =linint2_points(xm,ym,rainmod(0,:,:),True,x81,y81,0) forrest@0: forrest@0: rainmod933=linint2_points(xm,ym,rainmod(0,:,:),True,x933,y933,0) forrest@0: forrest@0: ;************************************************ forrest@0: ; Units for these four variables are: forrest@0: ; forrest@0: ; rain81 : mm/year forrest@0: ; rainmod : mm/s forrest@0: ; npp81 : g C/m^2/year forrest@0: ; nppmod81: g C/m^2/s forrest@0: ; forrest@0: ; We want to convert these to "m/year" and "g C/m^2/year". forrest@0: ; forrest@0: nsec_per_year = 60*60*24*365 forrest@0: forrest@0: rain81 = rain81 / 1000. forrest@0: rainmod81 = (rainmod81/ 1000.) * nsec_per_year forrest@0: nppmod81 = nppmod81 * nsec_per_year forrest@0: forrest@0: rain933 = rain933 / 1000. forrest@0: rainmod933 = (rainmod933/ 1000.) * nsec_per_year forrest@0: nppmod933 = nppmod933 * nsec_per_year forrest@0: forrest@0: npp81@units = "gC/m^2/yr" forrest@0: nppmod81@units = "gC/m^2/yr" forrest@0: npp933@units = "gC/m^2/yr" forrest@0: nppmod933@units = "gC/m^2/yr" forrest@0: rain81@units = "m/yr" forrest@0: rainmod81@units = "m/yr" forrest@0: rain933@units = "m/yr" forrest@0: rainmod933@units = "m/yr" forrest@0: forrest@0: npp81@long_name = "NPP (gC/m2/year)" forrest@0: npp933@long_name = "NPP (gC/m2/year)" forrest@0: nppmod81@long_name = "NPP (gC/m2/year)" forrest@0: nppmod933@long_name = "NPP (gC/m2/year)" forrest@0: rain81@long_name = "PREC (m/year)" forrest@0: rain933@long_name = "PREC (m/year)" forrest@0: rainmod81@long_name = "PREC (m/year)" forrest@0: rainmod933@long_name = "PREC (m/year)" forrest@0: forrest@0: ;(A) plot scatter ob 81 forrest@0: forrest@0: plot_name = "scatter_ob_81" forrest@0: title = "Observed 81 sites" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: res = True ; plot mods desired forrest@0: res@tiMainString = title ; add title forrest@0: res@xyMarkLineModes = "Markers" ; choose which have markers forrest@0: res@xyMarkers = 16 ; choose type of marker forrest@0: res@xyMarkerColor = "red" ; Marker color forrest@0: res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) forrest@0: res@tmLabelAutoStride = True ; nice tick mark labels forrest@0: forrest@0: plot = gsn_csm_xy (wks,id81,npp81,res) ; create plot forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: system("rm "+plot_name+"-1."+plot_type_new) forrest@0: system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: ;(B) plot scatter ob 933 forrest@0: forrest@0: plot_name = "scatter_ob_933" forrest@0: title = "Observed 933 sites" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: res = True ; plot mods desired forrest@0: res@tiMainString = title ; add title forrest@0: res@xyMarkLineModes = "Markers" ; choose which have markers forrest@0: res@xyMarkers = 16 ; choose type of marker forrest@0: res@xyMarkerColor = "red" ; Marker color forrest@0: res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) forrest@0: res@tmLabelAutoStride = True ; nice tick mark labels forrest@0: forrest@0: plot = gsn_csm_xy (wks,id933,npp933,res) ; create plot forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: system("rm "+plot_name+"-1."+plot_type_new) forrest@0: system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: ;(C) plot scatter model 81 forrest@0: forrest@0: plot_name = "scatter_model_81" forrest@0: title = "Model "+ model_name +" 81 sites" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: res = True ; plot mods desired forrest@0: res@tiMainString = title ; add title forrest@0: res@xyMarkLineModes = "Markers" ; choose which have markers forrest@0: res@xyMarkers = 16 ; choose type of marker forrest@0: res@xyMarkerColor = "red" ; Marker color forrest@0: res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) forrest@0: res@tmLabelAutoStride = True ; nice tick mark labels forrest@0: forrest@0: plot = gsn_csm_xy (wks,id81,nppmod81,res) ; create plot forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: system("rm "+plot_name+"-1."+plot_type_new) forrest@0: system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: ;(D) plot scatter model 933 forrest@0: forrest@0: plot_name = "scatter_model_933" forrest@0: title = "Model "+ model_name +" 933 sites" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: res = True ; plot mods desired forrest@0: res@tiMainString = title ; add title forrest@0: res@xyMarkLineModes = "Markers" ; choose which have markers forrest@0: res@xyMarkers = 16 ; choose type of marker forrest@0: res@xyMarkerColor = "red" ; Marker color forrest@0: res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) forrest@0: res@tmLabelAutoStride = True ; nice tick mark labels forrest@0: forrest@0: plot = gsn_csm_xy (wks,id933,nppmod933,res) ; create plot forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: system("rm "+plot_name+"-1."+plot_type_new) forrest@0: system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: ;(E) scatter model vs ob 81 forrest@0: forrest@0: plot_name = "scatter_model_vs_ob_81" forrest@0: title = model_name +" vs ob 81 sites" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: res = True ; plot mods desired forrest@0: res@tiMainString = title ; add title forrest@0: res@xyMarkLineModes = "Markers" ; choose which have markers forrest@0: res@xyMarkers = 16 ; choose type of marker forrest@0: res@xyMarkerColor = "red" ; Marker color forrest@0: res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) forrest@0: res@tmLabelAutoStride = True ; nice tick mark labels forrest@0: ;------------------------------- forrest@0: ;compute correlation coef. and M forrest@0: ccr81 = esccr(nppmod81,npp81,0) forrest@0: print (ccr81) forrest@0: forrest@0: ;new eq forrest@0: bias = sum(abs(nppmod81-npp81)/(nppmod81+npp81)) forrest@0: M81 = (1. - (bias/dimsizes(y81)))*5. forrest@0: print (M81) forrest@0: forrest@0: tRes = True forrest@0: tRes@txFontHeightF = 0.025 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes) forrest@0: ;-------------------------------- forrest@0: plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: system("rm "+plot_name+"-1."+plot_type_new) forrest@0: system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: ;(F) scatter model vs ob 933 forrest@0: forrest@0: plot_name = "scatter_model_vs_ob_933" forrest@0: title = model_name +" vs ob 933 sites" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: res = True ; plot mods desired forrest@0: res@tiMainString = title ; add title forrest@0: res@xyMarkLineModes = "Markers" ; choose which have markers forrest@0: res@xyMarkers = 16 ; choose type of marker forrest@0: res@xyMarkerColor = "red" ; Marker color forrest@0: res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01) forrest@0: res@tmLabelAutoStride = True ; nice tick mark labels forrest@0: ;------------------------------- forrest@0: ;compute correlation coef. and M forrest@0: ccr933 = esccr(nppmod933,npp933,0) forrest@0: print (ccr933) forrest@0: forrest@0: ;new eq forrest@0: bias = sum(abs(nppmod933-npp933)/(nppmod933+npp933)) forrest@0: M933 = (1. - (bias/dimsizes(y933)))*5. forrest@0: print (M933) forrest@0: forrest@0: tRes = True forrest@0: tRes@txFontHeightF = 0.025 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes) forrest@0: ;-------------------------------- forrest@0: plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: system("rm "+plot_name+"."+plot_type) forrest@0: system("rm "+plot_name+"-1."+plot_type_new) forrest@0: system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: ;************************************************************************** forrest@0: ;(G) histogram 81 forrest@0: forrest@0: ;-------------------------------------------------------------------- forrest@0: ;get data forrest@0: forrest@0: RAIN1_1D = ndtooned(rain81) forrest@0: RAIN2_1D = ndtooned(rainmod81) forrest@0: NPP1_1D = ndtooned(npp81) forrest@0: NPP2_1D = ndtooned(nppmod81) forrest@0: ; forrest@0: ; Calculate some "nice" bins for binning the data in equally spaced forrest@0: ; ranges. forrest@0: ; forrest@0: nbins = 15 ; Number of bins to use. forrest@0: forrest@0: nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True) forrest@0: nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1) forrest@0: range = fspan(nicevals(0),nicevals(1),nvals) forrest@0: ; forrest@0: ; Use this range information to grab all the values in a forrest@0: ; particular range, and then take an average. forrest@0: ; forrest@0: nr = dimsizes(range) forrest@0: nx = nr-1 forrest@0: xvalues = new((/2,nx/),typeof(RAIN1_1D)) forrest@0: xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2. forrest@0: dx = xvalues(0,1) - xvalues(0,0) ; range width forrest@0: dx4 = dx/4 ; 1/4 of the range forrest@0: xvalues(1,:) = xvalues(0,:) - dx/5. forrest@0: yvalues = new((/2,nx/),typeof(RAIN1_1D)) forrest@0: mn_yvalues = new((/2,nx/),typeof(RAIN1_1D)) forrest@0: mx_yvalues = new((/2,nx/),typeof(RAIN1_1D)) forrest@0: forrest@0: do nd=0,1 forrest@0: ; forrest@0: ; See if we are doing model or observational data. forrest@0: ; forrest@0: if(nd.eq.0) then forrest@0: data = RAIN1_1D forrest@0: npp_data = NPP1_1D forrest@0: else forrest@0: data = RAIN2_1D forrest@0: npp_data = NPP2_1D forrest@0: end if forrest@0: ; forrest@0: ; Loop through each range and check for values. forrest@0: ; forrest@0: do i=0,nr-2 forrest@0: if (i.ne.(nr-2)) then forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+","+range(i+1)+")") forrest@0: idx = ind((range(i).le.data).and.(data.lt.range(i+1))) forrest@0: else forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+",)") forrest@0: idx = ind(range(i).le.data) forrest@0: end if forrest@0: ; forrest@0: ; Calculate average, and get min and max. forrest@0: ; forrest@0: if(.not.any(ismissing(idx))) then forrest@0: yvalues(nd,i) = avg(npp_data(idx)) forrest@0: mn_yvalues(nd,i) = min(npp_data(idx)) forrest@0: mx_yvalues(nd,i) = max(npp_data(idx)) forrest@0: count = dimsizes(idx) forrest@0: else forrest@0: count = 0 forrest@0: yvalues(nd,i) = yvalues@_FillValue forrest@0: mn_yvalues(nd,i) = yvalues@_FillValue forrest@0: mx_yvalues(nd,i) = yvalues@_FillValue forrest@0: end if forrest@0: ; forrest@0: ; Print out information. forrest@0: ; forrest@0: ; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i)) forrest@0: ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) forrest@0: forrest@0: ; forrest@0: ; Clean up for next time in loop. forrest@0: ; forrest@0: delete(idx) forrest@0: end do forrest@0: delete(data) forrest@0: delete(npp_data) forrest@0: end do forrest@0: ;---------------------------------------- forrest@0: ;compute correlation coeff and M score forrest@0: u = yvalues(0,:) forrest@0: v = yvalues(1,:) forrest@0: forrest@0: good = ind(.not.ismissing(u) .and. .not.ismissing(v)) forrest@0: uu = u(good) forrest@0: vv = v(good) forrest@0: forrest@0: ccr81h = esccr(uu,vv,0) forrest@0: print (ccr81h) forrest@0: forrest@0: ;new eq forrest@0: bias = sum(abs(vv-uu)/(vv+uu)) forrest@0: M81h = (1.- (bias/dimsizes(uu)))*5. forrest@0: print (M81h) forrest@0: forrest@0: delete (u) forrest@0: delete (v) forrest@0: delete (uu) forrest@0: delete (vv) forrest@0: ;---------------------------------------------------------------------- forrest@0: ; histogram res forrest@0: forrest@0: res = True forrest@0: res@gsnMaximize = True forrest@0: res@gsnDraw = False forrest@0: res@gsnFrame = False forrest@0: res@xyMarkLineMode = "Markers" forrest@0: res@xyMarkerSizeF = 0.014 forrest@0: res@xyMarker = 16 forrest@0: res@xyMarkerColors = (/"Brown","Blue"/) forrest@0: res@trYMinF = min(mn_yvalues) - 10. forrest@0: res@trYMaxF = max(mx_yvalues) + 10. forrest@0: forrest@0: res@tiYAxisString = "NPP (g C/m2/year)" forrest@0: res@tiXAxisString = "Precipitation (m/year)" forrest@0: forrest@0: max_bar = new((/2,nx/),graphic) forrest@0: min_bar = new((/2,nx/),graphic) forrest@0: max_cap = new((/2,nx/),graphic) forrest@0: min_cap = new((/2,nx/),graphic) forrest@0: forrest@0: lnres = True forrest@0: line_colors = (/"brown","blue"/) forrest@0: ;================================================================= forrest@0: ; histogram ob 81 site only forrest@0: ; forrest@0: plot_name = "histogram_ob_81" forrest@0: title = "Observed 81 site" forrest@0: res@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: forrest@0: xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),res) forrest@0: forrest@0: ;------------------------------- forrest@0: ;Attach the vertical bar and the horizontal cap line forrest@0: forrest@0: do nd=0,0 forrest@0: lnres@gsLineColor = line_colors(nd) forrest@0: do i=0,nx-1 forrest@0: forrest@0: if(.not.ismissing(mn_yvalues(nd,i)).and. \ forrest@0: .not.ismissing(mx_yvalues(nd,i))) then forrest@0: ; forrest@0: ; Attach the vertical bar, both above and below the marker. forrest@0: ; forrest@0: x1 = xvalues(nd,i) forrest@0: y1 = yvalues(nd,i) forrest@0: y2 = mn_yvalues(nd,i) forrest@0: min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: forrest@0: y2 = mx_yvalues(nd,i) forrest@0: max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: ; forrest@0: ; Attach the horizontal cap line, both above and below the marker. forrest@0: ; forrest@0: x1 = xvalues(nd,i) - dx4 forrest@0: x2 = xvalues(nd,i) + dx4 forrest@0: y1 = mn_yvalues(nd,i) forrest@0: min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: forrest@0: y1 = mx_yvalues(nd,i) forrest@0: max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: end if forrest@0: end do forrest@0: end do forrest@0: forrest@0: draw(xy) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: ;=========================================================================== forrest@0: ; histogram model vs ob 81 site forrest@0: forrest@0: plot_name = "histogram_mod-ob_81" forrest@0: title = model_name+ " vs Observed 81 site" forrest@0: res@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: forrest@0: ;----------------------------- forrest@0: ; Add a boxed legend using the more simple method, which won't have forrest@0: ; vertical lines going through the markers. forrest@0: forrest@0: res@pmLegendDisplayMode = "Always" forrest@0: ; res@pmLegendWidthF = 0.1 forrest@0: res@pmLegendWidthF = 0.08 forrest@0: res@pmLegendHeightF = 0.05 forrest@0: res@pmLegendOrthogonalPosF = -1.17 forrest@0: ; res@pmLegendOrthogonalPosF = -1.00 ;(downward) forrest@0: ; res@pmLegendParallelPosF = 0.18 forrest@0: res@pmLegendParallelPosF = 0.23 ;(rightward) forrest@0: forrest@0: ; res@lgPerimOn = False forrest@0: res@lgLabelFontHeightF = 0.015 forrest@0: res@xyExplicitLegendLabels = (/"observed",model_name/) forrest@0: ;----------------------------- forrest@0: tRes = True forrest@0: tRes@txFontHeightF = 0.025 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.5,0.8,tRes) forrest@0: forrest@0: xy = gsn_csm_xy(wks,xvalues,yvalues,res) forrest@0: ;------------------------------- forrest@0: ;Attach the vertical bar and the horizontal cap line forrest@0: forrest@0: do nd=0,1 forrest@0: lnres@gsLineColor = line_colors(nd) forrest@0: do i=0,nx-1 forrest@0: forrest@0: if(.not.ismissing(mn_yvalues(nd,i)).and. \ forrest@0: .not.ismissing(mx_yvalues(nd,i))) then forrest@0: ; forrest@0: ; Attach the vertical bar, both above and below the marker. forrest@0: ; forrest@0: x1 = xvalues(nd,i) forrest@0: y1 = yvalues(nd,i) forrest@0: y2 = mn_yvalues(nd,i) forrest@0: min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: forrest@0: y2 = mx_yvalues(nd,i) forrest@0: max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: ; forrest@0: ; Attach the horizontal cap line, both above and below the marker. forrest@0: ; forrest@0: x1 = xvalues(nd,i) - dx4 forrest@0: x2 = xvalues(nd,i) + dx4 forrest@0: y1 = mn_yvalues(nd,i) forrest@0: min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: forrest@0: y1 = mx_yvalues(nd,i) forrest@0: max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: end if forrest@0: end do forrest@0: end do forrest@0: draw(xy) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: ;system("rm "+plot_name+"."+plot_type) forrest@0: ;system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ;system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: delete (RAIN1_1D) forrest@0: delete (RAIN2_1D) forrest@0: delete (NPP1_1D) forrest@0: delete (NPP2_1D) forrest@0: delete (range) forrest@0: delete (xvalues) forrest@0: delete (yvalues) forrest@0: delete (mn_yvalues) forrest@0: delete (mx_yvalues) forrest@0: delete (good) forrest@0: delete (max_bar) forrest@0: delete (min_bar) forrest@0: delete (max_cap) forrest@0: delete (min_cap) forrest@0: ;************************************************************************** forrest@0: ;(H) histogram 933 forrest@0: forrest@0: ;-------------------------------------------------------------------- forrest@0: ;get data forrest@0: forrest@0: RAIN1_1D = ndtooned(rain933) forrest@0: RAIN2_1D = ndtooned(rainmod933) forrest@0: NPP1_1D = ndtooned(npp933) forrest@0: NPP2_1D = ndtooned(nppmod933) forrest@0: ; forrest@0: ; Calculate some "nice" bins for binning the data in equally spaced forrest@0: ; ranges. forrest@0: ; forrest@0: nbins = 15 ; Number of bins to use. forrest@0: forrest@0: nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True) forrest@0: nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1) forrest@0: range = fspan(nicevals(0),nicevals(1),nvals) forrest@0: ; forrest@0: ; Use this range information to grab all the values in a forrest@0: ; particular range, and then take an average. forrest@0: ; forrest@0: nr = dimsizes(range) forrest@0: nx = nr-1 forrest@0: xvalues = new((/2,nx/),typeof(RAIN1_1D)) forrest@0: xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2. forrest@0: dx = xvalues(0,1) - xvalues(0,0) ; range width forrest@0: dx4 = dx/4 ; 1/4 of the range forrest@0: xvalues(1,:) = xvalues(0,:) - dx/5. forrest@0: yvalues = new((/2,nx/),typeof(RAIN1_1D)) forrest@0: mn_yvalues = new((/2,nx/),typeof(RAIN1_1D)) forrest@0: mx_yvalues = new((/2,nx/),typeof(RAIN1_1D)) forrest@0: forrest@0: do nd=0,1 forrest@0: ; forrest@0: ; See if we are doing model or observational data. forrest@0: ; forrest@0: if(nd.eq.0) then forrest@0: data = RAIN1_1D forrest@0: npp_data = NPP1_1D forrest@0: else forrest@0: data = RAIN2_1D forrest@0: npp_data = NPP2_1D forrest@0: end if forrest@0: ; forrest@0: ; Loop through each range and check for values. forrest@0: ; forrest@0: do i=0,nr-2 forrest@0: if (i.ne.(nr-2)) then forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+","+range(i+1)+")") forrest@0: idx = ind((range(i).le.data).and.(data.lt.range(i+1))) forrest@0: else forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+",)") forrest@0: idx = ind(range(i).le.data) forrest@0: end if forrest@0: ; forrest@0: ; Calculate average, and get min and max. forrest@0: ; forrest@0: if(.not.any(ismissing(idx))) then forrest@0: yvalues(nd,i) = avg(npp_data(idx)) forrest@0: mn_yvalues(nd,i) = min(npp_data(idx)) forrest@0: mx_yvalues(nd,i) = max(npp_data(idx)) forrest@0: count = dimsizes(idx) forrest@0: else forrest@0: count = 0 forrest@0: yvalues(nd,i) = yvalues@_FillValue forrest@0: mn_yvalues(nd,i) = yvalues@_FillValue forrest@0: mx_yvalues(nd,i) = yvalues@_FillValue forrest@0: end if forrest@0: ; forrest@0: ; Print out information. forrest@0: ; forrest@0: ; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i)) forrest@0: ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) forrest@0: forrest@0: ; forrest@0: ; Clean up for next time in loop. forrest@0: ; forrest@0: delete(idx) forrest@0: end do forrest@0: delete(data) forrest@0: delete(npp_data) forrest@0: end do forrest@0: ;---------------------------------------- forrest@0: ;compute correlation coeff and M score forrest@0: u = yvalues(0,:) forrest@0: v = yvalues(1,:) forrest@0: forrest@0: good = ind(.not.ismissing(u) .and. .not.ismissing(v)) forrest@0: uu = u(good) forrest@0: vv = v(good) forrest@0: forrest@0: ccr933h = esccr(uu,vv,0) forrest@0: print (ccr933h) forrest@0: forrest@0: ;new eq forrest@0: bias = sum(abs(vv-uu)/(vv+uu)) forrest@0: M933h = (1.- (bias/dimsizes(uu)))*5. forrest@0: print (M933h) forrest@0: forrest@0: delete (u) forrest@0: delete (v) forrest@0: delete (uu) forrest@0: delete (vv) forrest@0: ;---------------------------------------------------------------------- forrest@0: ; histogram res forrest@0: forrest@0: res = True forrest@0: res@gsnMaximize = True forrest@0: res@gsnDraw = False forrest@0: res@gsnFrame = False forrest@0: res@xyMarkLineMode = "Markers" forrest@0: res@xyMarkerSizeF = 0.014 forrest@0: res@xyMarker = 16 forrest@0: res@xyMarkerColors = (/"Brown","Blue"/) forrest@0: res@trYMinF = min(mn_yvalues) - 10. forrest@0: res@trYMaxF = max(mx_yvalues) + 10. forrest@0: forrest@0: res@tiYAxisString = "NPP (g C/m2/year)" forrest@0: res@tiXAxisString = "Precipitation (m/year)" forrest@0: forrest@0: max_bar = new((/2,nx/),graphic) forrest@0: min_bar = new((/2,nx/),graphic) forrest@0: max_cap = new((/2,nx/),graphic) forrest@0: min_cap = new((/2,nx/),graphic) forrest@0: forrest@0: lnres = True forrest@0: line_colors = (/"brown","blue"/) forrest@0: ;================================================================= forrest@0: ; histogram ob 933 site only forrest@0: ; forrest@0: plot_name = "histogram_ob_933" forrest@0: title = "Observed 933 site" forrest@0: res@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: forrest@0: xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),res) forrest@0: forrest@0: ;------------------------------- forrest@0: ;Attach the vertical bar and the horizontal cap line forrest@0: forrest@0: do nd=0,0 forrest@0: lnres@gsLineColor = line_colors(nd) forrest@0: do i=0,nx-1 forrest@0: forrest@0: if(.not.ismissing(mn_yvalues(nd,i)).and. \ forrest@0: .not.ismissing(mx_yvalues(nd,i))) then forrest@0: ; forrest@0: ; Attach the vertical bar, both above and below the marker. forrest@0: ; forrest@0: x1 = xvalues(nd,i) forrest@0: y1 = yvalues(nd,i) forrest@0: y2 = mn_yvalues(nd,i) forrest@0: min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: forrest@0: y2 = mx_yvalues(nd,i) forrest@0: max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: ; forrest@0: ; Attach the horizontal cap line, both above and below the marker. forrest@0: ; forrest@0: x1 = xvalues(nd,i) - dx4 forrest@0: x2 = xvalues(nd,i) + dx4 forrest@0: y1 = mn_yvalues(nd,i) forrest@0: min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: forrest@0: y1 = mx_yvalues(nd,i) forrest@0: max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: end if forrest@0: end do forrest@0: end do forrest@0: forrest@0: draw(xy) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: ;=========================================================================== forrest@0: ; histogram model vs ob 933 site forrest@0: forrest@0: plot_name = "histogram_mod-ob_933" forrest@0: title = model_name+ " vs Observed 933 site" forrest@0: res@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: forrest@0: ;----------------------------- forrest@0: ; Add a boxed legend using the more simple method, which won't have forrest@0: ; vertical lines going through the markers. forrest@0: forrest@0: res@pmLegendDisplayMode = "Always" forrest@0: ; res@pmLegendWidthF = 0.1 forrest@0: res@pmLegendWidthF = 0.08 forrest@0: res@pmLegendHeightF = 0.05 forrest@0: res@pmLegendOrthogonalPosF = -1.17 forrest@0: ; res@pmLegendOrthogonalPosF = -1.00 ;(downward) forrest@0: ; res@pmLegendParallelPosF = 0.18 forrest@0: res@pmLegendParallelPosF = 0.23 ;(rightward) forrest@0: forrest@0: ; res@lgPerimOn = False forrest@0: res@lgLabelFontHeightF = 0.015 forrest@0: res@xyExplicitLegendLabels = (/"observed",model_name/) forrest@0: ;----------------------------- forrest@0: tRes = True forrest@0: tRes@txFontHeightF = 0.025 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.5,0.8,tRes) forrest@0: forrest@0: xy = gsn_csm_xy(wks,xvalues,yvalues,res) forrest@0: ;------------------------------- forrest@0: ;Attach the vertical bar and the horizontal cap line forrest@0: forrest@0: do nd=0,1 forrest@0: lnres@gsLineColor = line_colors(nd) forrest@0: do i=0,nx-1 forrest@0: forrest@0: if(.not.ismissing(mn_yvalues(nd,i)).and. \ forrest@0: .not.ismissing(mx_yvalues(nd,i))) then forrest@0: ; forrest@0: ; Attach the vertical bar, both above and below the marker. forrest@0: ; forrest@0: x1 = xvalues(nd,i) forrest@0: y1 = yvalues(nd,i) forrest@0: y2 = mn_yvalues(nd,i) forrest@0: min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: forrest@0: y2 = mx_yvalues(nd,i) forrest@0: max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: ; forrest@0: ; Attach the horizontal cap line, both above and below the marker. forrest@0: ; forrest@0: x1 = xvalues(nd,i) - dx4 forrest@0: x2 = xvalues(nd,i) + dx4 forrest@0: y1 = mn_yvalues(nd,i) forrest@0: min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: forrest@0: y1 = mx_yvalues(nd,i) forrest@0: max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: end if forrest@0: end do forrest@0: end do forrest@0: draw(xy) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) forrest@0: ;system("rm "+plot_name+"."+plot_type) forrest@0: ;system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ;system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: ;------------------------------------------------------------------------ forrest@0: forrest@0: end