Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
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)
20 data_types = (/ "Obs", "Model" /)
21 data_names = (/ "data.81.nc" , "i01.03cn_1545-1569_ANN_climo.nc" /)
22 ; data_names = (/ "data.933.nc", "i01.03cn_1545-1569_ANN_climo.nc" /)
23 ; data_names = (/ "data.81.nc" , "i01.04casa_1605-1629_ANN_climo.nc" /)
24 ; data_names = (/ "data.933.nc", "i01.04casa_1605-1629_ANN_climo.nc" /)
25 filevar_names = (/ (/"PREC_ANN","TNPP_C"/), (/"RAIN","NPP"/) /)
26 ndata_types = dimsizes(data_types)
28 data_file_obs = addfile(data_names(0),"r") ; Open obs file
29 data_file_mod = addfile(data_names(1),"r") ; Open model file
31 ;************************************************
32 ; read in data: observed
33 ;************************************************
34 RAIN1 = tofloat(data_file_obs->PREC_ANN) ; for data.81
35 ;RAIN1 = data_file_obs->PREC ; for data.933
36 NPP1 = data_file_obs->TNPP_C
37 xo = data_file_obs->LONG_DD
38 yo = data_file_obs->LAT_DD
40 ; change longitude from (-180,180) to (0,360)
44 if (xo(i) .lt. 0.) then
50 ;************************************************
52 ;************************************************
53 ai = data_file_mod->RAIN
54 bi = data_file_mod->NPP
55 xi = data_file_mod->lon
56 yi = data_file_mod->lat
58 ;************************************************
59 ; interpolate from model grid to observed grid
60 ;************************************************
61 RAIN2 = linint2_points(xi,yi,ai,True,xo,yo,0)
62 NPP2 = linint2_points(xi,yi,bi,True,xo,yo,0)
64 ;************************************************
66 ;************************************************
67 ; Units for these four variables are:
74 ; We want to convert these to "m/year" and "g C/m^2/year".
76 nsec_per_year = 60*60*24*365 ; # seconds per year
78 ; Do the necessary conversions.
80 RAIN2 = (RAIN2/ 1000.) * nsec_per_year
81 NPP2 = NPP2 * nsec_per_year
86 NPP1@units = "gC/m^2/yr"
87 NPP2@units = "gC/m^2/yr"
89 ;************************************************
90 ; print min/max and unit
91 ;************************************************
92 pminmax(RAIN1,"RAIN1")
93 pminmax(RAIN2,"RAIN2")
97 RAIN1_1D = ndtooned(RAIN1)
98 RAIN2_1D = ndtooned(RAIN2)
99 NPP1_1D = ndtooned(NPP1)
100 NPP2_1D = ndtooned(NPP2)
102 ; Calculate some "nice" bins for binning the data in equally spaced
105 nbins = 15 ; Number of bins to use.
107 nicevals = nice_mnmxintvl(min(RAIN2_1D),max(RAIN2_1D),nbins,True)
108 nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
109 range = fspan(nicevals(0),nicevals(1),nvals)
111 ; Use this range information to grab all the values in a
112 ; particular range, and then take an average.
116 xvalues = new((/2,nx/),typeof(RAIN2_1D))
117 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
118 dx = xvalues(0,1) - xvalues(0,0) ; range width
119 dx4 = dx/4 ; 1/4 of the range
120 xvalues(1,:) = xvalues(0,:) - dx/5.
121 yvalues = new((/2,nx/),typeof(RAIN2_1D))
122 mn_yvalues = new((/2,nx/),typeof(RAIN2_1D))
123 mx_yvalues = new((/2,nx/),typeof(RAIN2_1D))
127 ; See if we are doing model or observational data.
137 ; Loop through each range and check for values.
140 if (i.ne.(nr-2)) then
142 ; print("In range ["+range(i)+","+range(i+1)+")")
143 idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
146 ; print("In range ["+range(i)+",)")
147 idx = ind(range(i).le.data)
150 ; Calculate average, and get min and max.
152 if(.not.any(ismissing(idx))) then
153 yvalues(nd,i) = avg(npp_data(idx))
154 mn_yvalues(nd,i) = min(npp_data(idx))
155 mx_yvalues(nd,i) = max(npp_data(idx))
156 count = dimsizes(idx)
159 yvalues(nd,i) = yvalues@_FillValue
160 mn_yvalues(nd,i) = yvalues@_FillValue
161 mx_yvalues(nd,i) = yvalues@_FillValue
164 ; Print out information.
166 ; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
167 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
170 ; Clean up for next time in loop.
179 ; Start the graphics.
181 wks = gsn_open_wks("png","xy")
184 res@gsnMaximize = True
187 res@xyMarkLineMode = "Markers"
188 res@xyMarkerSizeF = 0.014
190 res@xyMarkerColors = (/"Brown","Blue"/)
191 res@trYMinF = min(mn_yvalues) - 10.
192 res@trYMaxF = max(mx_yvalues) + 10.
194 ; res@tiMainString = "Observed vs i01.03cn 81 site"
195 ; res@tiMainString = "Observed vs i01.03cn 933 site"
196 res@tiMainString = "Observed vs i01.04casa 81 site"
197 ; res@tiMainString = "Observed vs i01.04casa 933 site"
198 res@tiYAxisString = "NPP (g C/m2/year)"
199 res@tiXAxisString = "Precipitation (m/year)"
202 ; Add a boxed legend using the more simple method, which won't have
203 ; vertical lines going through the markers.
205 res@pmLegendDisplayMode = "Always"
206 ; res@pmLegendWidthF = 0.1
207 res@pmLegendWidthF = 0.08
208 res@pmLegendHeightF = 0.05
209 res@pmLegendOrthogonalPosF = -1.17
210 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
211 ; res@pmLegendParallelPosF = 0.18
212 res@pmLegendParallelPosF = 0.23 ;(rightward)
214 ; res@lgPerimOn = False
215 res@lgLabelFontHeightF = 0.015
216 ; res@xyExplicitLegendLabels = (/"observed","model_i01.03cn"/)
217 res@xyExplicitLegendLabels = (/"observed","model_i01.04casa"/)
219 xy = gsn_csm_xy(wks,xvalues,yvalues,res)
221 max_bar = new((/2,nx/),graphic)
222 min_bar = new((/2,nx/),graphic)
223 max_cap = new((/2,nx/),graphic)
224 min_cap = new((/2,nx/),graphic)
228 line_colors = (/"brown","blue"/)
230 lnres@gsLineColor = line_colors(nd)
233 if(.not.ismissing(mn_yvalues(nd,i)).and. \
234 .not.ismissing(mx_yvalues(nd,i))) then
236 ; Attach the vertical bar, both above and below the marker.
240 y2 = mn_yvalues(nd,i)
241 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
243 y2 = mx_yvalues(nd,i)
244 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
246 ; Attach the horizontal cap line, both above and below the marker.
248 x1 = xvalues(nd,i) - dx4
249 x2 = xvalues(nd,i) + dx4
250 y1 = mn_yvalues(nd,i)
251 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
253 y1 = mx_yvalues(nd,i)
254 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
260 ; Here's how to do the legend by hand.
262 ; mkres = True ; Marker resources
263 ; txres = True ; Text resources
264 ; mkres@gsMarkerIndex = 16
265 ; mkres@gsMarkerSizeF = 0.02
266 ; txres@txFontHeightF = 0.02
267 ; txres@txJust = "CenterLeft"
269 ; Change these values if you want to move the marker legend location.
270 ; These values are in the same data space as the plot.
278 ; mkres@gsMarkerColor = "brown"
279 ; lnres@gsLineColor = "brown"
281 ; lg_mark_legend1 = gsn_add_polymarker(wks,xy,xlg1_cen,ylg1_cen,mkres)
282 ; lg_line_legend1 = gsn_add_polyline(wks,xy,(/xlg1_cen,xlg1_cen/), \
283 ; (/ylg1_cen-60,ylg1_cen+60/),lnres)
284 ; lg_cap_legend11 = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \
285 ; (/ylg1_cen-60,ylg1_cen-60/),lnres)
286 ; lg_cap_legend12 = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \
287 ; (/ylg1_cen+60,ylg1_cen+60/),lnres)
289 ; tx_legend1 = gsn_add_text(wks,xy,"observed",xlg1_cen+0.15,ylg1_cen,txres)
291 ; mkres@gsMarkerColor = "blue"
292 ; lnres@gsLineColor = "blue"
294 ; lg_mark_legend2 = gsn_add_polymarker(wks,xy,xlg2_cen,ylg2_cen,mkres)
295 ; lg_line_legend2 = gsn_add_polyline(wks,xy,(/xlg2_cen,xlg2_cen/), \
296 ; (/ylg2_cen-60,ylg2_cen+60/),lnres)
297 ; lg_cap_legend21 = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \
298 ; (/ylg2_cen-60,ylg2_cen-60/),lnres)
299 ; lg_cap_legend22 = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \
300 ; (/ylg2_cen+60,ylg2_cen+60/),lnres)
301 ; tx_legend2 = gsn_add_text(wks,xy,"model_i01.03cn",xlg2_cen+0.15,ylg2_cen,txres)
309 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
319 ;bias = sum(((vv-uu)/uu)^2)
320 ;M = (1.- sqrt(bias/nz))*5.
323 bias = sum(abs(vv-uu)/(vv+uu))
324 M = (1.- (bias/nz))*5.