Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
5 procedure pminmax(data:numeric,name:string)
7 print ("min/max " + name + " = " + min(data) + "/" + max(data))
8 if(isatt(data,"units")) then
9 print (name + " units = " + data@units)
17 data_types = (/ "Obs", "Model" /)
18 data_names = (/ "data.81.nc", "i01.03cn_1545-1569_ANN_climo.nc" /)
19 ; data_names = (/ "data.81.nc", "i01.04casa_1605-1629_ANN_climo.nc" /)
20 ; filevar_names = (/ (/"PREC_ANN","TNPP_C"/), (/"RAIN","NPP"/) /)
21 filevar_names = (/ (/"PREC_ANN","ANPP_C"/), (/"RAIN","AGNPP"/) /)
22 ndata_types = dimsizes(data_types)
24 data_file_obs = addfile(data_names(0),"r") ; Open obs file
25 data_file_mod = addfile(data_names(1),"r") ; Open model file
27 ;************************************************
28 ; read in data: observed
29 ;************************************************
30 PREC_ANN = tofloat(data_file_obs->PREC_ANN)
31 TNPP_C = data_file_obs->TNPP_C
32 xo = data_file_obs->LONG_DD
33 yo = data_file_obs->LAT_DD
35 ; change longitude from (-180,180) to (0,360)
39 if (xo(i) .lt. 0.) then
45 ;************************************************
47 ;************************************************
48 ai = data_file_mod->RAIN
49 bi = data_file_mod->NPP
50 xi = data_file_mod->lon
51 yi = data_file_mod->lat
53 ;************************************************
54 ; interpolate from model grid to observed grid
55 ;************************************************
56 RAIN = linint2_points(xi,yi,ai,True,xo,yo,0)
57 NPP = linint2_points(xi,yi,bi,True,xo,yo,0)
59 ;************************************************
61 ;************************************************
62 ; Units for these four variables are:
65 ; TNPP_C : g C/m^2/year
69 ; We want to convert these to "m/year" and "g C/m^2/year".
71 nsec_per_year = 60*60*24*365 ; # seconds per year
73 ; Do the necessary conversions.
74 PREC_ANN = PREC_ANN / 1000.
75 RAIN = (RAIN / 1000.) * nsec_per_year
76 NPP = NPP * nsec_per_year
79 PREC_ANN@units = "m/yr"
81 NPP@units = "gC/m^2/yr"
82 TNPP_C@units = "gC/m^2/yr"
84 ;************************************************
85 pminmax(PREC_ANN,"PREC_ANN")
86 pminmax(TNPP_C,"TNPP_C")
90 RAIN_1D = ndtooned(RAIN)
91 NPP_1D = ndtooned(NPP)
92 PREC_ANN_1D = ndtooned(PREC_ANN)
93 TNPP_C_1D = ndtooned(TNPP_C)
96 ; Calculate some "nice" bins for binning the data in equally spaced
99 nbins = 15 ; Number of bins to use.
101 nicevals = nice_mnmxintvl(min(RAIN_1D),max(RAIN_1D),nbins,True)
102 nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
103 range = fspan(nicevals(0),nicevals(1),nvals)
105 ; Use this range information to grab all the values in a
106 ; particular range, and then take an average.
110 xvalues = new((/2,nx/),typeof(RAIN_1D))
111 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
112 dx = xvalues(0,1) - xvalues(0,0) ; range width
113 dx4 = dx/4 ; 1/4 of the range
114 xvalues(1,:) = xvalues(0,:) - dx/5.
115 yvalues = new((/2,nx/),typeof(RAIN_1D))
116 mn_yvalues = new((/2,nx/),typeof(RAIN_1D))
117 mx_yvalues = new((/2,nx/),typeof(RAIN_1D))
121 ; See if we are doing model or observational data.
131 ; Loop through each range and check for values.
134 if (i.ne.(nr-2)) then
136 print("In range ["+range(i)+","+range(i+1)+")")
137 idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
140 print("In range ["+range(i)+",)")
141 idx = ind(range(i).le.data)
144 ; Calculate average, and get min and max.
146 if(.not.any(ismissing(idx))) then
147 yvalues(nd,i) = avg(npp_data(idx))
148 mn_yvalues(nd,i) = min(npp_data(idx))
149 mx_yvalues(nd,i) = max(npp_data(idx))
150 count = dimsizes(idx)
153 yvalues(nd,i) = yvalues@_FillValue
154 mn_yvalues(nd,i) = yvalues@_FillValue
155 mx_yvalues(nd,i) = yvalues@_FillValue
158 ; Print out information.
160 print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
161 print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
164 ; Clean up for next time in loop.
172 xvalues@long_name = "Mean Annual precipitation (m/year)"
173 yvalues@long_name = "NPP (g C/m2/year)"
176 ; Start the graphics.
178 wks = gsn_open_wks("png","npp")
181 res@gsnMaximize = False
184 res@xyMarkLineMode = "Markers"
185 res@xyMarkerSizeF = 0.014
187 ; res@xyMarkerColors = (/"Gray25","Gray50"/)
188 res@xyMarkerColors = (/"brown","blue"/)
189 res@trYMinF = min(mn_yvalues) - 10.
190 res@trYMaxF = max(mx_yvalues) + 10.
191 ; res@tiMainString = "Observed vs i01.03cn_81site"
192 res@tiMainString = "Observed vs i01.04casa_81site"
194 xy = gsn_csm_xy(wks,xvalues,yvalues,res)
196 max_bar = new((/2,nx/),graphic)
197 min_bar = new((/2,nx/),graphic)
198 max_cap = new((/2,nx/),graphic)
199 min_cap = new((/2,nx/),graphic)
203 line_colors = (/"brown","blue"/)
205 lnres@gsLineColor = line_colors(nd)
208 if(.not.ismissing(mn_yvalues(nd,i)).and. \
209 .not.ismissing(mx_yvalues(nd,i))) then
211 ; Attach the vertical bar, both above and below the marker.
215 y2 = mn_yvalues(nd,i)
216 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
218 y2 = mx_yvalues(nd,i)
219 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
221 ; Attach the horizontal cap line, both above and below the marker.
223 x1 = xvalues(nd,i) - dx4
224 x2 = xvalues(nd,i) + dx4
225 y1 = mn_yvalues(nd,i)
226 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
228 y1 = mx_yvalues(nd,i)
229 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)