1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lai/99.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,1185 @@
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 +; Main code.
1.20 +begin
1.21 +
1.22 +;nclass = 18
1.23 + nclass = 20
1.24 +
1.25 + plot_type = "ps"
1.26 + plot_type_new = "png"
1.27 +
1.28 +;************************************************
1.29 +; read in data: model
1.30 +;************************************************
1.31 +
1.32 + model_name = "b30.061n"
1.33 + model_grid = "T31"
1.34 +
1.35 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.36 + film = "b30.061n_1995-2004_MONS_climo_lnd.nc"
1.37 +;film = "i01.03cn_1545-1569_MONS_climo.nc"
1.38 + fm = addfile(dirm+film,"r")
1.39 +
1.40 + laimod = fm->TLAI
1.41 +
1.42 +;************************************************
1.43 +; read in data: observed
1.44 +;************************************************
1.45 +
1.46 + ob_name = "MODIS MOD 15A2 2000-2005"
1.47 +
1.48 + diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
1.49 + filo1 = "land_class_"+model_grid+".nc"
1.50 + filo2 = "LAI_2000-2005_MONS_"+model_grid+".nc"
1.51 +
1.52 + fo1 = addfile(diro+filo1,"r")
1.53 + fo2 = addfile(diro+filo2,"r")
1.54 +
1.55 + classob = tofloat(fo1->LAND_CLASS)
1.56 + laiob = fo2->LAI
1.57 +;*******************************************************************
1.58 +; Calculate "nice" bins for binning the data in equally spaced ranges
1.59 +;********************************************************************
1.60 + nclassn = nclass + 1
1.61 + range = fspan(0,nclassn-1,nclassn)
1.62 +; print (range)
1.63 +
1.64 +; Use this range information to grab all the values in a
1.65 +; particular range, and then take an average.
1.66 +
1.67 + nr = dimsizes(range)
1.68 + nx = nr-1
1.69 + xvalues = new((/2,nx/),float)
1.70 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.71 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.72 + dx4 = dx/4 ; 1/4 of the range
1.73 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.74 +;-----------------------------------------------------------------
1.75 +
1.76 +;-----------------------------------------------------------------
1.77 +;(B) max
1.78 +;--------------------------------------------------------------------
1.79 +; get data
1.80 +
1.81 +; observed
1.82 + laiob_max = laiob(0,:,:)
1.83 + s = laiob(:,0,0)
1.84 + laiob_max@long_name = "Leaf Area Index Max"
1.85 +
1.86 + dsizes_z = dimsizes(laiob)
1.87 + nlat = dsizes_z(1)
1.88 + nlon = dsizes_z(2)
1.89 +
1.90 + do j = 0,nlat-1
1.91 + do i = 0,nlon-1
1.92 + s = laiob(:,j,i)
1.93 + laiob_max(j,i) = max(s)
1.94 + end do
1.95 + end do
1.96 +
1.97 +; print (min(y)+"/"+max(y))
1.98 + delete (s)
1.99 + delete (dsizes_z)
1.100 +;-------------------------
1.101 +; model
1.102 + laimod_max = laimod(0,:,:)
1.103 + s = laimod(:,0,0)
1.104 + laimod_max@long_name = "Leaf Area Index Max"
1.105 +
1.106 + dsizes_z = dimsizes(laimod)
1.107 + nlat = dsizes_z(1)
1.108 + nlon = dsizes_z(2)
1.109 +
1.110 + do j = 0,nlat-1
1.111 + do i = 0,nlon-1
1.112 + s = laimod(:,j,i)
1.113 + laimod_max(j,i) = max(s)
1.114 + end do
1.115 + end do
1.116 +
1.117 +; print (min(laimod_max)+"/"+max(laimod_max))
1.118 + delete (s)
1.119 + delete (dsizes_z)
1.120 +;------------------------
1.121 + DATA11_1D = ndtooned(classob)
1.122 + DATA12_1D = ndtooned(laiob_max)
1.123 + DATA22_1D = ndtooned(laimod_max)
1.124 +
1.125 + yvalues = new((/2,nx/),float)
1.126 + mn_yvalues = new((/2,nx/),float)
1.127 + mx_yvalues = new((/2,nx/),float)
1.128 +
1.129 + do nd=0,1
1.130 +
1.131 +; See if we are doing model or observational data.
1.132 +
1.133 + if(nd.eq.0) then
1.134 + data_ob = DATA11_1D
1.135 + data_mod = DATA12_1D
1.136 + else
1.137 + data_ob = DATA11_1D
1.138 + data_mod = DATA22_1D
1.139 + end if
1.140 +
1.141 +; Loop through each range and check for values.
1.142 +
1.143 + do i=0,nr-2
1.144 + if (i.ne.(nr-2)) then
1.145 +; print("")
1.146 +; print("In range ["+range(i)+","+range(i+1)+")")
1.147 + idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))
1.148 + else
1.149 +; print("")
1.150 +; print("In range ["+range(i)+",)")
1.151 + idx = ind(range(i).le.data_ob)
1.152 + end if
1.153 +
1.154 +; Calculate average, and get min and max.
1.155 +
1.156 + if(.not.any(ismissing(idx))) then
1.157 + yvalues(nd,i) = avg(data_mod(idx))
1.158 + mn_yvalues(nd,i) = min(data_mod(idx))
1.159 + mx_yvalues(nd,i) = max(data_mod(idx))
1.160 + count = dimsizes(idx)
1.161 + else
1.162 + count = 0
1.163 + yvalues(nd,i) = yvalues@_FillValue
1.164 + mn_yvalues(nd,i) = yvalues@_FillValue
1.165 + mx_yvalues(nd,i) = yvalues@_FillValue
1.166 + end if
1.167 +
1.168 +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1.169 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.170 +
1.171 +; Clean up for next time in loop.
1.172 +
1.173 + delete(idx)
1.174 + end do
1.175 + delete(data_ob)
1.176 + delete(data_mod)
1.177 + end do
1.178 +;-----------------------------------------------------------------
1.179 +; compute correlation coef and M score
1.180 +
1.181 + u = yvalues(0,:)
1.182 + v = yvalues(1,:)
1.183 +
1.184 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.185 + uu = u(good)
1.186 + vv = v(good)
1.187 +
1.188 + ccrMax = esccr(uu,vv,0)
1.189 +; print (ccrMax)
1.190 +
1.191 +; new eq
1.192 + bias = sum(abs(vv-uu)/(vv+uu))
1.193 + Mmax = (1.- (bias/dimsizes(uu)))*5.
1.194 +
1.195 + print (Mmax)
1.196 +
1.197 + delete (u)
1.198 + delete (v)
1.199 + delete (uu)
1.200 + delete (vv)
1.201 +;--------------------------------------------------------------------
1.202 +; histogram res
1.203 +
1.204 + resm = True
1.205 + resm@gsnMaximize = True
1.206 + resm@gsnDraw = False
1.207 + resm@gsnFrame = False
1.208 + resm@xyMarkLineMode = "Markers"
1.209 + resm@xyMarkerSizeF = 0.014
1.210 + resm@xyMarker = 16
1.211 + resm@xyMarkerColors = (/"Brown","Blue"/)
1.212 +; resm@trYMinF = min(mn_yvalues) - 10.
1.213 +; resm@trYMaxF = max(mx_yvalues) + 10.
1.214 + resm@trYMinF = min(mn_yvalues) - 2
1.215 + resm@trYMaxF = max(mx_yvalues) + 4
1.216 +
1.217 + resm@tiYAxisString = "Max LAI (Leaf Area Index)"
1.218 + resm@tiXAxisString = "Land Cover Type"
1.219 +
1.220 + max_bar = new((/2,nx/),graphic)
1.221 + min_bar = new((/2,nx/),graphic)
1.222 + max_cap = new((/2,nx/),graphic)
1.223 + min_cap = new((/2,nx/),graphic)
1.224 +
1.225 + lnres = True
1.226 + line_colors = (/"brown","blue"/)
1.227 +;------------------------------------------------------------------
1.228 +; Start the graphics.
1.229 +
1.230 + plot_name = "histogram_max"
1.231 + title = model_name + " vs Observed"
1.232 + resm@tiMainString = title
1.233 +
1.234 + wks = gsn_open_wks (plot_type,plot_name)
1.235 +;-----------------------------
1.236 +; Add a boxed legend using the more simple method
1.237 +
1.238 + resm@pmLegendDisplayMode = "Always"
1.239 +; resm@pmLegendWidthF = 0.1
1.240 + resm@pmLegendWidthF = 0.08
1.241 + resm@pmLegendHeightF = 0.05
1.242 + resm@pmLegendOrthogonalPosF = -1.17
1.243 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.244 +; resm@pmLegendParallelPosF = 0.18
1.245 + resm@pmLegendParallelPosF = 0.88 ;(rightward)
1.246 +
1.247 +; resm@lgPerimOn = False
1.248 + resm@lgLabelFontHeightF = 0.015
1.249 + resm@xyExplicitLegendLabels = (/"observed",model_name/)
1.250 +;-----------------------------
1.251 + tRes = True
1.252 + tRes@txFontHeightF = 0.025
1.253 +
1.254 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")"
1.255 +
1.256 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.257 +
1.258 + xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
1.259 +;-------------------------------
1.260 +;Attach the vertical bar and the horizontal cap line
1.261 +
1.262 + do nd=0,1
1.263 + lnres@gsLineColor = line_colors(nd)
1.264 + do i=0,nx-1
1.265 +
1.266 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.267 + .not.ismissing(mx_yvalues(nd,i))) then
1.268 +
1.269 +; Attach the vertical bar, both above and below the marker.
1.270 +
1.271 + x1 = xvalues(nd,i)
1.272 + y1 = yvalues(nd,i)
1.273 + y2 = mn_yvalues(nd,i)
1.274 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.275 +
1.276 + y2 = mx_yvalues(nd,i)
1.277 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.278 +
1.279 +; Attach the horizontal cap line, both above and below the marker.
1.280 +
1.281 + x1 = xvalues(nd,i) - dx4
1.282 + x2 = xvalues(nd,i) + dx4
1.283 + y1 = mn_yvalues(nd,i)
1.284 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.285 +
1.286 + y1 = mx_yvalues(nd,i)
1.287 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.288 + end if
1.289 + end do
1.290 + end do
1.291 +
1.292 + draw(xy)
1.293 + frame(wks)
1.294 +
1.295 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.296 +; system("rm "+plot_name+"."+plot_type)
1.297 +; system("rm "+plot_name+"-1."+plot_type_new)
1.298 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.299 +
1.300 + clear (wks)
1.301 +
1.302 + delete (DATA11_1D)
1.303 + delete (DATA12_1D)
1.304 + delete (DATA22_1D)
1.305 +;delete (range)
1.306 +;delete (xvalues)
1.307 + delete (yvalues)
1.308 + delete (mn_yvalues)
1.309 + delete (mx_yvalues)
1.310 + delete (good)
1.311 + delete (max_bar)
1.312 + delete (min_bar)
1.313 + delete (max_cap)
1.314 + delete (min_cap)
1.315 +;-----------------------------------------------------------------
1.316 +;global res
1.317 +
1.318 + resg = True ; Use plot options
1.319 + resg@cnFillOn = True ; Turn on color fill
1.320 + resg@gsnSpreadColors = True ; use full colormap
1.321 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.322 +; resg@lbLabelAutoStride = True
1.323 + resg@cnLinesOn = False ; Turn off contourn lines
1.324 + resg@mpFillOn = False ; Turn off map fill
1.325 +
1.326 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.327 + resg@cnMinLevelValF = 0. ; Min level
1.328 + resg@cnMaxLevelValF = 10. ; Max level
1.329 + resg@cnLevelSpacingF = 1. ; interval
1.330 +
1.331 +;global contour ob
1.332 +
1.333 + delta = 0.00001
1.334 + laiob_max = where(ismissing(laiob_max).and.(ismissing(laimod_max).or.(laimod_max.lt.delta)),0.,laiob_max)
1.335 +
1.336 + plot_name = "global_max_ob"
1.337 + title = ob_name
1.338 + resg@tiMainString = title
1.339 +
1.340 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.341 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.342 +
1.343 + plot = gsn_csm_contour_map_ce(wks,laiob_max,resg)
1.344 + frame(wks)
1.345 +
1.346 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.347 +; system("rm "+plot_name+"."+plot_type)
1.348 +; system("rm "+plot_name+"-1."+plot_type_new)
1.349 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.350 +
1.351 + clear (wks)
1.352 +;------------------------------------------------------------------------
1.353 +;global contour model
1.354 +
1.355 + plot_name = "global_max_model"
1.356 + title = "Model " + model_name
1.357 + resg@tiMainString = title
1.358 +
1.359 + wks = gsn_open_wks (plot_type,plot_name)
1.360 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.361 +
1.362 + delete (plot)
1.363 + plot = gsn_csm_contour_map_ce(wks,laimod_max,resg)
1.364 + frame(wks)
1.365 +
1.366 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.367 +; system("rm "+plot_name+"."+plot_type)
1.368 +; system("rm "+plot_name+"-1."+plot_type_new)
1.369 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.370 +
1.371 + clear (wks)
1.372 +;------------------------------------------------------------------------
1.373 +;global contour model vs ob
1.374 +
1.375 + plot_name = "global_max_model_vs_ob"
1.376 +
1.377 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.378 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.379 +
1.380 + delete (plot)
1.381 + plot=new(3,graphic) ; create graphic array
1.382 +
1.383 + resg@gsnFrame = False ; Do not draw plot
1.384 + resg@gsnDraw = False ; Do not advance frame
1.385 +
1.386 +; plot correlation coef
1.387 +
1.388 + gRes = True
1.389 + gRes@txFontHeightF = 0.02
1.390 + gRes@txAngleF = 90
1.391 +
1.392 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")"
1.393 +
1.394 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.395 +;--------------------------------------------------------------------
1.396 +
1.397 +;(a) ob
1.398 +
1.399 + title = ob_name
1.400 + resg@tiMainString = title
1.401 +
1.402 + plot(0) = gsn_csm_contour_map_ce(wks,laiob_max,resg)
1.403 +
1.404 +;(b) model
1.405 +
1.406 + title = "Model "+ model_name
1.407 + resg@tiMainString = title
1.408 +
1.409 + plot(1) = gsn_csm_contour_map_ce(wks,laimod_max,resg)
1.410 +
1.411 +;(c) model-ob
1.412 +
1.413 + zz = laimod_max
1.414 + zz = laimod_max - laiob_max
1.415 + title = "Model_"+model_name+" - Observed"
1.416 + resg@tiMainString = title
1.417 +
1.418 + resg@cnMinLevelValF = -6. ; Min level
1.419 + resg@cnMaxLevelValF = 6. ; Max level
1.420 + resg@cnLevelSpacingF = 1. ; interval
1.421 +
1.422 +
1.423 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.424 +
1.425 + pres = True ; panel plot mods desired
1.426 + pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.427 + ; indiv. plots in panel
1.428 + pres@gsnMaximize = True ; fill the page
1.429 +
1.430 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.431 +
1.432 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.433 +; system("rm "+plot_name+"."+plot_type)
1.434 +; system("rm "+plot_name+"-1."+plot_type_new)
1.435 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.436 +
1.437 + frame (wks)
1.438 + clear (wks)
1.439 +
1.440 + delete (plot)
1.441 +;-----------------------------------------------------------------
1.442 +;(C) phase
1.443 +;--------------------------------------------------------------------
1.444 +; get data
1.445 +
1.446 +; observed
1.447 + laiob_phase = laiob(0,:,:)
1.448 + s = laiob(:,0,0)
1.449 + laiob_phase@long_name = "Leaf Area Index Max Month"
1.450 +
1.451 + dsizes_z = dimsizes(laiob)
1.452 + nlat = dsizes_z(1)
1.453 + nlon = dsizes_z(2)
1.454 +
1.455 + do j = 0,nlat-1
1.456 + do i = 0,nlon-1
1.457 + s = laiob(:,j,i)
1.458 + laiob_phase(j,i) = maxind(s) + 1
1.459 + end do
1.460 + end do
1.461 +
1.462 +; print (min(laiob_phase)+"/"+max(laiob_phase))
1.463 + delete (s)
1.464 + delete (dsizes_z)
1.465 +;-------------------------
1.466 +; model
1.467 + laimod_phase = laimod(0,:,:)
1.468 + s = laimod(:,0,0)
1.469 + laimod_phase@long_name = "Leaf Area Index Max Month"
1.470 +
1.471 + dsizes_z = dimsizes(laimod)
1.472 + nlat = dsizes_z(1)
1.473 + nlon = dsizes_z(2)
1.474 +
1.475 + do j = 0,nlat-1
1.476 + do i = 0,nlon-1
1.477 + s = laimod(:,j,i)
1.478 + laimod_phase(j,i) = maxind(s) + 1
1.479 + end do
1.480 + end do
1.481 +
1.482 +; print (min(laimod_phase)+"/"+max(laimod_phase))
1.483 + delete (s)
1.484 + delete (dsizes_z)
1.485 +;------------------------
1.486 + DATA11_1D = ndtooned(classob)
1.487 + DATA12_1D = ndtooned(laiob_phase)
1.488 + DATA22_1D = ndtooned(laimod_phase)
1.489 +
1.490 + yvalues = new((/2,nx/),float)
1.491 + mn_yvalues = new((/2,nx/),float)
1.492 + mx_yvalues = new((/2,nx/),float)
1.493 +
1.494 + do nd=0,1
1.495 +
1.496 +; See if we are doing model or observational data.
1.497 +
1.498 + if(nd.eq.0) then
1.499 + data_ob = DATA11_1D
1.500 + data_mod = DATA12_1D
1.501 + else
1.502 + data_ob = DATA11_1D
1.503 + data_mod = DATA22_1D
1.504 + end if
1.505 +
1.506 +; Loop through each range and check for values.
1.507 +
1.508 + do i=0,nr-2
1.509 + if (i.ne.(nr-2)) then
1.510 +; print("")
1.511 +; print("In range ["+range(i)+","+range(i+1)+")")
1.512 + idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))
1.513 + else
1.514 +; print("")
1.515 +; print("In range ["+range(i)+",)")
1.516 + idx = ind(range(i).le.data_ob)
1.517 + end if
1.518 +
1.519 +; Calculate average, and get min and max.
1.520 +
1.521 + if(.not.any(ismissing(idx))) then
1.522 + yvalues(nd,i) = avg(data_mod(idx))
1.523 + mn_yvalues(nd,i) = min(data_mod(idx))
1.524 + mx_yvalues(nd,i) = max(data_mod(idx))
1.525 + count = dimsizes(idx)
1.526 + else
1.527 + count = 0
1.528 + yvalues(nd,i) = yvalues@_FillValue
1.529 + mn_yvalues(nd,i) = yvalues@_FillValue
1.530 + mx_yvalues(nd,i) = yvalues@_FillValue
1.531 + end if
1.532 +
1.533 +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1.534 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.535 +
1.536 +; Clean up for next time in loop.
1.537 +
1.538 + delete(idx)
1.539 + end do
1.540 + delete(data_ob)
1.541 + delete(data_mod)
1.542 + end do
1.543 +;-----------------------------------------------------------------
1.544 +; compute correlation coef and M score
1.545 +
1.546 + u = yvalues(0,:)
1.547 + v = yvalues(1,:)
1.548 +
1.549 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.550 + uu = u(good)
1.551 + vv = v(good)
1.552 +
1.553 + ccrPhase = esccr(uu,vv,0)
1.554 +; print (ccrPhase)
1.555 +
1.556 +; old eq
1.557 +; bias = abs(avg(vv)-avg(uu))
1.558 +; new eq
1.559 + bias = avg(abs(vv-uu))
1.560 +
1.561 + bias = where((bias.gt. 6.),12.-bias,bias)
1.562 + Mphase = ((6. - bias)/6.)*5.
1.563 +
1.564 + print (Mphase)
1.565 +
1.566 + delete (u)
1.567 + delete (v)
1.568 + delete (uu)
1.569 + delete (vv)
1.570 +;--------------------------------------------------------------------
1.571 +; histogram res
1.572 +
1.573 + resm = True
1.574 + resm@gsnMaximize = True
1.575 + resm@gsnDraw = False
1.576 + resm@gsnFrame = False
1.577 + resm@xyMarkLineMode = "Markers"
1.578 + resm@xyMarkerSizeF = 0.014
1.579 + resm@xyMarker = 16
1.580 + resm@xyMarkerColors = (/"Brown","Blue"/)
1.581 +; resm@trYMinF = min(mn_yvalues) - 10.
1.582 +; resm@trYMaxF = max(mx_yvalues) + 10.
1.583 + resm@trYMinF = min(mn_yvalues) - 2
1.584 + resm@trYMaxF = max(mx_yvalues) + 4
1.585 +
1.586 + resm@tiYAxisString = "Max LAI (Leaf Area Index) Month"
1.587 + resm@tiXAxisString = "Land Cover Type"
1.588 +
1.589 + max_bar = new((/2,nx/),graphic)
1.590 + min_bar = new((/2,nx/),graphic)
1.591 + max_cap = new((/2,nx/),graphic)
1.592 + min_cap = new((/2,nx/),graphic)
1.593 +
1.594 + lnres = True
1.595 + line_colors = (/"brown","blue"/)
1.596 +;------------------------------------------------------------------
1.597 +; Start the graphics.
1.598 +
1.599 + plot_name = "histogram_phase"
1.600 + title = model_name + " vs Observed"
1.601 + resm@tiMainString = title
1.602 +
1.603 + wks = gsn_open_wks (plot_type,plot_name)
1.604 +;-----------------------------
1.605 +; Add a boxed legend using the more simple method
1.606 +
1.607 + resm@pmLegendDisplayMode = "Always"
1.608 +; resm@pmLegendWidthF = 0.1
1.609 + resm@pmLegendWidthF = 0.08
1.610 + resm@pmLegendHeightF = 0.05
1.611 + resm@pmLegendOrthogonalPosF = -1.17
1.612 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.613 +; resm@pmLegendParallelPosF = 0.18
1.614 + resm@pmLegendParallelPosF = 0.88 ;(rightward)
1.615 +
1.616 +; resm@lgPerimOn = False
1.617 + resm@lgLabelFontHeightF = 0.015
1.618 + resm@xyExplicitLegendLabels = (/"observed",model_name/)
1.619 +;-----------------------------
1.620 + tRes = True
1.621 + tRes@txFontHeightF = 0.025
1.622 +
1.623 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")"
1.624 +
1.625 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.626 +
1.627 + xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
1.628 +;-------------------------------
1.629 +;Attach the vertical bar and the horizontal cap line
1.630 +
1.631 + do nd=0,1
1.632 + lnres@gsLineColor = line_colors(nd)
1.633 + do i=0,nx-1
1.634 +
1.635 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.636 + .not.ismissing(mx_yvalues(nd,i))) then
1.637 +
1.638 +; Attach the vertical bar, both above and below the marker.
1.639 +
1.640 + x1 = xvalues(nd,i)
1.641 + y1 = yvalues(nd,i)
1.642 + y2 = mn_yvalues(nd,i)
1.643 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.644 +
1.645 + y2 = mx_yvalues(nd,i)
1.646 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.647 +
1.648 +; Attach the horizontal cap line, both above and below the marker.
1.649 +
1.650 + x1 = xvalues(nd,i) - dx4
1.651 + x2 = xvalues(nd,i) + dx4
1.652 + y1 = mn_yvalues(nd,i)
1.653 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.654 +
1.655 + y1 = mx_yvalues(nd,i)
1.656 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.657 + end if
1.658 + end do
1.659 + end do
1.660 +
1.661 + draw(xy)
1.662 + frame(wks)
1.663 +
1.664 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.665 +; system("rm "+plot_name+"."+plot_type)
1.666 +; system("rm "+plot_name+"-1."+plot_type_new)
1.667 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.668 +
1.669 + clear (wks)
1.670 +
1.671 + delete (DATA11_1D)
1.672 + delete (DATA12_1D)
1.673 + delete (DATA22_1D)
1.674 +;delete (range)
1.675 +;delete (xvalues)
1.676 + delete (yvalues)
1.677 + delete (mn_yvalues)
1.678 + delete (mx_yvalues)
1.679 + delete (good)
1.680 + delete (max_bar)
1.681 + delete (min_bar)
1.682 + delete (max_cap)
1.683 + delete (min_cap)
1.684 +;-----------------------------------------------------------------
1.685 +;global res
1.686 +
1.687 + resg = True ; Use plot options
1.688 + resg@cnFillOn = True ; Turn on color fill
1.689 + resg@gsnSpreadColors = True ; use full colormap
1.690 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.691 +; resg@lbLabelAutoStride = True
1.692 + resg@cnLinesOn = False ; Turn off contourn lines
1.693 + resg@mpFillOn = False ; Turn off map fill
1.694 +
1.695 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.696 + resg@cnMinLevelValF = 1. ; Min level
1.697 + resg@cnMaxLevelValF = 12. ; Max level
1.698 + resg@cnLevelSpacingF = 1. ; interval
1.699 +
1.700 +;global contour ob
1.701 +
1.702 + delta = 0.00001
1.703 + laiob_phase = where(ismissing(laiob_phase).and.(ismissing(laimod_phase).or.(laimod_phase.lt.delta)),0.,laiob_phase)
1.704 +
1.705 + plot_name = "global_phase_ob"
1.706 + title = ob_name
1.707 + resg@tiMainString = title
1.708 +
1.709 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.710 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.711 +
1.712 + plot = gsn_csm_contour_map_ce(wks,laiob_phase,resg)
1.713 + frame(wks)
1.714 +
1.715 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.716 +; system("rm "+plot_name+"."+plot_type)
1.717 +; system("rm "+plot_name+"-1."+plot_type_new)
1.718 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.719 +
1.720 + clear (wks)
1.721 +;------------------------------------------------------------------------
1.722 +;global contour model
1.723 +
1.724 + plot_name = "global_phase_model"
1.725 + title = "Model " + model_name
1.726 + resg@tiMainString = title
1.727 +
1.728 + wks = gsn_open_wks (plot_type,plot_name)
1.729 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.730 +
1.731 + delete (plot)
1.732 + plot = gsn_csm_contour_map_ce(wks,laimod_phase,resg)
1.733 + frame(wks)
1.734 +
1.735 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.736 +; system("rm "+plot_name+"."+plot_type)
1.737 +; system("rm "+plot_name+"-1."+plot_type_new)
1.738 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.739 +
1.740 + clear (wks)
1.741 +;------------------------------------------------------------------------
1.742 +;global contour model vs ob
1.743 +
1.744 + plot_name = "global_phase_model_vs_ob"
1.745 +
1.746 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.747 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.748 +
1.749 + delete (plot)
1.750 + plot=new(3,graphic) ; create graphic array
1.751 +
1.752 + resg@gsnFrame = False ; Do not draw plot
1.753 + resg@gsnDraw = False ; Do not advance frame
1.754 +
1.755 +; plot correlation coef
1.756 +
1.757 + gRes = True
1.758 + gRes@txFontHeightF = 0.02
1.759 + gRes@txAngleF = 90
1.760 +
1.761 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")"
1.762 +
1.763 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.764 +;--------------------------------------------------------------------
1.765 +
1.766 +;(a) ob
1.767 +
1.768 + title = ob_name
1.769 + resg@tiMainString = title
1.770 +
1.771 + plot(0) = gsn_csm_contour_map_ce(wks,laiob_phase,resg)
1.772 +
1.773 +;(b) model
1.774 +
1.775 + title = "Model "+ model_name
1.776 + resg@tiMainString = title
1.777 +
1.778 + plot(1) = gsn_csm_contour_map_ce(wks,laimod_phase,resg)
1.779 +
1.780 +;(c) model-ob
1.781 +
1.782 + delete (zz)
1.783 + zz = laimod_phase
1.784 + zz = laimod_phase - laiob_phase
1.785 + title = "Model_"+model_name+" - Observed"
1.786 + resg@tiMainString = title
1.787 +
1.788 + resg@cnMinLevelValF = -6. ; Min level
1.789 + resg@cnMaxLevelValF = 6. ; Max level
1.790 + resg@cnLevelSpacingF = 1. ; interval
1.791 +
1.792 +
1.793 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.794 +
1.795 +; pres = True ; panel plot mods desired
1.796 +; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.797 + ; indiv. plots in panel
1.798 +; pres@gsnMaximize = True ; fill the page
1.799 +
1.800 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.801 +
1.802 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.803 +; system("rm "+plot_name+"."+plot_type)
1.804 +; system("rm "+plot_name+"-1."+plot_type_new)
1.805 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.806 +
1.807 + frame (wks)
1.808 + clear (wks)
1.809 +
1.810 + delete (plot)
1.811 +;-----------------------------------------------------------------
1.812 +;(D) grow day
1.813 +;--------------------------------------------------------------------
1.814 +; get data
1.815 +
1.816 + day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/)
1.817 +
1.818 +; observed
1.819 + laiob_grow = laiob(0,:,:)
1.820 + laiob_grow@long_name = "Days of Growing Season"
1.821 +
1.822 + dsizes_z = dimsizes(laiob)
1.823 + ntime = dsizes_z(0)
1.824 + nlat = dsizes_z(1)
1.825 + nlon = dsizes_z(2)
1.826 +
1.827 + do j = 0,nlat-1
1.828 + do i = 0,nlon-1
1.829 + nday = 0.
1.830 + do k = 0,ntime-1
1.831 + if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then
1.832 + nday = nday + day_of_data(k)
1.833 + end if
1.834 + end do
1.835 +
1.836 + laiob_grow(j,i) = nday
1.837 + end do
1.838 + end do
1.839 +
1.840 +; print (min(laiob_grow)+"/"+max(laiob_grow))
1.841 +;-------------------------
1.842 +; model
1.843 + laimod_grow = laimod(0,:,:)
1.844 + laimod_grow@long_name = "Days of Growing Season"
1.845 +
1.846 + dsizes_z = dimsizes(laimod)
1.847 + ntime = dsizes_z(0)
1.848 + nlat = dsizes_z(1)
1.849 + nlon = dsizes_z(2)
1.850 +
1.851 + do j = 0,nlat-1
1.852 + do i = 0,nlon-1
1.853 + nday = 0.
1.854 + do k = 0,ntime-1
1.855 + if (.not. ismissing(laimod(k,j,i)) .and. laimod(k,j,i) .gt. 1.0) then
1.856 + nday = nday + day_of_data(k)
1.857 + end if
1.858 + end do
1.859 +
1.860 + laimod_grow(j,i) = nday
1.861 + end do
1.862 + end do
1.863 +
1.864 +; print (min(laimod_grow)+"/"+max(laimod_grow))
1.865 +;------------------------
1.866 + DATA11_1D = ndtooned(classob)
1.867 + DATA12_1D = ndtooned(laiob_grow)
1.868 + DATA22_1D = ndtooned(laimod_grow)
1.869 +
1.870 + yvalues = new((/2,nx/),float)
1.871 + mn_yvalues = new((/2,nx/),float)
1.872 + mx_yvalues = new((/2,nx/),float)
1.873 +
1.874 + do nd=0,1
1.875 +
1.876 +; See if we are doing model or observational data.
1.877 +
1.878 + if(nd.eq.0) then
1.879 + data_ob = DATA11_1D
1.880 + data_mod = DATA12_1D
1.881 + else
1.882 + data_ob = DATA11_1D
1.883 + data_mod = DATA22_1D
1.884 + end if
1.885 +
1.886 +; Loop through each range and check for values.
1.887 +
1.888 + do i=0,nr-2
1.889 + if (i.ne.(nr-2)) then
1.890 +; print("")
1.891 +; print("In range ["+range(i)+","+range(i+1)+")")
1.892 + idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))
1.893 + else
1.894 +; print("")
1.895 +; print("In range ["+range(i)+",)")
1.896 + idx = ind(range(i).le.data_ob)
1.897 + end if
1.898 +
1.899 +; Calculate average, and get min and max.
1.900 +
1.901 + if(.not.any(ismissing(idx))) then
1.902 + yvalues(nd,i) = avg(data_mod(idx))
1.903 + mn_yvalues(nd,i) = min(data_mod(idx))
1.904 + mx_yvalues(nd,i) = max(data_mod(idx))
1.905 + count = dimsizes(idx)
1.906 + else
1.907 + count = 0
1.908 + yvalues(nd,i) = yvalues@_FillValue
1.909 + mn_yvalues(nd,i) = yvalues@_FillValue
1.910 + mx_yvalues(nd,i) = yvalues@_FillValue
1.911 + end if
1.912 +
1.913 +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1.914 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.915 +
1.916 +; Clean up for next time in loop.
1.917 +
1.918 + delete(idx)
1.919 + end do
1.920 + delete(data_ob)
1.921 + delete(data_mod)
1.922 + end do
1.923 +;-----------------------------------------------------------------
1.924 +; compute correlation coef and M score
1.925 +
1.926 + u = yvalues(0,:)
1.927 + v = yvalues(1,:)
1.928 +
1.929 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.930 + uu = u(good)
1.931 + vv = v(good)
1.932 +
1.933 + ccrGrow = esccr(uu,vv,0)
1.934 +; print (ccrGrow)
1.935 +
1.936 +; new eq
1.937 + bias = sum(abs(vv-uu)/(vv+uu))
1.938 + Mgrow = (1.- (bias/dimsizes(uu)))*5.
1.939 +
1.940 + print (Mgrow)
1.941 +
1.942 + delete (u)
1.943 + delete (v)
1.944 + delete (uu)
1.945 + delete (vv)
1.946 +;--------------------------------------------------------------------
1.947 +; histogram res
1.948 +
1.949 + resm = True
1.950 + resm@gsnMaximize = True
1.951 + resm@gsnDraw = False
1.952 + resm@gsnFrame = False
1.953 + resm@xyMarkLineMode = "Markers"
1.954 + resm@xyMarkerSizeF = 0.014
1.955 + resm@xyMarker = 16
1.956 + resm@xyMarkerColors = (/"Brown","Blue"/)
1.957 +; resm@trYMinF = min(mn_yvalues) - 10.
1.958 +; resm@trYMaxF = max(mx_yvalues) + 10.
1.959 + resm@trYMinF = min(mn_yvalues) - 2
1.960 + resm@trYMaxF = max(mx_yvalues) + 4
1.961 +
1.962 + resm@tiYAxisString = "Days of Growing season"
1.963 + resm@tiXAxisString = "Land Cover Type"
1.964 +
1.965 + max_bar = new((/2,nx/),graphic)
1.966 + min_bar = new((/2,nx/),graphic)
1.967 + max_cap = new((/2,nx/),graphic)
1.968 + min_cap = new((/2,nx/),graphic)
1.969 +
1.970 + lnres = True
1.971 + line_colors = (/"brown","blue"/)
1.972 +;------------------------------------------------------------------
1.973 +; Start the graphics.
1.974 +
1.975 + plot_name = "histogram_grow"
1.976 + title = model_name + " vs Observed"
1.977 + resm@tiMainString = title
1.978 +
1.979 + wks = gsn_open_wks (plot_type,plot_name)
1.980 +;-----------------------------
1.981 +; Add a boxed legend using the more simple method
1.982 +
1.983 + resm@pmLegendDisplayMode = "Always"
1.984 +; resm@pmLegendWidthF = 0.1
1.985 + resm@pmLegendWidthF = 0.08
1.986 + resm@pmLegendHeightF = 0.05
1.987 + resm@pmLegendOrthogonalPosF = -1.17
1.988 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.989 +; resm@pmLegendParallelPosF = 0.18
1.990 + resm@pmLegendParallelPosF = 0.88 ;(rightward)
1.991 +
1.992 +; resm@lgPerimOn = False
1.993 + resm@lgLabelFontHeightF = 0.015
1.994 + resm@xyExplicitLegendLabels = (/"observed",model_name/)
1.995 +;-----------------------------
1.996 + tRes = True
1.997 + tRes@txFontHeightF = 0.025
1.998 +
1.999 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")"
1.1000 +
1.1001 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.1002 +
1.1003 + xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
1.1004 +;-------------------------------
1.1005 +;Attach the vertical bar and the horizontal cap line
1.1006 +
1.1007 + do nd=0,1
1.1008 + lnres@gsLineColor = line_colors(nd)
1.1009 + do i=0,nx-1
1.1010 +
1.1011 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1012 + .not.ismissing(mx_yvalues(nd,i))) then
1.1013 +
1.1014 +; Attach the vertical bar, both above and below the marker.
1.1015 +
1.1016 + x1 = xvalues(nd,i)
1.1017 + y1 = yvalues(nd,i)
1.1018 + y2 = mn_yvalues(nd,i)
1.1019 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1020 +
1.1021 + y2 = mx_yvalues(nd,i)
1.1022 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1023 +
1.1024 +; Attach the horizontal cap line, both above and below the marker.
1.1025 +
1.1026 + x1 = xvalues(nd,i) - dx4
1.1027 + x2 = xvalues(nd,i) + dx4
1.1028 + y1 = mn_yvalues(nd,i)
1.1029 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1030 +
1.1031 + y1 = mx_yvalues(nd,i)
1.1032 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1033 + end if
1.1034 + end do
1.1035 + end do
1.1036 +
1.1037 + draw(xy)
1.1038 + frame(wks)
1.1039 +
1.1040 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1041 +; system("rm "+plot_name+"."+plot_type)
1.1042 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1043 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1044 +
1.1045 + clear (wks)
1.1046 +
1.1047 + delete (DATA11_1D)
1.1048 + delete (DATA12_1D)
1.1049 + delete (DATA22_1D)
1.1050 +;delete (range)
1.1051 +;delete (xvalues)
1.1052 + delete (yvalues)
1.1053 + delete (mn_yvalues)
1.1054 + delete (mx_yvalues)
1.1055 + delete (good)
1.1056 + delete (max_bar)
1.1057 + delete (min_bar)
1.1058 + delete (max_cap)
1.1059 + delete (min_cap)
1.1060 +;-----------------------------------------------------------------
1.1061 +;global res
1.1062 +
1.1063 + resg = True ; Use plot options
1.1064 + resg@cnFillOn = True ; Turn on color fill
1.1065 + resg@gsnSpreadColors = True ; use full colormap
1.1066 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.1067 +; resg@lbLabelAutoStride = True
1.1068 + resg@cnLinesOn = False ; Turn off contourn lines
1.1069 + resg@mpFillOn = False ; Turn off map fill
1.1070 +
1.1071 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.1072 + resg@cnMinLevelValF = 0. ; Min level
1.1073 + resg@cnMaxLevelValF = 360. ; Max level
1.1074 + resg@cnLevelSpacingF = 30. ; interval
1.1075 +
1.1076 +;global contour ob
1.1077 +
1.1078 + delta = 0.00001
1.1079 + laiob_grow = where(ismissing(laiob_grow).and.(ismissing(laimod_grow).or.(laimod_grow.lt.delta)),0.,laiob_grow)
1.1080 +
1.1081 + plot_name = "global_grow_ob"
1.1082 + title = ob_name
1.1083 + resg@tiMainString = title
1.1084 +
1.1085 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1086 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1087 +
1.1088 + plot = gsn_csm_contour_map_ce(wks,laiob_grow,resg)
1.1089 + frame(wks)
1.1090 +
1.1091 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1092 +; system("rm "+plot_name+"."+plot_type)
1.1093 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1094 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1095 +
1.1096 + clear (wks)
1.1097 +;------------------------------------------------------------------------
1.1098 +;global contour model
1.1099 +
1.1100 + plot_name = "global_grow_model"
1.1101 + title = "Model " + model_name
1.1102 + resg@tiMainString = title
1.1103 +
1.1104 + wks = gsn_open_wks (plot_type,plot_name)
1.1105 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1106 +
1.1107 + delete (plot)
1.1108 + plot = gsn_csm_contour_map_ce(wks,laimod_grow,resg)
1.1109 + frame(wks)
1.1110 +
1.1111 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1112 +; system("rm "+plot_name+"."+plot_type)
1.1113 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1114 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1115 +
1.1116 + clear (wks)
1.1117 +;------------------------------------------------------------------------
1.1118 +;global contour model vs ob
1.1119 +
1.1120 + plot_name = "global_grow_model_vs_ob"
1.1121 +
1.1122 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1123 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1124 +
1.1125 + delete (plot)
1.1126 + plot=new(3,graphic) ; create graphic array
1.1127 +
1.1128 + resg@gsnFrame = False ; Do not draw plot
1.1129 + resg@gsnDraw = False ; Do not advance frame
1.1130 +
1.1131 +; plot correlation coef
1.1132 +
1.1133 + gRes = True
1.1134 + gRes@txFontHeightF = 0.02
1.1135 + gRes@txAngleF = 90
1.1136 +
1.1137 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")"
1.1138 +
1.1139 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.1140 +;--------------------------------------------------------------------
1.1141 +
1.1142 +;(a) ob
1.1143 +
1.1144 + title = ob_name
1.1145 + resg@tiMainString = title
1.1146 +
1.1147 + plot(0) = gsn_csm_contour_map_ce(wks,laiob_grow,resg)
1.1148 +
1.1149 +;(b) model
1.1150 +
1.1151 + title = "Model "+ model_name
1.1152 + resg@tiMainString = title
1.1153 +
1.1154 + plot(1) = gsn_csm_contour_map_ce(wks,laimod_grow,resg)
1.1155 +
1.1156 +;(c) model-ob
1.1157 +
1.1158 + delete (zz)
1.1159 + zz = laimod_grow
1.1160 + zz = laimod_grow - laiob_grow
1.1161 + title = "Model_"+model_name+" - Observed"
1.1162 + resg@tiMainString = title
1.1163 +
1.1164 + resg@cnMinLevelValF = -120. ; Min level
1.1165 + resg@cnMaxLevelValF = 120. ; Max level
1.1166 + resg@cnLevelSpacingF = 20. ; interval
1.1167 +
1.1168 +
1.1169 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.1170 +
1.1171 + pres = True ; panel plot mods desired
1.1172 + pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.1173 + ; indiv. plots in panel
1.1174 + pres@gsnMaximize = True ; fill the page
1.1175 +
1.1176 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.1177 +
1.1178 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1179 +; system("rm "+plot_name+"."+plot_type)
1.1180 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1181 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1182 +
1.1183 + frame (wks)
1.1184 + clear (wks)
1.1185 +
1.1186 + delete (plot)
1.1187 +end
1.1188 +