1 ;********************************************************
2 ; histogram normalized by rain and compute correleration
3 ;********************************************************
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
8 procedure pminmax(data:numeric,name:string)
10 print ("min/max " + name + " = " + min(data) + "/" + max(data))
11 if(isatt(data,"units")) then
12 print (name + " units = " + data@units)
23 day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/)
25 ;************************************************
26 ; read in data: observed
27 ;************************************************
28 diri1 = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
29 ;fili1 = "land_class_T42.nc"
30 fili1 = "land_class_T42_new.nc"
31 fili2 = "LAI_2000-2005_ensemble_T42.nc"
32 data_file_ob1 = addfile(diri1+fili1,"r")
33 data_file_ob2 = addfile(diri1+fili2,"r")
35 RAIN1 = tofloat(data_file_ob1->LAND_CLASS)
37 z = data_file_ob2->LAI
39 y@long_name = "Days of Growing Season"
41 dsizes_z = dimsizes(z)
50 if (.not. ismissing(z(k,j,i)) .and. z(k,j,i) .gt. 1.0) then
51 nday = nday + day_of_data(k)
58 print (min(y)+"/"+max(y))
64 ;************************************************
66 ;************************************************
67 diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
68 ;fili3 = "i01.03cn_1545-1569_MONS_climo.nc"
69 fili3 = "i01.04casa_1605-1629_MONS_climo.nc"
70 data_file_model = addfile(diri2+fili3,"r")
72 z = data_file_model->TLAI
74 y@long_name = "Days of Growing Season"
76 dsizes_z = dimsizes(z)
85 if (.not. ismissing(z(k,j,i)) .and. z(k,j,i) .gt. 1.0) then
86 nday = nday + day_of_data(k)
93 print (min(y)+"/"+max(y))
99 ;************************************************
100 ; print min/max and unit
101 ;************************************************
102 pminmax(RAIN1,"RAIN1")
106 RAIN1_1D = ndtooned(RAIN1)
107 NPP1_1D = ndtooned(NPP1)
108 NPP2_1D = ndtooned(NPP2)
110 ; Calculate some "nice" bins for binning the data in equally spaced
114 ; nbins = nclass + 1 ; Number of bins to use.
115 ; nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,False)
116 ; nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
117 ; range = fspan(nicevals(0),nicevals(1),nvals)
120 range = fspan(0,nclass-1,nclass)
128 ; Use this range information to grab all the values in a
129 ; particular range, and then take an average.
133 xvalues = new((/2,nx/),typeof(RAIN1_1D))
134 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
135 dx = xvalues(0,1) - xvalues(0,0) ; range width
136 dx4 = dx/4 ; 1/4 of the range
137 xvalues(1,:) = xvalues(0,:) - dx/5.
138 yvalues = new((/2,nx/),typeof(RAIN1_1D))
139 mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
140 mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
144 ; See if we are doing model or observational data.
154 ; Loop through each range and check for values.
157 if (i.ne.(nr-2)) then
159 print("In range ["+range(i)+","+range(i+1)+")")
160 idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
163 print("In range ["+range(i)+",)")
164 idx = ind(range(i).le.data)
167 ; Calculate average, and get min and max.
169 if(.not.any(ismissing(idx))) then
170 yvalues(nd,i) = avg(npp_data(idx))
171 mn_yvalues(nd,i) = min(npp_data(idx))
172 mx_yvalues(nd,i) = max(npp_data(idx))
173 count = dimsizes(idx)
176 yvalues(nd,i) = yvalues@_FillValue
177 mn_yvalues(nd,i) = yvalues@_FillValue
178 mx_yvalues(nd,i) = yvalues@_FillValue
181 ; Print out information.
183 print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
184 print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
187 ; Clean up for next time in loop.
196 ; Start the graphics.
198 wks = gsn_open_wks("ps","xy")
201 res@gsnMaximize = True
204 res@xyMarkLineMode = "Markers"
205 res@xyMarkerSizeF = 0.014
207 res@xyMarkerColors = (/"Brown","Blue"/)
208 ; res@trYMinF = min(mn_yvalues) - 10.
209 ; res@trYMaxF = max(mx_yvalues) + 10.
210 res@trYMinF = min(mn_yvalues) - 20.
211 res@trYMaxF = max(mx_yvalues) + 100.
213 ; res@tiMainString = "Observed vs i01.03cn"
214 res@tiMainString = "Observed vs i01.04casa"
216 res@tiYAxisString = "Days of Growing season"
217 res@tiXAxisString = "Land Cover Type"
219 ; Add a boxed legend using the more simple method, which won't have
220 ; vertical lines going through the markers.
222 res@pmLegendDisplayMode = "Always"
223 ; res@pmLegendWidthF = 0.1
224 res@pmLegendWidthF = 0.08
225 res@pmLegendHeightF = 0.05
226 res@pmLegendOrthogonalPosF = -1.17
227 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
228 ; res@pmLegendParallelPosF = 0.18
229 res@pmLegendParallelPosF = 0.23 ;(rightward)
231 ; res@lgPerimOn = False
232 res@lgLabelFontHeightF = 0.015
233 ; res@xyExplicitLegendLabels = (/"observed","model_i01.03cn"/)
234 res@xyExplicitLegendLabels = (/"observed","model_i01.04casa"/)
236 xy = gsn_csm_xy(wks,xvalues,yvalues,res)
238 max_bar = new((/2,nx/),graphic)
239 min_bar = new((/2,nx/),graphic)
240 max_cap = new((/2,nx/),graphic)
241 min_cap = new((/2,nx/),graphic)
245 line_colors = (/"brown","blue"/)
247 lnres@gsLineColor = line_colors(nd)
250 if(.not.ismissing(mn_yvalues(nd,i)).and. \
251 .not.ismissing(mx_yvalues(nd,i))) then
253 ; Attach the vertical bar, both above and below the marker.
257 y2 = mn_yvalues(nd,i)
258 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
260 y2 = mx_yvalues(nd,i)
261 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
263 ; Attach the horizontal cap line, both above and below the marker.
265 x1 = xvalues(nd,i) - dx4
266 x2 = xvalues(nd,i) + dx4
267 y1 = mn_yvalues(nd,i)
268 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
270 y1 = mx_yvalues(nd,i)
271 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
277 ; Here's how to do the legend by hand.
279 ; mkres = True ; Marker resources
280 ; txres = True ; Text resources
281 ; mkres@gsMarkerIndex = 16
282 ; mkres@gsMarkerSizeF = 0.02
283 ; txres@txFontHeightF = 0.02
284 ; txres@txJust = "CenterLeft"
286 ; Change these values if you want to move the marker legend location.
287 ; These values are in the same data space as the plot.
295 ; mkres@gsMarkerColor = "brown"
296 ; lnres@gsLineColor = "brown"
298 ; lg_mark_legend1 = gsn_add_polymarker(wks,xy,xlg1_cen,ylg1_cen,mkres)
299 ; lg_line_legend1 = gsn_add_polyline(wks,xy,(/xlg1_cen,xlg1_cen/), \
300 ; (/ylg1_cen-60,ylg1_cen+60/),lnres)
301 ; lg_cap_legend11 = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \
302 ; (/ylg1_cen-60,ylg1_cen-60/),lnres)
303 ; lg_cap_legend12 = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \
304 ; (/ylg1_cen+60,ylg1_cen+60/),lnres)
306 ; tx_legend1 = gsn_add_text(wks,xy,"observed",xlg1_cen+0.15,ylg1_cen,txres)
308 ; mkres@gsMarkerColor = "blue"
309 ; lnres@gsLineColor = "blue"
311 ; lg_mark_legend2 = gsn_add_polymarker(wks,xy,xlg2_cen,ylg2_cen,mkres)
312 ; lg_line_legend2 = gsn_add_polyline(wks,xy,(/xlg2_cen,xlg2_cen/), \
313 ; (/ylg2_cen-60,ylg2_cen+60/),lnres)
314 ; lg_cap_legend21 = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \
315 ; (/ylg2_cen-60,ylg2_cen-60/),lnres)
316 ; lg_cap_legend22 = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \
317 ; (/ylg2_cen+60,ylg2_cen+60/),lnres)
318 ; tx_legend2 = gsn_add_text(wks,xy,"model_i01.03cn",xlg2_cen+0.15,ylg2_cen,txres)
322 system("convert xy.ps xy.png")
327 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
336 ;bias = sum(((vv-uu)/uu)^2)
337 ;M = (1.- sqrt(bias/nz))*5.
340 bias = sum(abs(vv-uu)/(vv+uu))
341 M = (1.- (bias/nz))*5.