1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/npp/21.histogram.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,209 @@
1.4 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.5 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.6 +
1.7 +procedure pminmax(data:numeric,name:string)
1.8 +begin
1.9 + print ("min/max " + name + " = " + min(data) + "/" + max(data))
1.10 + if(isatt(data,"units")) then
1.11 + print (name + " units = " + data@units)
1.12 + end if
1.13 +end
1.14 +
1.15 +;
1.16 +; Main code.
1.17 +;
1.18 +begin
1.19 + data_types = (/ "Obs", "Model" /)
1.20 + data_names = (/ "data.81.nc", "i01.03cn_1545-1569_ANN_climo.nc" /)
1.21 +; data_names = (/ "data.81.nc", "i01.04casa_1605-1629_ANN_climo.nc" /)
1.22 + filevar_names = (/ (/"PREC_ANN","TNPP_C"/), (/"RAIN","NPP"/) /)
1.23 + ndata_types = dimsizes(data_types)
1.24 +
1.25 + data_file_obs = addfile(data_names(0),"r") ; Open obs file
1.26 + data_file_mod = addfile(data_names(1),"r") ; Open model file
1.27 +
1.28 +;
1.29 +; Read four variables from files.
1.30 +;
1.31 + PREC_ANN = tofloat(data_file_obs->PREC_ANN)
1.32 + TNPP_C = data_file_obs->TNPP_C
1.33 + RAIN = data_file_mod->RAIN
1.34 + NPP = data_file_mod->NPP
1.35 +;
1.36 +; Units for these four variables are:
1.37 +;
1.38 +; PREC_ANN : mm/year
1.39 +; TNPP_C : g C/m^2/year
1.40 +; RAIN : mm/s
1.41 +; NPP : g C/m^2/s
1.42 +;
1.43 +; We want to convert these to "m/year" and "g C/m^2/year".
1.44 +;
1.45 + nsec_per_year = 60*60*24*365 ; # seconds per year
1.46 +
1.47 +; Do the necessary conversions.
1.48 + PREC_ANN = PREC_ANN / 1000.
1.49 + RAIN = (RAIN / 1000.) * nsec_per_year
1.50 + NPP = NPP * nsec_per_year
1.51 +
1.52 +; Redo the units.
1.53 + PREC_ANN@units = "m/yr"
1.54 + RAIN@units = "m/yr"
1.55 + NPP@units = "gC/m^2/yr"
1.56 + TNPP_C@units = "gC/m^2/yr"
1.57 +
1.58 + pminmax(PREC_ANN,"PREC_ANN")
1.59 + pminmax(TNPP_C,"TNPP_C")
1.60 + pminmax(RAIN,"RAIN")
1.61 + pminmax(NPP,"NPP")
1.62 +
1.63 + RAIN_1D = ndtooned(RAIN)
1.64 + NPP_1D = ndtooned(NPP)
1.65 + PREC_ANN_1D = ndtooned(PREC_ANN)
1.66 + TNPP_C_1D = ndtooned(TNPP_C)
1.67 +
1.68 +;
1.69 +; Calculate some "nice" bins for binning the data in equally spaced
1.70 +; ranges.
1.71 +;
1.72 + nbins = 15 ; Number of bins to use.
1.73 +
1.74 + nicevals = nice_mnmxintvl(min(RAIN_1D),max(RAIN_1D),nbins,True)
1.75 + nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
1.76 + range = fspan(nicevals(0),nicevals(1),nvals)
1.77 +;
1.78 +; Use this range information to grab all the values in a
1.79 +; particular range, and then take an average.
1.80 +;
1.81 + nr = dimsizes(range)
1.82 + nx = nr-1
1.83 + xvalues = new((/2,nx/),typeof(RAIN_1D))
1.84 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.85 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.86 + dx4 = dx/4 ; 1/4 of the range
1.87 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.88 + yvalues = new((/2,nx/),typeof(RAIN_1D))
1.89 + mn_yvalues = new((/2,nx/),typeof(RAIN_1D))
1.90 + mx_yvalues = new((/2,nx/),typeof(RAIN_1D))
1.91 +
1.92 + do nd=0,1
1.93 +;
1.94 +; See if we are doing model or observational data.
1.95 +;
1.96 + if(nd.eq.0) then
1.97 + data = PREC_ANN_1D
1.98 + npp_data = TNPP_C_1D
1.99 + else
1.100 + data = RAIN_1D
1.101 + npp_data = NPP_1D
1.102 + end if
1.103 +;
1.104 +; Loop through each range and check for values.
1.105 +;
1.106 + do i=0,nr-2
1.107 + if (i.ne.(nr-2)) then
1.108 + print("")
1.109 + print("In range ["+range(i)+","+range(i+1)+")")
1.110 + idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
1.111 + else
1.112 + print("")
1.113 + print("In range ["+range(i)+",)")
1.114 + idx = ind(range(i).le.data)
1.115 + end if
1.116 +;
1.117 +; Calculate average, and get min and max.
1.118 +;
1.119 + if(.not.any(ismissing(idx))) then
1.120 + yvalues(nd,i) = avg(npp_data(idx))
1.121 + mn_yvalues(nd,i) = min(npp_data(idx))
1.122 + mx_yvalues(nd,i) = max(npp_data(idx))
1.123 + count = dimsizes(idx)
1.124 + else
1.125 + count = 0
1.126 + yvalues(nd,i) = yvalues@_FillValue
1.127 + mn_yvalues(nd,i) = yvalues@_FillValue
1.128 + mx_yvalues(nd,i) = yvalues@_FillValue
1.129 + end if
1.130 +;
1.131 +; Print out information.
1.132 +;
1.133 + print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
1.134 + print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.135 +
1.136 +;
1.137 +; Clean up for next time in loop.
1.138 +;
1.139 + delete(idx)
1.140 + end do
1.141 + delete(data)
1.142 + delete(npp_data)
1.143 + end do
1.144 +
1.145 + xvalues@long_name = "Mean Annual precipitation (m/year)"
1.146 + yvalues@long_name = "NPP (g C/m2/year)"
1.147 +
1.148 +;
1.149 +; Start the graphics.
1.150 +;
1.151 +; wks = gsn_open_wks("x11","npp")
1.152 + wks = gsn_open_wks("png","npp")
1.153 +
1.154 + res = True
1.155 + res@tiMainString = "Observed vs i01.03cn"
1.156 +; res@tiMainString = "Observed vs i01.04casa"
1.157 + res@gsnMaximize = False
1.158 + res@gsnDraw = False
1.159 + res@gsnFrame = False
1.160 + res@xyMarkLineMode = "Markers"
1.161 + res@xyMarkerSizeF = 0.014
1.162 + res@xyMarker = 16
1.163 +; res@xyMarkerColors = (/"Gray25","Gray50"/)
1.164 + res@xyMarkerColors = (/"brown","blue"/)
1.165 + res@trYMinF = min(mn_yvalues) - 10.
1.166 + res@trYMaxF = max(mx_yvalues) + 10.
1.167 +
1.168 + xy = gsn_csm_xy(wks,xvalues,yvalues,res)
1.169 +
1.170 + max_bar = new((/2,nx/),graphic)
1.171 + min_bar = new((/2,nx/),graphic)
1.172 + max_cap = new((/2,nx/),graphic)
1.173 + min_cap = new((/2,nx/),graphic)
1.174 +
1.175 + lnres = True
1.176 +
1.177 + line_colors = (/"brown","blue"/)
1.178 + do nd=0,1
1.179 + lnres@gsLineColor = line_colors(nd)
1.180 + do i=0,nx-1
1.181 +
1.182 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.183 + .not.ismissing(mx_yvalues(nd,i))) then
1.184 +;
1.185 +; Attach the vertical bar, both above and below the marker.
1.186 +;
1.187 + x1 = xvalues(nd,i)
1.188 + y1 = yvalues(nd,i)
1.189 + y2 = mn_yvalues(nd,i)
1.190 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.191 +
1.192 + y2 = mx_yvalues(nd,i)
1.193 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.194 +;
1.195 +; Attach the horizontal cap line, both above and below the marker.
1.196 +;
1.197 + x1 = xvalues(nd,i) - dx4
1.198 + x2 = xvalues(nd,i) + dx4
1.199 + y1 = mn_yvalues(nd,i)
1.200 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.201 +
1.202 + y1 = mx_yvalues(nd,i)
1.203 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.204 + end if
1.205 + end do
1.206 + end do
1.207 +
1.208 + draw(xy)
1.209 + frame(wks)
1.210 +
1.211 +end
1.212 +