|
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" |
|
7 |
|
8 procedure pminmax(data:numeric,name:string) |
|
9 begin |
|
10 print ("min/max " + name + " = " + min(data) + "/" + max(data)) |
|
11 if(isatt(data,"units")) then |
|
12 print (name + " units = " + data@units) |
|
13 end if |
|
14 end |
|
15 |
|
16 ; |
|
17 ; Main code. |
|
18 ; |
|
19 begin |
|
20 |
|
21 ;nclass = 18 |
|
22 nclass = 20 |
|
23 day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/) |
|
24 |
|
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") |
|
34 |
|
35 RAIN1 = tofloat(data_file_ob1->LAND_CLASS) |
|
36 |
|
37 z = data_file_ob2->LAI |
|
38 y = z(0,:,:) |
|
39 y@long_name = "Days of Growing Season" |
|
40 |
|
41 dsizes_z = dimsizes(z) |
|
42 ntime = dsizes_z(0) |
|
43 nlat = dsizes_z(1) |
|
44 nlon = dsizes_z(2) |
|
45 |
|
46 do j = 0,nlat-1 |
|
47 do i = 0,nlon-1 |
|
48 nday = 0. |
|
49 do k = 0,ntime-1 |
|
50 if (.not. ismissing(z(k,j,i)) .and. z(k,j,i) .gt. 1.0) then |
|
51 nday = nday + day_of_data(k) |
|
52 end if |
|
53 end do |
|
54 y(j,i) = nday |
|
55 end do |
|
56 end do |
|
57 |
|
58 print (min(y)+"/"+max(y)) |
|
59 |
|
60 NPP1 = y |
|
61 |
|
62 delete (z) |
|
63 delete (y) |
|
64 ;************************************************ |
|
65 ; read in data: model |
|
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") |
|
71 |
|
72 z = data_file_model->TLAI |
|
73 y = z(0,:,:) |
|
74 y@long_name = "Days of Growing Season" |
|
75 |
|
76 dsizes_z = dimsizes(z) |
|
77 ntime = dsizes_z(0) |
|
78 nlat = dsizes_z(1) |
|
79 nlon = dsizes_z(2) |
|
80 |
|
81 do j = 0,nlat-1 |
|
82 do i = 0,nlon-1 |
|
83 nday = 0. |
|
84 do k = 0,ntime-1 |
|
85 if (.not. ismissing(z(k,j,i)) .and. z(k,j,i) .gt. 1.0) then |
|
86 nday = nday + day_of_data(k) |
|
87 end if |
|
88 end do |
|
89 y(j,i) = nday |
|
90 end do |
|
91 end do |
|
92 |
|
93 print (min(y)+"/"+max(y)) |
|
94 |
|
95 NPP2 = y |
|
96 |
|
97 delete (z) |
|
98 delete (y) |
|
99 ;************************************************ |
|
100 ; print min/max and unit |
|
101 ;************************************************ |
|
102 pminmax(RAIN1,"RAIN1") |
|
103 pminmax(NPP1,"NPP1") |
|
104 pminmax(NPP2,"NPP2") |
|
105 |
|
106 RAIN1_1D = ndtooned(RAIN1) |
|
107 NPP1_1D = ndtooned(NPP1) |
|
108 NPP2_1D = ndtooned(NPP2) |
|
109 ; |
|
110 ; Calculate some "nice" bins for binning the data in equally spaced |
|
111 ; ranges. |
|
112 ; |
|
113 |
|
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) |
|
118 |
|
119 nclass = nclass + 1 |
|
120 range = fspan(0,nclass-1,nclass) |
|
121 |
|
122 ; print (nicevals) |
|
123 ; print (nvals) |
|
124 print (range) |
|
125 ; exit |
|
126 |
|
127 ; |
|
128 ; Use this range information to grab all the values in a |
|
129 ; particular range, and then take an average. |
|
130 ; |
|
131 nr = dimsizes(range) |
|
132 nx = nr-1 |
|
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)) |
|
141 |
|
142 do nd=0,1 |
|
143 ; |
|
144 ; See if we are doing model or observational data. |
|
145 ; |
|
146 if(nd.eq.0) then |
|
147 data = RAIN1_1D |
|
148 npp_data = NPP1_1D |
|
149 else |
|
150 data = RAIN1_1D |
|
151 npp_data = NPP2_1D |
|
152 end if |
|
153 ; |
|
154 ; Loop through each range and check for values. |
|
155 ; |
|
156 do i=0,nr-2 |
|
157 if (i.ne.(nr-2)) then |
|
158 print("") |
|
159 print("In range ["+range(i)+","+range(i+1)+")") |
|
160 idx = ind((range(i).le.data).and.(data.lt.range(i+1))) |
|
161 else |
|
162 print("") |
|
163 print("In range ["+range(i)+",)") |
|
164 idx = ind(range(i).le.data) |
|
165 end if |
|
166 ; |
|
167 ; Calculate average, and get min and max. |
|
168 ; |
|
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) |
|
174 else |
|
175 count = 0 |
|
176 yvalues(nd,i) = yvalues@_FillValue |
|
177 mn_yvalues(nd,i) = yvalues@_FillValue |
|
178 mx_yvalues(nd,i) = yvalues@_FillValue |
|
179 end if |
|
180 ; |
|
181 ; Print out information. |
|
182 ; |
|
183 print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) |
|
184 print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) |
|
185 |
|
186 ; |
|
187 ; Clean up for next time in loop. |
|
188 ; |
|
189 delete(idx) |
|
190 end do |
|
191 delete(data) |
|
192 delete(npp_data) |
|
193 end do |
|
194 |
|
195 ; |
|
196 ; Start the graphics. |
|
197 ; |
|
198 wks = gsn_open_wks("ps","xy") |
|
199 |
|
200 res = True |
|
201 res@gsnMaximize = True |
|
202 res@gsnDraw = False |
|
203 res@gsnFrame = False |
|
204 res@xyMarkLineMode = "Markers" |
|
205 res@xyMarkerSizeF = 0.014 |
|
206 res@xyMarker = 16 |
|
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. |
|
212 |
|
213 ; res@tiMainString = "Observed vs i01.03cn" |
|
214 res@tiMainString = "Observed vs i01.04casa" |
|
215 |
|
216 res@tiYAxisString = "Days of Growing season" |
|
217 res@tiXAxisString = "Land Cover Type" |
|
218 ; |
|
219 ; Add a boxed legend using the more simple method, which won't have |
|
220 ; vertical lines going through the markers. |
|
221 ; |
|
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) |
|
230 |
|
231 ; res@lgPerimOn = False |
|
232 res@lgLabelFontHeightF = 0.015 |
|
233 ; res@xyExplicitLegendLabels = (/"observed","model_i01.03cn"/) |
|
234 res@xyExplicitLegendLabels = (/"observed","model_i01.04casa"/) |
|
235 |
|
236 xy = gsn_csm_xy(wks,xvalues,yvalues,res) |
|
237 |
|
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) |
|
242 |
|
243 lnres = True |
|
244 |
|
245 line_colors = (/"brown","blue"/) |
|
246 do nd=0,1 |
|
247 lnres@gsLineColor = line_colors(nd) |
|
248 do i=0,nx-1 |
|
249 |
|
250 if(.not.ismissing(mn_yvalues(nd,i)).and. \ |
|
251 .not.ismissing(mx_yvalues(nd,i))) then |
|
252 ; |
|
253 ; Attach the vertical bar, both above and below the marker. |
|
254 ; |
|
255 x1 = xvalues(nd,i) |
|
256 y1 = yvalues(nd,i) |
|
257 y2 = mn_yvalues(nd,i) |
|
258 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
259 |
|
260 y2 = mx_yvalues(nd,i) |
|
261 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) |
|
262 ; |
|
263 ; Attach the horizontal cap line, both above and below the marker. |
|
264 ; |
|
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) |
|
269 |
|
270 y1 = mx_yvalues(nd,i) |
|
271 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) |
|
272 end if |
|
273 end do |
|
274 end do |
|
275 |
|
276 ; |
|
277 ; Here's how to do the legend by hand. |
|
278 ; |
|
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" |
|
285 ; |
|
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. |
|
288 ; |
|
289 ; xlg1_cen = 0.2 |
|
290 ; ylg1_cen = 900. |
|
291 |
|
292 ; xlg2_cen = 0.2 |
|
293 ; ylg2_cen = 760. |
|
294 |
|
295 ; mkres@gsMarkerColor = "brown" |
|
296 ; lnres@gsLineColor = "brown" |
|
297 |
|
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) |
|
305 |
|
306 ; tx_legend1 = gsn_add_text(wks,xy,"observed",xlg1_cen+0.15,ylg1_cen,txres) |
|
307 |
|
308 ; mkres@gsMarkerColor = "blue" |
|
309 ; lnres@gsLineColor = "blue" |
|
310 |
|
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) |
|
319 |
|
320 draw(xy) |
|
321 frame(wks) |
|
322 system("convert xy.ps xy.png") |
|
323 |
|
324 u = yvalues(0,:) |
|
325 v = yvalues(1,:) |
|
326 |
|
327 good = ind(.not.ismissing(u) .and. .not.ismissing(v)) |
|
328 uu = u(good) |
|
329 vv = v(good) |
|
330 nz = dimsizes(uu) |
|
331 print (nz) |
|
332 |
|
333 ccr = esccr(uu,vv,0) |
|
334 print (ccr) |
|
335 ;old eq |
|
336 ;bias = sum(((vv-uu)/uu)^2) |
|
337 ;M = (1.- sqrt(bias/nz))*5. |
|
338 |
|
339 ;new eq |
|
340 bias = sum(abs(vv-uu)/(vv+uu)) |
|
341 M = (1.- (bias/nz))*5. |
|
342 |
|
343 print (bias) |
|
344 print (M) |
|
345 |
|
346 end |
|
347 |