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" forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 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 = 18 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: forrest@0: model_name = "b30.061n" forrest@0: model_grid = "T31" forrest@0: forrest@0: dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" forrest@0: film = "b30.061n_1995-2004_MONS_climo_lnd.nc" forrest@0: ;film = "i01.03cn_1545-1569_MONS_climo.nc" 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: ; 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: 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: print (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: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) 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 = 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: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) 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: delete (plot) forrest@0: plot = 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: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: ;------------------------------------------------------------------------ forrest@0: ;global contour model vs ob 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: 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: 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: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: frame (wks) forrest@0: clear (wks) forrest@0: forrest@0: delete (plot) 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: print (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: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) 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: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) 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: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: ;------------------------------------------------------------------------ forrest@0: ;global contour model vs ob 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: 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: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: frame (wks) forrest@0: clear (wks) forrest@0: forrest@0: delete (plot) 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: print (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: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) 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 = 360. ; Max level forrest@0: resg@cnLevelSpacingF = 30. ; interval forrest@0: forrest@0: ;global contour ob forrest@0: forrest@0: delta = 0.00001 forrest@0: laiob_grow = where(ismissing(laiob_grow).and.(ismissing(laimod_grow).or.(laimod_grow.lt.delta)),0.,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: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: ;------------------------------------------------------------------------ forrest@0: ;global contour model 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: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: clear (wks) forrest@0: ;------------------------------------------------------------------------ forrest@0: ;global contour model vs ob 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: 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 = -120. ; Min level forrest@0: resg@cnMaxLevelValF = 120. ; Max level forrest@0: resg@cnLevelSpacingF = 20. ; 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: ; system("rm "+plot_name+"."+plot_type) forrest@0: ; system("rm "+plot_name+"-1."+plot_type_new) forrest@0: ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) forrest@0: forrest@0: frame (wks) forrest@0: clear (wks) forrest@0: forrest@0: delete (plot) forrest@0: end forrest@0: