diff -r 000000000000 -r 0c6405ab2ff4 lai/27.histogram+bias_grow.ncl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lai/27.histogram+bias_grow.ncl Mon Jan 26 22:08:20 2009 -0500 @@ -0,0 +1,347 @@ +;******************************************************** +; 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 + +;nclass = 18 + nclass = 20 + day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/) + +;************************************************ +; read in data: observed +;************************************************ + diri1 = "/fis/cgd/cseg/people/jeff/clamp_data/lai/" +;fili1 = "land_class_T42.nc" + fili1 = "land_class_T42_new.nc" + fili2 = "LAI_2000-2005_ensemble_T42.nc" + data_file_ob1 = addfile(diri1+fili1,"r") + data_file_ob2 = addfile(diri1+fili2,"r") + + RAIN1 = tofloat(data_file_ob1->LAND_CLASS) + + z = data_file_ob2->LAI + y = z(0,:,:) + y@long_name = "Days of Growing Season" + + dsizes_z = dimsizes(z) + ntime = dsizes_z(0) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + nday = 0. + do k = 0,ntime-1 + if (.not. ismissing(z(k,j,i)) .and. z(k,j,i) .gt. 1.0) then + nday = nday + day_of_data(k) + end if + end do + y(j,i) = nday + end do + end do + + print (min(y)+"/"+max(y)) + + NPP1 = y + + delete (z) + delete (y) +;************************************************ +; read in data: model +;************************************************ + diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/" +;fili3 = "i01.03cn_1545-1569_MONS_climo.nc" + fili3 = "i01.04casa_1605-1629_MONS_climo.nc" + data_file_model = addfile(diri2+fili3,"r") + + z = data_file_model->TLAI + y = z(0,:,:) + y@long_name = "Days of Growing Season" + + dsizes_z = dimsizes(z) + ntime = dsizes_z(0) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + nday = 0. + do k = 0,ntime-1 + if (.not. ismissing(z(k,j,i)) .and. z(k,j,i) .gt. 1.0) then + nday = nday + day_of_data(k) + end if + end do + y(j,i) = nday + end do + end do + + print (min(y)+"/"+max(y)) + + NPP2 = y + + delete (z) + delete (y) +;************************************************ +; print min/max and unit +;************************************************ + pminmax(RAIN1,"RAIN1") + pminmax(NPP1,"NPP1") + pminmax(NPP2,"NPP2") + + RAIN1_1D = ndtooned(RAIN1) + NPP1_1D = ndtooned(NPP1) + NPP2_1D = ndtooned(NPP2) +; +; Calculate some "nice" bins for binning the data in equally spaced +; ranges. +; + +; nbins = nclass + 1 ; Number of bins to use. +; nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,False) +; nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1) +; range = fspan(nicevals(0),nicevals(1),nvals) + + nclass = nclass + 1 + range = fspan(0,nclass-1,nclass) + +; print (nicevals) +; print (nvals) + print (range) +; exit + +; +; 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 = RAIN1_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(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("ps","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@trYMinF = min(mn_yvalues) - 20. + res@trYMaxF = max(mx_yvalues) + 100. + +; res@tiMainString = "Observed vs i01.03cn" + res@tiMainString = "Observed vs i01.04casa" + + res@tiYAxisString = "Days of Growing season" + res@tiXAxisString = "Land Cover Type" +; +; 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) + system("convert xy.ps xy.png") + + 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 +