diff -r 000000000000 -r 0c6405ab2ff4 lai/99.ncl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lai/99.ncl Mon Jan 26 22:08:20 2009 -0500 @@ -0,0 +1,1185 @@ +;******************************************************** +; histogram normalized by rain and compute correleration +;******************************************************** +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" + +procedure pminmax(data:numeric,name:string) +begin + print ("min/max " + name + " = " + min(data) + "/" + max(data)) + if(isatt(data,"units")) then + print (name + " units = " + data@units) + end if +end + +; Main code. +begin + +;nclass = 18 + nclass = 20 + + plot_type = "ps" + plot_type_new = "png" + +;************************************************ +; read in data: model +;************************************************ + + model_name = "b30.061n" + model_grid = "T31" + + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" + film = "b30.061n_1995-2004_MONS_climo_lnd.nc" +;film = "i01.03cn_1545-1569_MONS_climo.nc" + fm = addfile(dirm+film,"r") + + laimod = fm->TLAI + +;************************************************ +; read in data: observed +;************************************************ + + ob_name = "MODIS MOD 15A2 2000-2005" + + diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/" + filo1 = "land_class_"+model_grid+".nc" + filo2 = "LAI_2000-2005_MONS_"+model_grid+".nc" + + fo1 = addfile(diro+filo1,"r") + fo2 = addfile(diro+filo2,"r") + + classob = tofloat(fo1->LAND_CLASS) + laiob = fo2->LAI +;******************************************************************* +; Calculate "nice" bins for binning the data in equally spaced ranges +;******************************************************************** + nclassn = nclass + 1 + range = fspan(0,nclassn-1,nclassn) +; print (range) + +; Use this range information to grab all the values in a +; particular range, and then take an average. + + nr = dimsizes(range) + nx = nr-1 + xvalues = new((/2,nx/),float) + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2. + dx = xvalues(0,1) - xvalues(0,0) ; range width + dx4 = dx/4 ; 1/4 of the range + xvalues(1,:) = xvalues(0,:) - dx/5. +;----------------------------------------------------------------- + +;----------------------------------------------------------------- +;(B) max +;-------------------------------------------------------------------- +; get data + +; observed + laiob_max = laiob(0,:,:) + s = laiob(:,0,0) + laiob_max@long_name = "Leaf Area Index Max" + + dsizes_z = dimsizes(laiob) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + s = laiob(:,j,i) + laiob_max(j,i) = max(s) + end do + end do + +; print (min(y)+"/"+max(y)) + delete (s) + delete (dsizes_z) +;------------------------- +; model + laimod_max = laimod(0,:,:) + s = laimod(:,0,0) + laimod_max@long_name = "Leaf Area Index Max" + + dsizes_z = dimsizes(laimod) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + s = laimod(:,j,i) + laimod_max(j,i) = max(s) + end do + end do + +; print (min(laimod_max)+"/"+max(laimod_max)) + delete (s) + delete (dsizes_z) +;------------------------ + DATA11_1D = ndtooned(classob) + DATA12_1D = ndtooned(laiob_max) + DATA22_1D = ndtooned(laimod_max) + + yvalues = new((/2,nx/),float) + mn_yvalues = new((/2,nx/),float) + mx_yvalues = new((/2,nx/),float) + + do nd=0,1 + +; See if we are doing model or observational data. + + if(nd.eq.0) then + data_ob = DATA11_1D + data_mod = DATA12_1D + else + data_ob = DATA11_1D + data_mod = DATA22_1D + end if + +; Loop through each range and check for values. + + do i=0,nr-2 + if (i.ne.(nr-2)) then +; print("") +; print("In range ["+range(i)+","+range(i+1)+")") + idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) + else +; print("") +; print("In range ["+range(i)+",)") + idx = ind(range(i).le.data_ob) + end if + +; Calculate average, and get min and max. + + if(.not.any(ismissing(idx))) then + yvalues(nd,i) = avg(data_mod(idx)) + mn_yvalues(nd,i) = min(data_mod(idx)) + mx_yvalues(nd,i) = max(data_mod(idx)) + count = dimsizes(idx) + else + count = 0 + yvalues(nd,i) = yvalues@_FillValue + mn_yvalues(nd,i) = yvalues@_FillValue + mx_yvalues(nd,i) = yvalues@_FillValue + end if + +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) + +; Clean up for next time in loop. + + delete(idx) + end do + delete(data_ob) + delete(data_mod) + end do +;----------------------------------------------------------------- +; compute correlation coef and M score + + u = yvalues(0,:) + v = yvalues(1,:) + + good = ind(.not.ismissing(u) .and. .not.ismissing(v)) + uu = u(good) + vv = v(good) + + ccrMax = esccr(uu,vv,0) +; print (ccrMax) + +; new eq + bias = sum(abs(vv-uu)/(vv+uu)) + Mmax = (1.- (bias/dimsizes(uu)))*5. + + print (Mmax) + + delete (u) + delete (v) + delete (uu) + delete (vv) +;-------------------------------------------------------------------- +; histogram res + + resm = True + resm@gsnMaximize = True + resm@gsnDraw = False + resm@gsnFrame = False + resm@xyMarkLineMode = "Markers" + resm@xyMarkerSizeF = 0.014 + resm@xyMarker = 16 + resm@xyMarkerColors = (/"Brown","Blue"/) +; resm@trYMinF = min(mn_yvalues) - 10. +; resm@trYMaxF = max(mx_yvalues) + 10. + resm@trYMinF = min(mn_yvalues) - 2 + resm@trYMaxF = max(mx_yvalues) + 4 + + resm@tiYAxisString = "Max LAI (Leaf Area Index)" + resm@tiXAxisString = "Land Cover Type" + + max_bar = new((/2,nx/),graphic) + min_bar = new((/2,nx/),graphic) + max_cap = new((/2,nx/),graphic) + min_cap = new((/2,nx/),graphic) + + lnres = True + line_colors = (/"brown","blue"/) +;------------------------------------------------------------------ +; Start the graphics. + + plot_name = "histogram_max" + title = model_name + " vs Observed" + resm@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) +;----------------------------- +; Add a boxed legend using the more simple method + + resm@pmLegendDisplayMode = "Always" +; resm@pmLegendWidthF = 0.1 + resm@pmLegendWidthF = 0.08 + resm@pmLegendHeightF = 0.05 + resm@pmLegendOrthogonalPosF = -1.17 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) +; resm@pmLegendParallelPosF = 0.18 + resm@pmLegendParallelPosF = 0.88 ;(rightward) + +; resm@lgPerimOn = False + resm@lgLabelFontHeightF = 0.015 + resm@xyExplicitLegendLabels = (/"observed",model_name/) +;----------------------------- + tRes = True + tRes@txFontHeightF = 0.025 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")" + + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) + + xy = gsn_csm_xy(wks,xvalues,yvalues,resm) +;------------------------------- +;Attach the vertical bar and the horizontal cap line + + do nd=0,1 + lnres@gsLineColor = line_colors(nd) + do i=0,nx-1 + + if(.not.ismissing(mn_yvalues(nd,i)).and. \ + .not.ismissing(mx_yvalues(nd,i))) then + +; Attach the vertical bar, both above and below the marker. + + x1 = xvalues(nd,i) + y1 = yvalues(nd,i) + y2 = mn_yvalues(nd,i) + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + + y2 = mx_yvalues(nd,i) + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + +; Attach the horizontal cap line, both above and below the marker. + + x1 = xvalues(nd,i) - dx4 + x2 = xvalues(nd,i) + dx4 + y1 = mn_yvalues(nd,i) + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + + y1 = mx_yvalues(nd,i) + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + end if + end do + end do + + draw(xy) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +; system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) + + delete (DATA11_1D) + delete (DATA12_1D) + delete (DATA22_1D) +;delete (range) +;delete (xvalues) + delete (yvalues) + delete (mn_yvalues) + delete (mx_yvalues) + delete (good) + delete (max_bar) + delete (min_bar) + delete (max_cap) + delete (min_cap) +;----------------------------------------------------------------- +;global res + + resg = True ; Use plot options + resg@cnFillOn = True ; Turn on color fill + resg@gsnSpreadColors = True ; use full colormap +; resg@cnFillMode = "RasterFill" ; Turn on raster color +; resg@lbLabelAutoStride = True + resg@cnLinesOn = False ; Turn off contourn lines + resg@mpFillOn = False ; Turn off map fill + + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 0. ; Min level + resg@cnMaxLevelValF = 10. ; Max level + resg@cnLevelSpacingF = 1. ; interval + +;global contour ob + + delta = 0.00001 + laiob_max = where(ismissing(laiob_max).and.(ismissing(laimod_max).or.(laimod_max.lt.delta)),0.,laiob_max) + + plot_name = "global_max_ob" + title = ob_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot = gsn_csm_contour_map_ce(wks,laiob_max,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +; system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;global contour model + + plot_name = "global_max_model" + title = "Model " + model_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) + gsn_define_colormap(wks,"gui_default") ; choose colormap + + delete (plot) + plot = gsn_csm_contour_map_ce(wks,laimod_max,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +; system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;global contour model vs ob + + plot_name = "global_max_model_vs_ob" + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + delete (plot) + plot=new(3,graphic) ; create graphic array + + resg@gsnFrame = False ; Do not draw plot + resg@gsnDraw = False ; Do not advance frame + +; plot correlation coef + + gRes = True + gRes@txFontHeightF = 0.02 + gRes@txAngleF = 90 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")" + + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) +;-------------------------------------------------------------------- + +;(a) ob + + title = ob_name + resg@tiMainString = title + + plot(0) = gsn_csm_contour_map_ce(wks,laiob_max,resg) + +;(b) model + + title = "Model "+ model_name + resg@tiMainString = title + + plot(1) = gsn_csm_contour_map_ce(wks,laimod_max,resg) + +;(c) model-ob + + zz = laimod_max + zz = laimod_max - laiob_max + title = "Model_"+model_name+" - Observed" + resg@tiMainString = title + + resg@cnMinLevelValF = -6. ; Min level + resg@cnMaxLevelValF = 6. ; Max level + resg@cnLevelSpacingF = 1. ; interval + + + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) + + pres = True ; panel plot mods desired + pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around + ; indiv. plots in panel + pres@gsnMaximize = True ; fill the page + + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +; system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + frame (wks) + clear (wks) + + delete (plot) +;----------------------------------------------------------------- +;(C) phase +;-------------------------------------------------------------------- +; get data + +; observed + laiob_phase = laiob(0,:,:) + s = laiob(:,0,0) + laiob_phase@long_name = "Leaf Area Index Max Month" + + dsizes_z = dimsizes(laiob) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + s = laiob(:,j,i) + laiob_phase(j,i) = maxind(s) + 1 + end do + end do + +; print (min(laiob_phase)+"/"+max(laiob_phase)) + delete (s) + delete (dsizes_z) +;------------------------- +; model + laimod_phase = laimod(0,:,:) + s = laimod(:,0,0) + laimod_phase@long_name = "Leaf Area Index Max Month" + + dsizes_z = dimsizes(laimod) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + s = laimod(:,j,i) + laimod_phase(j,i) = maxind(s) + 1 + end do + end do + +; print (min(laimod_phase)+"/"+max(laimod_phase)) + delete (s) + delete (dsizes_z) +;------------------------ + DATA11_1D = ndtooned(classob) + DATA12_1D = ndtooned(laiob_phase) + DATA22_1D = ndtooned(laimod_phase) + + yvalues = new((/2,nx/),float) + mn_yvalues = new((/2,nx/),float) + mx_yvalues = new((/2,nx/),float) + + do nd=0,1 + +; See if we are doing model or observational data. + + if(nd.eq.0) then + data_ob = DATA11_1D + data_mod = DATA12_1D + else + data_ob = DATA11_1D + data_mod = DATA22_1D + end if + +; Loop through each range and check for values. + + do i=0,nr-2 + if (i.ne.(nr-2)) then +; print("") +; print("In range ["+range(i)+","+range(i+1)+")") + idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) + else +; print("") +; print("In range ["+range(i)+",)") + idx = ind(range(i).le.data_ob) + end if + +; Calculate average, and get min and max. + + if(.not.any(ismissing(idx))) then + yvalues(nd,i) = avg(data_mod(idx)) + mn_yvalues(nd,i) = min(data_mod(idx)) + mx_yvalues(nd,i) = max(data_mod(idx)) + count = dimsizes(idx) + else + count = 0 + yvalues(nd,i) = yvalues@_FillValue + mn_yvalues(nd,i) = yvalues@_FillValue + mx_yvalues(nd,i) = yvalues@_FillValue + end if + +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) + +; Clean up for next time in loop. + + delete(idx) + end do + delete(data_ob) + delete(data_mod) + end do +;----------------------------------------------------------------- +; compute correlation coef and M score + + u = yvalues(0,:) + v = yvalues(1,:) + + good = ind(.not.ismissing(u) .and. .not.ismissing(v)) + uu = u(good) + vv = v(good) + + ccrPhase = esccr(uu,vv,0) +; print (ccrPhase) + +; old eq +; bias = abs(avg(vv)-avg(uu)) +; new eq + bias = avg(abs(vv-uu)) + + bias = where((bias.gt. 6.),12.-bias,bias) + Mphase = ((6. - bias)/6.)*5. + + print (Mphase) + + delete (u) + delete (v) + delete (uu) + delete (vv) +;-------------------------------------------------------------------- +; histogram res + + resm = True + resm@gsnMaximize = True + resm@gsnDraw = False + resm@gsnFrame = False + resm@xyMarkLineMode = "Markers" + resm@xyMarkerSizeF = 0.014 + resm@xyMarker = 16 + resm@xyMarkerColors = (/"Brown","Blue"/) +; resm@trYMinF = min(mn_yvalues) - 10. +; resm@trYMaxF = max(mx_yvalues) + 10. + resm@trYMinF = min(mn_yvalues) - 2 + resm@trYMaxF = max(mx_yvalues) + 4 + + resm@tiYAxisString = "Max LAI (Leaf Area Index) Month" + resm@tiXAxisString = "Land Cover Type" + + max_bar = new((/2,nx/),graphic) + min_bar = new((/2,nx/),graphic) + max_cap = new((/2,nx/),graphic) + min_cap = new((/2,nx/),graphic) + + lnres = True + line_colors = (/"brown","blue"/) +;------------------------------------------------------------------ +; Start the graphics. + + plot_name = "histogram_phase" + title = model_name + " vs Observed" + resm@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) +;----------------------------- +; Add a boxed legend using the more simple method + + resm@pmLegendDisplayMode = "Always" +; resm@pmLegendWidthF = 0.1 + resm@pmLegendWidthF = 0.08 + resm@pmLegendHeightF = 0.05 + resm@pmLegendOrthogonalPosF = -1.17 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) +; resm@pmLegendParallelPosF = 0.18 + resm@pmLegendParallelPosF = 0.88 ;(rightward) + +; resm@lgPerimOn = False + resm@lgLabelFontHeightF = 0.015 + resm@xyExplicitLegendLabels = (/"observed",model_name/) +;----------------------------- + tRes = True + tRes@txFontHeightF = 0.025 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")" + + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) + + xy = gsn_csm_xy(wks,xvalues,yvalues,resm) +;------------------------------- +;Attach the vertical bar and the horizontal cap line + + do nd=0,1 + lnres@gsLineColor = line_colors(nd) + do i=0,nx-1 + + if(.not.ismissing(mn_yvalues(nd,i)).and. \ + .not.ismissing(mx_yvalues(nd,i))) then + +; Attach the vertical bar, both above and below the marker. + + x1 = xvalues(nd,i) + y1 = yvalues(nd,i) + y2 = mn_yvalues(nd,i) + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + + y2 = mx_yvalues(nd,i) + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + +; Attach the horizontal cap line, both above and below the marker. + + x1 = xvalues(nd,i) - dx4 + x2 = xvalues(nd,i) + dx4 + y1 = mn_yvalues(nd,i) + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + + y1 = mx_yvalues(nd,i) + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + end if + end do + end do + + draw(xy) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +; system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) + + delete (DATA11_1D) + delete (DATA12_1D) + delete (DATA22_1D) +;delete (range) +;delete (xvalues) + delete (yvalues) + delete (mn_yvalues) + delete (mx_yvalues) + delete (good) + delete (max_bar) + delete (min_bar) + delete (max_cap) + delete (min_cap) +;----------------------------------------------------------------- +;global res + + resg = True ; Use plot options + resg@cnFillOn = True ; Turn on color fill + resg@gsnSpreadColors = True ; use full colormap +; resg@cnFillMode = "RasterFill" ; Turn on raster color +; resg@lbLabelAutoStride = True + resg@cnLinesOn = False ; Turn off contourn lines + resg@mpFillOn = False ; Turn off map fill + + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 1. ; Min level + resg@cnMaxLevelValF = 12. ; Max level + resg@cnLevelSpacingF = 1. ; interval + +;global contour ob + + delta = 0.00001 + laiob_phase = where(ismissing(laiob_phase).and.(ismissing(laimod_phase).or.(laimod_phase.lt.delta)),0.,laiob_phase) + + plot_name = "global_phase_ob" + title = ob_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot = gsn_csm_contour_map_ce(wks,laiob_phase,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +; system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;global contour model + + plot_name = "global_phase_model" + title = "Model " + model_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) + gsn_define_colormap(wks,"gui_default") ; choose colormap + + delete (plot) + plot = gsn_csm_contour_map_ce(wks,laimod_phase,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +; system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;global contour model vs ob + + plot_name = "global_phase_model_vs_ob" + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + delete (plot) + plot=new(3,graphic) ; create graphic array + + resg@gsnFrame = False ; Do not draw plot + resg@gsnDraw = False ; Do not advance frame + +; plot correlation coef + + gRes = True + gRes@txFontHeightF = 0.02 + gRes@txAngleF = 90 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")" + + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) +;-------------------------------------------------------------------- + +;(a) ob + + title = ob_name + resg@tiMainString = title + + plot(0) = gsn_csm_contour_map_ce(wks,laiob_phase,resg) + +;(b) model + + title = "Model "+ model_name + resg@tiMainString = title + + plot(1) = gsn_csm_contour_map_ce(wks,laimod_phase,resg) + +;(c) model-ob + + delete (zz) + zz = laimod_phase + zz = laimod_phase - laiob_phase + title = "Model_"+model_name+" - Observed" + resg@tiMainString = title + + resg@cnMinLevelValF = -6. ; Min level + resg@cnMaxLevelValF = 6. ; Max level + resg@cnLevelSpacingF = 1. ; interval + + + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) + +; pres = True ; panel plot mods desired +; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around + ; indiv. plots in panel +; pres@gsnMaximize = True ; fill the page + + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +; system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + frame (wks) + clear (wks) + + delete (plot) +;----------------------------------------------------------------- +;(D) grow day +;-------------------------------------------------------------------- +; get data + + day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/) + +; observed + laiob_grow = laiob(0,:,:) + laiob_grow@long_name = "Days of Growing Season" + + dsizes_z = dimsizes(laiob) + ntime = dsizes_z(0) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + nday = 0. + do k = 0,ntime-1 + if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then + nday = nday + day_of_data(k) + end if + end do + + laiob_grow(j,i) = nday + end do + end do + +; print (min(laiob_grow)+"/"+max(laiob_grow)) +;------------------------- +; model + laimod_grow = laimod(0,:,:) + laimod_grow@long_name = "Days of Growing Season" + + dsizes_z = dimsizes(laimod) + ntime = dsizes_z(0) + nlat = dsizes_z(1) + nlon = dsizes_z(2) + + do j = 0,nlat-1 + do i = 0,nlon-1 + nday = 0. + do k = 0,ntime-1 + if (.not. ismissing(laimod(k,j,i)) .and. laimod(k,j,i) .gt. 1.0) then + nday = nday + day_of_data(k) + end if + end do + + laimod_grow(j,i) = nday + end do + end do + +; print (min(laimod_grow)+"/"+max(laimod_grow)) +;------------------------ + DATA11_1D = ndtooned(classob) + DATA12_1D = ndtooned(laiob_grow) + DATA22_1D = ndtooned(laimod_grow) + + yvalues = new((/2,nx/),float) + mn_yvalues = new((/2,nx/),float) + mx_yvalues = new((/2,nx/),float) + + do nd=0,1 + +; See if we are doing model or observational data. + + if(nd.eq.0) then + data_ob = DATA11_1D + data_mod = DATA12_1D + else + data_ob = DATA11_1D + data_mod = DATA22_1D + end if + +; Loop through each range and check for values. + + do i=0,nr-2 + if (i.ne.(nr-2)) then +; print("") +; print("In range ["+range(i)+","+range(i+1)+")") + idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1))) + else +; print("") +; print("In range ["+range(i)+",)") + idx = ind(range(i).le.data_ob) + end if + +; Calculate average, and get min and max. + + if(.not.any(ismissing(idx))) then + yvalues(nd,i) = avg(data_mod(idx)) + mn_yvalues(nd,i) = min(data_mod(idx)) + mx_yvalues(nd,i) = max(data_mod(idx)) + count = dimsizes(idx) + else + count = 0 + yvalues(nd,i) = yvalues@_FillValue + mn_yvalues(nd,i) = yvalues@_FillValue + mx_yvalues(nd,i) = yvalues@_FillValue + end if + +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i)) +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i)) + +; Clean up for next time in loop. + + delete(idx) + end do + delete(data_ob) + delete(data_mod) + end do +;----------------------------------------------------------------- +; compute correlation coef and M score + + u = yvalues(0,:) + v = yvalues(1,:) + + good = ind(.not.ismissing(u) .and. .not.ismissing(v)) + uu = u(good) + vv = v(good) + + ccrGrow = esccr(uu,vv,0) +; print (ccrGrow) + +; new eq + bias = sum(abs(vv-uu)/(vv+uu)) + Mgrow = (1.- (bias/dimsizes(uu)))*5. + + print (Mgrow) + + delete (u) + delete (v) + delete (uu) + delete (vv) +;-------------------------------------------------------------------- +; histogram res + + resm = True + resm@gsnMaximize = True + resm@gsnDraw = False + resm@gsnFrame = False + resm@xyMarkLineMode = "Markers" + resm@xyMarkerSizeF = 0.014 + resm@xyMarker = 16 + resm@xyMarkerColors = (/"Brown","Blue"/) +; resm@trYMinF = min(mn_yvalues) - 10. +; resm@trYMaxF = max(mx_yvalues) + 10. + resm@trYMinF = min(mn_yvalues) - 2 + resm@trYMaxF = max(mx_yvalues) + 4 + + resm@tiYAxisString = "Days of Growing season" + resm@tiXAxisString = "Land Cover Type" + + max_bar = new((/2,nx/),graphic) + min_bar = new((/2,nx/),graphic) + max_cap = new((/2,nx/),graphic) + min_cap = new((/2,nx/),graphic) + + lnres = True + line_colors = (/"brown","blue"/) +;------------------------------------------------------------------ +; Start the graphics. + + plot_name = "histogram_grow" + title = model_name + " vs Observed" + resm@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) +;----------------------------- +; Add a boxed legend using the more simple method + + resm@pmLegendDisplayMode = "Always" +; resm@pmLegendWidthF = 0.1 + resm@pmLegendWidthF = 0.08 + resm@pmLegendHeightF = 0.05 + resm@pmLegendOrthogonalPosF = -1.17 +; resm@pmLegendOrthogonalPosF = -1.00 ;(downward) +; resm@pmLegendParallelPosF = 0.18 + resm@pmLegendParallelPosF = 0.88 ;(rightward) + +; resm@lgPerimOn = False + resm@lgLabelFontHeightF = 0.015 + resm@xyExplicitLegendLabels = (/"observed",model_name/) +;----------------------------- + tRes = True + tRes@txFontHeightF = 0.025 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")" + + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) + + xy = gsn_csm_xy(wks,xvalues,yvalues,resm) +;------------------------------- +;Attach the vertical bar and the horizontal cap line + + do nd=0,1 + lnres@gsLineColor = line_colors(nd) + do i=0,nx-1 + + if(.not.ismissing(mn_yvalues(nd,i)).and. \ + .not.ismissing(mx_yvalues(nd,i))) then + +; Attach the vertical bar, both above and below the marker. + + x1 = xvalues(nd,i) + y1 = yvalues(nd,i) + y2 = mn_yvalues(nd,i) + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + + y2 = mx_yvalues(nd,i) + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres) + +; Attach the horizontal cap line, both above and below the marker. + + x1 = xvalues(nd,i) - dx4 + x2 = xvalues(nd,i) + dx4 + y1 = mn_yvalues(nd,i) + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + + y1 = mx_yvalues(nd,i) + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres) + end if + end do + end do + + draw(xy) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +; system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) + + delete (DATA11_1D) + delete (DATA12_1D) + delete (DATA22_1D) +;delete (range) +;delete (xvalues) + delete (yvalues) + delete (mn_yvalues) + delete (mx_yvalues) + delete (good) + delete (max_bar) + delete (min_bar) + delete (max_cap) + delete (min_cap) +;----------------------------------------------------------------- +;global res + + resg = True ; Use plot options + resg@cnFillOn = True ; Turn on color fill + resg@gsnSpreadColors = True ; use full colormap +; resg@cnFillMode = "RasterFill" ; Turn on raster color +; resg@lbLabelAutoStride = True + resg@cnLinesOn = False ; Turn off contourn lines + resg@mpFillOn = False ; Turn off map fill + + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 0. ; Min level + resg@cnMaxLevelValF = 360. ; Max level + resg@cnLevelSpacingF = 30. ; interval + +;global contour ob + + delta = 0.00001 + laiob_grow = where(ismissing(laiob_grow).and.(ismissing(laimod_grow).or.(laimod_grow.lt.delta)),0.,laiob_grow) + + plot_name = "global_grow_ob" + title = ob_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + plot = gsn_csm_contour_map_ce(wks,laiob_grow,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +; system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;global contour model + + plot_name = "global_grow_model" + title = "Model " + model_name + resg@tiMainString = title + + wks = gsn_open_wks (plot_type,plot_name) + gsn_define_colormap(wks,"gui_default") ; choose colormap + + delete (plot) + plot = gsn_csm_contour_map_ce(wks,laimod_grow,resg) + frame(wks) + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +; system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + clear (wks) +;------------------------------------------------------------------------ +;global contour model vs ob + + plot_name = "global_grow_model_vs_ob" + + wks = gsn_open_wks (plot_type,plot_name) ; open workstation + gsn_define_colormap(wks,"gui_default") ; choose colormap + + delete (plot) + plot=new(3,graphic) ; create graphic array + + resg@gsnFrame = False ; Do not draw plot + resg@gsnDraw = False ; Do not advance frame + +; plot correlation coef + + gRes = True + gRes@txFontHeightF = 0.02 + gRes@txAngleF = 90 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")" + + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) +;-------------------------------------------------------------------- + +;(a) ob + + title = ob_name + resg@tiMainString = title + + plot(0) = gsn_csm_contour_map_ce(wks,laiob_grow,resg) + +;(b) model + + title = "Model "+ model_name + resg@tiMainString = title + + plot(1) = gsn_csm_contour_map_ce(wks,laimod_grow,resg) + +;(c) model-ob + + delete (zz) + zz = laimod_grow + zz = laimod_grow - laiob_grow + title = "Model_"+model_name+" - Observed" + resg@tiMainString = title + + resg@cnMinLevelValF = -120. ; Min level + resg@cnMaxLevelValF = 120. ; Max level + resg@cnLevelSpacingF = 20. ; interval + + + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) + + pres = True ; panel plot mods desired + pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around + ; indiv. plots in panel + pres@gsnMaximize = True ; fill the page + + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot + + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) +; system("rm "+plot_name+"."+plot_type) +; system("rm "+plot_name+"-1."+plot_type_new) +; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new) + + frame (wks) + clear (wks) + + delete (plot) +end +