1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/npp/16.scatter_bias.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,840 @@
1.4 +; ***********************************************
1.5 +; combine all scatter
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 +; Main code.
1.20 +begin
1.21 +
1.22 + plot_type = "ps"
1.23 + plot_type_new = "png"
1.24 +
1.25 +;************************************************
1.26 +; read in ob data
1.27 +;************************************************
1.28 + dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.29 + fil81 = "data.81.nc"
1.30 + f81 = addfile (dir81+fil81,"r")
1.31 +
1.32 + id81 = f81->SITE_ID
1.33 + npp81 = f81->TNPP_C
1.34 + rain81 = tofloat(f81->PREC_ANN)
1.35 + x81 = f81->LONG_DD
1.36 + y81 = f81->LAT_DD
1.37 +
1.38 + id81@long_name = "SITE_ID"
1.39 + npp81@long_name = "NPP (gC/m2/year)"
1.40 +
1.41 +;change longitude from (-180,180) to (0,360)
1.42 +;for model data interpolation
1.43 +
1.44 + do i= 0,dimsizes(x81)-1
1.45 + if (x81(i) .lt. 0.) then
1.46 + x81(i) = x81(i)+ 360.
1.47 + end if
1.48 + end do
1.49 +;print (x81)
1.50 +;-------------------------------------
1.51 + dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.52 + fil933 = "data.933.nc"
1.53 + f933 = addfile (dir933+fil933,"r")
1.54 +
1.55 + id933 = f933->SITE_ID
1.56 + npp933 = f933->TNPP_C
1.57 + rain933 = f933->PREC
1.58 + x933 = f933->LONG_DD
1.59 + y933 = f933->LAT_DD
1.60 +
1.61 + id933@long_name = "SITE_ID"
1.62 + npp933@long_name = "NPP (gC/m2/year)"
1.63 +
1.64 +;change longitude from (-180,180) to (0,360)
1.65 +;for model data interpolation
1.66 +
1.67 + do i= 0,dimsizes(x933)-1
1.68 + if (x933(i) .lt. 0.) then
1.69 + x933(i) = x933(i)+ 360.
1.70 + end if
1.71 + end do
1.72 +;print (x933)
1.73 +;************************************************
1.74 +; read in model data
1.75 +;************************************************
1.76 + model_name = "b30.061n"
1.77 +
1.78 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.79 + film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
1.80 + fm = addfile (dirm+film,"r")
1.81 +
1.82 + nppmod = fm->NPP
1.83 + rainmod = fm->RAIN
1.84 + xm = fm->lon
1.85 + ym = fm->lat
1.86 +
1.87 + nppmod81 =linint2_points(xm,ym,nppmod(0,:,:),True,x81,y81,0)
1.88 +
1.89 + nppmod933 =linint2_points(xm,ym,nppmod(0,:,:),True,x933,y933,0)
1.90 +
1.91 + rainmod81 =linint2_points(xm,ym,rainmod(0,:,:),True,x81,y81,0)
1.92 +
1.93 + rainmod933=linint2_points(xm,ym,rainmod(0,:,:),True,x933,y933,0)
1.94 +
1.95 +;************************************************
1.96 +; Units for these four variables are:
1.97 +;
1.98 +; rain81 : mm/year
1.99 +; rainmod : mm/s
1.100 +; npp81 : g C/m^2/year
1.101 +; nppmod81: g C/m^2/s
1.102 +;
1.103 +; We want to convert these to "m/year" and "g C/m^2/year".
1.104 +;
1.105 + nsec_per_year = 60*60*24*365
1.106 +
1.107 + rain81 = rain81 / 1000.
1.108 + rainmod81 = (rainmod81/ 1000.) * nsec_per_year
1.109 + nppmod81 = nppmod81 * nsec_per_year
1.110 +
1.111 + rain933 = rain933 / 1000.
1.112 + rainmod933 = (rainmod933/ 1000.) * nsec_per_year
1.113 + nppmod933 = nppmod933 * nsec_per_year
1.114 +
1.115 + npp81@units = "gC/m^2/yr"
1.116 + nppmod81@units = "gC/m^2/yr"
1.117 + npp933@units = "gC/m^2/yr"
1.118 + nppmod933@units = "gC/m^2/yr"
1.119 + rain81@units = "m/yr"
1.120 + rainmod81@units = "m/yr"
1.121 + rain933@units = "m/yr"
1.122 + rainmod933@units = "m/yr"
1.123 +
1.124 + npp81@long_name = "NPP (gC/m2/year)"
1.125 + npp933@long_name = "NPP (gC/m2/year)"
1.126 + nppmod81@long_name = "NPP (gC/m2/year)"
1.127 + nppmod933@long_name = "NPP (gC/m2/year)"
1.128 + rain81@long_name = "PREC (m/year)"
1.129 + rain933@long_name = "PREC (m/year)"
1.130 + rainmod81@long_name = "PREC (m/year)"
1.131 + rainmod933@long_name = "PREC (m/year)"
1.132 +
1.133 +;(A) plot scatter ob 81
1.134 +
1.135 + plot_name = "scatter_ob_81"
1.136 + title = "Observed 81 sites"
1.137 +
1.138 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.139 + res = True ; plot mods desired
1.140 + res@tiMainString = title ; add title
1.141 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.142 + res@xyMarkers = 16 ; choose type of marker
1.143 + res@xyMarkerColor = "red" ; Marker color
1.144 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.145 + res@tmLabelAutoStride = True ; nice tick mark labels
1.146 +
1.147 + plot = gsn_csm_xy (wks,id81,npp81,res) ; create plot
1.148 + frame(wks)
1.149 +
1.150 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.151 + system("rm "+plot_name+"."+plot_type)
1.152 + system("rm "+plot_name+"-1."+plot_type_new)
1.153 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.154 +
1.155 + clear (wks)
1.156 +
1.157 +;(B) plot scatter ob 933
1.158 +
1.159 + plot_name = "scatter_ob_933"
1.160 + title = "Observed 933 sites"
1.161 +
1.162 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.163 + res = True ; plot mods desired
1.164 + res@tiMainString = title ; add title
1.165 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.166 + res@xyMarkers = 16 ; choose type of marker
1.167 + res@xyMarkerColor = "red" ; Marker color
1.168 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.169 + res@tmLabelAutoStride = True ; nice tick mark labels
1.170 +
1.171 + plot = gsn_csm_xy (wks,id933,npp933,res) ; create plot
1.172 + frame(wks)
1.173 +
1.174 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.175 + system("rm "+plot_name+"."+plot_type)
1.176 + system("rm "+plot_name+"-1."+plot_type_new)
1.177 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.178 +
1.179 + clear (wks)
1.180 +
1.181 +;(C) plot scatter model 81
1.182 +
1.183 + plot_name = "scatter_model_81"
1.184 + title = "Model "+ model_name +" 81 sites"
1.185 +
1.186 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.187 + res = True ; plot mods desired
1.188 + res@tiMainString = title ; add title
1.189 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.190 + res@xyMarkers = 16 ; choose type of marker
1.191 + res@xyMarkerColor = "red" ; Marker color
1.192 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.193 + res@tmLabelAutoStride = True ; nice tick mark labels
1.194 +
1.195 + plot = gsn_csm_xy (wks,id81,nppmod81,res) ; create plot
1.196 + frame(wks)
1.197 +
1.198 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.199 + system("rm "+plot_name+"."+plot_type)
1.200 + system("rm "+plot_name+"-1."+plot_type_new)
1.201 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.202 +
1.203 + clear (wks)
1.204 +
1.205 +;(D) plot scatter model 933
1.206 +
1.207 + plot_name = "scatter_model_933"
1.208 + title = "Model "+ model_name +" 933 sites"
1.209 +
1.210 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.211 + res = True ; plot mods desired
1.212 + res@tiMainString = title ; add title
1.213 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.214 + res@xyMarkers = 16 ; choose type of marker
1.215 + res@xyMarkerColor = "red" ; Marker color
1.216 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.217 + res@tmLabelAutoStride = True ; nice tick mark labels
1.218 +
1.219 + plot = gsn_csm_xy (wks,id933,nppmod933,res) ; create plot
1.220 + frame(wks)
1.221 +
1.222 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.223 + system("rm "+plot_name+"."+plot_type)
1.224 + system("rm "+plot_name+"-1."+plot_type_new)
1.225 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.226 +
1.227 + clear (wks)
1.228 +
1.229 +;(E) scatter model vs ob 81
1.230 +
1.231 + plot_name = "scatter_model_vs_ob_81"
1.232 + title = model_name +" vs ob 81 sites"
1.233 +
1.234 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.235 + res = True ; plot mods desired
1.236 + res@tiMainString = title ; add title
1.237 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.238 + res@xyMarkers = 16 ; choose type of marker
1.239 + res@xyMarkerColor = "red" ; Marker color
1.240 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.241 + res@tmLabelAutoStride = True ; nice tick mark labels
1.242 +;-------------------------------
1.243 +;compute correlation coef. and M
1.244 + ccr81 = esccr(nppmod81,npp81,0)
1.245 + print (ccr81)
1.246 +
1.247 +;new eq
1.248 + bias = sum(abs(nppmod81-npp81)/(nppmod81+npp81))
1.249 + M81 = (1. - (bias/dimsizes(y81)))*5.
1.250 + print (M81)
1.251 +
1.252 + tRes = True
1.253 + tRes@txFontHeightF = 0.025
1.254 +
1.255 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
1.256 +
1.257 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.258 +;--------------------------------
1.259 + plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot
1.260 + frame(wks)
1.261 +
1.262 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.263 + system("rm "+plot_name+"."+plot_type)
1.264 + system("rm "+plot_name+"-1."+plot_type_new)
1.265 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.266 +
1.267 + clear (wks)
1.268 +
1.269 +;(F) scatter model vs ob 933
1.270 +
1.271 + plot_name = "scatter_model_vs_ob_933"
1.272 + title = model_name +" vs ob 933 sites"
1.273 +
1.274 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.275 + res = True ; plot mods desired
1.276 + res@tiMainString = title ; add title
1.277 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.278 + res@xyMarkers = 16 ; choose type of marker
1.279 + res@xyMarkerColor = "red" ; Marker color
1.280 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.281 + res@tmLabelAutoStride = True ; nice tick mark labels
1.282 +;-------------------------------
1.283 +;compute correlation coef. and M
1.284 + ccr933 = esccr(nppmod933,npp933,0)
1.285 + print (ccr933)
1.286 +
1.287 +;new eq
1.288 + bias = sum(abs(nppmod933-npp933)/(nppmod933+npp933))
1.289 + M933 = (1. - (bias/dimsizes(y933)))*5.
1.290 + print (M933)
1.291 +
1.292 + tRes = True
1.293 + tRes@txFontHeightF = 0.025
1.294 +
1.295 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
1.296 +
1.297 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.298 +;--------------------------------
1.299 + plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot
1.300 + frame(wks)
1.301 +
1.302 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.303 + system("rm "+plot_name+"."+plot_type)
1.304 + system("rm "+plot_name+"-1."+plot_type_new)
1.305 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.306 +
1.307 + clear (wks)
1.308 +;**************************************************************************
1.309 +;(G) histogram 81
1.310 +
1.311 +;--------------------------------------------------------------------
1.312 +;get data
1.313 +
1.314 + RAIN1_1D = ndtooned(rain81)
1.315 + RAIN2_1D = ndtooned(rainmod81)
1.316 + NPP1_1D = ndtooned(npp81)
1.317 + NPP2_1D = ndtooned(nppmod81)
1.318 +;
1.319 +; Calculate some "nice" bins for binning the data in equally spaced
1.320 +; ranges.
1.321 +;
1.322 + nbins = 15 ; Number of bins to use.
1.323 +
1.324 + nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
1.325 + nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
1.326 + range = fspan(nicevals(0),nicevals(1),nvals)
1.327 +;
1.328 +; Use this range information to grab all the values in a
1.329 +; particular range, and then take an average.
1.330 +;
1.331 + nr = dimsizes(range)
1.332 + nx = nr-1
1.333 + xvalues = new((/2,nx/),typeof(RAIN1_1D))
1.334 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.335 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.336 + dx4 = dx/4 ; 1/4 of the range
1.337 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.338 + yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.339 + mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.340 + mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.341 +
1.342 + do nd=0,1
1.343 +;
1.344 +; See if we are doing model or observational data.
1.345 +;
1.346 + if(nd.eq.0) then
1.347 + data = RAIN1_1D
1.348 + npp_data = NPP1_1D
1.349 + else
1.350 + data = RAIN2_1D
1.351 + npp_data = NPP2_1D
1.352 + end if
1.353 +;
1.354 +; Loop through each range and check for values.
1.355 +;
1.356 + do i=0,nr-2
1.357 + if (i.ne.(nr-2)) then
1.358 +; print("")
1.359 +; print("In range ["+range(i)+","+range(i+1)+")")
1.360 + idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
1.361 + else
1.362 +; print("")
1.363 +; print("In range ["+range(i)+",)")
1.364 + idx = ind(range(i).le.data)
1.365 + end if
1.366 +;
1.367 +; Calculate average, and get min and max.
1.368 +;
1.369 + if(.not.any(ismissing(idx))) then
1.370 + yvalues(nd,i) = avg(npp_data(idx))
1.371 + mn_yvalues(nd,i) = min(npp_data(idx))
1.372 + mx_yvalues(nd,i) = max(npp_data(idx))
1.373 + count = dimsizes(idx)
1.374 + else
1.375 + count = 0
1.376 + yvalues(nd,i) = yvalues@_FillValue
1.377 + mn_yvalues(nd,i) = yvalues@_FillValue
1.378 + mx_yvalues(nd,i) = yvalues@_FillValue
1.379 + end if
1.380 +;
1.381 +; Print out information.
1.382 +;
1.383 +; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
1.384 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.385 +
1.386 +;
1.387 +; Clean up for next time in loop.
1.388 +;
1.389 + delete(idx)
1.390 + end do
1.391 + delete(data)
1.392 + delete(npp_data)
1.393 + end do
1.394 +;----------------------------------------
1.395 +;compute correlation coeff and M score
1.396 + u = yvalues(0,:)
1.397 + v = yvalues(1,:)
1.398 +
1.399 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.400 + uu = u(good)
1.401 + vv = v(good)
1.402 +
1.403 + ccr81h = esccr(uu,vv,0)
1.404 + print (ccr81h)
1.405 +
1.406 +;new eq
1.407 + bias = sum(abs(vv-uu)/(vv+uu))
1.408 + M81h = (1.- (bias/dimsizes(uu)))*5.
1.409 + print (M81h)
1.410 +
1.411 + delete (u)
1.412 + delete (v)
1.413 + delete (uu)
1.414 + delete (vv)
1.415 +;----------------------------------------------------------------------
1.416 +; histogram res
1.417 +
1.418 + res = True
1.419 + res@gsnMaximize = True
1.420 + res@gsnDraw = False
1.421 + res@gsnFrame = False
1.422 + res@xyMarkLineMode = "Markers"
1.423 + res@xyMarkerSizeF = 0.014
1.424 + res@xyMarker = 16
1.425 + res@xyMarkerColors = (/"Brown","Blue"/)
1.426 + res@trYMinF = min(mn_yvalues) - 10.
1.427 + res@trYMaxF = max(mx_yvalues) + 10.
1.428 +
1.429 + res@tiYAxisString = "NPP (g C/m2/year)"
1.430 + res@tiXAxisString = "Precipitation (m/year)"
1.431 +
1.432 + max_bar = new((/2,nx/),graphic)
1.433 + min_bar = new((/2,nx/),graphic)
1.434 + max_cap = new((/2,nx/),graphic)
1.435 + min_cap = new((/2,nx/),graphic)
1.436 +
1.437 + lnres = True
1.438 + line_colors = (/"brown","blue"/)
1.439 +;=================================================================
1.440 +; histogram ob 81 site only
1.441 +;
1.442 + plot_name = "histogram_ob_81"
1.443 + title = "Observed 81 site"
1.444 + res@tiMainString = title
1.445 +
1.446 + wks = gsn_open_wks (plot_type,plot_name)
1.447 +
1.448 + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),res)
1.449 +
1.450 +;-------------------------------
1.451 +;Attach the vertical bar and the horizontal cap line
1.452 +
1.453 + do nd=0,0
1.454 + lnres@gsLineColor = line_colors(nd)
1.455 + do i=0,nx-1
1.456 +
1.457 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.458 + .not.ismissing(mx_yvalues(nd,i))) then
1.459 +;
1.460 +; Attach the vertical bar, both above and below the marker.
1.461 +;
1.462 + x1 = xvalues(nd,i)
1.463 + y1 = yvalues(nd,i)
1.464 + y2 = mn_yvalues(nd,i)
1.465 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.466 +
1.467 + y2 = mx_yvalues(nd,i)
1.468 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.469 +;
1.470 +; Attach the horizontal cap line, both above and below the marker.
1.471 +;
1.472 + x1 = xvalues(nd,i) - dx4
1.473 + x2 = xvalues(nd,i) + dx4
1.474 + y1 = mn_yvalues(nd,i)
1.475 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.476 +
1.477 + y1 = mx_yvalues(nd,i)
1.478 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.479 + end if
1.480 + end do
1.481 + end do
1.482 +
1.483 + draw(xy)
1.484 + frame(wks)
1.485 +
1.486 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.487 +; system("rm "+plot_name+"."+plot_type)
1.488 +; system("rm "+plot_name+"-1."+plot_type_new)
1.489 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.490 +
1.491 + clear (wks)
1.492 +;===========================================================================
1.493 +; histogram model vs ob 81 site
1.494 +
1.495 + plot_name = "histogram_mod-ob_81"
1.496 + title = model_name+ " vs Observed 81 site"
1.497 + res@tiMainString = title
1.498 +
1.499 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.500 +
1.501 +;-----------------------------
1.502 +; Add a boxed legend using the more simple method, which won't have
1.503 +; vertical lines going through the markers.
1.504 +
1.505 + res@pmLegendDisplayMode = "Always"
1.506 +; res@pmLegendWidthF = 0.1
1.507 + res@pmLegendWidthF = 0.08
1.508 + res@pmLegendHeightF = 0.05
1.509 + res@pmLegendOrthogonalPosF = -1.17
1.510 +; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.511 +; res@pmLegendParallelPosF = 0.18
1.512 + res@pmLegendParallelPosF = 0.23 ;(rightward)
1.513 +
1.514 +; res@lgPerimOn = False
1.515 + res@lgLabelFontHeightF = 0.015
1.516 + res@xyExplicitLegendLabels = (/"observed",model_name/)
1.517 +;-----------------------------
1.518 + tRes = True
1.519 + tRes@txFontHeightF = 0.025
1.520 +
1.521 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"
1.522 +
1.523 + gsn_text_ndc(wks,correlation_text,0.5,0.8,tRes)
1.524 +
1.525 + xy = gsn_csm_xy(wks,xvalues,yvalues,res)
1.526 +;-------------------------------
1.527 +;Attach the vertical bar and the horizontal cap line
1.528 +
1.529 + do nd=0,1
1.530 + lnres@gsLineColor = line_colors(nd)
1.531 + do i=0,nx-1
1.532 +
1.533 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.534 + .not.ismissing(mx_yvalues(nd,i))) then
1.535 +;
1.536 +; Attach the vertical bar, both above and below the marker.
1.537 +;
1.538 + x1 = xvalues(nd,i)
1.539 + y1 = yvalues(nd,i)
1.540 + y2 = mn_yvalues(nd,i)
1.541 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.542 +
1.543 + y2 = mx_yvalues(nd,i)
1.544 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.545 +;
1.546 +; Attach the horizontal cap line, both above and below the marker.
1.547 +;
1.548 + x1 = xvalues(nd,i) - dx4
1.549 + x2 = xvalues(nd,i) + dx4
1.550 + y1 = mn_yvalues(nd,i)
1.551 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.552 +
1.553 + y1 = mx_yvalues(nd,i)
1.554 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.555 + end if
1.556 + end do
1.557 + end do
1.558 + draw(xy)
1.559 + frame(wks)
1.560 +
1.561 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.562 +;system("rm "+plot_name+"."+plot_type)
1.563 +;system("rm "+plot_name+"-1."+plot_type_new)
1.564 +;system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.565 +
1.566 + clear (wks)
1.567 +
1.568 + delete (RAIN1_1D)
1.569 + delete (RAIN2_1D)
1.570 + delete (NPP1_1D)
1.571 + delete (NPP2_1D)
1.572 + delete (range)
1.573 + delete (xvalues)
1.574 + delete (yvalues)
1.575 + delete (mn_yvalues)
1.576 + delete (mx_yvalues)
1.577 + delete (good)
1.578 + delete (max_bar)
1.579 + delete (min_bar)
1.580 + delete (max_cap)
1.581 + delete (min_cap)
1.582 +;**************************************************************************
1.583 +;(H) histogram 933
1.584 +
1.585 +;--------------------------------------------------------------------
1.586 +;get data
1.587 +
1.588 + RAIN1_1D = ndtooned(rain933)
1.589 + RAIN2_1D = ndtooned(rainmod933)
1.590 + NPP1_1D = ndtooned(npp933)
1.591 + NPP2_1D = ndtooned(nppmod933)
1.592 +;
1.593 +; Calculate some "nice" bins for binning the data in equally spaced
1.594 +; ranges.
1.595 +;
1.596 + nbins = 15 ; Number of bins to use.
1.597 +
1.598 + nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
1.599 + nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
1.600 + range = fspan(nicevals(0),nicevals(1),nvals)
1.601 +;
1.602 +; Use this range information to grab all the values in a
1.603 +; particular range, and then take an average.
1.604 +;
1.605 + nr = dimsizes(range)
1.606 + nx = nr-1
1.607 + xvalues = new((/2,nx/),typeof(RAIN1_1D))
1.608 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.609 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.610 + dx4 = dx/4 ; 1/4 of the range
1.611 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.612 + yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.613 + mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.614 + mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.615 +
1.616 + do nd=0,1
1.617 +;
1.618 +; See if we are doing model or observational data.
1.619 +;
1.620 + if(nd.eq.0) then
1.621 + data = RAIN1_1D
1.622 + npp_data = NPP1_1D
1.623 + else
1.624 + data = RAIN2_1D
1.625 + npp_data = NPP2_1D
1.626 + end if
1.627 +;
1.628 +; Loop through each range and check for values.
1.629 +;
1.630 + do i=0,nr-2
1.631 + if (i.ne.(nr-2)) then
1.632 +; print("")
1.633 +; print("In range ["+range(i)+","+range(i+1)+")")
1.634 + idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
1.635 + else
1.636 +; print("")
1.637 +; print("In range ["+range(i)+",)")
1.638 + idx = ind(range(i).le.data)
1.639 + end if
1.640 +;
1.641 +; Calculate average, and get min and max.
1.642 +;
1.643 + if(.not.any(ismissing(idx))) then
1.644 + yvalues(nd,i) = avg(npp_data(idx))
1.645 + mn_yvalues(nd,i) = min(npp_data(idx))
1.646 + mx_yvalues(nd,i) = max(npp_data(idx))
1.647 + count = dimsizes(idx)
1.648 + else
1.649 + count = 0
1.650 + yvalues(nd,i) = yvalues@_FillValue
1.651 + mn_yvalues(nd,i) = yvalues@_FillValue
1.652 + mx_yvalues(nd,i) = yvalues@_FillValue
1.653 + end if
1.654 +;
1.655 +; Print out information.
1.656 +;
1.657 +; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
1.658 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.659 +
1.660 +;
1.661 +; Clean up for next time in loop.
1.662 +;
1.663 + delete(idx)
1.664 + end do
1.665 + delete(data)
1.666 + delete(npp_data)
1.667 + end do
1.668 +;----------------------------------------
1.669 +;compute correlation coeff and M score
1.670 + u = yvalues(0,:)
1.671 + v = yvalues(1,:)
1.672 +
1.673 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.674 + uu = u(good)
1.675 + vv = v(good)
1.676 +
1.677 + ccr933h = esccr(uu,vv,0)
1.678 + print (ccr933h)
1.679 +
1.680 +;new eq
1.681 + bias = sum(abs(vv-uu)/(vv+uu))
1.682 + M933h = (1.- (bias/dimsizes(uu)))*5.
1.683 + print (M933h)
1.684 +
1.685 + delete (u)
1.686 + delete (v)
1.687 + delete (uu)
1.688 + delete (vv)
1.689 +;----------------------------------------------------------------------
1.690 +; histogram res
1.691 +
1.692 + res = True
1.693 + res@gsnMaximize = True
1.694 + res@gsnDraw = False
1.695 + res@gsnFrame = False
1.696 + res@xyMarkLineMode = "Markers"
1.697 + res@xyMarkerSizeF = 0.014
1.698 + res@xyMarker = 16
1.699 + res@xyMarkerColors = (/"Brown","Blue"/)
1.700 + res@trYMinF = min(mn_yvalues) - 10.
1.701 + res@trYMaxF = max(mx_yvalues) + 10.
1.702 +
1.703 + res@tiYAxisString = "NPP (g C/m2/year)"
1.704 + res@tiXAxisString = "Precipitation (m/year)"
1.705 +
1.706 + max_bar = new((/2,nx/),graphic)
1.707 + min_bar = new((/2,nx/),graphic)
1.708 + max_cap = new((/2,nx/),graphic)
1.709 + min_cap = new((/2,nx/),graphic)
1.710 +
1.711 + lnres = True
1.712 + line_colors = (/"brown","blue"/)
1.713 +;=================================================================
1.714 +; histogram ob 933 site only
1.715 +;
1.716 + plot_name = "histogram_ob_933"
1.717 + title = "Observed 933 site"
1.718 + res@tiMainString = title
1.719 +
1.720 + wks = gsn_open_wks (plot_type,plot_name)
1.721 +
1.722 + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),res)
1.723 +
1.724 +;-------------------------------
1.725 +;Attach the vertical bar and the horizontal cap line
1.726 +
1.727 + do nd=0,0
1.728 + lnres@gsLineColor = line_colors(nd)
1.729 + do i=0,nx-1
1.730 +
1.731 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.732 + .not.ismissing(mx_yvalues(nd,i))) then
1.733 +;
1.734 +; Attach the vertical bar, both above and below the marker.
1.735 +;
1.736 + x1 = xvalues(nd,i)
1.737 + y1 = yvalues(nd,i)
1.738 + y2 = mn_yvalues(nd,i)
1.739 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.740 +
1.741 + y2 = mx_yvalues(nd,i)
1.742 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.743 +;
1.744 +; Attach the horizontal cap line, both above and below the marker.
1.745 +;
1.746 + x1 = xvalues(nd,i) - dx4
1.747 + x2 = xvalues(nd,i) + dx4
1.748 + y1 = mn_yvalues(nd,i)
1.749 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.750 +
1.751 + y1 = mx_yvalues(nd,i)
1.752 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.753 + end if
1.754 + end do
1.755 + end do
1.756 +
1.757 + draw(xy)
1.758 + frame(wks)
1.759 +
1.760 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.761 +; system("rm "+plot_name+"."+plot_type)
1.762 +; system("rm "+plot_name+"-1."+plot_type_new)
1.763 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.764 +
1.765 + clear (wks)
1.766 +;===========================================================================
1.767 +; histogram model vs ob 933 site
1.768 +
1.769 + plot_name = "histogram_mod-ob_933"
1.770 + title = model_name+ " vs Observed 933 site"
1.771 + res@tiMainString = title
1.772 +
1.773 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.774 +
1.775 +;-----------------------------
1.776 +; Add a boxed legend using the more simple method, which won't have
1.777 +; vertical lines going through the markers.
1.778 +
1.779 + res@pmLegendDisplayMode = "Always"
1.780 +; res@pmLegendWidthF = 0.1
1.781 + res@pmLegendWidthF = 0.08
1.782 + res@pmLegendHeightF = 0.05
1.783 + res@pmLegendOrthogonalPosF = -1.17
1.784 +; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.785 +; res@pmLegendParallelPosF = 0.18
1.786 + res@pmLegendParallelPosF = 0.23 ;(rightward)
1.787 +
1.788 +; res@lgPerimOn = False
1.789 + res@lgLabelFontHeightF = 0.015
1.790 + res@xyExplicitLegendLabels = (/"observed",model_name/)
1.791 +;-----------------------------
1.792 + tRes = True
1.793 + tRes@txFontHeightF = 0.025
1.794 +
1.795 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"
1.796 +
1.797 + gsn_text_ndc(wks,correlation_text,0.5,0.8,tRes)
1.798 +
1.799 + xy = gsn_csm_xy(wks,xvalues,yvalues,res)
1.800 +;-------------------------------
1.801 +;Attach the vertical bar and the horizontal cap line
1.802 +
1.803 + do nd=0,1
1.804 + lnres@gsLineColor = line_colors(nd)
1.805 + do i=0,nx-1
1.806 +
1.807 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.808 + .not.ismissing(mx_yvalues(nd,i))) then
1.809 +;
1.810 +; Attach the vertical bar, both above and below the marker.
1.811 +;
1.812 + x1 = xvalues(nd,i)
1.813 + y1 = yvalues(nd,i)
1.814 + y2 = mn_yvalues(nd,i)
1.815 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.816 +
1.817 + y2 = mx_yvalues(nd,i)
1.818 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.819 +;
1.820 +; Attach the horizontal cap line, both above and below the marker.
1.821 +;
1.822 + x1 = xvalues(nd,i) - dx4
1.823 + x2 = xvalues(nd,i) + dx4
1.824 + y1 = mn_yvalues(nd,i)
1.825 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.826 +
1.827 + y1 = mx_yvalues(nd,i)
1.828 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.829 + end if
1.830 + end do
1.831 + end do
1.832 + draw(xy)
1.833 + frame(wks)
1.834 +
1.835 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.836 +;system("rm "+plot_name+"."+plot_type)
1.837 +;system("rm "+plot_name+"-1."+plot_type_new)
1.838 +;system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.839 +
1.840 + clear (wks)
1.841 +;------------------------------------------------------------------------
1.842 +
1.843 +end