1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lai/99.all.ncl.0 Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,1747 @@
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.test"
1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
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 = 20
1.23 +
1.24 + plot_type = "ps"
1.25 + plot_type_new = "png"
1.26 +
1.27 +;************************************************
1.28 +; read in data: model
1.29 +;************************************************
1.30 +;film = "b30.061n_1995-2004_MONS_climo_lnd.nc"
1.31 +;model_name = "b30.061n"
1.32 +;model_grid = "T31"
1.33 +
1.34 +;film = "newcn05_ncep_1i_MONS_climo_lnd.nc"
1.35 +;model_name = "newcn"
1.36 +;model_grid = "1.9"
1.37 +
1.38 +;film = "i01.06cn_1798-2004_MONS_climo.nc"
1.39 +;model_name = "06cn"
1.40 +;model_grid = "T42"
1.41 +
1.42 +;film = "i01.06casa_1798-2004_MONS_climo.nc"
1.43 +;model_name = "06casa"
1.44 +;model_grid = "T42"
1.45 +
1.46 +;film = "i01.10cn_1948-2004_MONS_climo.nc"
1.47 +;model_name = "10cn"
1.48 +;model_grid = "T42"
1.49 +
1.50 + film = "i01.10casa_1948-2004_MONS_climo.nc"
1.51 + model_name = "10casa"
1.52 + model_grid = "T42"
1.53 +;------------------------------------------------
1.54 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.55 + fm = addfile(dirm+film,"r")
1.56 +
1.57 + laimod = fm->TLAI
1.58 +
1.59 +;************************************************
1.60 +; read in data: observed
1.61 +;************************************************
1.62 +
1.63 + ob_name = "MODIS MOD 15A2 2000-2005"
1.64 +
1.65 + diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
1.66 + filo1 = "land_class_"+model_grid+".nc"
1.67 + filo2 = "LAI_2000-2005_MONS_"+model_grid+".nc"
1.68 +
1.69 + fo1 = addfile(diro+filo1,"r")
1.70 + fo2 = addfile(diro+filo2,"r")
1.71 +
1.72 + classob = tofloat(fo1->LAND_CLASS)
1.73 + laiob = fo2->LAI
1.74 +
1.75 +;*******************************************************************
1.76 +; for plotting table
1.77 +;*******************************************************************
1.78 +; table header name
1.79 + table_header_name = "LAI"
1.80 +
1.81 +; column (not including header column)
1.82 +
1.83 + col_header1 = (/"Mean","Max","Phase","Growth"/)
1.84 + col_header2 = (/"ob","model","M" \
1.85 + ,"ob","model","M" \
1.86 + ,"ob","model","M" \
1.87 + ,"ob","model","M" \
1.88 + /)
1.89 +
1.90 + ncol1 = dimsizes(col_header1)
1.91 + ncol2 = dimsizes(col_header2)
1.92 + ncol = ncol2
1.93 +
1.94 +; row (not including header row)
1.95 + row_header = (/"Water Bodies" \
1.96 + ,"Evergreen Needleleaf Forests" \
1.97 + ,"Evergreen Broadleaf Forests" \
1.98 + ,"Deciduous Needleleaf Forest" \
1.99 + ,"Deciduous Broadleaf Forests" \
1.100 + ,"Mixed Forests" \
1.101 + ,"Closed Bushlands" \
1.102 + ,"Open Bushlands" \
1.103 + ,"Woody Savannas (S. Hem.)" \
1.104 + ,"Savannas (S. Hem.)" \
1.105 + ,"Grasslands" \
1.106 + ,"Permanent Wetlands" \
1.107 + ,"Croplands" \
1.108 + ,"Urban and Built-Up" \
1.109 + ,"Cropland/Natural Vegetation Mosaic" \
1.110 + ,"Permanent Snow and Ice" \
1.111 + ,"Barren or Sparsely Vegetated" \
1.112 + ,"Unclassified" \
1.113 + ,"Woody Savannas (N. Hem.)" \
1.114 + ,"Savannas (N. Hem.)" \
1.115 + ,"All biome average" \
1.116 + /)
1.117 + nrow = dimsizes(row_header)
1.118 +
1.119 +; arrays to be passed to table.
1.120 + value = new ((/nrow, ncol/),string )
1.121 +
1.122 + table_length = 0.995
1.123 +
1.124 +; Table header
1.125 + ncr1 = (/1,1/) ; 1 row, 1 column
1.126 + xx1 = (/0.005,0.25/) ; Start and end X
1.127 + yy1 = (/0.900,0.995/) ; Start and end Y
1.128 + text1 = table_header_name
1.129 + res1 = True
1.130 + res1@txFontHeightF = 0.03
1.131 + res1@gsFillColor = "CornFlowerBlue"
1.132 +
1.133 +; Column header (equally space in x2)
1.134 + ncr21 = (/1,ncol1/) ; 1 rows, 4 columns
1.135 + xx21 = (/xx1(1),0.995/) ; start from end of x1
1.136 + yy21 = (/0.9475,0.995/) ; half of y1
1.137 + text21 = col_header1
1.138 + res21 = True
1.139 + res21@txFontHeightF = 0.015
1.140 + res21@gsFillColor = "Gray"
1.141 +
1.142 + ncr22 = (/1,ncol2/) ; 1 rows, 12 columns
1.143 + xx22 = xx21 ; start from end of x1
1.144 + yy22 = (/0.900,0.9475/) ; half of y1
1.145 + text22 = col_header2
1.146 + res22 = True
1.147 + res22@txFontHeightF = 0.015
1.148 + res22@gsFillColor = "Gray"
1.149 +
1.150 +; Row header (equally space in y2)
1.151 + ncr3 = (/nrow,1/) ; 20 rows, 1 columns
1.152 + xx3 = xx1 ; same as x1
1.153 + yy3 = (/1.0-table_length,0.900/) ; end at start of y1
1.154 + text3 = row_header
1.155 + res3 = True
1.156 + res3@txFontHeightF = 0.01
1.157 + res3@gsFillColor = "Gray"
1.158 +
1.159 +; Main table body
1.160 + ncr4 = (/nrow,ncol/) ; 5 rows, 5 columns
1.161 + xx4 = xx21 ; Start and end x
1.162 + yy4 = yy3 ; Start and end Y
1.163 + text4 = new((/nrow,ncol/),string)
1.164 +
1.165 + color_fill4 = new((/nrow,ncol/),string)
1.166 + color_fill4 = "white"
1.167 + color_fill4(nrow-1,:) = "CornFlowerBlue"
1.168 +
1.169 + res4 = True ; Set up resource list
1.170 +; res4@gsnDebug = True ; Useful to print NDC row,col values used.
1.171 + res4@txFontHeightF = 0.015
1.172 + res4@gsFillColor = color_fill4
1.173 +
1.174 + delete (color_fill4)
1.175 +
1.176 +;************************************************
1.177 +; plot global land class: observed
1.178 +;************************************************
1.179 +;global res
1.180 +
1.181 + resg = True ; Use plot options
1.182 + resg@cnFillOn = True ; Turn on color fill
1.183 + resg@gsnSpreadColors = True ; use full colormap
1.184 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.185 +; resg@lbLabelAutoStride = True
1.186 + resg@cnLinesOn = False ; Turn off contourn lines
1.187 + resg@mpFillOn = False ; Turn off map fill
1.188 +
1.189 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.190 + resg@cnMinLevelValF = 1. ; Min level
1.191 + resg@cnMaxLevelValF = 19. ; Max level
1.192 + resg@cnLevelSpacingF = 1. ; interval
1.193 +
1.194 +;global contour ob
1.195 + classob@_FillValue = 1.e+36
1.196 + classob = where(classob.eq.0,classob@_FillValue,classob)
1.197 +
1.198 + plot_name = "global_class_ob"
1.199 + title = ob_name
1.200 + resg@tiMainString = title
1.201 +
1.202 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.203 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.204 +
1.205 + plot = gsn_csm_contour_map_ce(wks,classob,resg)
1.206 + frame(wks)
1.207 +
1.208 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.209 + system("rm "+plot_name+"."+plot_type)
1.210 + system("rm "+plot_name+"-1."+plot_type_new)
1.211 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.212 +
1.213 + clear (wks)
1.214 +;*******************************************************************
1.215 +; Calculate "nice" bins for binning the data in equally spaced ranges
1.216 +;********************************************************************
1.217 + nclassn = nclass + 1
1.218 + range = fspan(0,nclassn-1,nclassn)
1.219 +; print (range)
1.220 +
1.221 +; Use this range information to grab all the values in a
1.222 +; particular range, and then take an average.
1.223 +
1.224 + nr = dimsizes(range)
1.225 + nx = nr-1
1.226 + xvalues = new((/2,nx/),float)
1.227 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.228 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.229 + dx4 = dx/4 ; 1/4 of the range
1.230 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.231 +;-----------------------------------------------------------------
1.232 +;(A) mean
1.233 +;--------------------------------------------------------------------
1.234 +; get data
1.235 +
1.236 + laiob_mean = dim_avg_Wrap(laiob(lat|:,lon|:,time|:))
1.237 + laimod_mean = dim_avg_Wrap(laimod(lat|:,lon|:,time|:))
1.238 +
1.239 + DATA11_1D = ndtooned(classob)
1.240 + DATA12_1D = ndtooned(laiob_mean)
1.241 + DATA22_1D = ndtooned(laimod_mean)
1.242 +
1.243 + yvalues = new((/2,nx/),float)
1.244 + mn_yvalues = new((/2,nx/),float)
1.245 + mx_yvalues = new((/2,nx/),float)
1.246 +
1.247 + do nd=0,1
1.248 +
1.249 +; See if we are doing model or observational data.
1.250 +
1.251 + if(nd.eq.0) then
1.252 + data_ob = DATA11_1D
1.253 + data_mod = DATA12_1D
1.254 + else
1.255 + data_ob = DATA11_1D
1.256 + data_mod = DATA22_1D
1.257 + end if
1.258 +
1.259 +; Loop through each range and check for values.
1.260 +
1.261 + do i=0,nr-2
1.262 + if (i.ne.(nr-2)) then
1.263 +; print("")
1.264 +; print("In range ["+range(i)+","+range(i+1)+")")
1.265 + idx = ind((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1)))
1.266 + else
1.267 +; print("")
1.268 +; print("In range ["+range(i)+",)")
1.269 + idx = ind(data_ob.ge.range(i))
1.270 + end if
1.271 +
1.272 +; Calculate average, and get min and max.
1.273 +
1.274 + if(.not.any(ismissing(idx))) then
1.275 + yvalues(nd,i) = avg(data_mod(idx))
1.276 + mn_yvalues(nd,i) = min(data_mod(idx))
1.277 + mx_yvalues(nd,i) = max(data_mod(idx))
1.278 + count = dimsizes(idx)
1.279 + else
1.280 + count = 0
1.281 + yvalues(nd,i) = yvalues@_FillValue
1.282 + mn_yvalues(nd,i) = yvalues@_FillValue
1.283 + mx_yvalues(nd,i) = yvalues@_FillValue
1.284 + end if
1.285 +
1.286 +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1.287 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.288 +
1.289 +; Clean up for next time in loop.
1.290 +
1.291 + delete(idx)
1.292 + end do
1.293 + delete(data_ob)
1.294 + delete(data_mod)
1.295 + end do
1.296 +;-----------------------------------------------------------------
1.297 +; compute correlation coef and M score
1.298 +
1.299 + u = yvalues(0,:)
1.300 + v = yvalues(1,:)
1.301 +
1.302 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.303 + uu = u(good)
1.304 + vv = v(good)
1.305 +
1.306 + ccrMean = esccr(uu,vv,0)
1.307 +; print (ccrMean)
1.308 +
1.309 +; new eq
1.310 + bias = sum(abs(vv-uu)/(vv+uu))
1.311 + Mmean = (1.- (bias/dimsizes(uu)))*5.
1.312 +
1.313 + print (Mmean)
1.314 +
1.315 + do i=0,nrow-2
1.316 + text4(i,0) = sprintf("%5.2f",u(i))
1.317 + text4(i,1) = sprintf("%5.2f",v(i))
1.318 + text4(i,2) = "-"
1.319 + end do
1.320 + text4(nrow-1,0) = sprintf("%5.2f",avg(u))
1.321 + text4(nrow-1,1) = sprintf("%5.2f",avg(v))
1.322 + text4(nrow-1,2) = sprintf("%5.2f",Mmean)
1.323 +
1.324 + delete (u)
1.325 + delete (v)
1.326 + delete (uu)
1.327 + delete (vv)
1.328 +;--------------------------------------------------------------------
1.329 +; histogram res
1.330 +
1.331 + resm = True
1.332 + resm@gsnMaximize = True
1.333 + resm@gsnDraw = False
1.334 + resm@gsnFrame = False
1.335 + resm@xyMarkLineMode = "Markers"
1.336 + resm@xyMarkerSizeF = 0.014
1.337 + resm@xyMarker = 16
1.338 + resm@xyMarkerColors = (/"Brown","Blue"/)
1.339 +; resm@trYMinF = min(mn_yvalues) - 10.
1.340 +; resm@trYMaxF = max(mx_yvalues) + 10.
1.341 + resm@trYMinF = min(mn_yvalues) - 2
1.342 + resm@trYMaxF = max(mx_yvalues) + 4
1.343 +
1.344 + resm@tiYAxisString = "Mean LAI (Leaf Area Index)"
1.345 + resm@tiXAxisString = "Land Cover Type"
1.346 +
1.347 + max_bar = new((/2,nx/),graphic)
1.348 + min_bar = new((/2,nx/),graphic)
1.349 + max_cap = new((/2,nx/),graphic)
1.350 + min_cap = new((/2,nx/),graphic)
1.351 +
1.352 + lnres = True
1.353 + line_colors = (/"brown","blue"/)
1.354 +;------------------------------------------------------------------
1.355 +; Start the graphics.
1.356 +
1.357 + plot_name = "histogram_mean"
1.358 + title = model_name + " vs Observed"
1.359 + resm@tiMainString = title
1.360 +
1.361 + wks = gsn_open_wks (plot_type,plot_name)
1.362 +;-----------------------------
1.363 +; Add a boxed legend using the more simple method
1.364 +
1.365 + resm@pmLegendDisplayMode = "Always"
1.366 +; resm@pmLegendWidthF = 0.1
1.367 + resm@pmLegendWidthF = 0.08
1.368 + resm@pmLegendHeightF = 0.05
1.369 + resm@pmLegendOrthogonalPosF = -1.17
1.370 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.371 +; resm@pmLegendParallelPosF = 0.18
1.372 + resm@pmLegendParallelPosF = 0.88 ;(rightward)
1.373 +
1.374 +; resm@lgPerimOn = False
1.375 + resm@lgLabelFontHeightF = 0.015
1.376 + resm@xyExplicitLegendLabels = (/"observed",model_name/)
1.377 +;-----------------------------
1.378 + tRes = True
1.379 + tRes@txFontHeightF = 0.025
1.380 +
1.381 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")"
1.382 +
1.383 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.384 +
1.385 + xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
1.386 +;-------------------------------
1.387 +;Attach the vertical bar and the horizontal cap line
1.388 +
1.389 + do nd=0,1
1.390 + lnres@gsLineColor = line_colors(nd)
1.391 + do i=0,nx-1
1.392 +
1.393 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.394 + .not.ismissing(mx_yvalues(nd,i))) then
1.395 +
1.396 +; Attach the vertical bar, both above and below the marker.
1.397 +
1.398 + x1 = xvalues(nd,i)
1.399 + y1 = yvalues(nd,i)
1.400 + y2 = mn_yvalues(nd,i)
1.401 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.402 +
1.403 + y2 = mx_yvalues(nd,i)
1.404 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.405 +
1.406 +; Attach the horizontal cap line, both above and below the marker.
1.407 +
1.408 + x1 = xvalues(nd,i) - dx4
1.409 + x2 = xvalues(nd,i) + dx4
1.410 + y1 = mn_yvalues(nd,i)
1.411 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.412 +
1.413 + y1 = mx_yvalues(nd,i)
1.414 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.415 + end if
1.416 + end do
1.417 + end do
1.418 +
1.419 + draw(xy)
1.420 + frame(wks)
1.421 +
1.422 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.423 + system("rm "+plot_name+"."+plot_type)
1.424 +; system("rm "+plot_name+"-1."+plot_type_new)
1.425 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.426 +
1.427 + clear (wks)
1.428 +
1.429 + delete (DATA11_1D)
1.430 + delete (DATA12_1D)
1.431 + delete (DATA22_1D)
1.432 +;delete (range)
1.433 +;delete (xvalues)
1.434 + delete (yvalues)
1.435 + delete (mn_yvalues)
1.436 + delete (mx_yvalues)
1.437 + delete (good)
1.438 + delete (max_bar)
1.439 + delete (min_bar)
1.440 + delete (max_cap)
1.441 + delete (min_cap)
1.442 +;-----------------------------------------------------------------
1.443 +;global res
1.444 +
1.445 + resg = True ; Use plot options
1.446 + resg@cnFillOn = True ; Turn on color fill
1.447 + resg@gsnSpreadColors = True ; use full colormap
1.448 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.449 +; resg@lbLabelAutoStride = True
1.450 + resg@cnLinesOn = False ; Turn off contourn lines
1.451 + resg@mpFillOn = False ; Turn off map fill
1.452 +
1.453 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.454 + resg@cnMinLevelValF = 0. ; Min level
1.455 + resg@cnMaxLevelValF = 10. ; Max level
1.456 + resg@cnLevelSpacingF = 1. ; interval
1.457 +
1.458 +;global contour ob
1.459 +
1.460 + delta = 0.00001
1.461 + laiob_mean = where(ismissing(laiob_mean).and.(ismissing(laimod_mean).or.(laimod_mean.lt.delta)),0.,laiob_mean)
1.462 +
1.463 + plot_name = "global_mean_ob"
1.464 + title = ob_name
1.465 + resg@tiMainString = title
1.466 +
1.467 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.468 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.469 +
1.470 + plot = gsn_csm_contour_map_ce(wks,laiob_mean,resg)
1.471 + frame(wks)
1.472 +
1.473 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.474 + system("rm "+plot_name+"."+plot_type)
1.475 + system("rm "+plot_name+"-1."+plot_type_new)
1.476 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.477 +
1.478 + clear (wks)
1.479 +;------------------------------------------------------------------------
1.480 +;global contour model
1.481 +
1.482 + plot_name = "global_mean_model"
1.483 + title = "Model " + model_name
1.484 + resg@tiMainString = title
1.485 +
1.486 + wks = gsn_open_wks (plot_type,plot_name)
1.487 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.488 +
1.489 + plot = gsn_csm_contour_map_ce(wks,laimod_mean,resg)
1.490 + frame(wks)
1.491 +
1.492 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.493 + system("rm "+plot_name+"."+plot_type)
1.494 + system("rm "+plot_name+"-1."+plot_type_new)
1.495 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.496 +
1.497 + clear (wks)
1.498 +;-----------------------------------------------------------------
1.499 +;(B) max
1.500 +;--------------------------------------------------------------------
1.501 +; get data
1.502 +
1.503 +; observed
1.504 + laiob_max = laiob(0,:,:)
1.505 + s = laiob(:,0,0)
1.506 + laiob_max@long_name = "Leaf Area Index Max"
1.507 +
1.508 + dsizes_z = dimsizes(laiob)
1.509 + nlat = dsizes_z(1)
1.510 + nlon = dsizes_z(2)
1.511 +
1.512 + do j = 0,nlat-1
1.513 + do i = 0,nlon-1
1.514 + s = laiob(:,j,i)
1.515 + laiob_max(j,i) = max(s)
1.516 + end do
1.517 + end do
1.518 +
1.519 +; print (min(y)+"/"+max(y))
1.520 + delete (s)
1.521 + delete (dsizes_z)
1.522 +;-------------------------
1.523 +; model
1.524 + laimod_max = laimod(0,:,:)
1.525 + s = laimod(:,0,0)
1.526 + laimod_max@long_name = "Leaf Area Index Max"
1.527 +
1.528 + dsizes_z = dimsizes(laimod)
1.529 + nlat = dsizes_z(1)
1.530 + nlon = dsizes_z(2)
1.531 +
1.532 + do j = 0,nlat-1
1.533 + do i = 0,nlon-1
1.534 + s = laimod(:,j,i)
1.535 + laimod_max(j,i) = max(s)
1.536 + end do
1.537 + end do
1.538 +
1.539 +; print (min(laimod_max)+"/"+max(laimod_max))
1.540 + delete (s)
1.541 + delete (dsizes_z)
1.542 +;------------------------
1.543 + DATA11_1D = ndtooned(classob)
1.544 + DATA12_1D = ndtooned(laiob_max)
1.545 + DATA22_1D = ndtooned(laimod_max)
1.546 +
1.547 + yvalues = new((/2,nx/),float)
1.548 + mn_yvalues = new((/2,nx/),float)
1.549 + mx_yvalues = new((/2,nx/),float)
1.550 +
1.551 + do nd=0,1
1.552 +
1.553 +; See if we are doing model or observational data.
1.554 +
1.555 + if(nd.eq.0) then
1.556 + data_ob = DATA11_1D
1.557 + data_mod = DATA12_1D
1.558 + else
1.559 + data_ob = DATA11_1D
1.560 + data_mod = DATA22_1D
1.561 + end if
1.562 +
1.563 +; Loop through each range and check for values.
1.564 +
1.565 + do i=0,nr-2
1.566 + if (i.ne.(nr-2)) then
1.567 +; print("")
1.568 +; print("In range ["+range(i)+","+range(i+1)+")")
1.569 + idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))
1.570 + else
1.571 +; print("")
1.572 +; print("In range ["+range(i)+",)")
1.573 + idx = ind(range(i).le.data_ob)
1.574 + end if
1.575 +
1.576 +; Calculate average, and get min and max.
1.577 +
1.578 + if(.not.any(ismissing(idx))) then
1.579 + yvalues(nd,i) = avg(data_mod(idx))
1.580 + mn_yvalues(nd,i) = min(data_mod(idx))
1.581 + mx_yvalues(nd,i) = max(data_mod(idx))
1.582 + count = dimsizes(idx)
1.583 + else
1.584 + count = 0
1.585 + yvalues(nd,i) = yvalues@_FillValue
1.586 + mn_yvalues(nd,i) = yvalues@_FillValue
1.587 + mx_yvalues(nd,i) = yvalues@_FillValue
1.588 + end if
1.589 +
1.590 +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1.591 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.592 +
1.593 +; Clean up for next time in loop.
1.594 +
1.595 + delete(idx)
1.596 + end do
1.597 + delete(data_ob)
1.598 + delete(data_mod)
1.599 + end do
1.600 +;-----------------------------------------------------------------
1.601 +; compute correlation coef and M score
1.602 +
1.603 + u = yvalues(0,:)
1.604 + v = yvalues(1,:)
1.605 +
1.606 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.607 + uu = u(good)
1.608 + vv = v(good)
1.609 +
1.610 + ccrMax = esccr(uu,vv,0)
1.611 +; print (ccrMax)
1.612 +
1.613 +; new eq
1.614 + bias = sum(abs(vv-uu)/(vv+uu))
1.615 + Mmax = (1.- (bias/dimsizes(uu)))*5.
1.616 +
1.617 + print (Mmax)
1.618 +
1.619 + do i=0,nrow-2
1.620 + text4(i,3) = sprintf("%5.2f",u(i))
1.621 + text4(i,4) = sprintf("%5.2f",v(i))
1.622 + text4(i,5) = "-"
1.623 + end do
1.624 + text4(nrow-1,3) = sprintf("%5.2f",avg(u))
1.625 + text4(nrow-1,4) = sprintf("%5.2f",avg(v))
1.626 + text4(nrow-1,5) = sprintf("%5.2f",Mmax)
1.627 +
1.628 + delete (u)
1.629 + delete (v)
1.630 + delete (uu)
1.631 + delete (vv)
1.632 +;--------------------------------------------------------------------
1.633 +; histogram res
1.634 +
1.635 + resm = True
1.636 + resm@gsnMaximize = True
1.637 + resm@gsnDraw = False
1.638 + resm@gsnFrame = False
1.639 + resm@xyMarkLineMode = "Markers"
1.640 + resm@xyMarkerSizeF = 0.014
1.641 + resm@xyMarker = 16
1.642 + resm@xyMarkerColors = (/"Brown","Blue"/)
1.643 +; resm@trYMinF = min(mn_yvalues) - 10.
1.644 +; resm@trYMaxF = max(mx_yvalues) + 10.
1.645 + resm@trYMinF = min(mn_yvalues) - 2
1.646 + resm@trYMaxF = max(mx_yvalues) + 4
1.647 +
1.648 + resm@tiYAxisString = "Max LAI (Leaf Area Index)"
1.649 + resm@tiXAxisString = "Land Cover Type"
1.650 +
1.651 + max_bar = new((/2,nx/),graphic)
1.652 + min_bar = new((/2,nx/),graphic)
1.653 + max_cap = new((/2,nx/),graphic)
1.654 + min_cap = new((/2,nx/),graphic)
1.655 +
1.656 + lnres = True
1.657 + line_colors = (/"brown","blue"/)
1.658 +;------------------------------------------------------------------
1.659 +; Start the graphics.
1.660 +
1.661 + plot_name = "histogram_max"
1.662 + title = model_name + " vs Observed"
1.663 + resm@tiMainString = title
1.664 +
1.665 + wks = gsn_open_wks (plot_type,plot_name)
1.666 +;-----------------------------
1.667 +; Add a boxed legend using the more simple method
1.668 +
1.669 + resm@pmLegendDisplayMode = "Always"
1.670 +; resm@pmLegendWidthF = 0.1
1.671 + resm@pmLegendWidthF = 0.08
1.672 + resm@pmLegendHeightF = 0.05
1.673 + resm@pmLegendOrthogonalPosF = -1.17
1.674 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.675 +; resm@pmLegendParallelPosF = 0.18
1.676 + resm@pmLegendParallelPosF = 0.88 ;(rightward)
1.677 +
1.678 +; resm@lgPerimOn = False
1.679 + resm@lgLabelFontHeightF = 0.015
1.680 + resm@xyExplicitLegendLabels = (/"observed",model_name/)
1.681 +;-----------------------------
1.682 + tRes = True
1.683 + tRes@txFontHeightF = 0.025
1.684 +
1.685 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")"
1.686 +
1.687 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.688 +
1.689 + xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
1.690 +;-------------------------------
1.691 +;Attach the vertical bar and the horizontal cap line
1.692 +
1.693 + do nd=0,1
1.694 + lnres@gsLineColor = line_colors(nd)
1.695 + do i=0,nx-1
1.696 +
1.697 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.698 + .not.ismissing(mx_yvalues(nd,i))) then
1.699 +
1.700 +; Attach the vertical bar, both above and below the marker.
1.701 +
1.702 + x1 = xvalues(nd,i)
1.703 + y1 = yvalues(nd,i)
1.704 + y2 = mn_yvalues(nd,i)
1.705 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.706 +
1.707 + y2 = mx_yvalues(nd,i)
1.708 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.709 +
1.710 +; Attach the horizontal cap line, both above and below the marker.
1.711 +
1.712 + x1 = xvalues(nd,i) - dx4
1.713 + x2 = xvalues(nd,i) + dx4
1.714 + y1 = mn_yvalues(nd,i)
1.715 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.716 +
1.717 + y1 = mx_yvalues(nd,i)
1.718 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.719 + end if
1.720 + end do
1.721 + end do
1.722 +
1.723 + draw(xy)
1.724 + frame(wks)
1.725 +
1.726 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.727 + system("rm "+plot_name+"."+plot_type)
1.728 +; system("rm "+plot_name+"-1."+plot_type_new)
1.729 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.730 +
1.731 + clear (wks)
1.732 +
1.733 + delete (DATA11_1D)
1.734 + delete (DATA12_1D)
1.735 + delete (DATA22_1D)
1.736 +;delete (range)
1.737 +;delete (xvalues)
1.738 + delete (yvalues)
1.739 + delete (mn_yvalues)
1.740 + delete (mx_yvalues)
1.741 + delete (good)
1.742 + delete (max_bar)
1.743 + delete (min_bar)
1.744 + delete (max_cap)
1.745 + delete (min_cap)
1.746 +;-----------------------------------------------------------------
1.747 +;global res
1.748 +
1.749 + resg = True ; Use plot options
1.750 + resg@cnFillOn = True ; Turn on color fill
1.751 + resg@gsnSpreadColors = True ; use full colormap
1.752 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.753 +; resg@lbLabelAutoStride = True
1.754 + resg@cnLinesOn = False ; Turn off contourn lines
1.755 + resg@mpFillOn = False ; Turn off map fill
1.756 +
1.757 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.758 + resg@cnMinLevelValF = 0. ; Min level
1.759 + resg@cnMaxLevelValF = 10. ; Max level
1.760 + resg@cnLevelSpacingF = 1. ; interval
1.761 +
1.762 +;global contour ob
1.763 +
1.764 + delta = 0.00001
1.765 + laiob_max = where(ismissing(laiob_max).and.(ismissing(laimod_max).or.(laimod_max.lt.delta)),0.,laiob_max)
1.766 +
1.767 + plot_name = "global_max_ob"
1.768 + title = ob_name
1.769 + resg@tiMainString = title
1.770 +
1.771 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.772 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.773 +
1.774 + plot_max = gsn_csm_contour_map_ce(wks,laiob_max,resg)
1.775 + frame(wks)
1.776 +
1.777 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.778 + system("rm "+plot_name+"."+plot_type)
1.779 + system("rm "+plot_name+"-1."+plot_type_new)
1.780 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.781 +
1.782 + clear (wks)
1.783 +;------------------------------------------------------------------------
1.784 +;global contour model
1.785 +
1.786 + plot_name = "global_max_model"
1.787 + title = "Model " + model_name
1.788 + resg@tiMainString = title
1.789 +
1.790 + wks = gsn_open_wks (plot_type,plot_name)
1.791 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.792 +
1.793 + plot_max = gsn_csm_contour_map_ce(wks,laimod_max,resg)
1.794 + frame(wks)
1.795 +
1.796 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.797 + system("rm "+plot_name+"."+plot_type)
1.798 + system("rm "+plot_name+"-1."+plot_type_new)
1.799 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.800 +
1.801 + clear (wks)
1.802 +;------------------------------------------------------------------------
1.803 +;(C) phase
1.804 +;--------------------------------------------------------------------
1.805 +; get data
1.806 +
1.807 +; observed
1.808 + laiob_phase = laiob(0,:,:)
1.809 + s = laiob(:,0,0)
1.810 + laiob_phase@long_name = "Leaf Area Index Max Month"
1.811 +
1.812 + dsizes_z = dimsizes(laiob)
1.813 + nlat = dsizes_z(1)
1.814 + nlon = dsizes_z(2)
1.815 +
1.816 + do j = 0,nlat-1
1.817 + do i = 0,nlon-1
1.818 + s = laiob(:,j,i)
1.819 + laiob_phase(j,i) = maxind(s) + 1
1.820 + end do
1.821 + end do
1.822 +
1.823 +; print (min(laiob_phase)+"/"+max(laiob_phase))
1.824 + delete (s)
1.825 + delete (dsizes_z)
1.826 +;-------------------------
1.827 +; model
1.828 + laimod_phase = laimod(0,:,:)
1.829 + s = laimod(:,0,0)
1.830 + laimod_phase@long_name = "Leaf Area Index Max Month"
1.831 +
1.832 + dsizes_z = dimsizes(laimod)
1.833 + nlat = dsizes_z(1)
1.834 + nlon = dsizes_z(2)
1.835 +
1.836 + do j = 0,nlat-1
1.837 + do i = 0,nlon-1
1.838 + s = laimod(:,j,i)
1.839 + laimod_phase(j,i) = maxind(s) + 1
1.840 + end do
1.841 + end do
1.842 +
1.843 +; print (min(laimod_phase)+"/"+max(laimod_phase))
1.844 + delete (s)
1.845 + delete (dsizes_z)
1.846 +;------------------------
1.847 + DATA11_1D = ndtooned(classob)
1.848 + DATA12_1D = ndtooned(laiob_phase)
1.849 + DATA22_1D = ndtooned(laimod_phase)
1.850 +
1.851 + yvalues = new((/2,nx/),float)
1.852 + mn_yvalues = new((/2,nx/),float)
1.853 + mx_yvalues = new((/2,nx/),float)
1.854 +
1.855 + do nd=0,1
1.856 +
1.857 +; See if we are doing model or observational data.
1.858 +
1.859 + if(nd.eq.0) then
1.860 + data_ob = DATA11_1D
1.861 + data_mod = DATA12_1D
1.862 + else
1.863 + data_ob = DATA11_1D
1.864 + data_mod = DATA22_1D
1.865 + end if
1.866 +
1.867 +; Loop through each range and check for values.
1.868 +
1.869 + do i=0,nr-2
1.870 + if (i.ne.(nr-2)) then
1.871 +; print("")
1.872 +; print("In range ["+range(i)+","+range(i+1)+")")
1.873 + idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))
1.874 + else
1.875 +; print("")
1.876 +; print("In range ["+range(i)+",)")
1.877 + idx = ind(range(i).le.data_ob)
1.878 + end if
1.879 +
1.880 +; Calculate average, and get min and max.
1.881 +
1.882 + if(.not.any(ismissing(idx))) then
1.883 + yvalues(nd,i) = avg(data_mod(idx))
1.884 + mn_yvalues(nd,i) = min(data_mod(idx))
1.885 + mx_yvalues(nd,i) = max(data_mod(idx))
1.886 + count = dimsizes(idx)
1.887 + else
1.888 + count = 0
1.889 + yvalues(nd,i) = yvalues@_FillValue
1.890 + mn_yvalues(nd,i) = yvalues@_FillValue
1.891 + mx_yvalues(nd,i) = yvalues@_FillValue
1.892 + end if
1.893 +
1.894 +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1.895 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.896 +
1.897 +; Clean up for next time in loop.
1.898 +
1.899 + delete(idx)
1.900 + end do
1.901 + delete(data_ob)
1.902 + delete(data_mod)
1.903 + end do
1.904 +;-----------------------------------------------------------------
1.905 +; compute correlation coef and M score
1.906 +
1.907 + u = yvalues(0,:)
1.908 + v = yvalues(1,:)
1.909 +
1.910 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.911 + uu = u(good)
1.912 + vv = v(good)
1.913 +
1.914 + ccrPhase = esccr(uu,vv,0)
1.915 +; print (ccrPhase)
1.916 +
1.917 +; old eq
1.918 +; bias = abs(avg(vv)-avg(uu))
1.919 +; new eq
1.920 + bias = avg(abs(vv-uu))
1.921 +
1.922 + bias = where((bias.gt. 6.),12.-bias,bias)
1.923 + Mphase = ((6. - bias)/6.)*5.
1.924 +
1.925 + print (Mphase)
1.926 +
1.927 + do i=0,nrow-2
1.928 + text4(i,6) = sprintf("%5.2f",u(i))
1.929 + text4(i,7) = sprintf("%5.2f",v(i))
1.930 + text4(i,8) = "-"
1.931 + end do
1.932 + text4(nrow-1,6) = sprintf("%5.2f",avg(u))
1.933 + text4(nrow-1,7) = sprintf("%5.2f",avg(v))
1.934 + text4(nrow-1,8) = sprintf("%5.2f",Mphase)
1.935 +
1.936 + delete (u)
1.937 + delete (v)
1.938 + delete (uu)
1.939 + delete (vv)
1.940 +;--------------------------------------------------------------------
1.941 +; histogram res
1.942 +
1.943 + resm = True
1.944 + resm@gsnMaximize = True
1.945 + resm@gsnDraw = False
1.946 + resm@gsnFrame = False
1.947 + resm@xyMarkLineMode = "Markers"
1.948 + resm@xyMarkerSizeF = 0.014
1.949 + resm@xyMarker = 16
1.950 + resm@xyMarkerColors = (/"Brown","Blue"/)
1.951 +; resm@trYMinF = min(mn_yvalues) - 10.
1.952 +; resm@trYMaxF = max(mx_yvalues) + 10.
1.953 + resm@trYMinF = min(mn_yvalues) - 2
1.954 + resm@trYMaxF = max(mx_yvalues) + 4
1.955 +
1.956 + resm@tiYAxisString = "Max LAI (Leaf Area Index) Month"
1.957 + resm@tiXAxisString = "Land Cover Type"
1.958 +
1.959 + max_bar = new((/2,nx/),graphic)
1.960 + min_bar = new((/2,nx/),graphic)
1.961 + max_cap = new((/2,nx/),graphic)
1.962 + min_cap = new((/2,nx/),graphic)
1.963 +
1.964 + lnres = True
1.965 + line_colors = (/"brown","blue"/)
1.966 +;------------------------------------------------------------------
1.967 +; Start the graphics.
1.968 +
1.969 + plot_name = "histogram_phase"
1.970 + title = model_name + " vs Observed"
1.971 + resm@tiMainString = title
1.972 +
1.973 + wks = gsn_open_wks (plot_type,plot_name)
1.974 +;-----------------------------
1.975 +; Add a boxed legend using the more simple method
1.976 +
1.977 + resm@pmLegendDisplayMode = "Always"
1.978 +; resm@pmLegendWidthF = 0.1
1.979 + resm@pmLegendWidthF = 0.08
1.980 + resm@pmLegendHeightF = 0.05
1.981 + resm@pmLegendOrthogonalPosF = -1.17
1.982 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.983 +; resm@pmLegendParallelPosF = 0.18
1.984 + resm@pmLegendParallelPosF = 0.88 ;(rightward)
1.985 +
1.986 +; resm@lgPerimOn = False
1.987 + resm@lgLabelFontHeightF = 0.015
1.988 + resm@xyExplicitLegendLabels = (/"observed",model_name/)
1.989 +;-----------------------------
1.990 + tRes = True
1.991 + tRes@txFontHeightF = 0.025
1.992 +
1.993 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")"
1.994 +
1.995 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.996 +
1.997 + xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
1.998 +;-------------------------------
1.999 +;Attach the vertical bar and the horizontal cap line
1.1000 +
1.1001 + do nd=0,1
1.1002 + lnres@gsLineColor = line_colors(nd)
1.1003 + do i=0,nx-1
1.1004 +
1.1005 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1006 + .not.ismissing(mx_yvalues(nd,i))) then
1.1007 +
1.1008 +; Attach the vertical bar, both above and below the marker.
1.1009 +
1.1010 + x1 = xvalues(nd,i)
1.1011 + y1 = yvalues(nd,i)
1.1012 + y2 = mn_yvalues(nd,i)
1.1013 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1014 +
1.1015 + y2 = mx_yvalues(nd,i)
1.1016 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1017 +
1.1018 +; Attach the horizontal cap line, both above and below the marker.
1.1019 +
1.1020 + x1 = xvalues(nd,i) - dx4
1.1021 + x2 = xvalues(nd,i) + dx4
1.1022 + y1 = mn_yvalues(nd,i)
1.1023 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1024 +
1.1025 + y1 = mx_yvalues(nd,i)
1.1026 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1027 + end if
1.1028 + end do
1.1029 + end do
1.1030 +
1.1031 + draw(xy)
1.1032 + frame(wks)
1.1033 +
1.1034 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1035 + system("rm "+plot_name+"."+plot_type)
1.1036 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1037 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1038 +
1.1039 + clear (wks)
1.1040 +
1.1041 + delete (DATA11_1D)
1.1042 + delete (DATA12_1D)
1.1043 + delete (DATA22_1D)
1.1044 +;delete (range)
1.1045 +;delete (xvalues)
1.1046 + delete (yvalues)
1.1047 + delete (mn_yvalues)
1.1048 + delete (mx_yvalues)
1.1049 + delete (good)
1.1050 + delete (max_bar)
1.1051 + delete (min_bar)
1.1052 + delete (max_cap)
1.1053 + delete (min_cap)
1.1054 +;-----------------------------------------------------------------
1.1055 +;global res
1.1056 +
1.1057 + resg = True ; Use plot options
1.1058 + resg@cnFillOn = True ; Turn on color fill
1.1059 + resg@gsnSpreadColors = True ; use full colormap
1.1060 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.1061 +; resg@lbLabelAutoStride = True
1.1062 + resg@cnLinesOn = False ; Turn off contourn lines
1.1063 + resg@mpFillOn = False ; Turn off map fill
1.1064 +
1.1065 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.1066 + resg@cnMinLevelValF = 1. ; Min level
1.1067 + resg@cnMaxLevelValF = 12. ; Max level
1.1068 + resg@cnLevelSpacingF = 1. ; interval
1.1069 +
1.1070 +;global contour ob
1.1071 +
1.1072 + delta = 0.00001
1.1073 + laiob_phase = where(ismissing(laiob_phase).and.(ismissing(laimod_phase).or.(laimod_phase.lt.delta)),0.,laiob_phase)
1.1074 +
1.1075 + plot_name = "global_phase_ob"
1.1076 + title = ob_name
1.1077 + resg@tiMainString = title
1.1078 +
1.1079 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1080 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1081 +
1.1082 + plot = gsn_csm_contour_map_ce(wks,laiob_phase,resg)
1.1083 + frame(wks)
1.1084 +
1.1085 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1086 + system("rm "+plot_name+"."+plot_type)
1.1087 + system("rm "+plot_name+"-1."+plot_type_new)
1.1088 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1089 +
1.1090 + clear (wks)
1.1091 +;------------------------------------------------------------------------
1.1092 +;global contour model
1.1093 +
1.1094 + plot_name = "global_phase_model"
1.1095 + title = "Model " + model_name
1.1096 + resg@tiMainString = title
1.1097 +
1.1098 + wks = gsn_open_wks (plot_type,plot_name)
1.1099 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1100 +
1.1101 + delete (plot)
1.1102 + plot = gsn_csm_contour_map_ce(wks,laimod_phase,resg)
1.1103 + frame(wks)
1.1104 +
1.1105 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1106 + system("rm "+plot_name+"."+plot_type)
1.1107 + system("rm "+plot_name+"-1."+plot_type_new)
1.1108 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1109 +
1.1110 + clear (wks)
1.1111 +;-----------------------------------------------------------------
1.1112 +;(D) grow day
1.1113 +;--------------------------------------------------------------------
1.1114 +; get data
1.1115 +
1.1116 + day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/)
1.1117 +
1.1118 +; observed
1.1119 + laiob_grow = laiob(0,:,:)
1.1120 + laiob_grow@long_name = "Days of Growing Season"
1.1121 +
1.1122 + dsizes_z = dimsizes(laiob)
1.1123 + ntime = dsizes_z(0)
1.1124 + nlat = dsizes_z(1)
1.1125 + nlon = dsizes_z(2)
1.1126 +
1.1127 + do j = 0,nlat-1
1.1128 + do i = 0,nlon-1
1.1129 + nday = 0.
1.1130 + do k = 0,ntime-1
1.1131 + if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then
1.1132 + nday = nday + day_of_data(k)
1.1133 + end if
1.1134 + end do
1.1135 +
1.1136 + laiob_grow(j,i) = nday
1.1137 + end do
1.1138 + end do
1.1139 +
1.1140 +; print (min(laiob_grow)+"/"+max(laiob_grow))
1.1141 +;-------------------------
1.1142 +; model
1.1143 + laimod_grow = laimod(0,:,:)
1.1144 + laimod_grow@long_name = "Days of Growing Season"
1.1145 +
1.1146 + dsizes_z = dimsizes(laimod)
1.1147 + ntime = dsizes_z(0)
1.1148 + nlat = dsizes_z(1)
1.1149 + nlon = dsizes_z(2)
1.1150 +
1.1151 + do j = 0,nlat-1
1.1152 + do i = 0,nlon-1
1.1153 + nday = 0.
1.1154 + do k = 0,ntime-1
1.1155 + if (.not. ismissing(laimod(k,j,i)) .and. laimod(k,j,i) .gt. 1.0) then
1.1156 + nday = nday + day_of_data(k)
1.1157 + end if
1.1158 + end do
1.1159 +
1.1160 + laimod_grow(j,i) = nday
1.1161 + end do
1.1162 + end do
1.1163 +
1.1164 +; print (min(laimod_grow)+"/"+max(laimod_grow))
1.1165 +;------------------------
1.1166 + DATA11_1D = ndtooned(classob)
1.1167 + DATA12_1D = ndtooned(laiob_grow)
1.1168 + DATA22_1D = ndtooned(laimod_grow)
1.1169 +
1.1170 + yvalues = new((/2,nx/),float)
1.1171 + mn_yvalues = new((/2,nx/),float)
1.1172 + mx_yvalues = new((/2,nx/),float)
1.1173 +
1.1174 + do nd=0,1
1.1175 +
1.1176 +; See if we are doing model or observational data.
1.1177 +
1.1178 + if(nd.eq.0) then
1.1179 + data_ob = DATA11_1D
1.1180 + data_mod = DATA12_1D
1.1181 + else
1.1182 + data_ob = DATA11_1D
1.1183 + data_mod = DATA22_1D
1.1184 + end if
1.1185 +
1.1186 +; Loop through each range and check for values.
1.1187 +
1.1188 + do i=0,nr-2
1.1189 + if (i.ne.(nr-2)) then
1.1190 +; print("")
1.1191 +; print("In range ["+range(i)+","+range(i+1)+")")
1.1192 + idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))
1.1193 + else
1.1194 +; print("")
1.1195 +; print("In range ["+range(i)+",)")
1.1196 + idx = ind(range(i).le.data_ob)
1.1197 + end if
1.1198 +
1.1199 +; Calculate average, and get min and max.
1.1200 +
1.1201 + if(.not.any(ismissing(idx))) then
1.1202 + yvalues(nd,i) = avg(data_mod(idx))
1.1203 + mn_yvalues(nd,i) = min(data_mod(idx))
1.1204 + mx_yvalues(nd,i) = max(data_mod(idx))
1.1205 + count = dimsizes(idx)
1.1206 + else
1.1207 + count = 0
1.1208 + yvalues(nd,i) = yvalues@_FillValue
1.1209 + mn_yvalues(nd,i) = yvalues@_FillValue
1.1210 + mx_yvalues(nd,i) = yvalues@_FillValue
1.1211 + end if
1.1212 +
1.1213 +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1.1214 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.1215 +
1.1216 +; Clean up for next time in loop.
1.1217 +
1.1218 + delete(idx)
1.1219 + end do
1.1220 + delete(data_ob)
1.1221 + delete(data_mod)
1.1222 + end do
1.1223 +;-----------------------------------------------------------------
1.1224 +; compute correlation coef and M score
1.1225 +
1.1226 + u = yvalues(0,:)
1.1227 + v = yvalues(1,:)
1.1228 +
1.1229 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.1230 + uu = u(good)
1.1231 + vv = v(good)
1.1232 +
1.1233 + ccrGrow = esccr(uu,vv,0)
1.1234 +; print (ccrGrow)
1.1235 +
1.1236 +; new eq
1.1237 + bias = sum(abs(vv-uu)/(vv+uu))
1.1238 + Mgrow = (1.- (bias/dimsizes(uu)))*5.
1.1239 +
1.1240 + print (Mgrow)
1.1241 +
1.1242 + do i=0,nrow-2
1.1243 + text4(i,9) = sprintf("%5.2f",u(i))
1.1244 + text4(i,10) = sprintf("%5.2f",v(i))
1.1245 + text4(i,11) = "-"
1.1246 + end do
1.1247 + text4(nrow-1,9) = sprintf("%5.2f",avg(u))
1.1248 + text4(nrow-1,10) = sprintf("%5.2f",avg(v))
1.1249 + text4(nrow-1,11) = sprintf("%5.2f",Mgrow)
1.1250 +
1.1251 + delete (u)
1.1252 + delete (v)
1.1253 + delete (uu)
1.1254 + delete (vv)
1.1255 +;--------------------------------------------------------------------
1.1256 +; histogram res
1.1257 +
1.1258 + resm = True
1.1259 + resm@gsnMaximize = True
1.1260 + resm@gsnDraw = False
1.1261 + resm@gsnFrame = False
1.1262 + resm@xyMarkLineMode = "Markers"
1.1263 + resm@xyMarkerSizeF = 0.014
1.1264 + resm@xyMarker = 16
1.1265 + resm@xyMarkerColors = (/"Brown","Blue"/)
1.1266 +; resm@trYMinF = min(mn_yvalues) - 10.
1.1267 +; resm@trYMaxF = max(mx_yvalues) + 10.
1.1268 + resm@trYMinF = min(mn_yvalues) - 2
1.1269 + resm@trYMaxF = max(mx_yvalues) + 4
1.1270 +
1.1271 + resm@tiYAxisString = "Days of Growing season"
1.1272 + resm@tiXAxisString = "Land Cover Type"
1.1273 +
1.1274 + max_bar = new((/2,nx/),graphic)
1.1275 + min_bar = new((/2,nx/),graphic)
1.1276 + max_cap = new((/2,nx/),graphic)
1.1277 + min_cap = new((/2,nx/),graphic)
1.1278 +
1.1279 + lnres = True
1.1280 + line_colors = (/"brown","blue"/)
1.1281 +;------------------------------------------------------------------
1.1282 +; Start the graphics.
1.1283 +
1.1284 + plot_name = "histogram_grow"
1.1285 + title = model_name + " vs Observed"
1.1286 + resm@tiMainString = title
1.1287 +
1.1288 + wks = gsn_open_wks (plot_type,plot_name)
1.1289 +;-----------------------------
1.1290 +; Add a boxed legend using the more simple method
1.1291 +
1.1292 + resm@pmLegendDisplayMode = "Always"
1.1293 +; resm@pmLegendWidthF = 0.1
1.1294 + resm@pmLegendWidthF = 0.08
1.1295 + resm@pmLegendHeightF = 0.05
1.1296 + resm@pmLegendOrthogonalPosF = -1.17
1.1297 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.1298 +; resm@pmLegendParallelPosF = 0.18
1.1299 + resm@pmLegendParallelPosF = 0.88 ;(rightward)
1.1300 +
1.1301 +; resm@lgPerimOn = False
1.1302 + resm@lgLabelFontHeightF = 0.015
1.1303 + resm@xyExplicitLegendLabels = (/"observed",model_name/)
1.1304 +;-----------------------------
1.1305 + tRes = True
1.1306 + tRes@txFontHeightF = 0.025
1.1307 +
1.1308 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")"
1.1309 +
1.1310 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.1311 +
1.1312 + xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
1.1313 +;-------------------------------
1.1314 +;Attach the vertical bar and the horizontal cap line
1.1315 +
1.1316 + do nd=0,1
1.1317 + lnres@gsLineColor = line_colors(nd)
1.1318 + do i=0,nx-1
1.1319 +
1.1320 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1321 + .not.ismissing(mx_yvalues(nd,i))) then
1.1322 +
1.1323 +; Attach the vertical bar, both above and below the marker.
1.1324 +
1.1325 + x1 = xvalues(nd,i)
1.1326 + y1 = yvalues(nd,i)
1.1327 + y2 = mn_yvalues(nd,i)
1.1328 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1329 +
1.1330 + y2 = mx_yvalues(nd,i)
1.1331 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1332 +
1.1333 +; Attach the horizontal cap line, both above and below the marker.
1.1334 +
1.1335 + x1 = xvalues(nd,i) - dx4
1.1336 + x2 = xvalues(nd,i) + dx4
1.1337 + y1 = mn_yvalues(nd,i)
1.1338 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1339 +
1.1340 + y1 = mx_yvalues(nd,i)
1.1341 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1342 + end if
1.1343 + end do
1.1344 + end do
1.1345 +
1.1346 + draw(xy)
1.1347 + frame(wks)
1.1348 +
1.1349 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1350 + system("rm "+plot_name+"."+plot_type)
1.1351 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1352 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1353 +
1.1354 + clear (wks)
1.1355 +
1.1356 + delete (DATA11_1D)
1.1357 + delete (DATA12_1D)
1.1358 + delete (DATA22_1D)
1.1359 +;delete (range)
1.1360 +;delete (xvalues)
1.1361 + delete (yvalues)
1.1362 + delete (mn_yvalues)
1.1363 + delete (mx_yvalues)
1.1364 + delete (good)
1.1365 + delete (max_bar)
1.1366 + delete (min_bar)
1.1367 + delete (max_cap)
1.1368 + delete (min_cap)
1.1369 +;-----------------------------------------------------------------
1.1370 +;global res
1.1371 +
1.1372 + resg = True ; Use plot options
1.1373 + resg@cnFillOn = True ; Turn on color fill
1.1374 + resg@gsnSpreadColors = True ; use full colormap
1.1375 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.1376 +; resg@lbLabelAutoStride = True
1.1377 + resg@cnLinesOn = False ; Turn off contourn lines
1.1378 + resg@mpFillOn = False ; Turn off map fill
1.1379 +
1.1380 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.1381 + resg@cnMinLevelValF = 60. ; Min level
1.1382 + resg@cnMaxLevelValF = 360. ; Max level
1.1383 + resg@cnLevelSpacingF = 20. ; interval
1.1384 +
1.1385 +;global contour ob
1.1386 +
1.1387 + laiob_grow@_FillValue = 1.e+36
1.1388 + laiob_grow = where(laiob_grow .lt. 10.,laiob_grow@_FillValue,laiob_grow)
1.1389 +
1.1390 + plot_name = "global_grow_ob"
1.1391 + title = ob_name
1.1392 + resg@tiMainString = title
1.1393 +
1.1394 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1395 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1396 +
1.1397 + plot = gsn_csm_contour_map_ce(wks,laiob_grow,resg)
1.1398 + frame(wks)
1.1399 +
1.1400 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1401 + system("rm "+plot_name+"."+plot_type)
1.1402 + system("rm "+plot_name+"-1."+plot_type_new)
1.1403 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1404 +
1.1405 + clear (wks)
1.1406 +;------------------------------------------------------------------------
1.1407 +;global contour model
1.1408 +
1.1409 + laimod_grow@_FillValue = 1.e+36
1.1410 + laimod_grow = where(laimod_grow .lt. 10.,laimod_grow@_FillValue,laimod_grow)
1.1411 +
1.1412 + plot_name = "global_grow_model"
1.1413 + title = "Model " + model_name
1.1414 + resg@tiMainString = title
1.1415 +
1.1416 + wks = gsn_open_wks (plot_type,plot_name)
1.1417 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1418 +
1.1419 + delete (plot)
1.1420 + plot = gsn_csm_contour_map_ce(wks,laimod_grow,resg)
1.1421 + frame(wks)
1.1422 +
1.1423 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1424 + system("rm "+plot_name+"."+plot_type)
1.1425 + system("rm "+plot_name+"-1."+plot_type_new)
1.1426 + system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1427 +
1.1428 + clear (wks)
1.1429 +;------------------------------------------------------------------------
1.1430 +;global contour model vs ob
1.1431 +
1.1432 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.1433 + resg@cnMinLevelValF = 0. ; Min level
1.1434 + resg@cnMaxLevelValF = 10. ; Max level
1.1435 + resg@cnLevelSpacingF = 1. ; interval
1.1436 +
1.1437 + plot_name = "global_mean_model_vs_ob"
1.1438 +
1.1439 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1440 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1441 +
1.1442 + delete (plot)
1.1443 + plot=new(3,graphic) ; create graphic array
1.1444 +
1.1445 + resg@gsnFrame = False ; Do not draw plot
1.1446 + resg@gsnDraw = False ; Do not advance frame
1.1447 +
1.1448 +; plot correlation coef
1.1449 +
1.1450 + gRes = True
1.1451 + gRes@txFontHeightF = 0.02
1.1452 + gRes@txAngleF = 90
1.1453 +
1.1454 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")"
1.1455 +
1.1456 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.1457 +;--------------------------------------------------------------------
1.1458 +
1.1459 +;(a) ob
1.1460 +
1.1461 + title = ob_name
1.1462 + resg@tiMainString = title
1.1463 +
1.1464 + plot(0) = gsn_csm_contour_map_ce(wks,laiob_mean,resg)
1.1465 +
1.1466 +;(b) model
1.1467 +
1.1468 + title = "Model "+ model_name
1.1469 + resg@tiMainString = title
1.1470 +
1.1471 + plot(1) = gsn_csm_contour_map_ce(wks,laimod_mean,resg)
1.1472 +
1.1473 +;(c) model-ob
1.1474 +
1.1475 + zz = laimod_mean
1.1476 + zz = laimod_mean - laiob_mean
1.1477 + title = "Model_"+model_name+" - Observed"
1.1478 + resg@tiMainString = title
1.1479 +
1.1480 + resg@cnMinLevelValF = -2. ; Min level
1.1481 + resg@cnMaxLevelValF = 2. ; Max level
1.1482 + resg@cnLevelSpacingF = 0.4 ; interval
1.1483 +
1.1484 +
1.1485 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.1486 +
1.1487 + pres = True ; panel plot mods desired
1.1488 + pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.1489 + ; indiv. plots in panel
1.1490 + pres@gsnMaximize = True ; fill the page
1.1491 +
1.1492 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.1493 +
1.1494 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1495 + system("rm "+plot_name+"."+plot_type)
1.1496 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1497 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1498 +
1.1499 + frame (wks)
1.1500 + clear (wks)
1.1501 +
1.1502 + delete (plot)
1.1503 +;-----------------------------------------------------------------
1.1504 +;global contour model vs ob
1.1505 +
1.1506 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.1507 + resg@cnMinLevelValF = 0. ; Min level
1.1508 + resg@cnMaxLevelValF = 10. ; Max level
1.1509 + resg@cnLevelSpacingF = 1. ; interval
1.1510 +
1.1511 + plot_name = "global_max_model_vs_ob"
1.1512 +
1.1513 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1514 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1515 +
1.1516 +; delete(plot)
1.1517 + plot=new(3,graphic) ; create graphic array
1.1518 +
1.1519 + resg@gsnFrame = False ; Do not draw plot
1.1520 + resg@gsnDraw = False ; Do not advance frame
1.1521 +
1.1522 +; plot correlation coef
1.1523 +
1.1524 + gRes = True
1.1525 + gRes@txFontHeightF = 0.02
1.1526 + gRes@txAngleF = 90
1.1527 +
1.1528 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")"
1.1529 +
1.1530 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.1531 +;--------------------------------------------------------------------
1.1532 +;(a) ob
1.1533 +
1.1534 + title = ob_name
1.1535 + resg@tiMainString = title
1.1536 +
1.1537 + plot(0) = gsn_csm_contour_map_ce(wks,laiob_max,resg)
1.1538 +
1.1539 +;(b) model
1.1540 +
1.1541 + title = "Model "+ model_name
1.1542 + resg@tiMainString = title
1.1543 +
1.1544 + plot(1) = gsn_csm_contour_map_ce(wks,laimod_max,resg)
1.1545 +
1.1546 +;(c) model-ob
1.1547 +
1.1548 + delete (zz)
1.1549 + zz = laimod_max
1.1550 + zz = laimod_max - laiob_max
1.1551 + title = "Model_"+model_name+" - Observed"
1.1552 + resg@tiMainString = title
1.1553 +
1.1554 + resg@cnMinLevelValF = -6. ; Min level
1.1555 + resg@cnMaxLevelValF = 6. ; Max level
1.1556 + resg@cnLevelSpacingF = 1. ; interval
1.1557 +
1.1558 +
1.1559 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.1560 +
1.1561 + pres = True ; panel plot mods desired
1.1562 + pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.1563 + ; indiv. plots in panel
1.1564 + pres@gsnMaximize = True ; fill the page
1.1565 +
1.1566 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.1567 +
1.1568 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1569 + system("rm "+plot_name+"."+plot_type)
1.1570 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1571 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1572 +
1.1573 + frame (wks)
1.1574 + clear (wks)
1.1575 +
1.1576 + delete (plot)
1.1577 +;-----------------------------------------------------------------
1.1578 +;global contour model vs ob
1.1579 +
1.1580 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.1581 + resg@cnMinLevelValF = 1. ; Min level
1.1582 + resg@cnMaxLevelValF = 12. ; Max level
1.1583 + resg@cnLevelSpacingF = 1. ; interval
1.1584 +
1.1585 + plot_name = "global_phase_model_vs_ob"
1.1586 +
1.1587 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1588 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1589 +
1.1590 +; delete (plot)
1.1591 + plot=new(3,graphic) ; create graphic array
1.1592 +
1.1593 + resg@gsnFrame = False ; Do not draw plot
1.1594 + resg@gsnDraw = False ; Do not advance frame
1.1595 +
1.1596 +; plot correlation coef
1.1597 +
1.1598 + gRes = True
1.1599 + gRes@txFontHeightF = 0.02
1.1600 + gRes@txAngleF = 90
1.1601 +
1.1602 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")"
1.1603 +
1.1604 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.1605 +;--------------------------------------------------------------------
1.1606 +;(a) ob
1.1607 +
1.1608 + title = ob_name
1.1609 + resg@tiMainString = title
1.1610 +
1.1611 + plot(0) = gsn_csm_contour_map_ce(wks,laiob_phase,resg)
1.1612 +
1.1613 +;(b) model
1.1614 +
1.1615 + title = "Model "+ model_name
1.1616 + resg@tiMainString = title
1.1617 +
1.1618 + plot(1) = gsn_csm_contour_map_ce(wks,laimod_phase,resg)
1.1619 +
1.1620 +;(c) model-ob
1.1621 +
1.1622 + delete (zz)
1.1623 + zz = laimod_phase
1.1624 + zz = laimod_phase - laiob_phase
1.1625 + title = "Model_"+model_name+" - Observed"
1.1626 + resg@tiMainString = title
1.1627 +
1.1628 + resg@cnMinLevelValF = -6. ; Min level
1.1629 + resg@cnMaxLevelValF = 6. ; Max level
1.1630 + resg@cnLevelSpacingF = 1. ; interval
1.1631 +
1.1632 +
1.1633 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.1634 +
1.1635 +; pres = True ; panel plot mods desired
1.1636 +; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.1637 + ; indiv. plots in panel
1.1638 +; pres@gsnMaximize = True ; fill the page
1.1639 +
1.1640 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.1641 +
1.1642 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1643 + system("rm "+plot_name+"."+plot_type)
1.1644 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1645 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1646 +
1.1647 + frame (wks)
1.1648 + clear (wks)
1.1649 +
1.1650 + delete (plot)
1.1651 +;------------------------------------------------------------------
1.1652 +;global contour model vs ob
1.1653 +
1.1654 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.1655 + resg@cnMinLevelValF = 60. ; Min level
1.1656 + resg@cnMaxLevelValF = 360. ; Max level
1.1657 + resg@cnLevelSpacingF = 20. ; interval
1.1658 +
1.1659 + plot_name = "global_grow_model_vs_ob"
1.1660 +
1.1661 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1662 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1663 +
1.1664 +; delete (plot)
1.1665 + plot=new(3,graphic) ; create graphic array
1.1666 +
1.1667 + resg@gsnFrame = False ; Do not draw plot
1.1668 + resg@gsnDraw = False ; Do not advance frame
1.1669 +
1.1670 +; plot correlation coef
1.1671 +
1.1672 + gRes = True
1.1673 + gRes@txFontHeightF = 0.02
1.1674 + gRes@txAngleF = 90
1.1675 +
1.1676 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")"
1.1677 +
1.1678 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.1679 +;--------------------------------------------------------------------
1.1680 +;(a) ob
1.1681 +
1.1682 + title = ob_name
1.1683 + resg@tiMainString = title
1.1684 +
1.1685 + plot(0) = gsn_csm_contour_map_ce(wks,laiob_grow,resg)
1.1686 +
1.1687 +;(b) model
1.1688 +
1.1689 + title = "Model "+ model_name
1.1690 + resg@tiMainString = title
1.1691 +
1.1692 + plot(1) = gsn_csm_contour_map_ce(wks,laimod_grow,resg)
1.1693 +
1.1694 +;(c) model-ob
1.1695 +
1.1696 + delete (zz)
1.1697 + zz = laimod_grow
1.1698 + zz = laimod_grow - laiob_grow
1.1699 + title = "Model_"+model_name+" - Observed"
1.1700 + resg@tiMainString = title
1.1701 +
1.1702 + resg@cnMinLevelValF = -100. ; Min level
1.1703 + resg@cnMaxLevelValF = 100. ; Max level
1.1704 + resg@cnLevelSpacingF = 20. ; interval
1.1705 +
1.1706 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.1707 +
1.1708 + pres = True ; panel plot mods desired
1.1709 + pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.1710 + ; indiv. plots in panel
1.1711 + pres@gsnMaximize = True ; fill the page
1.1712 +
1.1713 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.1714 +
1.1715 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1716 + system("rm "+plot_name+"."+plot_type)
1.1717 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1718 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1719 +
1.1720 + frame (wks)
1.1721 + clear (wks)
1.1722 +
1.1723 + delete (plot)
1.1724 +;**************************************************
1.1725 +; plot lai table
1.1726 +;**************************************************
1.1727 +
1.1728 + plot_name = "table_lai"
1.1729 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1730 +
1.1731 + gsn_table(wks,ncr1,xx1,yy1,text1,res1)
1.1732 + gsn_table(wks,ncr21,xx21,yy21,text21,res21)
1.1733 + gsn_table(wks,ncr22,xx22,yy22,text22,res22)
1.1734 + gsn_table(wks,ncr3,xx3,yy3,text3,res3)
1.1735 + gsn_table(wks,ncr4,xx4,yy4,text4,res4)
1.1736 +
1.1737 + frame(wks)
1.1738 +
1.1739 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
1.1740 +; system("rm "+plot_name+"."+plot_type)
1.1741 +; system("rm "+plot_name+"-1."+plot_type_new)
1.1742 +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
1.1743 +;-------------------------------------------------------------------
1.1744 + temp_name = "temp." + model_name
1.1745 + system("mkdir -p " + temp_name)
1.1746 + system("mv *.png " + temp_name)
1.1747 + system("tar cf "+ temp_name +".tar " + temp_name)
1.1748 +;-------------------------------------------------------------------
1.1749 +end
1.1750 +