1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/npp/18.scatter_bias.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,1096 @@
1.4 +; ****************************************************
1.5 +; combine scatter, histogram, global and zonal plots
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 +;(1)
1.29 + dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.30 + fil81 = "data.81.nc"
1.31 + f81 = addfile (dir81+fil81,"r")
1.32 +
1.33 + id81 = f81->SITE_ID
1.34 + npp81 = f81->TNPP_C
1.35 + rain81 = tofloat(f81->PREC_ANN)
1.36 + x81 = f81->LONG_DD
1.37 + y81 = f81->LAT_DD
1.38 +
1.39 + id81@long_name = "SITE_ID"
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 +;(2)
1.52 + dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.53 + fil933 = "data.933.nc"
1.54 + f933 = addfile (dir933+fil933,"r")
1.55 +
1.56 + id933 = f933->SITE_ID
1.57 + npp933 = f933->TNPP_C
1.58 + rain933 = f933->PREC
1.59 + x933 = f933->LONG_DD
1.60 + y933 = f933->LAT_DD
1.61 +
1.62 + id933@long_name = "SITE_ID"
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 +;(3)
1.75 + dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
1.76 + filglobe = "Npp_T31_mean.nc"
1.77 + fglobe = addfile (dirglobe+filglobe,"r")
1.78 +
1.79 + nppglobe0 = fglobe->NPP
1.80 + nppglobe = nppglobe0
1.81 +;************************************************
1.82 +; read in model data
1.83 +;************************************************
1.84 + model_name = "b30.061n"
1.85 +
1.86 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.87 + film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
1.88 + fm = addfile (dirm+film,"r")
1.89 +
1.90 + nppmod0 = fm->NPP
1.91 + rainmod0 = fm->RAIN
1.92 + xm = fm->lon
1.93 + ym = fm->lat
1.94 +
1.95 + nppmod = nppmod0(0,:,:)
1.96 + rainmod = rainmod0(0,:,:)
1.97 + delete (nppmod0)
1.98 + delete (rainmod0)
1.99 +
1.100 + nppmod81 =linint2_points(xm,ym,nppmod,True,x81,y81,0)
1.101 +
1.102 + nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)
1.103 +
1.104 + rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)
1.105 +
1.106 + rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)
1.107 +
1.108 +;************************************************
1.109 +; Units for these variables are:
1.110 +;
1.111 +; rain81 : mm/year
1.112 +; rainmod : mm/s
1.113 +; npp81 : g C/m^2/year
1.114 +; nppmod81: g C/m^2/s
1.115 +; nppglobe: g C/m^2/year
1.116 +;
1.117 +; We want to convert these to "m/year" and "g C/m^2/year".
1.118 +;
1.119 + nsec_per_year = 60*60*24*365
1.120 +
1.121 + rain81 = rain81 / 1000.
1.122 + rainmod81 = (rainmod81/ 1000.) * nsec_per_year
1.123 + nppmod81 = nppmod81 * nsec_per_year
1.124 +
1.125 + rain933 = rain933 / 1000.
1.126 + rainmod933 = (rainmod933/ 1000.) * nsec_per_year
1.127 + nppmod933 = nppmod933 * nsec_per_year
1.128 +
1.129 + nppmod = nppmod * nsec_per_year
1.130 +
1.131 + npp81@units = "gC/m^2/yr"
1.132 + nppmod81@units = "gC/m^2/yr"
1.133 + npp933@units = "gC/m^2/yr"
1.134 + nppmod933@units = "gC/m^2/yr"
1.135 + nppmod@units = "gC/m^2/yr"
1.136 + nppglobe@units = "gC/m^2/yr"
1.137 + rain81@units = "m/yr"
1.138 + rainmod81@units = "m/yr"
1.139 + rain933@units = "m/yr"
1.140 + rainmod933@units = "m/yr"
1.141 +
1.142 + npp81@long_name = "NPP (gC/m2/year)"
1.143 + npp933@long_name = "NPP (gC/m2/year)"
1.144 + nppmod81@long_name = "NPP (gC/m2/year)"
1.145 + nppmod933@long_name = "NPP (gC/m2/year)"
1.146 + nppmod@long_name = "NPP (gC/m2/year)"
1.147 + nppglobe@long_name = "NPP (gC/m2/year)"
1.148 + rain81@long_name = "PREC (m/year)"
1.149 + rain933@long_name = "PREC (m/year)"
1.150 + rainmod81@long_name = "PREC (m/year)"
1.151 + rainmod933@long_name = "PREC (m/year)"
1.152 +
1.153 +;(A) plot scatter ob 81
1.154 +
1.155 + plot_name = "scatter_ob_81"
1.156 + title = "Observed 81 sites"
1.157 +
1.158 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.159 + res = True ; plot mods desired
1.160 + res@tiMainString = title ; add title
1.161 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.162 + res@xyMarkers = 16 ; choose type of marker
1.163 + res@xyMarkerColor = "red" ; Marker color
1.164 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.165 + res@tmLabelAutoStride = True ; nice tick mark labels
1.166 +
1.167 + plot = gsn_csm_xy (wks,id81,npp81,res) ; create plot
1.168 + frame(wks)
1.169 +
1.170 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.171 + system("rm "+plot_name+"."+plot_type)
1.172 + system("rm "+plot_name+"-1."+plot_type_new)
1.173 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.174 +
1.175 + clear (wks)
1.176 +
1.177 +;(B) plot scatter ob 933
1.178 +
1.179 + plot_name = "scatter_ob_933"
1.180 + title = "Observed 933 sites"
1.181 +
1.182 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.183 + res = True ; plot mods desired
1.184 + res@tiMainString = title ; add title
1.185 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.186 + res@xyMarkers = 16 ; choose type of marker
1.187 + res@xyMarkerColor = "red" ; Marker color
1.188 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.189 + res@tmLabelAutoStride = True ; nice tick mark labels
1.190 +
1.191 + plot = gsn_csm_xy (wks,id933,npp933,res) ; create plot
1.192 + frame(wks)
1.193 +
1.194 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.195 + system("rm "+plot_name+"."+plot_type)
1.196 + system("rm "+plot_name+"-1."+plot_type_new)
1.197 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.198 +
1.199 + clear (wks)
1.200 +
1.201 +;(C) plot scatter model 81
1.202 +
1.203 + plot_name = "scatter_model_81"
1.204 + title = "Model "+ model_name +" 81 sites"
1.205 +
1.206 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.207 + res = True ; plot mods desired
1.208 + res@tiMainString = title ; add title
1.209 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.210 + res@xyMarkers = 16 ; choose type of marker
1.211 + res@xyMarkerColor = "red" ; Marker color
1.212 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.213 + res@tmLabelAutoStride = True ; nice tick mark labels
1.214 +
1.215 + plot = gsn_csm_xy (wks,id81,nppmod81,res) ; create plot
1.216 + frame(wks)
1.217 +
1.218 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.219 + system("rm "+plot_name+"."+plot_type)
1.220 + system("rm "+plot_name+"-1."+plot_type_new)
1.221 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.222 +
1.223 + clear (wks)
1.224 +
1.225 +;(D) plot scatter model 933
1.226 +
1.227 + plot_name = "scatter_model_933"
1.228 + title = "Model "+ model_name +" 933 sites"
1.229 +
1.230 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.231 + res = True ; plot mods desired
1.232 + res@tiMainString = title ; add title
1.233 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.234 + res@xyMarkers = 16 ; choose type of marker
1.235 + res@xyMarkerColor = "red" ; Marker color
1.236 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.237 + res@tmLabelAutoStride = True ; nice tick mark labels
1.238 +
1.239 + plot = gsn_csm_xy (wks,id933,nppmod933,res) ; create plot
1.240 + frame(wks)
1.241 +
1.242 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.243 + system("rm "+plot_name+"."+plot_type)
1.244 + system("rm "+plot_name+"-1."+plot_type_new)
1.245 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.246 +
1.247 + clear (wks)
1.248 +
1.249 +;(E) scatter model vs ob 81
1.250 +
1.251 + plot_name = "scatter_model_vs_ob_81"
1.252 + title = model_name +" vs ob 81 sites"
1.253 +
1.254 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.255 + res = True ; plot mods desired
1.256 + res@tiMainString = title ; add title
1.257 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.258 + res@xyMarkers = 16 ; choose type of marker
1.259 + res@xyMarkerColor = "red" ; Marker color
1.260 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.261 + res@tmLabelAutoStride = True ; nice tick mark labels
1.262 +;-------------------------------
1.263 +;compute correlation coef. and M
1.264 + ccr81 = esccr(nppmod81,npp81,0)
1.265 + print (ccr81)
1.266 +
1.267 +;new eq
1.268 + bias = sum(abs(nppmod81-npp81)/(nppmod81+npp81))
1.269 + M81 = (1. - (bias/dimsizes(y81)))*5.
1.270 + print (M81)
1.271 +
1.272 + tRes = True
1.273 + tRes@txFontHeightF = 0.025
1.274 +
1.275 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
1.276 +
1.277 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.278 +;--------------------------------
1.279 + plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot
1.280 + frame(wks)
1.281 +
1.282 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.283 + system("rm "+plot_name+"."+plot_type)
1.284 + system("rm "+plot_name+"-1."+plot_type_new)
1.285 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.286 +
1.287 + clear (wks)
1.288 +
1.289 +;(F) scatter model vs ob 933
1.290 +
1.291 + plot_name = "scatter_model_vs_ob_933"
1.292 + title = model_name +" vs ob 933 sites"
1.293 +
1.294 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.295 + res = True ; plot mods desired
1.296 + res@tiMainString = title ; add title
1.297 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.298 + res@xyMarkers = 16 ; choose type of marker
1.299 + res@xyMarkerColor = "red" ; Marker color
1.300 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.301 + res@tmLabelAutoStride = True ; nice tick mark labels
1.302 +;-------------------------------
1.303 +;compute correlation coef. and M
1.304 + ccr933 = esccr(nppmod933,npp933,0)
1.305 + print (ccr933)
1.306 +
1.307 +;new eq
1.308 + bias = sum(abs(nppmod933-npp933)/(nppmod933+npp933))
1.309 + M933 = (1. - (bias/dimsizes(y933)))*5.
1.310 + print (M933)
1.311 +
1.312 + tRes = True
1.313 + tRes@txFontHeightF = 0.025
1.314 +
1.315 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
1.316 +
1.317 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.318 +;--------------------------------
1.319 + plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot
1.320 + frame(wks)
1.321 +
1.322 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.323 + system("rm "+plot_name+"."+plot_type)
1.324 + system("rm "+plot_name+"-1."+plot_type_new)
1.325 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.326 +
1.327 + clear (wks)
1.328 +;**************************************************************************
1.329 +;(G) histogram 81
1.330 +
1.331 +;--------------------------------------------------------------------
1.332 +;get data
1.333 +
1.334 + RAIN1_1D = ndtooned(rain81)
1.335 + RAIN2_1D = ndtooned(rainmod81)
1.336 + NPP1_1D = ndtooned(npp81)
1.337 + NPP2_1D = ndtooned(nppmod81)
1.338 +
1.339 +; Calculaee "nice" bins for binning the data in equally spaced ranges.
1.340 +
1.341 + nbins = 15 ; Number of bins to use.
1.342 +
1.343 + nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
1.344 + nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
1.345 + range = fspan(nicevals(0),nicevals(1),nvals)
1.346 +
1.347 +; Use this range information to grab all the values in a
1.348 +; particular range, and then take an average.
1.349 +
1.350 + nr = dimsizes(range)
1.351 + nx = nr-1
1.352 + xvalues = new((/2,nx/),typeof(RAIN1_1D))
1.353 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.354 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.355 + dx4 = dx/4 ; 1/4 of the range
1.356 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.357 + yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.358 + mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.359 + mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.360 +
1.361 + do nd=0,1
1.362 +
1.363 +; See if we are doing model or observational data.
1.364 +
1.365 + if(nd.eq.0) then
1.366 + data = RAIN1_1D
1.367 + npp_data = NPP1_1D
1.368 + else
1.369 + data = RAIN2_1D
1.370 + npp_data = NPP2_1D
1.371 + end if
1.372 +
1.373 +; Loop through each range and check for values.
1.374 +
1.375 + do i=0,nr-2
1.376 + if (i.ne.(nr-2)) then
1.377 +; print("")
1.378 +; print("In range ["+range(i)+","+range(i+1)+")")
1.379 + idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
1.380 + else
1.381 +; print("")
1.382 +; print("In range ["+range(i)+",)")
1.383 + idx = ind(range(i).le.data)
1.384 + end if
1.385 +
1.386 +; Calculate average, and get min and max.
1.387 +
1.388 + if(.not.any(ismissing(idx))) then
1.389 + yvalues(nd,i) = avg(npp_data(idx))
1.390 + mn_yvalues(nd,i) = min(npp_data(idx))
1.391 + mx_yvalues(nd,i) = max(npp_data(idx))
1.392 + count = dimsizes(idx)
1.393 + else
1.394 + count = 0
1.395 + yvalues(nd,i) = yvalues@_FillValue
1.396 + mn_yvalues(nd,i) = yvalues@_FillValue
1.397 + mx_yvalues(nd,i) = yvalues@_FillValue
1.398 + end if
1.399 +
1.400 +; Print out information.
1.401 +
1.402 +; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
1.403 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.404 +
1.405 +
1.406 +; Clean up for next time in loop.
1.407 +
1.408 + delete(idx)
1.409 + end do
1.410 + delete(data)
1.411 + delete(npp_data)
1.412 + end do
1.413 +;----------------------------------------
1.414 +;compute correlation coeff and M score
1.415 +
1.416 + u = yvalues(0,:)
1.417 + v = yvalues(1,:)
1.418 +
1.419 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.420 + uu = u(good)
1.421 + vv = v(good)
1.422 +
1.423 + ccr81h = esccr(uu,vv,0)
1.424 + print (ccr81h)
1.425 +
1.426 +;new eq
1.427 + bias = sum(abs(vv-uu)/(vv+uu))
1.428 + M81h = (1.- (bias/dimsizes(uu)))*5.
1.429 + print (M81h)
1.430 +
1.431 + delete (u)
1.432 + delete (v)
1.433 + delete (uu)
1.434 + delete (vv)
1.435 +;----------------------------------------------------------------------
1.436 +; histogram res
1.437 +
1.438 + resh = True
1.439 + resh@gsnMaximize = True
1.440 + resh@gsnDraw = False
1.441 + resh@gsnFrame = False
1.442 + resh@xyMarkLineMode = "Markers"
1.443 + resh@xyMarkerSizeF = 0.014
1.444 + resh@xyMarker = 16
1.445 + resh@xyMarkerColors = (/"Brown","Blue"/)
1.446 + resh@trYMinF = min(mn_yvalues) - 10.
1.447 + resh@trYMaxF = max(mx_yvalues) + 10.
1.448 +
1.449 + resh@tiYAxisString = "NPP (g C/m2/year)"
1.450 + resh@tiXAxisString = "Precipitation (m/year)"
1.451 +
1.452 + max_bar = new((/2,nx/),graphic)
1.453 + min_bar = new((/2,nx/),graphic)
1.454 + max_cap = new((/2,nx/),graphic)
1.455 + min_cap = new((/2,nx/),graphic)
1.456 +
1.457 + lnres = True
1.458 + line_colors = (/"brown","blue"/)
1.459 +;=================================================================
1.460 +; histogram ob 81 site only
1.461 +;
1.462 + plot_name = "histogram_ob_81"
1.463 + title = "Observed 81 site"
1.464 + resh@tiMainString = title
1.465 +
1.466 + wks = gsn_open_wks (plot_type,plot_name)
1.467 +
1.468 + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
1.469 +
1.470 +;-------------------------------
1.471 +;Attach the vertical bar and the horizontal cap line
1.472 +
1.473 + do nd=0,0
1.474 + lnres@gsLineColor = line_colors(nd)
1.475 + do i=0,nx-1
1.476 +
1.477 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.478 + .not.ismissing(mx_yvalues(nd,i))) then
1.479 +;
1.480 +; Attach the vertical bar, both above and below the marker.
1.481 +;
1.482 + x1 = xvalues(nd,i)
1.483 + y1 = yvalues(nd,i)
1.484 + y2 = mn_yvalues(nd,i)
1.485 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.486 +
1.487 + y2 = mx_yvalues(nd,i)
1.488 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.489 +;
1.490 +; Attach the horizontal cap line, both above and below the marker.
1.491 +;
1.492 + x1 = xvalues(nd,i) - dx4
1.493 + x2 = xvalues(nd,i) + dx4
1.494 + y1 = mn_yvalues(nd,i)
1.495 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.496 +
1.497 + y1 = mx_yvalues(nd,i)
1.498 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.499 + end if
1.500 + end do
1.501 + end do
1.502 +
1.503 + draw(xy)
1.504 + frame(wks)
1.505 +
1.506 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.507 + system("rm "+plot_name+"."+plot_type)
1.508 +; system("rm "+plot_name+"-1."+plot_type_new)
1.509 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.510 +
1.511 + clear (wks)
1.512 +;===========================================================================
1.513 +; histogram model vs ob 81 site
1.514 +
1.515 + plot_name = "histogram_mod-ob_81"
1.516 + title = model_name+ " vs Observed 81 site"
1.517 + resh@tiMainString = title
1.518 +
1.519 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.520 +
1.521 +;-----------------------------
1.522 +; Add a boxed legend using the more simple method, which won't have
1.523 +; vertical lines going through the markers.
1.524 +
1.525 + resh@pmLegendDisplayMode = "Always"
1.526 +; resh@pmLegendWidthF = 0.1
1.527 + resh@pmLegendWidthF = 0.08
1.528 + resh@pmLegendHeightF = 0.05
1.529 + resh@pmLegendOrthogonalPosF = -1.17
1.530 +; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.531 +; resh@pmLegendParallelPosF = 0.18
1.532 + resh@pmLegendParallelPosF = 0.88 ;(rightward)
1.533 +
1.534 +; resh@lgPerimOn = False
1.535 + resh@lgLabelFontHeightF = 0.015
1.536 + resh@xyExplicitLegendLabels = (/"observed",model_name/)
1.537 +;-----------------------------
1.538 + tRes = True
1.539 + tRes@txFontHeightF = 0.025
1.540 +
1.541 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"
1.542 +
1.543 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.544 +
1.545 + xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
1.546 +;-------------------------------
1.547 +;Attach the vertical bar and the horizontal cap line
1.548 +
1.549 + do nd=0,1
1.550 + lnres@gsLineColor = line_colors(nd)
1.551 + do i=0,nx-1
1.552 +
1.553 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.554 + .not.ismissing(mx_yvalues(nd,i))) then
1.555 +;
1.556 +; Attach the vertical bar, both above and below the marker.
1.557 +;
1.558 + x1 = xvalues(nd,i)
1.559 + y1 = yvalues(nd,i)
1.560 + y2 = mn_yvalues(nd,i)
1.561 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.562 +
1.563 + y2 = mx_yvalues(nd,i)
1.564 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.565 +;
1.566 +; Attach the horizontal cap line, both above and below the marker.
1.567 +;
1.568 + x1 = xvalues(nd,i) - dx4
1.569 + x2 = xvalues(nd,i) + dx4
1.570 + y1 = mn_yvalues(nd,i)
1.571 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.572 +
1.573 + y1 = mx_yvalues(nd,i)
1.574 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.575 + end if
1.576 + end do
1.577 + end do
1.578 + draw(xy)
1.579 + frame(wks)
1.580 +
1.581 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.582 + system("rm "+plot_name+"."+plot_type)
1.583 +; system("rm "+plot_name+"-1."+plot_type_new)
1.584 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.585 +
1.586 + clear (wks)
1.587 +
1.588 + delete (RAIN1_1D)
1.589 + delete (RAIN2_1D)
1.590 + delete (NPP1_1D)
1.591 + delete (NPP2_1D)
1.592 + delete (range)
1.593 + delete (xvalues)
1.594 + delete (yvalues)
1.595 + delete (mn_yvalues)
1.596 + delete (mx_yvalues)
1.597 + delete (good)
1.598 + delete (max_bar)
1.599 + delete (min_bar)
1.600 + delete (max_cap)
1.601 + delete (min_cap)
1.602 +;**************************************************************************
1.603 +;(H) histogram 933
1.604 +
1.605 +;--------------------------------------------------------------------
1.606 +;get data
1.607 +
1.608 + RAIN1_1D = ndtooned(rain933)
1.609 + RAIN2_1D = ndtooned(rainmod933)
1.610 + NPP1_1D = ndtooned(npp933)
1.611 + NPP2_1D = ndtooned(nppmod933)
1.612 +
1.613 +; Calculate "nice" bins for binning the data in equally spaced ranges.
1.614 +
1.615 + nbins = 15 ; Number of bins to use.
1.616 +
1.617 + nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
1.618 + nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
1.619 + range = fspan(nicevals(0),nicevals(1),nvals)
1.620 +
1.621 +; Use this range information to grab all the values in a
1.622 +; particular range, and then take an average.
1.623 +
1.624 + nr = dimsizes(range)
1.625 + nx = nr-1
1.626 + xvalues = new((/2,nx/),typeof(RAIN1_1D))
1.627 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.628 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.629 + dx4 = dx/4 ; 1/4 of the range
1.630 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.631 + yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.632 + mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.633 + mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.634 +
1.635 + do nd=0,1
1.636 +
1.637 +; See if we are doing model or observational data.
1.638 +
1.639 + if(nd.eq.0) then
1.640 + data = RAIN1_1D
1.641 + npp_data = NPP1_1D
1.642 + else
1.643 + data = RAIN2_1D
1.644 + npp_data = NPP2_1D
1.645 + end if
1.646 +
1.647 +; Loop through each range and check for values.
1.648 +
1.649 + do i=0,nr-2
1.650 + if (i.ne.(nr-2)) then
1.651 +; print("")
1.652 +; print("In range ["+range(i)+","+range(i+1)+")")
1.653 + idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
1.654 + else
1.655 +; print("")
1.656 +; print("In range ["+range(i)+",)")
1.657 + idx = ind(range(i).le.data)
1.658 + end if
1.659 +
1.660 +; Calculate average, and get min and max.
1.661 +
1.662 + if(.not.any(ismissing(idx))) then
1.663 + yvalues(nd,i) = avg(npp_data(idx))
1.664 + mn_yvalues(nd,i) = min(npp_data(idx))
1.665 + mx_yvalues(nd,i) = max(npp_data(idx))
1.666 + count = dimsizes(idx)
1.667 + else
1.668 + count = 0
1.669 + yvalues(nd,i) = yvalues@_FillValue
1.670 + mn_yvalues(nd,i) = yvalues@_FillValue
1.671 + mx_yvalues(nd,i) = yvalues@_FillValue
1.672 + end if
1.673 +
1.674 +; Print out information.
1.675 +
1.676 +; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
1.677 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.678 +
1.679 +; Clean up for next time in loop.
1.680 +
1.681 + delete(idx)
1.682 + end do
1.683 + delete(data)
1.684 + delete(npp_data)
1.685 + end do
1.686 +;----------------------------------------
1.687 +;compute correlation coeff and M score
1.688 +
1.689 + u = yvalues(0,:)
1.690 + v = yvalues(1,:)
1.691 +
1.692 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.693 + uu = u(good)
1.694 + vv = v(good)
1.695 +
1.696 + ccr933h = esccr(uu,vv,0)
1.697 + print (ccr933h)
1.698 +
1.699 +;new eq
1.700 + bias = sum(abs(vv-uu)/(vv+uu))
1.701 + M933h = (1.- (bias/dimsizes(uu)))*5.
1.702 + print (M933h)
1.703 +
1.704 + delete (u)
1.705 + delete (v)
1.706 + delete (uu)
1.707 + delete (vv)
1.708 +;----------------------------------------------------------------------
1.709 +; histogram res
1.710 +
1.711 + delete (resh)
1.712 + resh = True
1.713 + resh@gsnMaximize = True
1.714 + resh@gsnDraw = False
1.715 + resh@gsnFrame = False
1.716 + resh@xyMarkLineMode = "Markers"
1.717 + resh@xyMarkerSizeF = 0.014
1.718 + resh@xyMarker = 16
1.719 + resh@xyMarkerColors = (/"Brown","Blue"/)
1.720 + resh@trYMinF = min(mn_yvalues) - 10.
1.721 + resh@trYMaxF = max(mx_yvalues) + 10.
1.722 +
1.723 + resh@tiYAxisString = "NPP (g C/m2/year)"
1.724 + resh@tiXAxisString = "Precipitation (m/year)"
1.725 +
1.726 + max_bar = new((/2,nx/),graphic)
1.727 + min_bar = new((/2,nx/),graphic)
1.728 + max_cap = new((/2,nx/),graphic)
1.729 + min_cap = new((/2,nx/),graphic)
1.730 +
1.731 + lnres = True
1.732 + line_colors = (/"brown","blue"/)
1.733 +;=================================================================
1.734 +; histogram ob 933 site only
1.735 +
1.736 + plot_name = "histogram_ob_933"
1.737 + title = "Observed 933 site"
1.738 + resh@tiMainString = title
1.739 +
1.740 + wks = gsn_open_wks (plot_type,plot_name)
1.741 +
1.742 + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
1.743 +
1.744 +;-------------------------------
1.745 +;Attach the vertical bar and the horizontal cap line
1.746 +
1.747 + do nd=0,0
1.748 + lnres@gsLineColor = line_colors(nd)
1.749 + do i=0,nx-1
1.750 +
1.751 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.752 + .not.ismissing(mx_yvalues(nd,i))) then
1.753 +
1.754 +; Attach the vertical bar, both above and below the marker.
1.755 +
1.756 + x1 = xvalues(nd,i)
1.757 + y1 = yvalues(nd,i)
1.758 + y2 = mn_yvalues(nd,i)
1.759 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.760 +
1.761 + y2 = mx_yvalues(nd,i)
1.762 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.763 +
1.764 +; Attach the horizontal cap line, both above and below the marker.
1.765 +
1.766 + x1 = xvalues(nd,i) - dx4
1.767 + x2 = xvalues(nd,i) + dx4
1.768 + y1 = mn_yvalues(nd,i)
1.769 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.770 +
1.771 + y1 = mx_yvalues(nd,i)
1.772 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.773 + end if
1.774 + end do
1.775 + end do
1.776 +
1.777 + draw(xy)
1.778 + frame(wks)
1.779 +
1.780 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.781 + system("rm "+plot_name+"."+plot_type)
1.782 +; system("rm "+plot_name+"-1."+plot_type_new)
1.783 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.784 +
1.785 + delete (xy)
1.786 + clear (wks)
1.787 +
1.788 +;===========================================================================
1.789 +; histogram model vs ob 933 site
1.790 +
1.791 + plot_name = "histogram_mod-ob_933"
1.792 + title = model_name+ " vs Observed 933 site"
1.793 + resh@tiMainString = title
1.794 +
1.795 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.796 +
1.797 +;-----------------------------
1.798 +; Add a boxed legend using the more simple method, which won't have
1.799 +; vertical lines going through the markers.
1.800 +
1.801 + resh@pmLegendDisplayMode = "Always"
1.802 +; resh@pmLegendWidthF = 0.1
1.803 + resh@pmLegendWidthF = 0.08
1.804 + resh@pmLegendHeightF = 0.05
1.805 + resh@pmLegendOrthogonalPosF = -1.17
1.806 +; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.807 +; resh@pmLegendParallelPosF = 0.18
1.808 + resh@pmLegendParallelPosF = 0.88 ;(rightward)
1.809 +
1.810 +; resh@lgPerimOn = False
1.811 + resh@lgLabelFontHeightF = 0.015
1.812 + resh@xyExplicitLegendLabels = (/"observed",model_name/)
1.813 +;-----------------------------
1.814 + tRes = True
1.815 + tRes@txFontHeightF = 0.025
1.816 +
1.817 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"
1.818 +
1.819 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.820 +
1.821 + xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
1.822 +;-------------------------------
1.823 +;Attach the vertical bar and the horizontal cap line
1.824 +
1.825 + do nd=0,1
1.826 + lnres@gsLineColor = line_colors(nd)
1.827 + do i=0,nx-1
1.828 +
1.829 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.830 + .not.ismissing(mx_yvalues(nd,i))) then
1.831 +;
1.832 +; Attach the vertical bar, both above and below the marker.
1.833 +;
1.834 + x1 = xvalues(nd,i)
1.835 + y1 = yvalues(nd,i)
1.836 + y2 = mn_yvalues(nd,i)
1.837 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.838 +
1.839 + y2 = mx_yvalues(nd,i)
1.840 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.841 +;
1.842 +; Attach the horizontal cap line, both above and below the marker.
1.843 +;
1.844 + x1 = xvalues(nd,i) - dx4
1.845 + x2 = xvalues(nd,i) + dx4
1.846 + y1 = mn_yvalues(nd,i)
1.847 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.848 +
1.849 + y1 = mx_yvalues(nd,i)
1.850 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.851 + end if
1.852 + end do
1.853 + end do
1.854 + draw(xy)
1.855 + frame(wks)
1.856 +
1.857 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.858 + system("rm "+plot_name+"."+plot_type)
1.859 +; system("rm "+plot_name+"-1."+plot_type_new)
1.860 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.861 +
1.862 + delete (xy)
1.863 + clear (wks)
1.864 +;------------------------------------------------------------------------
1.865 +; global contour
1.866 +
1.867 +;res
1.868 + resg = True ; Use plot options
1.869 + resg@cnFillOn = True ; Turn on color fill
1.870 + resg@gsnSpreadColors = True ; use full colormap
1.871 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.872 +; resg@lbLabelAutoStride = True
1.873 + resg@cnLinesOn = False ; Turn off contourn lines
1.874 + resg@mpFillOn = False ; Turn off map fill
1.875 +
1.876 + resg@gsnSpreadColors = True ; use full colormap
1.877 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.878 + resg@cnMinLevelValF = 0. ; Min level
1.879 + resg@cnMaxLevelValF = 2200. ; Max level
1.880 + resg@cnLevelSpacingF = 200. ; interval
1.881 +;------------------------------------------------------------------------
1.882 +;global contour ob
1.883 +
1.884 + delta = 0.00000000001
1.885 + nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)
1.886 +
1.887 + plot_name = "global_ob"
1.888 + title = "Observed MODIS MOD 17"
1.889 + resg@tiMainString = title
1.890 +
1.891 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.892 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.893 +
1.894 + plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1.895 + frame(wks)
1.896 +
1.897 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.898 + system("rm "+plot_name+"."+plot_type)
1.899 + system("rm "+plot_name+"-1."+plot_type_new)
1.900 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.901 +
1.902 + clear (wks)
1.903 +;------------------------------------------------------------------------
1.904 +;global contour model
1.905 +
1.906 + plot_name = "global_model"
1.907 + title = "Model "+ model_name
1.908 + resg@tiMainString = title
1.909 +
1.910 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.911 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.912 +
1.913 + plot = gsn_csm_contour_map_ce(wks,nppmod,resg)
1.914 + frame(wks)
1.915 +
1.916 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.917 + system("rm "+plot_name+"."+plot_type)
1.918 + system("rm "+plot_name+"-1."+plot_type_new)
1.919 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.920 +
1.921 + clear (wks)
1.922 +;------------------------------------------------------------------------
1.923 +;global contour model vs ob
1.924 +
1.925 + plot_name = "global_model_vs_ob"
1.926 +
1.927 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.928 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.929 +
1.930 + delete (plot)
1.931 + plot=new(3,graphic) ; create graphic array
1.932 +
1.933 + resg@gsnFrame = False ; Do not draw plot
1.934 + resg@gsnDraw = False ; Do not advance frame
1.935 +
1.936 +;(d) compute correlation coef and M score
1.937 +
1.938 + uu1 = ndtooned(nppmod)
1.939 + vv1 = ndtooned(nppglobe)
1.940 +
1.941 + delete (good)
1.942 + good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))
1.943 +
1.944 + ug = uu1(good)
1.945 + vg = vv1(good)
1.946 +
1.947 + ccrG = esccr(ug,vg,0)
1.948 + print (ccrG)
1.949 +
1.950 + MG = (ccrG*ccrG)* 5.0
1.951 + print (MG)
1.952 +
1.953 +; plot correlation coef
1.954 +
1.955 + gRes = True
1.956 + gRes@txFontHeightF = 0.02
1.957 + gRes@txAngleF = 90
1.958 +
1.959 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
1.960 +
1.961 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.962 +;--------------------------------------------------------------------
1.963 +
1.964 +;(a) ob
1.965 +
1.966 + title = "Observed MODIS MOD 17"
1.967 + resg@tiMainString = title
1.968 +
1.969 + plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1.970 +
1.971 +;(b) model
1.972 +
1.973 + title = "Model "+ model_name
1.974 + resg@tiMainString = title
1.975 +
1.976 + plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)
1.977 +
1.978 +;(c) model-ob
1.979 +
1.980 + zz = nppmod
1.981 + zz = nppmod - nppglobe
1.982 + title = "Model_"+model_name+" - Observed"
1.983 +
1.984 + resg@cnMinLevelValF = -500 ; Min level
1.985 + resg@cnMaxLevelValF = 500. ; Max level
1.986 + resg@cnLevelSpacingF = 50. ; interval
1.987 + resg@tiMainString = title
1.988 +
1.989 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.990 +
1.991 + pres = True ; panel plot mods desired
1.992 + pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.993 + ; indiv. plots in panel
1.994 + pres@gsnMaximize = True ; fill the page
1.995 +
1.996 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.997 +
1.998 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.999 +; system("rm "+plot_name+"."+plot_type)
1.1000 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1001 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1002 +
1.1003 + frame (wks)
1.1004 + clear (wks)
1.1005 +
1.1006 + delete (plot)
1.1007 +;---------------------------------------------------------------------
1.1008 +; zonal line plot ob
1.1009 +
1.1010 + resz = True
1.1011 +
1.1012 + vv = zonalAve(nppglobe)
1.1013 + vv@long_name = nppglobe@long_name
1.1014 +
1.1015 + plot_name = "zonal_ob"
1.1016 + title = "Observed MODIS MOD 17"
1.1017 + resz@tiMainString = title
1.1018 +
1.1019 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1020 +
1.1021 + plot = gsn_csm_xy (wks,ym,vv,resz)
1.1022 + frame(wks)
1.1023 +
1.1024 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1025 + system("rm "+plot_name+"."+plot_type)
1.1026 + system("rm "+plot_name+"-1."+plot_type_new)
1.1027 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1028 +
1.1029 + clear (wks)
1.1030 +;---------------------------------------------------------------------
1.1031 +; zonal line plot model vs ob
1.1032 +
1.1033 + s = new ((/2,dimsizes(ym)/), float)
1.1034 +
1.1035 + s(0,:) = zonalAve(nppglobe)
1.1036 + s(1,:) = zonalAve(nppmod)
1.1037 +
1.1038 + s@long_name = nppglobe@long_name
1.1039 +;-------------------------------------------
1.1040 +;(d) compute correlation coef and M score
1.1041 +
1.1042 + ccrZ = esccr(s(1,:), s(0,:),0)
1.1043 + print (ccrZ)
1.1044 +
1.1045 + MZ = (ccrZ*ccrZ)* 5.0
1.1046 + print (MZ)
1.1047 +;-------------------------------------------
1.1048 + plot_name = "zonal_model_vs_ob"
1.1049 + title = "Zonal Average"
1.1050 + resz@tiMainString = title
1.1051 +
1.1052 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1053 +
1.1054 +; resz@vpHeightF = 0.4 ; change aspect ratio of plot
1.1055 +; resz@vpWidthF = 0.7
1.1056 +
1.1057 + resz@xyMonoLineColor = "False" ; want colored lines
1.1058 + resz@xyLineColors = (/"black","red"/) ; colors chosen
1.1059 +; resz@xyLineThicknesses = (/3.,3./) ; line thicknesses
1.1060 + resz@xyLineThicknesses = (/2.,2./) ; line thicknesses
1.1061 + resz@xyDashPatterns = (/0.,0./) ; make all lines solid
1.1062 +
1.1063 + resz@tiYAxisString = s@long_name ; add a axis title
1.1064 + resz@txFontHeightF = 0.0195 ; change title font heights
1.1065 +
1.1066 +; Legent
1.1067 + resz@pmLegendDisplayMode = "Always" ; turn on legend
1.1068 + resz@pmLegendSide = "Top" ; Change location of
1.1069 +; resz@pmLegendParallelPosF = .45 ; move units right
1.1070 + resz@pmLegendParallelPosF = .82 ; move units right
1.1071 + resz@pmLegendOrthogonalPosF = -0.4 ; move units down
1.1072 +
1.1073 + resz@pmLegendWidthF = 0.10 ; Change width and
1.1074 + resz@pmLegendHeightF = 0.10 ; height of legend.
1.1075 + resz@lgLabelFontHeightF = .02 ; change font height
1.1076 +; resz@lgTitleOn = True ; turn on legend title
1.1077 +; resz@lgTitleString = "Example" ; create legend title
1.1078 +; resz@lgTitleFontHeightF = .025 ; font of legend title
1.1079 + resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels
1.1080 +;--------------------------------------------------------------------
1.1081 + zRes = True
1.1082 + zRes@txFontHeightF = 0.025
1.1083 +
1.1084 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"
1.1085 +
1.1086 + gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)
1.1087 +;--------------------------------------------------------------------
1.1088 +
1.1089 + plot = gsn_csm_xy (wks,ym,s,resz) ; create plot
1.1090 +
1.1091 + frame(wks) ; advance frame
1.1092 +
1.1093 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1094 + system("rm "+plot_name+"."+plot_type)
1.1095 + system("rm "+plot_name+"-1."+plot_type_new)
1.1096 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1097 +
1.1098 + clear (wks)
1.1099 +end