forrest@0: ;******************************************************** forrest@0: ; histogram normalized by rain and compute correleration forrest@0: ;******************************************************** forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test" forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test" forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" forrest@0: forrest@0: procedure pminmax(data:numeric,name:string) forrest@0: begin forrest@0: print ("min/max " + name + " = " + min(data) + "/" + max(data)) forrest@0: if(isatt(data,"units")) then forrest@0: print (name + " units = " + data@units) forrest@0: end if forrest@0: end forrest@0: forrest@0: ; Main code. forrest@0: begin forrest@0: forrest@0: nclass = 20 forrest@0: forrest@0: plot_type = "ps" forrest@0: plot_type_new = "png" forrest@0: forrest@0: ;************************************************ forrest@0: ; read in data: model forrest@0: ;************************************************ forrest@0: ;film = "b30.061n_1995-2004_MONS_climo_lnd.nc" forrest@0: ;model_name = "b30.061n" forrest@0: ;model_grid = "T31" forrest@0: forrest@0: ;film = "newcn05_ncep_1i_MONS_climo_lnd.nc" forrest@0: ;model_name = "newcn" forrest@0: ;model_grid = "1.9" forrest@0: forrest@0: ;film = "i01.06cn_1798-2004_MONS_climo.nc" forrest@0: ;model_name = "06cn" forrest@0: ;model_grid = "T42" forrest@0: forrest@0: ;film = "i01.06casa_1798-2004_MONS_climo.nc" forrest@0: ;model_name = "06casa" forrest@0: ;model_grid = "T42" forrest@0: forrest@0: film = "i01.10cn_1948-2004_MONS_climo.nc" forrest@0: model_name = "10cn" forrest@0: model_grid = "T42" forrest@0: forrest@0: ;film = "i01.10casa_1948-2004_MONS_climo.nc" forrest@0: ;model_name = "10casa" forrest@0: ;model_grid = "T42" forrest@0: forrest@0: html_name = "table.html." + model_name forrest@0: html_new = html_name +".new" forrest@0: system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \ forrest@0: "mv -f "+html_new+" "+html_name) forrest@0: ;------------------------------------------------ forrest@0: dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" forrest@0: fm = addfile(dirm+film,"r") forrest@0: forrest@0: laimod = fm->TLAI forrest@0: forrest@0: ;************************************************ forrest@0: ; read in data: observed forrest@0: ;************************************************ forrest@0: forrest@0: ob_name = "MODIS MOD 15A2 2000-2005" forrest@0: forrest@0: diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/" forrest@0: filo1 = "land_class_"+model_grid+".nc" forrest@0: filo2 = "LAI_2000-2005_MONS_"+model_grid+".nc" forrest@0: forrest@0: fo1 = addfile(diro+filo1,"r") forrest@0: fo2 = addfile(diro+filo2,"r") forrest@0: forrest@0: classob = tofloat(fo1->LAND_CLASS) forrest@0: laiob = fo2->LAI forrest@0: forrest@0: ;******************************************************************* forrest@0: ; for plotting table forrest@0: ;******************************************************************* forrest@0: ; table header name forrest@0: table_header_name = "LAI" forrest@0: forrest@0: ; column (not including header column) forrest@0: forrest@0: col_header1 = (/"Mean","Max","Phase","Growth"/) forrest@0: col_header2 = (/"ob","model","M" \ forrest@0: ,"ob","model","M" \ forrest@0: ,"ob","model","M" \ forrest@0: ,"ob","model","M" \ forrest@0: /) forrest@0: forrest@0: ncol1 = dimsizes(col_header1) forrest@0: ncol2 = dimsizes(col_header2) forrest@0: ncol = ncol2 forrest@0: forrest@0: ; row (not including header row) forrest@0: row_header = (/"Water Bodies" \ forrest@0: ,"Evergreen Needleleaf Forests" \ forrest@0: ,"Evergreen Broadleaf Forests" \ forrest@0: ,"Deciduous Needleleaf Forest" \ forrest@0: ,"Deciduous Broadleaf Forests" \ forrest@0: ,"Mixed Forests" \ forrest@0: ,"Closed Bushlands" \ forrest@0: ,"Open Bushlands" \ forrest@0: ,"Woody Savannas (S. Hem.)" \ forrest@0: ,"Savannas (S. Hem.)" \ forrest@0: ,"Grasslands" \ forrest@0: ,"Permanent Wetlands" \ forrest@0: ,"Croplands" \ forrest@0: ,"Urban and Built-Up" \ forrest@0: ,"Cropland/Natural Vegetation Mosaic" \ forrest@0: ,"Permanent Snow and Ice" \ forrest@0: ,"Barren or Sparsely Vegetated" \ forrest@0: ,"Unclassified" \ forrest@0: ,"Woody Savannas (N. Hem.)" \ forrest@0: ,"Savannas (N. Hem.)" \ forrest@0: ,"All biome average" \ forrest@0: /) forrest@0: nrow = dimsizes(row_header) forrest@0: forrest@0: ; arrays to be passed to table. forrest@0: value = new ((/nrow, ncol/),string ) forrest@0: forrest@0: table_length = 0.995 forrest@0: forrest@0: ; Table header forrest@0: ncr1 = (/1,1/) ; 1 row, 1 column forrest@0: xx1 = (/0.005,0.25/) ; Start and end X forrest@0: yy1 = (/0.900,0.995/) ; Start and end Y forrest@0: text1 = table_header_name forrest@0: res1 = True forrest@0: res1@txFontHeightF = 0.03 forrest@0: res1@gsFillColor = "CornFlowerBlue" forrest@0: forrest@0: ; Column header (equally space in x2) forrest@0: ncr21 = (/1,ncol1/) ; 1 rows, 4 columns forrest@0: xx21 = (/xx1(1),0.995/) ; start from end of x1 forrest@0: yy21 = (/0.9475,0.995/) ; half of y1 forrest@0: text21 = col_header1 forrest@0: res21 = True forrest@0: res21@txFontHeightF = 0.015 forrest@0: res21@gsFillColor = "Gray" forrest@0: forrest@0: ncr22 = (/1,ncol2/) ; 1 rows, 12 columns forrest@0: xx22 = xx21 ; start from end of x1 forrest@0: yy22 = (/0.900,0.9475/) ; half of y1 forrest@0: text22 = col_header2 forrest@0: res22 = True forrest@0: res22@txFontHeightF = 0.015 forrest@0: res22@gsFillColor = "Gray" forrest@0: forrest@0: ; Row header (equally space in y2) forrest@0: ncr3 = (/nrow,1/) ; 20 rows, 1 columns forrest@0: xx3 = xx1 ; same as x1 forrest@0: yy3 = (/1.0-table_length,0.900/) ; end at start of y1 forrest@0: text3 = row_header forrest@0: res3 = True forrest@0: res3@txFontHeightF = 0.01 forrest@0: res3@gsFillColor = "Gray" forrest@0: forrest@0: ; Main table body forrest@0: ncr4 = (/nrow,ncol/) ; 5 rows, 5 columns forrest@0: xx4 = xx21 ; Start and end x forrest@0: yy4 = yy3 ; Start and end Y forrest@0: text4 = new((/nrow,ncol/),string) forrest@0: forrest@0: color_fill4 = new((/nrow,ncol/),string) forrest@0: color_fill4 = "white" forrest@0: color_fill4(nrow-1,:) = "CornFlowerBlue" forrest@0: forrest@0: res4 = True ; Set up resource list forrest@0: ; res4@gsnDebug = True ; Useful to print NDC row,col values used. forrest@0: res4@txFontHeightF = 0.015 forrest@0: res4@gsFillColor = color_fill4 forrest@0: forrest@0: delete (color_fill4) forrest@0: forrest@0: ;************************************************ forrest@0: ; plot global land class: observed forrest@0: ;************************************************ forrest@0: ;global res forrest@0: forrest@0: resg = True ; Use plot options forrest@0: resg@cnFillOn = True ; Turn on color fill forrest@0: resg@gsnSpreadColors = True ; use full colormap forrest@0: ; resg@cnFillMode = "RasterFill" ; Turn on raster color forrest@0: ; resg@lbLabelAutoStride = True forrest@0: resg@cnLinesOn = False ; Turn off contourn lines forrest@0: resg@mpFillOn = False ; Turn off map fill forrest@0: forrest@0: resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals forrest@0: resg@cnMinLevelValF = 1. ; Min level forrest@0: resg@cnMaxLevelValF = 19. ; Max level forrest@0: resg@cnLevelSpacingF = 1. ; interval forrest@0: forrest@0: ;global contour ob forrest@0: classob@_FillValue = 1.e+36 forrest@0: classob = where(classob.eq.0,classob@_FillValue,classob) forrest@0: forrest@0: plot_name = "global_class_ob" forrest@0: title = ob_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: plot = gsn_csm_contour_map_ce(wks,classob,resg) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"-*."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: ;******************************************************************* forrest@0: ; Calculate "nice" bins for binning the data in equally spaced ranges forrest@0: ;******************************************************************** forrest@0: nclassn = nclass + 1 forrest@0: range = fspan(0,nclassn-1,nclassn) forrest@0: ; print (range) forrest@0: forrest@0: ; Use this range information to grab all the values in a forrest@0: ; particular range, and then take an average. forrest@0: forrest@0: nr = dimsizes(range) forrest@0: nx = nr-1 forrest@0: xvalues = new((/2,nx/),float) forrest@0: xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2. forrest@0: dx = xvalues(0,1) - xvalues(0,0) ; range width forrest@0: dx4 = dx/4 ; 1/4 of the range forrest@0: xvalues(1,:) = xvalues(0,:) - dx/5. forrest@0: ;----------------------------------------------------------------- forrest@0: ;(A) mean forrest@0: ;-------------------------------------------------------------------- forrest@0: ; get data forrest@0: forrest@0: laiob_mean = dim_avg_Wrap(laiob(lat|:,lon|:,time|:)) forrest@0: laimod_mean = dim_avg_Wrap(laimod(lat|:,lon|:,time|:)) forrest@0: forrest@0: DATA11_1D = ndtooned(classob) forrest@0: DATA12_1D = ndtooned(laiob_mean) forrest@0: DATA22_1D = ndtooned(laimod_mean) forrest@0: forrest@0: yvalues = new((/2,nx/),float) forrest@0: mn_yvalues = new((/2,nx/),float) forrest@0: mx_yvalues = new((/2,nx/),float) forrest@0: forrest@0: do nd=0,1 forrest@0: forrest@0: ; See if we are doing model or observational data. forrest@0: forrest@0: if(nd.eq.0) then forrest@0: data_ob = DATA11_1D forrest@0: data_mod = DATA12_1D forrest@0: else forrest@0: data_ob = DATA11_1D forrest@0: data_mod = DATA22_1D forrest@0: end if forrest@0: forrest@0: ; Loop through each range and check for values. forrest@0: forrest@0: do i=0,nr-2 forrest@0: if (i.ne.(nr-2)) then forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+","+range(i+1)+")") forrest@0: idx = ind((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1))) forrest@0: else forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+",)") forrest@0: idx = ind(data_ob.ge.range(i)) forrest@0: end if forrest@0: forrest@0: ; Calculate average, and get min and max. forrest@0: forrest@0: if(.not.any(ismissing(idx))) then forrest@0: yvalues(nd,i) = avg(data_mod(idx)) forrest@0: mn_yvalues(nd,i) = min(data_mod(idx)) forrest@0: mx_yvalues(nd,i) = max(data_mod(idx)) forrest@0: count = dimsizes(idx) forrest@0: else forrest@0: count = 0 forrest@0: yvalues(nd,i) = yvalues@_FillValue forrest@0: mn_yvalues(nd,i) = yvalues@_FillValue forrest@0: mx_yvalues(nd,i) = yvalues@_FillValue forrest@0: end if forrest@0: forrest@0: ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) forrest@0: ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) forrest@0: forrest@0: ; Clean up for next time in loop. forrest@0: forrest@0: delete(idx) forrest@0: end do forrest@0: delete(data_ob) forrest@0: delete(data_mod) forrest@0: end do forrest@0: ;----------------------------------------------------------------- forrest@0: ; compute correlation coef and M score forrest@0: forrest@0: u = yvalues(0,:) forrest@0: v = yvalues(1,:) forrest@0: forrest@0: good = ind(.not.ismissing(u) .and. .not.ismissing(v)) forrest@0: uu = u(good) forrest@0: vv = v(good) forrest@0: forrest@0: ccrMean = esccr(uu,vv,0) forrest@0: ; print (ccrMean) forrest@0: forrest@0: ; new eq forrest@0: bias = sum(abs(vv-uu)/(vv+uu)) forrest@0: Mmean = (1.- (bias/dimsizes(uu)))*5. forrest@0: forrest@0: M_lai_mean = sprintf("%.2f", Mmean) forrest@0: system("sed s#M_lai_mean#"+M_lai_mean+"# "+html_name+" > "+html_new+";"+ \ forrest@0: "mv -f "+html_new+" "+html_name) forrest@0: print (M_lai_mean) forrest@0: forrest@0: do i=0,nrow-2 forrest@0: text4(i,0) = sprintf("%5.2f",u(i)) forrest@0: text4(i,1) = sprintf("%5.2f",v(i)) forrest@0: text4(i,2) = "-" forrest@0: end do forrest@0: text4(nrow-1,0) = sprintf("%5.2f",avg(u)) forrest@0: text4(nrow-1,1) = sprintf("%5.2f",avg(v)) forrest@0: text4(nrow-1,2) = sprintf("%5.2f",Mmean) forrest@0: forrest@0: delete (u) forrest@0: delete (v) forrest@0: delete (uu) forrest@0: delete (vv) forrest@0: ;-------------------------------------------------------------------- forrest@0: ; histogram res forrest@0: forrest@0: resm = True forrest@0: resm@gsnMaximize = True forrest@0: resm@gsnDraw = False forrest@0: resm@gsnFrame = False forrest@0: resm@xyMarkLineMode = "Markers" forrest@0: resm@xyMarkerSizeF = 0.014 forrest@0: resm@xyMarker = 16 forrest@0: resm@xyMarkerColors = (/"Brown","Blue"/) forrest@0: ; resm@trYMinF = min(mn_yvalues) - 10. forrest@0: ; resm@trYMaxF = max(mx_yvalues) + 10. forrest@0: resm@trYMinF = min(mn_yvalues) - 2 forrest@0: resm@trYMaxF = max(mx_yvalues) + 4 forrest@0: forrest@0: resm@tiYAxisString = "Mean LAI (Leaf Area Index)" forrest@0: resm@tiXAxisString = "Land Cover Type" forrest@0: forrest@0: max_bar = new((/2,nx/),graphic) forrest@0: min_bar = new((/2,nx/),graphic) forrest@0: max_cap = new((/2,nx/),graphic) forrest@0: min_cap = new((/2,nx/),graphic) forrest@0: forrest@0: lnres = True forrest@0: line_colors = (/"brown","blue"/) forrest@0: ;------------------------------------------------------------------ forrest@0: ; Start the graphics. forrest@0: forrest@0: plot_name = "histogram_mean" forrest@0: title = model_name + " vs Observed" forrest@0: resm@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: ;----------------------------- forrest@0: ; Add a boxed legend using the more simple method forrest@0: forrest@0: resm@pmLegendDisplayMode = "Always" forrest@0: ; resm@pmLegendWidthF = 0.1 forrest@0: resm@pmLegendWidthF = 0.08 forrest@0: resm@pmLegendHeightF = 0.05 forrest@0: resm@pmLegendOrthogonalPosF = -1.17 forrest@0: ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) forrest@0: ; resm@pmLegendParallelPosF = 0.18 forrest@0: resm@pmLegendParallelPosF = 0.88 ;(rightward) forrest@0: forrest@0: ; resm@lgPerimOn = False forrest@0: resm@lgLabelFontHeightF = 0.015 forrest@0: resm@xyExplicitLegendLabels = (/"observed",model_name/) forrest@0: ;----------------------------- forrest@0: tRes = True forrest@0: tRes@txFontHeightF = 0.025 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) forrest@0: forrest@0: xy = gsn_csm_xy(wks,xvalues,yvalues,resm) forrest@0: ;------------------------------- forrest@0: ;Attach the vertical bar and the horizontal cap line forrest@0: forrest@0: do nd=0,1 forrest@0: lnres@gsLineColor = line_colors(nd) forrest@0: do i=0,nx-1 forrest@0: forrest@0: if(.not.ismissing(mn_yvalues(nd,i)).and. \ forrest@0: .not.ismissing(mx_yvalues(nd,i))) then forrest@0: forrest@0: ; Attach the vertical bar, both above and below the marker. forrest@0: forrest@0: x1 = xvalues(nd,i) forrest@0: y1 = yvalues(nd,i) forrest@0: y2 = mn_yvalues(nd,i) forrest@0: min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: forrest@0: y2 = mx_yvalues(nd,i) forrest@0: max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: forrest@0: ; Attach the horizontal cap line, both above and below the marker. forrest@0: forrest@0: x1 = xvalues(nd,i) - dx4 forrest@0: x2 = xvalues(nd,i) + dx4 forrest@0: y1 = mn_yvalues(nd,i) forrest@0: min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: forrest@0: y1 = mx_yvalues(nd,i) forrest@0: max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: end if forrest@0: end do forrest@0: end do forrest@0: forrest@0: draw(xy) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: delete (DATA11_1D) forrest@0: delete (DATA12_1D) forrest@0: delete (DATA22_1D) forrest@0: ;delete (range) forrest@0: ;delete (xvalues) forrest@0: delete (yvalues) forrest@0: delete (mn_yvalues) forrest@0: delete (mx_yvalues) forrest@0: delete (good) forrest@0: delete (max_bar) forrest@0: delete (min_bar) forrest@0: delete (max_cap) forrest@0: delete (min_cap) forrest@0: ;----------------------------------------------------------------- forrest@0: ;global res forrest@0: forrest@0: resg = True ; Use plot options forrest@0: resg@cnFillOn = True ; Turn on color fill forrest@0: resg@gsnSpreadColors = True ; use full colormap forrest@0: ; resg@cnFillMode = "RasterFill" ; Turn on raster color forrest@0: ; resg@lbLabelAutoStride = True forrest@0: resg@cnLinesOn = False ; Turn off contourn lines forrest@0: resg@mpFillOn = False ; Turn off map fill forrest@0: forrest@0: resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals forrest@0: resg@cnMinLevelValF = 0. ; Min level forrest@0: resg@cnMaxLevelValF = 10. ; Max level forrest@0: resg@cnLevelSpacingF = 1. ; interval forrest@0: forrest@0: ;global contour ob forrest@0: forrest@0: delta = 0.00001 forrest@0: laiob_mean = where(ismissing(laiob_mean).and.(ismissing(laimod_mean).or.(laimod_mean.lt.delta)),0.,laiob_mean) forrest@0: forrest@0: plot_name = "global_mean_ob" forrest@0: title = ob_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: plot = gsn_csm_contour_map_ce(wks,laiob_mean,resg) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"-*."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: ;------------------------------------------------------------------------ forrest@0: ;global contour model forrest@0: forrest@0: plot_name = "global_mean_model" forrest@0: title = "Model " + model_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: plot = gsn_csm_contour_map_ce(wks,laimod_mean,resg) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"-*."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: ;----------------------------------------------------------------- forrest@0: ;(B) max forrest@0: ;-------------------------------------------------------------------- forrest@0: ; get data forrest@0: forrest@0: ; observed forrest@0: laiob_max = laiob(0,:,:) forrest@0: s = laiob(:,0,0) forrest@0: laiob_max@long_name = "Leaf Area Index Max" forrest@0: forrest@0: dsizes_z = dimsizes(laiob) forrest@0: nlat = dsizes_z(1) forrest@0: nlon = dsizes_z(2) forrest@0: forrest@0: do j = 0,nlat-1 forrest@0: do i = 0,nlon-1 forrest@0: s = laiob(:,j,i) forrest@0: laiob_max(j,i) = max(s) forrest@0: end do forrest@0: end do forrest@0: forrest@0: ; print (min(y)+"/"+max(y)) forrest@0: delete (s) forrest@0: delete (dsizes_z) forrest@0: ;------------------------- forrest@0: ; model forrest@0: laimod_max = laimod(0,:,:) forrest@0: s = laimod(:,0,0) forrest@0: laimod_max@long_name = "Leaf Area Index Max" forrest@0: forrest@0: dsizes_z = dimsizes(laimod) forrest@0: nlat = dsizes_z(1) forrest@0: nlon = dsizes_z(2) forrest@0: forrest@0: do j = 0,nlat-1 forrest@0: do i = 0,nlon-1 forrest@0: s = laimod(:,j,i) forrest@0: laimod_max(j,i) = max(s) forrest@0: end do forrest@0: end do forrest@0: forrest@0: ; print (min(laimod_max)+"/"+max(laimod_max)) forrest@0: delete (s) forrest@0: delete (dsizes_z) forrest@0: ;------------------------ forrest@0: DATA11_1D = ndtooned(classob) forrest@0: DATA12_1D = ndtooned(laiob_max) forrest@0: DATA22_1D = ndtooned(laimod_max) forrest@0: forrest@0: yvalues = new((/2,nx/),float) forrest@0: mn_yvalues = new((/2,nx/),float) forrest@0: mx_yvalues = new((/2,nx/),float) forrest@0: forrest@0: do nd=0,1 forrest@0: forrest@0: ; See if we are doing model or observational data. forrest@0: forrest@0: if(nd.eq.0) then forrest@0: data_ob = DATA11_1D forrest@0: data_mod = DATA12_1D forrest@0: else forrest@0: data_ob = DATA11_1D forrest@0: data_mod = DATA22_1D forrest@0: end if forrest@0: forrest@0: ; Loop through each range and check for values. forrest@0: forrest@0: do i=0,nr-2 forrest@0: if (i.ne.(nr-2)) then forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+","+range(i+1)+")") forrest@0: idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) forrest@0: else forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+",)") forrest@0: idx = ind(range(i).le.data_ob) forrest@0: end if forrest@0: forrest@0: ; Calculate average, and get min and max. forrest@0: forrest@0: if(.not.any(ismissing(idx))) then forrest@0: yvalues(nd,i) = avg(data_mod(idx)) forrest@0: mn_yvalues(nd,i) = min(data_mod(idx)) forrest@0: mx_yvalues(nd,i) = max(data_mod(idx)) forrest@0: count = dimsizes(idx) forrest@0: else forrest@0: count = 0 forrest@0: yvalues(nd,i) = yvalues@_FillValue forrest@0: mn_yvalues(nd,i) = yvalues@_FillValue forrest@0: mx_yvalues(nd,i) = yvalues@_FillValue forrest@0: end if forrest@0: forrest@0: ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) forrest@0: ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) forrest@0: forrest@0: ; Clean up for next time in loop. forrest@0: forrest@0: delete(idx) forrest@0: end do forrest@0: delete(data_ob) forrest@0: delete(data_mod) forrest@0: end do forrest@0: ;----------------------------------------------------------------- forrest@0: ; compute correlation coef and M score forrest@0: forrest@0: u = yvalues(0,:) forrest@0: v = yvalues(1,:) forrest@0: forrest@0: good = ind(.not.ismissing(u) .and. .not.ismissing(v)) forrest@0: uu = u(good) forrest@0: vv = v(good) forrest@0: forrest@0: ccrMax = esccr(uu,vv,0) forrest@0: ; print (ccrMax) forrest@0: forrest@0: ; new eq forrest@0: bias = sum(abs(vv-uu)/(vv+uu)) forrest@0: Mmax = (1.- (bias/dimsizes(uu)))*5. forrest@0: forrest@0: M_lai_max = sprintf("%.2f", Mmax) forrest@0: system("sed s#M_lai_max#"+M_lai_max+"# "+html_name+" > "+html_new+";"+ \ forrest@0: "mv -f "+html_new+" "+html_name) forrest@0: print (M_lai_max) forrest@0: forrest@0: do i=0,nrow-2 forrest@0: text4(i,3) = sprintf("%5.2f",u(i)) forrest@0: text4(i,4) = sprintf("%5.2f",v(i)) forrest@0: text4(i,5) = "-" forrest@0: end do forrest@0: text4(nrow-1,3) = sprintf("%5.2f",avg(u)) forrest@0: text4(nrow-1,4) = sprintf("%5.2f",avg(v)) forrest@0: text4(nrow-1,5) = sprintf("%5.2f",Mmax) forrest@0: forrest@0: delete (u) forrest@0: delete (v) forrest@0: delete (uu) forrest@0: delete (vv) forrest@0: ;-------------------------------------------------------------------- forrest@0: ; histogram res forrest@0: forrest@0: resm = True forrest@0: resm@gsnMaximize = True forrest@0: resm@gsnDraw = False forrest@0: resm@gsnFrame = False forrest@0: resm@xyMarkLineMode = "Markers" forrest@0: resm@xyMarkerSizeF = 0.014 forrest@0: resm@xyMarker = 16 forrest@0: resm@xyMarkerColors = (/"Brown","Blue"/) forrest@0: ; resm@trYMinF = min(mn_yvalues) - 10. forrest@0: ; resm@trYMaxF = max(mx_yvalues) + 10. forrest@0: resm@trYMinF = min(mn_yvalues) - 2 forrest@0: resm@trYMaxF = max(mx_yvalues) + 4 forrest@0: forrest@0: resm@tiYAxisString = "Max LAI (Leaf Area Index)" forrest@0: resm@tiXAxisString = "Land Cover Type" forrest@0: forrest@0: max_bar = new((/2,nx/),graphic) forrest@0: min_bar = new((/2,nx/),graphic) forrest@0: max_cap = new((/2,nx/),graphic) forrest@0: min_cap = new((/2,nx/),graphic) forrest@0: forrest@0: lnres = True forrest@0: line_colors = (/"brown","blue"/) forrest@0: ;------------------------------------------------------------------ forrest@0: ; Start the graphics. forrest@0: forrest@0: plot_name = "histogram_max" forrest@0: title = model_name + " vs Observed" forrest@0: resm@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: ;----------------------------- forrest@0: ; Add a boxed legend using the more simple method forrest@0: forrest@0: resm@pmLegendDisplayMode = "Always" forrest@0: ; resm@pmLegendWidthF = 0.1 forrest@0: resm@pmLegendWidthF = 0.08 forrest@0: resm@pmLegendHeightF = 0.05 forrest@0: resm@pmLegendOrthogonalPosF = -1.17 forrest@0: ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) forrest@0: ; resm@pmLegendParallelPosF = 0.18 forrest@0: resm@pmLegendParallelPosF = 0.88 ;(rightward) forrest@0: forrest@0: ; resm@lgPerimOn = False forrest@0: resm@lgLabelFontHeightF = 0.015 forrest@0: resm@xyExplicitLegendLabels = (/"observed",model_name/) forrest@0: ;----------------------------- forrest@0: tRes = True forrest@0: tRes@txFontHeightF = 0.025 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) forrest@0: forrest@0: xy = gsn_csm_xy(wks,xvalues,yvalues,resm) forrest@0: ;------------------------------- forrest@0: ;Attach the vertical bar and the horizontal cap line forrest@0: forrest@0: do nd=0,1 forrest@0: lnres@gsLineColor = line_colors(nd) forrest@0: do i=0,nx-1 forrest@0: forrest@0: if(.not.ismissing(mn_yvalues(nd,i)).and. \ forrest@0: .not.ismissing(mx_yvalues(nd,i))) then forrest@0: forrest@0: ; Attach the vertical bar, both above and below the marker. forrest@0: forrest@0: x1 = xvalues(nd,i) forrest@0: y1 = yvalues(nd,i) forrest@0: y2 = mn_yvalues(nd,i) forrest@0: min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: forrest@0: y2 = mx_yvalues(nd,i) forrest@0: max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: forrest@0: ; Attach the horizontal cap line, both above and below the marker. forrest@0: forrest@0: x1 = xvalues(nd,i) - dx4 forrest@0: x2 = xvalues(nd,i) + dx4 forrest@0: y1 = mn_yvalues(nd,i) forrest@0: min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: forrest@0: y1 = mx_yvalues(nd,i) forrest@0: max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: end if forrest@0: end do forrest@0: end do forrest@0: forrest@0: draw(xy) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: delete (DATA11_1D) forrest@0: delete (DATA12_1D) forrest@0: delete (DATA22_1D) forrest@0: ;delete (range) forrest@0: ;delete (xvalues) forrest@0: delete (yvalues) forrest@0: delete (mn_yvalues) forrest@0: delete (mx_yvalues) forrest@0: delete (good) forrest@0: delete (max_bar) forrest@0: delete (min_bar) forrest@0: delete (max_cap) forrest@0: delete (min_cap) forrest@0: ;----------------------------------------------------------------- forrest@0: ;global res forrest@0: forrest@0: resg = True ; Use plot options forrest@0: resg@cnFillOn = True ; Turn on color fill forrest@0: resg@gsnSpreadColors = True ; use full colormap forrest@0: ; resg@cnFillMode = "RasterFill" ; Turn on raster color forrest@0: ; resg@lbLabelAutoStride = True forrest@0: resg@cnLinesOn = False ; Turn off contourn lines forrest@0: resg@mpFillOn = False ; Turn off map fill forrest@0: forrest@0: resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals forrest@0: resg@cnMinLevelValF = 0. ; Min level forrest@0: resg@cnMaxLevelValF = 10. ; Max level forrest@0: resg@cnLevelSpacingF = 1. ; interval forrest@0: forrest@0: ;global contour ob forrest@0: forrest@0: delta = 0.00001 forrest@0: laiob_max = where(ismissing(laiob_max).and.(ismissing(laimod_max).or.(laimod_max.lt.delta)),0.,laiob_max) forrest@0: forrest@0: plot_name = "global_max_ob" forrest@0: title = ob_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: plot_max = gsn_csm_contour_map_ce(wks,laiob_max,resg) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"-*."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: ;------------------------------------------------------------------------ forrest@0: ;global contour model forrest@0: forrest@0: plot_name = "global_max_model" forrest@0: title = "Model " + model_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: plot_max = gsn_csm_contour_map_ce(wks,laimod_max,resg) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"-*."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: ;------------------------------------------------------------------------ forrest@0: ;(C) phase forrest@0: ;-------------------------------------------------------------------- forrest@0: ; get data forrest@0: forrest@0: ; observed forrest@0: laiob_phase = laiob(0,:,:) forrest@0: s = laiob(:,0,0) forrest@0: laiob_phase@long_name = "Leaf Area Index Max Month" forrest@0: forrest@0: dsizes_z = dimsizes(laiob) forrest@0: nlat = dsizes_z(1) forrest@0: nlon = dsizes_z(2) forrest@0: forrest@0: do j = 0,nlat-1 forrest@0: do i = 0,nlon-1 forrest@0: s = laiob(:,j,i) forrest@0: laiob_phase(j,i) = maxind(s) + 1 forrest@0: end do forrest@0: end do forrest@0: forrest@0: ; print (min(laiob_phase)+"/"+max(laiob_phase)) forrest@0: delete (s) forrest@0: delete (dsizes_z) forrest@0: ;------------------------- forrest@0: ; model forrest@0: laimod_phase = laimod(0,:,:) forrest@0: s = laimod(:,0,0) forrest@0: laimod_phase@long_name = "Leaf Area Index Max Month" forrest@0: forrest@0: dsizes_z = dimsizes(laimod) forrest@0: nlat = dsizes_z(1) forrest@0: nlon = dsizes_z(2) forrest@0: forrest@0: do j = 0,nlat-1 forrest@0: do i = 0,nlon-1 forrest@0: s = laimod(:,j,i) forrest@0: laimod_phase(j,i) = maxind(s) + 1 forrest@0: end do forrest@0: end do forrest@0: forrest@0: ; print (min(laimod_phase)+"/"+max(laimod_phase)) forrest@0: delete (s) forrest@0: delete (dsizes_z) forrest@0: ;------------------------ forrest@0: DATA11_1D = ndtooned(classob) forrest@0: DATA12_1D = ndtooned(laiob_phase) forrest@0: DATA22_1D = ndtooned(laimod_phase) forrest@0: forrest@0: yvalues = new((/2,nx/),float) forrest@0: mn_yvalues = new((/2,nx/),float) forrest@0: mx_yvalues = new((/2,nx/),float) forrest@0: forrest@0: do nd=0,1 forrest@0: forrest@0: ; See if we are doing model or observational data. forrest@0: forrest@0: if(nd.eq.0) then forrest@0: data_ob = DATA11_1D forrest@0: data_mod = DATA12_1D forrest@0: else forrest@0: data_ob = DATA11_1D forrest@0: data_mod = DATA22_1D forrest@0: end if forrest@0: forrest@0: ; Loop through each range and check for values. forrest@0: forrest@0: do i=0,nr-2 forrest@0: if (i.ne.(nr-2)) then forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+","+range(i+1)+")") forrest@0: idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) forrest@0: else forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+",)") forrest@0: idx = ind(range(i).le.data_ob) forrest@0: end if forrest@0: forrest@0: ; Calculate average, and get min and max. forrest@0: forrest@0: if(.not.any(ismissing(idx))) then forrest@0: yvalues(nd,i) = avg(data_mod(idx)) forrest@0: mn_yvalues(nd,i) = min(data_mod(idx)) forrest@0: mx_yvalues(nd,i) = max(data_mod(idx)) forrest@0: count = dimsizes(idx) forrest@0: else forrest@0: count = 0 forrest@0: yvalues(nd,i) = yvalues@_FillValue forrest@0: mn_yvalues(nd,i) = yvalues@_FillValue forrest@0: mx_yvalues(nd,i) = yvalues@_FillValue forrest@0: end if forrest@0: forrest@0: ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) forrest@0: ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) forrest@0: forrest@0: ; Clean up for next time in loop. forrest@0: forrest@0: delete(idx) forrest@0: end do forrest@0: delete(data_ob) forrest@0: delete(data_mod) forrest@0: end do forrest@0: ;----------------------------------------------------------------- forrest@0: ; compute correlation coef and M score forrest@0: forrest@0: u = yvalues(0,:) forrest@0: v = yvalues(1,:) forrest@0: forrest@0: good = ind(.not.ismissing(u) .and. .not.ismissing(v)) forrest@0: uu = u(good) forrest@0: vv = v(good) forrest@0: forrest@0: ccrPhase = esccr(uu,vv,0) forrest@0: ; print (ccrPhase) forrest@0: forrest@0: ; old eq forrest@0: ; bias = abs(avg(vv)-avg(uu)) forrest@0: ; new eq forrest@0: bias = avg(abs(vv-uu)) forrest@0: forrest@0: bias = where((bias.gt. 6.),12.-bias,bias) forrest@0: Mphase = ((6. - bias)/6.)*5. forrest@0: forrest@0: M_lai_phase = sprintf("%.2f", Mphase) forrest@0: system("sed s#M_lai_phase#"+M_lai_phase+"# "+html_name+" > "+html_new+";"+ \ forrest@0: "mv -f "+html_new+" "+html_name) forrest@0: print (M_lai_phase) forrest@0: forrest@0: do i=0,nrow-2 forrest@0: text4(i,6) = sprintf("%5.2f",u(i)) forrest@0: text4(i,7) = sprintf("%5.2f",v(i)) forrest@0: text4(i,8) = "-" forrest@0: end do forrest@0: text4(nrow-1,6) = sprintf("%5.2f",avg(u)) forrest@0: text4(nrow-1,7) = sprintf("%5.2f",avg(v)) forrest@0: text4(nrow-1,8) = sprintf("%5.2f",Mphase) forrest@0: forrest@0: delete (u) forrest@0: delete (v) forrest@0: delete (uu) forrest@0: delete (vv) forrest@0: ;-------------------------------------------------------------------- forrest@0: ; histogram res forrest@0: forrest@0: resm = True forrest@0: resm@gsnMaximize = True forrest@0: resm@gsnDraw = False forrest@0: resm@gsnFrame = False forrest@0: resm@xyMarkLineMode = "Markers" forrest@0: resm@xyMarkerSizeF = 0.014 forrest@0: resm@xyMarker = 16 forrest@0: resm@xyMarkerColors = (/"Brown","Blue"/) forrest@0: ; resm@trYMinF = min(mn_yvalues) - 10. forrest@0: ; resm@trYMaxF = max(mx_yvalues) + 10. forrest@0: resm@trYMinF = min(mn_yvalues) - 2 forrest@0: resm@trYMaxF = max(mx_yvalues) + 4 forrest@0: forrest@0: resm@tiYAxisString = "Max LAI (Leaf Area Index) Month" forrest@0: resm@tiXAxisString = "Land Cover Type" forrest@0: forrest@0: max_bar = new((/2,nx/),graphic) forrest@0: min_bar = new((/2,nx/),graphic) forrest@0: max_cap = new((/2,nx/),graphic) forrest@0: min_cap = new((/2,nx/),graphic) forrest@0: forrest@0: lnres = True forrest@0: line_colors = (/"brown","blue"/) forrest@0: ;------------------------------------------------------------------ forrest@0: ; Start the graphics. forrest@0: forrest@0: plot_name = "histogram_phase" forrest@0: title = model_name + " vs Observed" forrest@0: resm@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: ;----------------------------- forrest@0: ; Add a boxed legend using the more simple method forrest@0: forrest@0: resm@pmLegendDisplayMode = "Always" forrest@0: ; resm@pmLegendWidthF = 0.1 forrest@0: resm@pmLegendWidthF = 0.08 forrest@0: resm@pmLegendHeightF = 0.05 forrest@0: resm@pmLegendOrthogonalPosF = -1.17 forrest@0: ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) forrest@0: ; resm@pmLegendParallelPosF = 0.18 forrest@0: resm@pmLegendParallelPosF = 0.88 ;(rightward) forrest@0: forrest@0: ; resm@lgPerimOn = False forrest@0: resm@lgLabelFontHeightF = 0.015 forrest@0: resm@xyExplicitLegendLabels = (/"observed",model_name/) forrest@0: ;----------------------------- forrest@0: tRes = True forrest@0: tRes@txFontHeightF = 0.025 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) forrest@0: forrest@0: xy = gsn_csm_xy(wks,xvalues,yvalues,resm) forrest@0: ;------------------------------- forrest@0: ;Attach the vertical bar and the horizontal cap line forrest@0: forrest@0: do nd=0,1 forrest@0: lnres@gsLineColor = line_colors(nd) forrest@0: do i=0,nx-1 forrest@0: forrest@0: if(.not.ismissing(mn_yvalues(nd,i)).and. \ forrest@0: .not.ismissing(mx_yvalues(nd,i))) then forrest@0: forrest@0: ; Attach the vertical bar, both above and below the marker. forrest@0: forrest@0: x1 = xvalues(nd,i) forrest@0: y1 = yvalues(nd,i) forrest@0: y2 = mn_yvalues(nd,i) forrest@0: min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: forrest@0: y2 = mx_yvalues(nd,i) forrest@0: max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: forrest@0: ; Attach the horizontal cap line, both above and below the marker. forrest@0: forrest@0: x1 = xvalues(nd,i) - dx4 forrest@0: x2 = xvalues(nd,i) + dx4 forrest@0: y1 = mn_yvalues(nd,i) forrest@0: min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: forrest@0: y1 = mx_yvalues(nd,i) forrest@0: max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: end if forrest@0: end do forrest@0: end do forrest@0: forrest@0: draw(xy) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: delete (DATA11_1D) forrest@0: delete (DATA12_1D) forrest@0: delete (DATA22_1D) forrest@0: ;delete (range) forrest@0: ;delete (xvalues) forrest@0: delete (yvalues) forrest@0: delete (mn_yvalues) forrest@0: delete (mx_yvalues) forrest@0: delete (good) forrest@0: delete (max_bar) forrest@0: delete (min_bar) forrest@0: delete (max_cap) forrest@0: delete (min_cap) forrest@0: ;----------------------------------------------------------------- forrest@0: ;global res forrest@0: forrest@0: resg = True ; Use plot options forrest@0: resg@cnFillOn = True ; Turn on color fill forrest@0: resg@gsnSpreadColors = True ; use full colormap forrest@0: ; resg@cnFillMode = "RasterFill" ; Turn on raster color forrest@0: ; resg@lbLabelAutoStride = True forrest@0: resg@cnLinesOn = False ; Turn off contourn lines forrest@0: resg@mpFillOn = False ; Turn off map fill forrest@0: forrest@0: resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals forrest@0: resg@cnMinLevelValF = 1. ; Min level forrest@0: resg@cnMaxLevelValF = 12. ; Max level forrest@0: resg@cnLevelSpacingF = 1. ; interval forrest@0: forrest@0: ;global contour ob forrest@0: forrest@0: delta = 0.00001 forrest@0: laiob_phase = where(ismissing(laiob_phase).and.(ismissing(laimod_phase).or.(laimod_phase.lt.delta)),0.,laiob_phase) forrest@0: forrest@0: plot_name = "global_phase_ob" forrest@0: title = ob_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: plot = gsn_csm_contour_map_ce(wks,laiob_phase,resg) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"-*."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: ;------------------------------------------------------------------------ forrest@0: ;global contour model forrest@0: forrest@0: plot_name = "global_phase_model" forrest@0: title = "Model " + model_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: delete (plot) forrest@0: plot = gsn_csm_contour_map_ce(wks,laimod_phase,resg) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"-*."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: ;----------------------------------------------------------------- forrest@0: ;(D) grow day forrest@0: ;-------------------------------------------------------------------- forrest@0: ; get data forrest@0: forrest@0: day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/) forrest@0: forrest@0: ; observed forrest@0: laiob_grow = laiob(0,:,:) forrest@0: laiob_grow@long_name = "Days of Growing Season" forrest@0: forrest@0: dsizes_z = dimsizes(laiob) forrest@0: ntime = dsizes_z(0) forrest@0: nlat = dsizes_z(1) forrest@0: nlon = dsizes_z(2) forrest@0: forrest@0: do j = 0,nlat-1 forrest@0: do i = 0,nlon-1 forrest@0: nday = 0. forrest@0: do k = 0,ntime-1 forrest@0: if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then forrest@0: nday = nday + day_of_data(k) forrest@0: end if forrest@0: end do forrest@0: forrest@0: laiob_grow(j,i) = nday forrest@0: end do forrest@0: end do forrest@0: forrest@0: ; print (min(laiob_grow)+"/"+max(laiob_grow)) forrest@0: ;------------------------- forrest@0: ; model forrest@0: laimod_grow = laimod(0,:,:) forrest@0: laimod_grow@long_name = "Days of Growing Season" forrest@0: forrest@0: dsizes_z = dimsizes(laimod) forrest@0: ntime = dsizes_z(0) forrest@0: nlat = dsizes_z(1) forrest@0: nlon = dsizes_z(2) forrest@0: forrest@0: do j = 0,nlat-1 forrest@0: do i = 0,nlon-1 forrest@0: nday = 0. forrest@0: do k = 0,ntime-1 forrest@0: if (.not. ismissing(laimod(k,j,i)) .and. laimod(k,j,i) .gt. 1.0) then forrest@0: nday = nday + day_of_data(k) forrest@0: end if forrest@0: end do forrest@0: forrest@0: laimod_grow(j,i) = nday forrest@0: end do forrest@0: end do forrest@0: forrest@0: ; print (min(laimod_grow)+"/"+max(laimod_grow)) forrest@0: ;------------------------ forrest@0: DATA11_1D = ndtooned(classob) forrest@0: DATA12_1D = ndtooned(laiob_grow) forrest@0: DATA22_1D = ndtooned(laimod_grow) forrest@0: forrest@0: yvalues = new((/2,nx/),float) forrest@0: mn_yvalues = new((/2,nx/),float) forrest@0: mx_yvalues = new((/2,nx/),float) forrest@0: forrest@0: do nd=0,1 forrest@0: forrest@0: ; See if we are doing model or observational data. forrest@0: forrest@0: if(nd.eq.0) then forrest@0: data_ob = DATA11_1D forrest@0: data_mod = DATA12_1D forrest@0: else forrest@0: data_ob = DATA11_1D forrest@0: data_mod = DATA22_1D forrest@0: end if forrest@0: forrest@0: ; Loop through each range and check for values. forrest@0: forrest@0: do i=0,nr-2 forrest@0: if (i.ne.(nr-2)) then forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+","+range(i+1)+")") forrest@0: idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) forrest@0: else forrest@0: ; print("") forrest@0: ; print("In range ["+range(i)+",)") forrest@0: idx = ind(range(i).le.data_ob) forrest@0: end if forrest@0: forrest@0: ; Calculate average, and get min and max. forrest@0: forrest@0: if(.not.any(ismissing(idx))) then forrest@0: yvalues(nd,i) = avg(data_mod(idx)) forrest@0: mn_yvalues(nd,i) = min(data_mod(idx)) forrest@0: mx_yvalues(nd,i) = max(data_mod(idx)) forrest@0: count = dimsizes(idx) forrest@0: else forrest@0: count = 0 forrest@0: yvalues(nd,i) = yvalues@_FillValue forrest@0: mn_yvalues(nd,i) = yvalues@_FillValue forrest@0: mx_yvalues(nd,i) = yvalues@_FillValue forrest@0: end if forrest@0: forrest@0: ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) forrest@0: ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) forrest@0: forrest@0: ; Clean up for next time in loop. forrest@0: forrest@0: delete(idx) forrest@0: end do forrest@0: delete(data_ob) forrest@0: delete(data_mod) forrest@0: end do forrest@0: ;----------------------------------------------------------------- forrest@0: ; compute correlation coef and M score forrest@0: forrest@0: u = yvalues(0,:) forrest@0: v = yvalues(1,:) forrest@0: forrest@0: good = ind(.not.ismissing(u) .and. .not.ismissing(v)) forrest@0: uu = u(good) forrest@0: vv = v(good) forrest@0: forrest@0: ccrGrow = esccr(uu,vv,0) forrest@0: ; print (ccrGrow) forrest@0: forrest@0: ; new eq forrest@0: bias = sum(abs(vv-uu)/(vv+uu)) forrest@0: Mgrow = (1.- (bias/dimsizes(uu)))*5. forrest@0: forrest@0: M_lai_grow = sprintf("%.2f", Mgrow) forrest@0: system("sed s#M_lai_grow#"+M_lai_grow+"# "+html_name+" > "+html_new+";"+ \ forrest@0: "mv -f "+html_new+" "+html_name) forrest@0: print (M_lai_grow) forrest@0: forrest@0: do i=0,nrow-2 forrest@0: text4(i,9) = sprintf("%5.2f",u(i)) forrest@0: text4(i,10) = sprintf("%5.2f",v(i)) forrest@0: text4(i,11) = "-" forrest@0: end do forrest@0: text4(nrow-1,9) = sprintf("%5.2f",avg(u)) forrest@0: text4(nrow-1,10) = sprintf("%5.2f",avg(v)) forrest@0: text4(nrow-1,11) = sprintf("%5.2f",Mgrow) forrest@0: forrest@0: delete (u) forrest@0: delete (v) forrest@0: delete (uu) forrest@0: delete (vv) forrest@0: ;-------------------------------------------------------------------- forrest@0: ; histogram res forrest@0: forrest@0: resm = True forrest@0: resm@gsnMaximize = True forrest@0: resm@gsnDraw = False forrest@0: resm@gsnFrame = False forrest@0: resm@xyMarkLineMode = "Markers" forrest@0: resm@xyMarkerSizeF = 0.014 forrest@0: resm@xyMarker = 16 forrest@0: resm@xyMarkerColors = (/"Brown","Blue"/) forrest@0: ; resm@trYMinF = min(mn_yvalues) - 10. forrest@0: ; resm@trYMaxF = max(mx_yvalues) + 10. forrest@0: resm@trYMinF = min(mn_yvalues) - 2 forrest@0: resm@trYMaxF = max(mx_yvalues) + 4 forrest@0: forrest@0: resm@tiYAxisString = "Days of Growing season" forrest@0: resm@tiXAxisString = "Land Cover Type" forrest@0: forrest@0: max_bar = new((/2,nx/),graphic) forrest@0: min_bar = new((/2,nx/),graphic) forrest@0: max_cap = new((/2,nx/),graphic) forrest@0: min_cap = new((/2,nx/),graphic) forrest@0: forrest@0: lnres = True forrest@0: line_colors = (/"brown","blue"/) forrest@0: ;------------------------------------------------------------------ forrest@0: ; Start the graphics. forrest@0: forrest@0: plot_name = "histogram_grow" forrest@0: title = model_name + " vs Observed" forrest@0: resm@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: ;----------------------------- forrest@0: ; Add a boxed legend using the more simple method forrest@0: forrest@0: resm@pmLegendDisplayMode = "Always" forrest@0: ; resm@pmLegendWidthF = 0.1 forrest@0: resm@pmLegendWidthF = 0.08 forrest@0: resm@pmLegendHeightF = 0.05 forrest@0: resm@pmLegendOrthogonalPosF = -1.17 forrest@0: ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) forrest@0: ; resm@pmLegendParallelPosF = 0.18 forrest@0: resm@pmLegendParallelPosF = 0.88 ;(rightward) forrest@0: forrest@0: ; resm@lgPerimOn = False forrest@0: resm@lgLabelFontHeightF = 0.015 forrest@0: resm@xyExplicitLegendLabels = (/"observed",model_name/) forrest@0: ;----------------------------- forrest@0: tRes = True forrest@0: tRes@txFontHeightF = 0.025 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) forrest@0: forrest@0: xy = gsn_csm_xy(wks,xvalues,yvalues,resm) forrest@0: ;------------------------------- forrest@0: ;Attach the vertical bar and the horizontal cap line forrest@0: forrest@0: do nd=0,1 forrest@0: lnres@gsLineColor = line_colors(nd) forrest@0: do i=0,nx-1 forrest@0: forrest@0: if(.not.ismissing(mn_yvalues(nd,i)).and. \ forrest@0: .not.ismissing(mx_yvalues(nd,i))) then forrest@0: forrest@0: ; Attach the vertical bar, both above and below the marker. forrest@0: forrest@0: x1 = xvalues(nd,i) forrest@0: y1 = yvalues(nd,i) forrest@0: y2 = mn_yvalues(nd,i) forrest@0: min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: forrest@0: y2 = mx_yvalues(nd,i) forrest@0: max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) forrest@0: forrest@0: ; Attach the horizontal cap line, both above and below the marker. forrest@0: forrest@0: x1 = xvalues(nd,i) - dx4 forrest@0: x2 = xvalues(nd,i) + dx4 forrest@0: y1 = mn_yvalues(nd,i) forrest@0: min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: forrest@0: y1 = mx_yvalues(nd,i) forrest@0: max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) forrest@0: end if forrest@0: end do forrest@0: end do forrest@0: forrest@0: draw(xy) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: forrest@0: delete (DATA11_1D) forrest@0: delete (DATA12_1D) forrest@0: delete (DATA22_1D) forrest@0: ;delete (range) forrest@0: ;delete (xvalues) forrest@0: delete (yvalues) forrest@0: delete (mn_yvalues) forrest@0: delete (mx_yvalues) forrest@0: delete (good) forrest@0: delete (max_bar) forrest@0: delete (min_bar) forrest@0: delete (max_cap) forrest@0: delete (min_cap) forrest@0: ;----------------------------------------------------------------- forrest@0: ;global res forrest@0: forrest@0: resg = True ; Use plot options forrest@0: resg@cnFillOn = True ; Turn on color fill forrest@0: resg@gsnSpreadColors = True ; use full colormap forrest@0: ; resg@cnFillMode = "RasterFill" ; Turn on raster color forrest@0: ; resg@lbLabelAutoStride = True forrest@0: resg@cnLinesOn = False ; Turn off contourn lines forrest@0: resg@mpFillOn = False ; Turn off map fill forrest@0: forrest@0: resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals forrest@0: resg@cnMinLevelValF = 60. ; Min level forrest@0: resg@cnMaxLevelValF = 360. ; Max level forrest@0: resg@cnLevelSpacingF = 20. ; interval forrest@0: forrest@0: ;global contour ob forrest@0: forrest@0: laiob_grow@_FillValue = 1.e+36 forrest@0: laiob_grow = where(laiob_grow .lt. 10.,laiob_grow@_FillValue,laiob_grow) forrest@0: forrest@0: plot_name = "global_grow_ob" forrest@0: title = ob_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: plot = gsn_csm_contour_map_ce(wks,laiob_grow,resg) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"-*."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: ;------------------------------------------------------------------------ forrest@0: ;global contour model forrest@0: forrest@0: laimod_grow@_FillValue = 1.e+36 forrest@0: laimod_grow = where(laimod_grow .lt. 10.,laimod_grow@_FillValue,laimod_grow) forrest@0: forrest@0: plot_name = "global_grow_model" forrest@0: title = "Model " + model_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: delete (plot) forrest@0: plot = gsn_csm_contour_map_ce(wks,laimod_grow,resg) forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"-*."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: ;------------------------------------------------------------------------ forrest@0: ;global contour model vs ob forrest@0: forrest@0: resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals forrest@0: resg@cnMinLevelValF = 0. ; Min level forrest@0: resg@cnMaxLevelValF = 10. ; Max level forrest@0: resg@cnLevelSpacingF = 1. ; interval forrest@0: forrest@0: plot_name = "global_mean_model_vs_ob" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: delete (plot) forrest@0: plot=new(3,graphic) ; create graphic array forrest@0: forrest@0: resg@gsnFrame = False ; Do not draw plot forrest@0: resg@gsnDraw = False ; Do not advance frame forrest@0: forrest@0: ; plot correlation coef forrest@0: forrest@0: gRes = True forrest@0: gRes@txFontHeightF = 0.02 forrest@0: gRes@txAngleF = 90 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) forrest@0: ;-------------------------------------------------------------------- forrest@0: forrest@0: ;(a) ob forrest@0: forrest@0: title = ob_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(0) = gsn_csm_contour_map_ce(wks,laiob_mean,resg) forrest@0: forrest@0: ;(b) model forrest@0: forrest@0: title = "Model "+ model_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(1) = gsn_csm_contour_map_ce(wks,laimod_mean,resg) forrest@0: forrest@0: ;(c) model-ob forrest@0: forrest@0: zz = laimod_mean forrest@0: zz = laimod_mean - laiob_mean forrest@0: title = "Model_"+model_name+" - Observed" forrest@0: resg@tiMainString = title forrest@0: forrest@0: resg@cnMinLevelValF = -2. ; Min level forrest@0: resg@cnMaxLevelValF = 2. ; Max level forrest@0: resg@cnLevelSpacingF = 0.4 ; interval forrest@0: forrest@0: forrest@0: plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) forrest@0: forrest@0: pres = True ; panel plot mods desired forrest@0: ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around forrest@0: ; indiv. plots in panel forrest@0: pres@gsnMaximize = True ; fill the page forrest@0: forrest@0: gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: delete (plot) forrest@0: ;----------------------------------------------------------------- forrest@0: ;global contour model vs ob forrest@0: forrest@0: resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals forrest@0: resg@cnMinLevelValF = 0. ; Min level forrest@0: resg@cnMaxLevelValF = 10. ; Max level forrest@0: resg@cnLevelSpacingF = 1. ; interval forrest@0: forrest@0: plot_name = "global_max_model_vs_ob" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: ; delete(plot) forrest@0: plot=new(3,graphic) ; create graphic array forrest@0: forrest@0: resg@gsnFrame = False ; Do not draw plot forrest@0: resg@gsnDraw = False ; Do not advance frame forrest@0: forrest@0: ; plot correlation coef forrest@0: forrest@0: gRes = True forrest@0: gRes@txFontHeightF = 0.02 forrest@0: gRes@txAngleF = 90 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) forrest@0: ;-------------------------------------------------------------------- forrest@0: ;(a) ob forrest@0: forrest@0: title = ob_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(0) = gsn_csm_contour_map_ce(wks,laiob_max,resg) forrest@0: forrest@0: ;(b) model forrest@0: forrest@0: title = "Model "+ model_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(1) = gsn_csm_contour_map_ce(wks,laimod_max,resg) forrest@0: forrest@0: ;(c) model-ob forrest@0: forrest@0: delete (zz) forrest@0: zz = laimod_max forrest@0: zz = laimod_max - laiob_max forrest@0: title = "Model_"+model_name+" - Observed" forrest@0: resg@tiMainString = title forrest@0: forrest@0: resg@cnMinLevelValF = -6. ; Min level forrest@0: resg@cnMaxLevelValF = 6. ; Max level forrest@0: resg@cnLevelSpacingF = 1. ; interval forrest@0: forrest@0: forrest@0: plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) forrest@0: forrest@0: pres = True ; panel plot mods desired forrest@0: ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around forrest@0: ; indiv. plots in panel forrest@0: pres@gsnMaximize = True ; fill the page forrest@0: forrest@0: gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: delete (plot) forrest@0: ;----------------------------------------------------------------- forrest@0: ;global contour model vs ob forrest@0: forrest@0: resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals forrest@0: resg@cnMinLevelValF = 1. ; Min level forrest@0: resg@cnMaxLevelValF = 12. ; Max level forrest@0: resg@cnLevelSpacingF = 1. ; interval forrest@0: forrest@0: plot_name = "global_phase_model_vs_ob" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: ; delete (plot) forrest@0: plot=new(3,graphic) ; create graphic array forrest@0: forrest@0: resg@gsnFrame = False ; Do not draw plot forrest@0: resg@gsnDraw = False ; Do not advance frame forrest@0: forrest@0: ; plot correlation coef forrest@0: forrest@0: gRes = True forrest@0: gRes@txFontHeightF = 0.02 forrest@0: gRes@txAngleF = 90 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) forrest@0: ;-------------------------------------------------------------------- forrest@0: ;(a) ob forrest@0: forrest@0: title = ob_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(0) = gsn_csm_contour_map_ce(wks,laiob_phase,resg) forrest@0: forrest@0: ;(b) model forrest@0: forrest@0: title = "Model "+ model_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(1) = gsn_csm_contour_map_ce(wks,laimod_phase,resg) forrest@0: forrest@0: ;(c) model-ob forrest@0: forrest@0: delete (zz) forrest@0: zz = laimod_phase forrest@0: zz = laimod_phase - laiob_phase forrest@0: title = "Model_"+model_name+" - Observed" forrest@0: resg@tiMainString = title forrest@0: forrest@0: resg@cnMinLevelValF = -6. ; Min level forrest@0: resg@cnMaxLevelValF = 6. ; Max level forrest@0: resg@cnLevelSpacingF = 1. ; interval forrest@0: forrest@0: forrest@0: plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) forrest@0: forrest@0: pres = True ; panel plot mods desired forrest@0: ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around forrest@0: ; indiv. plots in panel forrest@0: pres@gsnMaximize = True ; fill the page forrest@0: forrest@0: gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: delete (plot) forrest@0: ;------------------------------------------------------------------ forrest@0: ;global contour model vs ob forrest@0: forrest@0: resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals forrest@0: resg@cnMinLevelValF = 60. ; Min level forrest@0: resg@cnMaxLevelValF = 360. ; Max level forrest@0: resg@cnLevelSpacingF = 20. ; interval forrest@0: forrest@0: plot_name = "global_grow_model_vs_ob" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: ; delete (plot) forrest@0: plot=new(3,graphic) ; create graphic array forrest@0: forrest@0: resg@gsnFrame = False ; Do not draw plot forrest@0: resg@gsnDraw = False ; Do not advance frame forrest@0: forrest@0: ; plot correlation coef forrest@0: forrest@0: gRes = True forrest@0: gRes@txFontHeightF = 0.02 forrest@0: gRes@txAngleF = 90 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) forrest@0: ;-------------------------------------------------------------------- forrest@0: ;(a) ob forrest@0: forrest@0: title = ob_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(0) = gsn_csm_contour_map_ce(wks,laiob_grow,resg) forrest@0: forrest@0: ;(b) model forrest@0: forrest@0: title = "Model "+ model_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(1) = gsn_csm_contour_map_ce(wks,laimod_grow,resg) forrest@0: forrest@0: ;(c) model-ob forrest@0: forrest@0: delete (zz) forrest@0: zz = laimod_grow forrest@0: zz = laimod_grow - laiob_grow forrest@0: title = "Model_"+model_name+" - Observed" forrest@0: resg@tiMainString = title forrest@0: forrest@0: resg@cnMinLevelValF = -100. ; Min level forrest@0: resg@cnMaxLevelValF = 100. ; Max level forrest@0: resg@cnLevelSpacingF = 20. ; interval forrest@0: forrest@0: plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) forrest@0: forrest@0: pres = True ; panel plot mods desired forrest@0: ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around forrest@0: ; indiv. plots in panel forrest@0: pres@gsnMaximize = True ; fill the page forrest@0: forrest@0: gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: delete (plot) forrest@0: ;************************************************** forrest@0: ; plot lai table forrest@0: ;************************************************** forrest@0: forrest@0: plot_name = "table_lai" forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: forrest@0: gsn_table(wks,ncr1,xx1,yy1,text1,res1) forrest@0: gsn_table(wks,ncr21,xx21,yy21,text21,res21) forrest@0: gsn_table(wks,ncr22,xx22,yy22,text22,res22) forrest@0: gsn_table(wks,ncr3,xx3,yy3,text3,res3) forrest@0: gsn_table(wks,ncr4,xx4,yy4,text4,res4) forrest@0: forrest@0: frame(wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: ;------------------------------------------------------------------- forrest@0: temp_name = "lai." + model_name forrest@0: system("mkdir -p " + temp_name+";"+ \ forrest@0: "cp "+ html_name + " " +temp_name+"/table.html"+";"+ \ forrest@0: "mv *.png " + temp_name +";"+ \ forrest@0: "tar cf "+ temp_name +".tar " + temp_name) forrest@0: ;------------------------------------------------------------------- forrest@0: end forrest@0: