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