Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
1 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
4 procedure pminmax(data:numeric,name:string)
6 print ("min/max " + name + " = " + min(data) + "/" + max(data))
7 if(isatt(data,"units")) then
8 print (name + " units = " + data@units)
16 data_types = (/ "Obs", "Model" /)
17 data_names = (/ "data.81.nc", "i01.03cn_1545-1569_ANN_climo.nc" /)
18 ; data_names = (/ "data.81.nc", "i01.04casa_1605-1629_ANN_climo.nc" /)
19 filevar_names = (/ (/"PREC_ANN","TNPP_C"/), (/"RAIN","NPP"/) /)
20 ndata_types = dimsizes(data_types)
22 data_file_obs = addfile(data_names(0),"r") ; Open obs file
23 data_file_mod = addfile(data_names(1),"r") ; Open model file
26 ; Read four variables from files.
28 PREC_ANN = tofloat(data_file_obs->PREC_ANN)
29 TNPP_C = data_file_obs->TNPP_C
30 RAIN = data_file_mod->RAIN
31 NPP = data_file_mod->NPP
33 ; Units for these four variables are:
36 ; TNPP_C : g C/m^2/year
40 ; We want to convert these to "m/year" and "g C/m^2/year".
42 nsec_per_year = 60*60*24*365 ; # seconds per year
44 ; Do the necessary conversions.
45 PREC_ANN = PREC_ANN / 1000.
46 RAIN = (RAIN / 1000.) * nsec_per_year
47 NPP = NPP * nsec_per_year
50 PREC_ANN@units = "m/yr"
52 NPP@units = "gC/m^2/yr"
53 TNPP_C@units = "gC/m^2/yr"
55 pminmax(PREC_ANN,"PREC_ANN")
56 pminmax(TNPP_C,"TNPP_C")
60 RAIN_1D = ndtooned(RAIN)
61 NPP_1D = ndtooned(NPP)
62 PREC_ANN_1D = ndtooned(PREC_ANN)
63 TNPP_C_1D = ndtooned(TNPP_C)
66 ; Calculate some "nice" bins for binning the data in equally spaced
69 nbins = 15 ; Number of bins to use.
71 nicevals = nice_mnmxintvl(min(RAIN_1D),max(RAIN_1D),nbins,True)
72 nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
73 range = fspan(nicevals(0),nicevals(1),nvals)
75 ; Use this range information to grab all the values in a
76 ; particular range, and then take an average.
80 xvalues = new((/2,nx/),typeof(RAIN_1D))
81 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
82 dx = xvalues(0,1) - xvalues(0,0) ; range width
83 dx4 = dx/4 ; 1/4 of the range
84 xvalues(1,:) = xvalues(0,:) - dx/5.
85 yvalues = new((/2,nx/),typeof(RAIN_1D))
86 mn_yvalues = new((/2,nx/),typeof(RAIN_1D))
87 mx_yvalues = new((/2,nx/),typeof(RAIN_1D))
91 ; See if we are doing model or observational data.
101 ; Loop through each range and check for values.
104 if (i.ne.(nr-2)) then
106 print("In range ["+range(i)+","+range(i+1)+")")
107 idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
110 print("In range ["+range(i)+",)")
111 idx = ind(range(i).le.data)
114 ; Calculate average, and get min and max.
116 if(.not.any(ismissing(idx))) then
117 yvalues(nd,i) = avg(npp_data(idx))
118 mn_yvalues(nd,i) = min(npp_data(idx))
119 mx_yvalues(nd,i) = max(npp_data(idx))
120 count = dimsizes(idx)
123 yvalues(nd,i) = yvalues@_FillValue
124 mn_yvalues(nd,i) = yvalues@_FillValue
125 mx_yvalues(nd,i) = yvalues@_FillValue
128 ; Print out information.
130 print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
131 print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
134 ; Clean up for next time in loop.
142 xvalues@long_name = "Mean Annual precipitation (m/year)"
143 yvalues@long_name = "NPP (g C/m2/year)"
146 ; Start the graphics.
148 ; wks = gsn_open_wks("x11","npp")
149 wks = gsn_open_wks("png","npp")
152 res@tiMainString = "Observed vs i01.03cn"
153 ; res@tiMainString = "Observed vs i01.04casa"
154 res@gsnMaximize = False
157 res@xyMarkLineMode = "Markers"
158 res@xyMarkerSizeF = 0.014
160 ; res@xyMarkerColors = (/"Gray25","Gray50"/)
161 res@xyMarkerColors = (/"brown","blue"/)
162 res@trYMinF = min(mn_yvalues) - 10.
163 res@trYMaxF = max(mx_yvalues) + 10.
165 xy = gsn_csm_xy(wks,xvalues,yvalues,res)
167 max_bar = new((/2,nx/),graphic)
168 min_bar = new((/2,nx/),graphic)
169 max_cap = new((/2,nx/),graphic)
170 min_cap = new((/2,nx/),graphic)
174 line_colors = (/"brown","blue"/)
176 lnres@gsLineColor = line_colors(nd)
179 if(.not.ismissing(mn_yvalues(nd,i)).and. \
180 .not.ismissing(mx_yvalues(nd,i))) then
182 ; Attach the vertical bar, both above and below the marker.
186 y2 = mn_yvalues(nd,i)
187 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
189 y2 = mx_yvalues(nd,i)
190 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
192 ; Attach the horizontal cap line, both above and below the marker.
194 x1 = xvalues(nd,i) - dx4
195 x2 = xvalues(nd,i) + dx4
196 y1 = mn_yvalues(nd,i)
197 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
199 y1 = mx_yvalues(nd,i)
200 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)