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