1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/npp/24.histogram_ob.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,306 @@
1.4 +;********************************************************
1.5 +; histogram normalized by rain and compute correleration
1.6 +;********************************************************
1.7 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.10 +
1.11 +procedure pminmax(data:numeric,name:string)
1.12 +begin
1.13 + print ("min/max " + name + " = " + min(data) + "/" + max(data))
1.14 + if(isatt(data,"units")) then
1.15 + print (name + " units = " + data@units)
1.16 + end if
1.17 +end
1.18 +
1.19 +;
1.20 +; Main code.
1.21 +;
1.22 +begin
1.23 + data_types = (/ "Obs", "Model" /)
1.24 + data_names = (/ "data.81.nc" , "i01.03cn_1545-1569_ANN_climo.nc" /)
1.25 +; data_names = (/ "data.933.nc", "i01.03cn_1545-1569_ANN_climo.nc" /)
1.26 +; data_names = (/ "data.81.nc" , "i01.04casa_1605-1629_ANN_climo.nc" /)
1.27 +; data_names = (/ "data.933.nc", "i01.04casa_1605-1629_ANN_climo.nc" /)
1.28 + filevar_names = (/ (/"PREC_ANN","TNPP_C"/), (/"RAIN","NPP"/) /)
1.29 + ndata_types = dimsizes(data_types)
1.30 +
1.31 + data_file_obs = addfile(data_names(0),"r") ; Open obs file
1.32 + data_file_mod = addfile(data_names(1),"r") ; Open model file
1.33 +
1.34 +;************************************************
1.35 +; read in data: observed
1.36 +;************************************************
1.37 + RAIN1 = tofloat(data_file_obs->PREC_ANN) ; for data.81
1.38 +;RAIN1 = data_file_obs->PREC ; for data.933
1.39 + NPP1 = data_file_obs->TNPP_C
1.40 + xo = data_file_obs->LONG_DD
1.41 + yo = data_file_obs->LAT_DD
1.42 +
1.43 +; change longitude from (-180,180) to (0,360)
1.44 + print (xo)
1.45 + nx = dimsizes(xo)
1.46 + do i= 0,nx-1
1.47 + if (xo(i) .lt. 0.) then
1.48 + xo(i) = xo(i)+ 360.
1.49 + end if
1.50 + end do
1.51 + print (xo)
1.52 +
1.53 +;************************************************
1.54 +; read in data: model
1.55 +;************************************************
1.56 + ai = data_file_mod->RAIN
1.57 + bi = data_file_mod->NPP
1.58 + xi = data_file_mod->lon
1.59 + yi = data_file_mod->lat
1.60 +
1.61 +;************************************************
1.62 +; interpolate from model grid to observed grid
1.63 +;************************************************
1.64 + RAIN2 = linint2_points(xi,yi,ai,True,xo,yo,0)
1.65 + NPP2 = linint2_points(xi,yi,bi,True,xo,yo,0)
1.66 +
1.67 +;************************************************
1.68 +; convert unit
1.69 +;************************************************
1.70 +; Units for these four variables are:
1.71 +;
1.72 +; RAIN1 : mm/year
1.73 +; RAIN2 : mm/s
1.74 +; NPP1 : g C/m^2/year
1.75 +; NPP2 : g C/m^2/s
1.76 +;
1.77 +; We want to convert these to "m/year" and "g C/m^2/year".
1.78 +;
1.79 + nsec_per_year = 60*60*24*365 ; # seconds per year
1.80 +
1.81 +; Do the necessary conversions.
1.82 + RAIN1 = RAIN1 / 1000.
1.83 + RAIN2 = (RAIN2/ 1000.) * nsec_per_year
1.84 + NPP2 = NPP2 * nsec_per_year
1.85 +
1.86 +; Redo the units.
1.87 + RAIN1@units = "m/yr"
1.88 + RAIN2@units = "m/yr"
1.89 + NPP1@units = "gC/m^2/yr"
1.90 + NPP2@units = "gC/m^2/yr"
1.91 +
1.92 +;************************************************
1.93 +; print min/max and unit
1.94 +;************************************************
1.95 + pminmax(RAIN1,"RAIN1")
1.96 + pminmax(RAIN2,"RAIN2")
1.97 + pminmax(NPP1,"NPP1")
1.98 + pminmax(NPP2,"NPP2")
1.99 +
1.100 + RAIN1_1D = ndtooned(RAIN1)
1.101 + RAIN2_1D = ndtooned(RAIN2)
1.102 + NPP1_1D = ndtooned(NPP1)
1.103 + NPP2_1D = ndtooned(NPP2)
1.104 +;
1.105 +; Calculate some "nice" bins for binning the data in equally spaced
1.106 +; ranges.
1.107 +;
1.108 + nbins = 15 ; Number of bins to use.
1.109 +
1.110 + nicevals = nice_mnmxintvl(min(RAIN2_1D),max(RAIN2_1D),nbins,True)
1.111 + nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
1.112 + range = fspan(nicevals(0),nicevals(1),nvals)
1.113 +;
1.114 +; Use this range information to grab all the values in a
1.115 +; particular range, and then take an average.
1.116 +;
1.117 + nr = dimsizes(range)
1.118 + nx = nr-1
1.119 + xvalues = new((/2,nx/),typeof(RAIN2_1D))
1.120 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.121 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.122 + dx4 = dx/4 ; 1/4 of the range
1.123 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.124 + yvalues = new((/2,nx/),typeof(RAIN2_1D))
1.125 + mn_yvalues = new((/2,nx/),typeof(RAIN2_1D))
1.126 + mx_yvalues = new((/2,nx/),typeof(RAIN2_1D))
1.127 +
1.128 + do nd=0,1
1.129 +;
1.130 +; See if we are doing model or observational data.
1.131 +;
1.132 + if(nd.eq.0) then
1.133 + data = RAIN1_1D
1.134 + npp_data = NPP1_1D
1.135 + else
1.136 + data = RAIN2_1D
1.137 + npp_data = NPP2_1D
1.138 + end if
1.139 +;
1.140 +; Loop through each range and check for values.
1.141 +;
1.142 + do i=0,nr-2
1.143 + if (i.ne.(nr-2)) then
1.144 +; print("")
1.145 +; print("In range ["+range(i)+","+range(i+1)+")")
1.146 + idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
1.147 + else
1.148 +; print("")
1.149 +; print("In range ["+range(i)+",)")
1.150 + idx = ind(range(i).le.data)
1.151 + end if
1.152 +;
1.153 +; Calculate average, and get min and max.
1.154 +;
1.155 + if(.not.any(ismissing(idx))) then
1.156 + yvalues(nd,i) = avg(npp_data(idx))
1.157 + mn_yvalues(nd,i) = min(npp_data(idx))
1.158 + mx_yvalues(nd,i) = max(npp_data(idx))
1.159 + count = dimsizes(idx)
1.160 + else
1.161 + count = 0
1.162 + yvalues(nd,i) = yvalues@_FillValue
1.163 + mn_yvalues(nd,i) = yvalues@_FillValue
1.164 + mx_yvalues(nd,i) = yvalues@_FillValue
1.165 + end if
1.166 +;
1.167 +; Print out information.
1.168 +;
1.169 +; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
1.170 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.171 +
1.172 +;
1.173 +; Clean up for next time in loop.
1.174 +;
1.175 + delete(idx)
1.176 + end do
1.177 + delete(data)
1.178 + delete(npp_data)
1.179 + end do
1.180 +
1.181 +;
1.182 +; Start the graphics.
1.183 +;
1.184 + wks = gsn_open_wks("png","hist")
1.185 +
1.186 + res = True
1.187 + res@gsnMaximize = True
1.188 + res@gsnDraw = False
1.189 + res@gsnFrame = False
1.190 + res@xyMarkLineMode = "Markers"
1.191 + res@xyMarkerSizeF = 0.014
1.192 + res@xyMarker = 16
1.193 + res@xyMarkerColors = (/"Brown","Blue"/)
1.194 + res@trYMinF = min(mn_yvalues) - 10.
1.195 + res@trYMaxF = max(mx_yvalues) + 10.
1.196 +
1.197 + res@tiMainString = "Observed 81 sites"
1.198 +; res@tiMainString = "Observed 933 sites"
1.199 + res@tiYAxisString = "NPP (g C/m2/year)"
1.200 + res@tiXAxisString = "Precipitation (m/year)"
1.201 +
1.202 +;
1.203 +; Add a boxed legend using the more simple method, which won't have
1.204 +; vertical lines going through the markers.
1.205 +;
1.206 +; res@pmLegendDisplayMode = "Always"
1.207 +;;res@pmLegendWidthF = 0.1
1.208 +; res@pmLegendWidthF = 0.08
1.209 +; res@pmLegendHeightF = 0.05
1.210 +; res@pmLegendOrthogonalPosF = -1.17
1.211 +;;res@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.212 +;;res@pmLegendParallelPosF = 0.18
1.213 +; res@pmLegendParallelPosF = 0.23 ;(rightward)
1.214 +
1.215 +; res@lgPerimOn = False
1.216 +; res@lgLabelFontHeightF = 0.015
1.217 +; res@xyExplicitLegendLabels = (/"observed","model_i01.03cn"/)
1.218 +;;res@xyExplicitLegendLabels = (/"observed","model_i01.04casa"/)
1.219 +
1.220 + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),res)
1.221 +
1.222 + max_bar = new((/2,nx/),graphic)
1.223 + min_bar = new((/2,nx/),graphic)
1.224 + max_cap = new((/2,nx/),graphic)
1.225 + min_cap = new((/2,nx/),graphic)
1.226 +
1.227 + lnres = True
1.228 +
1.229 + line_colors = (/"brown","blue"/)
1.230 +; do nd=0,1
1.231 + do nd=0,0
1.232 + lnres@gsLineColor = line_colors(nd)
1.233 + do i=0,nx-1
1.234 +
1.235 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.236 + .not.ismissing(mx_yvalues(nd,i))) then
1.237 +;
1.238 +; Attach the vertical bar, both above and below the marker.
1.239 +;
1.240 + x1 = xvalues(nd,i)
1.241 + y1 = yvalues(nd,i)
1.242 + y2 = mn_yvalues(nd,i)
1.243 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.244 +
1.245 + y2 = mx_yvalues(nd,i)
1.246 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.247 +;
1.248 +; Attach the horizontal cap line, both above and below the marker.
1.249 +;
1.250 + x1 = xvalues(nd,i) - dx4
1.251 + x2 = xvalues(nd,i) + dx4
1.252 + y1 = mn_yvalues(nd,i)
1.253 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.254 +
1.255 + y1 = mx_yvalues(nd,i)
1.256 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.257 + end if
1.258 + end do
1.259 + end do
1.260 +
1.261 +;
1.262 +; Here's how to do the legend by hand.
1.263 +;
1.264 +; mkres = True ; Marker resources
1.265 +; txres = True ; Text resources
1.266 +; mkres@gsMarkerIndex = 16
1.267 +; mkres@gsMarkerSizeF = 0.02
1.268 +; txres@txFontHeightF = 0.02
1.269 +; txres@txJust = "CenterLeft"
1.270 +;
1.271 +; Change these values if you want to move the marker legend location.
1.272 +; These values are in the same data space as the plot.
1.273 +;
1.274 +; xlg1_cen = 0.2
1.275 +; ylg1_cen = 900.
1.276 +
1.277 +; xlg2_cen = 0.2
1.278 +; ylg2_cen = 760.
1.279 +
1.280 +; mkres@gsMarkerColor = "brown"
1.281 +; lnres@gsLineColor = "brown"
1.282 +
1.283 +; lg_mark_legend1 = gsn_add_polymarker(wks,xy,xlg1_cen,ylg1_cen,mkres)
1.284 +; lg_line_legend1 = gsn_add_polyline(wks,xy,(/xlg1_cen,xlg1_cen/), \
1.285 +; (/ylg1_cen-60,ylg1_cen+60/),lnres)
1.286 +; lg_cap_legend11 = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \
1.287 +; (/ylg1_cen-60,ylg1_cen-60/),lnres)
1.288 +; lg_cap_legend12 = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \
1.289 +; (/ylg1_cen+60,ylg1_cen+60/),lnres)
1.290 +
1.291 +; tx_legend1 = gsn_add_text(wks,xy,"observed",xlg1_cen+0.15,ylg1_cen,txres)
1.292 +
1.293 +; mkres@gsMarkerColor = "blue"
1.294 +; lnres@gsLineColor = "blue"
1.295 +
1.296 +; lg_mark_legend2 = gsn_add_polymarker(wks,xy,xlg2_cen,ylg2_cen,mkres)
1.297 +; lg_line_legend2 = gsn_add_polyline(wks,xy,(/xlg2_cen,xlg2_cen/), \
1.298 +; (/ylg2_cen-60,ylg2_cen+60/),lnres)
1.299 +; lg_cap_legend21 = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \
1.300 +; (/ylg2_cen-60,ylg2_cen-60/),lnres)
1.301 +; lg_cap_legend22 = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \
1.302 +; (/ylg2_cen+60,ylg2_cen+60/),lnres)
1.303 +; tx_legend2 = gsn_add_text(wks,xy,"model_i01.03cn",xlg2_cen+0.15,ylg2_cen,txres)
1.304 +
1.305 + draw(xy)
1.306 + frame(wks)
1.307 +
1.308 +end
1.309 +