1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lai/27.histogram+bias_grow.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,347 @@
1.4 +;********************************************************
1.5 +; histogram normalized by rain and compute correleration
1.6 +;********************************************************
1.7 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.10 +
1.11 +procedure pminmax(data:numeric,name:string)
1.12 +begin
1.13 + print ("min/max " + name + " = " + min(data) + "/" + max(data))
1.14 + if(isatt(data,"units")) then
1.15 + print (name + " units = " + data@units)
1.16 + end if
1.17 +end
1.18 +
1.19 +;
1.20 +; Main code.
1.21 +;
1.22 +begin
1.23 +
1.24 +;nclass = 18
1.25 + nclass = 20
1.26 + day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/)
1.27 +
1.28 +;************************************************
1.29 +; read in data: observed
1.30 +;************************************************
1.31 + diri1 = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
1.32 +;fili1 = "land_class_T42.nc"
1.33 + fili1 = "land_class_T42_new.nc"
1.34 + fili2 = "LAI_2000-2005_ensemble_T42.nc"
1.35 + data_file_ob1 = addfile(diri1+fili1,"r")
1.36 + data_file_ob2 = addfile(diri1+fili2,"r")
1.37 +
1.38 + RAIN1 = tofloat(data_file_ob1->LAND_CLASS)
1.39 +
1.40 + z = data_file_ob2->LAI
1.41 + y = z(0,:,:)
1.42 + y@long_name = "Days of Growing Season"
1.43 +
1.44 + dsizes_z = dimsizes(z)
1.45 + ntime = dsizes_z(0)
1.46 + nlat = dsizes_z(1)
1.47 + nlon = dsizes_z(2)
1.48 +
1.49 + do j = 0,nlat-1
1.50 + do i = 0,nlon-1
1.51 + nday = 0.
1.52 + do k = 0,ntime-1
1.53 + if (.not. ismissing(z(k,j,i)) .and. z(k,j,i) .gt. 1.0) then
1.54 + nday = nday + day_of_data(k)
1.55 + end if
1.56 + end do
1.57 + y(j,i) = nday
1.58 + end do
1.59 + end do
1.60 +
1.61 + print (min(y)+"/"+max(y))
1.62 +
1.63 + NPP1 = y
1.64 +
1.65 + delete (z)
1.66 + delete (y)
1.67 +;************************************************
1.68 +; read in data: model
1.69 +;************************************************
1.70 + diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.71 +;fili3 = "i01.03cn_1545-1569_MONS_climo.nc"
1.72 + fili3 = "i01.04casa_1605-1629_MONS_climo.nc"
1.73 + data_file_model = addfile(diri2+fili3,"r")
1.74 +
1.75 + z = data_file_model->TLAI
1.76 + y = z(0,:,:)
1.77 + y@long_name = "Days of Growing Season"
1.78 +
1.79 + dsizes_z = dimsizes(z)
1.80 + ntime = dsizes_z(0)
1.81 + nlat = dsizes_z(1)
1.82 + nlon = dsizes_z(2)
1.83 +
1.84 + do j = 0,nlat-1
1.85 + do i = 0,nlon-1
1.86 + nday = 0.
1.87 + do k = 0,ntime-1
1.88 + if (.not. ismissing(z(k,j,i)) .and. z(k,j,i) .gt. 1.0) then
1.89 + nday = nday + day_of_data(k)
1.90 + end if
1.91 + end do
1.92 + y(j,i) = nday
1.93 + end do
1.94 + end do
1.95 +
1.96 + print (min(y)+"/"+max(y))
1.97 +
1.98 + NPP2 = y
1.99 +
1.100 + delete (z)
1.101 + delete (y)
1.102 +;************************************************
1.103 +; print min/max and unit
1.104 +;************************************************
1.105 + pminmax(RAIN1,"RAIN1")
1.106 + pminmax(NPP1,"NPP1")
1.107 + pminmax(NPP2,"NPP2")
1.108 +
1.109 + RAIN1_1D = ndtooned(RAIN1)
1.110 + NPP1_1D = ndtooned(NPP1)
1.111 + NPP2_1D = ndtooned(NPP2)
1.112 +;
1.113 +; Calculate some "nice" bins for binning the data in equally spaced
1.114 +; ranges.
1.115 +;
1.116 +
1.117 +; nbins = nclass + 1 ; Number of bins to use.
1.118 +; nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,False)
1.119 +; nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
1.120 +; range = fspan(nicevals(0),nicevals(1),nvals)
1.121 +
1.122 + nclass = nclass + 1
1.123 + range = fspan(0,nclass-1,nclass)
1.124 +
1.125 +; print (nicevals)
1.126 +; print (nvals)
1.127 + print (range)
1.128 +; exit
1.129 +
1.130 +;
1.131 +; Use this range information to grab all the values in a
1.132 +; particular range, and then take an average.
1.133 +;
1.134 + nr = dimsizes(range)
1.135 + nx = nr-1
1.136 + xvalues = new((/2,nx/),typeof(RAIN1_1D))
1.137 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.138 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.139 + dx4 = dx/4 ; 1/4 of the range
1.140 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.141 + yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.142 + mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.143 + mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.144 +
1.145 + do nd=0,1
1.146 +;
1.147 +; See if we are doing model or observational data.
1.148 +;
1.149 + if(nd.eq.0) then
1.150 + data = RAIN1_1D
1.151 + npp_data = NPP1_1D
1.152 + else
1.153 + data = RAIN1_1D
1.154 + npp_data = NPP2_1D
1.155 + end if
1.156 +;
1.157 +; Loop through each range and check for values.
1.158 +;
1.159 + do i=0,nr-2
1.160 + if (i.ne.(nr-2)) then
1.161 + print("")
1.162 + print("In range ["+range(i)+","+range(i+1)+")")
1.163 + idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
1.164 + else
1.165 + print("")
1.166 + print("In range ["+range(i)+",)")
1.167 + idx = ind(range(i).le.data)
1.168 + end if
1.169 +;
1.170 +; Calculate average, and get min and max.
1.171 +;
1.172 + if(.not.any(ismissing(idx))) then
1.173 + yvalues(nd,i) = avg(npp_data(idx))
1.174 + mn_yvalues(nd,i) = min(npp_data(idx))
1.175 + mx_yvalues(nd,i) = max(npp_data(idx))
1.176 + count = dimsizes(idx)
1.177 + else
1.178 + count = 0
1.179 + yvalues(nd,i) = yvalues@_FillValue
1.180 + mn_yvalues(nd,i) = yvalues@_FillValue
1.181 + mx_yvalues(nd,i) = yvalues@_FillValue
1.182 + end if
1.183 +;
1.184 +; Print out information.
1.185 +;
1.186 + print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1.187 + print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.188 +
1.189 +;
1.190 +; Clean up for next time in loop.
1.191 +;
1.192 + delete(idx)
1.193 + end do
1.194 + delete(data)
1.195 + delete(npp_data)
1.196 + end do
1.197 +
1.198 +;
1.199 +; Start the graphics.
1.200 +;
1.201 + wks = gsn_open_wks("ps","xy")
1.202 +
1.203 + res = True
1.204 + res@gsnMaximize = True
1.205 + res@gsnDraw = False
1.206 + res@gsnFrame = False
1.207 + res@xyMarkLineMode = "Markers"
1.208 + res@xyMarkerSizeF = 0.014
1.209 + res@xyMarker = 16
1.210 + res@xyMarkerColors = (/"Brown","Blue"/)
1.211 +; res@trYMinF = min(mn_yvalues) - 10.
1.212 +; res@trYMaxF = max(mx_yvalues) + 10.
1.213 + res@trYMinF = min(mn_yvalues) - 20.
1.214 + res@trYMaxF = max(mx_yvalues) + 100.
1.215 +
1.216 +; res@tiMainString = "Observed vs i01.03cn"
1.217 + res@tiMainString = "Observed vs i01.04casa"
1.218 +
1.219 + res@tiYAxisString = "Days of Growing season"
1.220 + res@tiXAxisString = "Land Cover Type"
1.221 +;
1.222 +; Add a boxed legend using the more simple method, which won't have
1.223 +; vertical lines going through the markers.
1.224 +;
1.225 + res@pmLegendDisplayMode = "Always"
1.226 +; res@pmLegendWidthF = 0.1
1.227 + res@pmLegendWidthF = 0.08
1.228 + res@pmLegendHeightF = 0.05
1.229 + res@pmLegendOrthogonalPosF = -1.17
1.230 +; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.231 +; res@pmLegendParallelPosF = 0.18
1.232 + res@pmLegendParallelPosF = 0.23 ;(rightward)
1.233 +
1.234 +; res@lgPerimOn = False
1.235 + res@lgLabelFontHeightF = 0.015
1.236 +; res@xyExplicitLegendLabels = (/"observed","model_i01.03cn"/)
1.237 + res@xyExplicitLegendLabels = (/"observed","model_i01.04casa"/)
1.238 +
1.239 + xy = gsn_csm_xy(wks,xvalues,yvalues,res)
1.240 +
1.241 + max_bar = new((/2,nx/),graphic)
1.242 + min_bar = new((/2,nx/),graphic)
1.243 + max_cap = new((/2,nx/),graphic)
1.244 + min_cap = new((/2,nx/),graphic)
1.245 +
1.246 + lnres = True
1.247 +
1.248 + line_colors = (/"brown","blue"/)
1.249 + do nd=0,1
1.250 + lnres@gsLineColor = line_colors(nd)
1.251 + do i=0,nx-1
1.252 +
1.253 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.254 + .not.ismissing(mx_yvalues(nd,i))) then
1.255 +;
1.256 +; Attach the vertical bar, both above and below the marker.
1.257 +;
1.258 + x1 = xvalues(nd,i)
1.259 + y1 = yvalues(nd,i)
1.260 + y2 = mn_yvalues(nd,i)
1.261 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.262 +
1.263 + y2 = mx_yvalues(nd,i)
1.264 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.265 +;
1.266 +; Attach the horizontal cap line, both above and below the marker.
1.267 +;
1.268 + x1 = xvalues(nd,i) - dx4
1.269 + x2 = xvalues(nd,i) + dx4
1.270 + y1 = mn_yvalues(nd,i)
1.271 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.272 +
1.273 + y1 = mx_yvalues(nd,i)
1.274 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.275 + end if
1.276 + end do
1.277 + end do
1.278 +
1.279 +;
1.280 +; Here's how to do the legend by hand.
1.281 +;
1.282 +; mkres = True ; Marker resources
1.283 +; txres = True ; Text resources
1.284 +; mkres@gsMarkerIndex = 16
1.285 +; mkres@gsMarkerSizeF = 0.02
1.286 +; txres@txFontHeightF = 0.02
1.287 +; txres@txJust = "CenterLeft"
1.288 +;
1.289 +; Change these values if you want to move the marker legend location.
1.290 +; These values are in the same data space as the plot.
1.291 +;
1.292 +; xlg1_cen = 0.2
1.293 +; ylg1_cen = 900.
1.294 +
1.295 +; xlg2_cen = 0.2
1.296 +; ylg2_cen = 760.
1.297 +
1.298 +; mkres@gsMarkerColor = "brown"
1.299 +; lnres@gsLineColor = "brown"
1.300 +
1.301 +; lg_mark_legend1 = gsn_add_polymarker(wks,xy,xlg1_cen,ylg1_cen,mkres)
1.302 +; lg_line_legend1 = gsn_add_polyline(wks,xy,(/xlg1_cen,xlg1_cen/), \
1.303 +; (/ylg1_cen-60,ylg1_cen+60/),lnres)
1.304 +; lg_cap_legend11 = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \
1.305 +; (/ylg1_cen-60,ylg1_cen-60/),lnres)
1.306 +; lg_cap_legend12 = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \
1.307 +; (/ylg1_cen+60,ylg1_cen+60/),lnres)
1.308 +
1.309 +; tx_legend1 = gsn_add_text(wks,xy,"observed",xlg1_cen+0.15,ylg1_cen,txres)
1.310 +
1.311 +; mkres@gsMarkerColor = "blue"
1.312 +; lnres@gsLineColor = "blue"
1.313 +
1.314 +; lg_mark_legend2 = gsn_add_polymarker(wks,xy,xlg2_cen,ylg2_cen,mkres)
1.315 +; lg_line_legend2 = gsn_add_polyline(wks,xy,(/xlg2_cen,xlg2_cen/), \
1.316 +; (/ylg2_cen-60,ylg2_cen+60/),lnres)
1.317 +; lg_cap_legend21 = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \
1.318 +; (/ylg2_cen-60,ylg2_cen-60/),lnres)
1.319 +; lg_cap_legend22 = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \
1.320 +; (/ylg2_cen+60,ylg2_cen+60/),lnres)
1.321 +; tx_legend2 = gsn_add_text(wks,xy,"model_i01.03cn",xlg2_cen+0.15,ylg2_cen,txres)
1.322 +
1.323 + draw(xy)
1.324 + frame(wks)
1.325 + system("convert xy.ps xy.png")
1.326 +
1.327 + u = yvalues(0,:)
1.328 + v = yvalues(1,:)
1.329 +
1.330 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.331 + uu = u(good)
1.332 + vv = v(good)
1.333 + nz = dimsizes(uu)
1.334 + print (nz)
1.335 +
1.336 + ccr = esccr(uu,vv,0)
1.337 + print (ccr)
1.338 +;old eq
1.339 +;bias = sum(((vv-uu)/uu)^2)
1.340 +;M = (1.- sqrt(bias/nz))*5.
1.341 +
1.342 +;new eq
1.343 + bias = sum(abs(vv-uu)/(vv+uu))
1.344 + M = (1.- (bias/nz))*5.
1.345 +
1.346 + print (bias)
1.347 + print (M)
1.348 +
1.349 +end
1.350 +