diff -r 000000000000 -r 0c6405ab2ff4 npp/24.histogram+bias.ncl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/npp/24.histogram+bias.ncl Mon Jan 26 22:08:20 2009 -0500 @@ -0,0 +1,329 @@ +;******************************************************** +; histogram normalized by rain and compute correleration +;******************************************************** +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 + data_types = (/ "Obs", "Model" /) + data_names = (/ "data.81.nc" , "i01.03cn_1545-1569_ANN_climo.nc" /) +; data_names = (/ "data.933.nc", "i01.03cn_1545-1569_ANN_climo.nc" /) +; data_names = (/ "data.81.nc" , "i01.04casa_1605-1629_ANN_climo.nc" /) +; data_names = (/ "data.933.nc", "i01.04casa_1605-1629_ANN_climo.nc" /) + filevar_names = (/ (/"PREC_ANN","TNPP_C"/), (/"RAIN","NPP"/) /) + ndata_types = dimsizes(data_types) + + data_file_obs = addfile(data_names(0),"r") ; Open obs file + data_file_mod = addfile(data_names(1),"r") ; Open model file + +;************************************************ +; read in data: observed +;************************************************ + RAIN1 = tofloat(data_file_obs->PREC_ANN) ; for data.81 +;RAIN1 = data_file_obs->PREC ; for data.933 + NPP1 = data_file_obs->TNPP_C + xo = data_file_obs->LONG_DD + yo = data_file_obs->LAT_DD + +; change longitude from (-180,180) to (0,360) + print (xo) + nx = dimsizes(xo) + do i= 0,nx-1 + if (xo(i) .lt. 0.) then + xo(i) = xo(i)+ 360. + end if + end do + print (xo) + +;************************************************ +; read in data: model +;************************************************ + ai = data_file_mod->RAIN + bi = data_file_mod->NPP + xi = data_file_mod->lon + yi = data_file_mod->lat + +;************************************************ +; interpolate from model grid to observed grid +;************************************************ + RAIN2 = linint2_points(xi,yi,ai,True,xo,yo,0) + NPP2 = linint2_points(xi,yi,bi,True,xo,yo,0) + +;************************************************ +; convert unit +;************************************************ +; Units for these four variables are: +; +; RAIN1 : mm/year +; RAIN2 : mm/s +; NPP1 : g C/m^2/year +; NPP2 : g C/m^2/s +; +; We want to convert these to "m/year" and "g C/m^2/year". +; + nsec_per_year = 60*60*24*365 ; # seconds per year + +; Do the necessary conversions. + RAIN1 = RAIN1 / 1000. + RAIN2 = (RAIN2/ 1000.) * nsec_per_year + NPP2 = NPP2 * nsec_per_year + +; Redo the units. + RAIN1@units = "m/yr" + RAIN2@units = "m/yr" + NPP1@units = "gC/m^2/yr" + NPP2@units = "gC/m^2/yr" + +;************************************************ +; print min/max and unit +;************************************************ + pminmax(RAIN1,"RAIN1") + pminmax(RAIN2,"RAIN2") + pminmax(NPP1,"NPP1") + pminmax(NPP2,"NPP2") + + RAIN1_1D = ndtooned(RAIN1) + RAIN2_1D = ndtooned(RAIN2) + NPP1_1D = ndtooned(NPP1) + NPP2_1D = ndtooned(NPP2) +; +; Calculate some "nice" bins for binning the data in equally spaced +; ranges. +; + nbins = 15 ; Number of bins to use. + + nicevals = nice_mnmxintvl(min(RAIN2_1D),max(RAIN2_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(RAIN2_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(RAIN2_1D)) + mn_yvalues = new((/2,nx/),typeof(RAIN2_1D)) + mx_yvalues = new((/2,nx/),typeof(RAIN2_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 + +; +; Start the graphics. +; + wks = gsn_open_wks("png","xy") + + res = True + res@gsnMaximize = True + res@gsnDraw = False + res@gsnFrame = False + res@xyMarkLineMode = "Markers" + res@xyMarkerSizeF = 0.014 + res@xyMarker = 16 + res@xyMarkerColors = (/"Brown","Blue"/) + res@trYMinF = min(mn_yvalues) - 10. + res@trYMaxF = max(mx_yvalues) + 10. + +; res@tiMainString = "Observed vs i01.03cn 81 site" +; res@tiMainString = "Observed vs i01.03cn 933 site" + res@tiMainString = "Observed vs i01.04casa 81 site" +; res@tiMainString = "Observed vs i01.04casa 933 site" + res@tiYAxisString = "NPP (g C/m2/year)" + res@tiXAxisString = "Precipitation (m/year)" + +; +; Add a boxed legend using the more simple method, which won't have +; vertical lines going through the markers. +; + res@pmLegendDisplayMode = "Always" +; res@pmLegendWidthF = 0.1 + res@pmLegendWidthF = 0.08 + res@pmLegendHeightF = 0.05 + res@pmLegendOrthogonalPosF = -1.17 +; res@pmLegendOrthogonalPosF = -1.00 ;(downward) +; res@pmLegendParallelPosF = 0.18 + res@pmLegendParallelPosF = 0.23 ;(rightward) + +; res@lgPerimOn = False + res@lgLabelFontHeightF = 0.015 +; res@xyExplicitLegendLabels = (/"observed","model_i01.03cn"/) + res@xyExplicitLegendLabels = (/"observed","model_i01.04casa"/) + + xy = gsn_csm_xy(wks,xvalues,yvalues,res) + + 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"/) + 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 + +; +; Here's how to do the legend by hand. +; +; mkres = True ; Marker resources +; txres = True ; Text resources +; mkres@gsMarkerIndex = 16 +; mkres@gsMarkerSizeF = 0.02 +; txres@txFontHeightF = 0.02 +; txres@txJust = "CenterLeft" +; +; Change these values if you want to move the marker legend location. +; These values are in the same data space as the plot. +; +; xlg1_cen = 0.2 +; ylg1_cen = 900. + +; xlg2_cen = 0.2 +; ylg2_cen = 760. + +; mkres@gsMarkerColor = "brown" +; lnres@gsLineColor = "brown" + +; lg_mark_legend1 = gsn_add_polymarker(wks,xy,xlg1_cen,ylg1_cen,mkres) +; lg_line_legend1 = gsn_add_polyline(wks,xy,(/xlg1_cen,xlg1_cen/), \ +; (/ylg1_cen-60,ylg1_cen+60/),lnres) +; lg_cap_legend11 = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \ +; (/ylg1_cen-60,ylg1_cen-60/),lnres) +; lg_cap_legend12 = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \ +; (/ylg1_cen+60,ylg1_cen+60/),lnres) + +; tx_legend1 = gsn_add_text(wks,xy,"observed",xlg1_cen+0.15,ylg1_cen,txres) + +; mkres@gsMarkerColor = "blue" +; lnres@gsLineColor = "blue" + +; lg_mark_legend2 = gsn_add_polymarker(wks,xy,xlg2_cen,ylg2_cen,mkres) +; lg_line_legend2 = gsn_add_polyline(wks,xy,(/xlg2_cen,xlg2_cen/), \ +; (/ylg2_cen-60,ylg2_cen+60/),lnres) +; lg_cap_legend21 = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \ +; (/ylg2_cen-60,ylg2_cen-60/),lnres) +; lg_cap_legend22 = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \ +; (/ylg2_cen+60,ylg2_cen+60/),lnres) +; tx_legend2 = gsn_add_text(wks,xy,"model_i01.03cn",xlg2_cen+0.15,ylg2_cen,txres) + + draw(xy) + frame(wks) + + u = yvalues(0,:) + v = yvalues(1,:) + + good = ind(.not.ismissing(u) .and. .not.ismissing(v)) + uu = u(good) + vv = v(good) + nz = dimsizes(uu) + print (nz) + + ccr = esccr(uu,vv,0) + print (ccr) + +;old eq +;bias = sum(((vv-uu)/uu)^2) +;M = (1.- sqrt(bias/nz))*5. + +;new eq + bias = sum(abs(vv-uu)/(vv+uu)) + M = (1.- (bias/nz))*5. + print (bias) + print (M) + +end +