diff -r 000000000000 -r 0c6405ab2ff4 biomass/99.all.ncl.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/biomass/99.all.ncl.1 Mon Jan 26 22:08:20 2009 -0500 @@ -0,0 +1,264 @@ +; **************************************************** +; combine scatter, histogram, global and zonal plots +; compute all correlation coef and M score +; **************************************************** +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" +;************************************************ +begin + + plot_type = "ps" + plot_type_new = "png" +;************************************************ +; read in model data +;************************************************ +;film = "b30.061n_1995-2004_ANN_climo_lnd.nc" +;model_name = "b30.061n" +;model_grid = "T31" + +;film = "newcn05_ncep_1i_ANN_climo_lnd.nc" +;model_name = "newcn" +;model_grid = "1.9" +;-------------------------------------------- +;film = "i01.10cn_1948-2004_ANN_climo.nc" +;model_name = "10cn" +;model_grid = "T42" + +;dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" +;fm = addfile (dirm+film,"r") + +;unit: gC/m2 +;data1 = fm->LIVESTEMC +;data2 = fm->DEADSTEMC +;data3 = fm->LEAFC +;datamod0 = data1(0,:,:) +;datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:) +;------------------------------------------- + film = "i01.10casa_1948-2004_ANN_climo.nc" + model_name = "10casa" + model_grid = "T42" + + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" + fm = addfile (dirm+film,"r") + +;unit: gC/m2 + factor_WOODC = 0.7 + data1 = fm->WOODC + data2 = fm->LEAFC + datamod0 = data1(0,:,:) + datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:) +;---------------------------------------------------- + + system("sed s#model_name#"+model_name+"# table.html > table.html.new") + system("mv -f table.html.new table.html") + + xm = fm->lon + ym = fm->lat +;************************************************ +; read in ob data +;************************************************ + ob_name = "LC15_Amazon_Biomass" + + diro = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/" + filo = "amazon_biomass_"+model_grid+".nc" + fo = addfile (diro+filo,"r") + + dataob = fo->BIOMASS + xo = fo->lon + yo = fo->lat +;************************************************ +; Units for these variables are: +; dataob : MgC/ha +; datamod0 : gC/m2 +; We want to convert these to KgC/m2 +; ha = 100m*100m = 10,000 m2 +; MgC/ha*1000/10,000 = KgC/m2 + + factor_aboveground = 0.5 + factor_unit_ob = 0.1 + factor_unit_mod = 0.001 + + dataob = dataob * factor_aboveground * factor_unit_ob + datamod0 = datamod0 * factor_unit_mod + + dataob@units = "KgC/m2" + datamod0@units = "KgC/m2" + + dataob@long_name = "Amazon Biomass" + datamod0@long_name = "Amazon Biomass" +;******************************************************** +; get subset of datamod0 that match dataob + + nlon = dimsizes(xo) + nlat = dimsizes(yo) + + ind_lonL = ind(xm .eq. xo(0)) + ind_lonR = ind(xm .eq. xo(nlon-1)) + ind_latS = ind(ym .eq. yo(0)) + ind_latN = ind(ym .eq. yo(nlat-1)) + + datamod = dataob + datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR) + +; print (datamod) + +;---------------------------------------------------------------------- +; global contour + +;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@gsnAddCyclic = False + + resg@gsnSpreadColors = True ; use full colormap + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + resg@cnMinLevelValF = 0. ; Min level + resg@cnMaxLevelValF = 30. ; Max level + resg@cnLevelSpacingF = 2. ; interval + + resg@mpMinLatF = -21.1 ; range to zoom in on + resg@mpMaxLatF = 13.8 + resg@mpMinLonF = 277.28 + resg@mpMaxLonF = 326.43 +;------------------------------------------------------------------------ +;global contour ob + + plot_name = "global_ob" + title = ob_name+" "+ model_grid + 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,dataob,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_model" + title = "Model "+ model_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,datamod,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_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 + +;(d) compute correlation coef and M score + + uu = ndtooned(datamod) + vv = ndtooned(dataob) + + good = ind(.not.ismissing(uu) .and. .not.ismissing(vv)) + + ug = uu(good) + vg = vv(good) + + ccrG = esccr(ug,vg,0) +; print (ccrG) + + MG = (ccrG*ccrG)* 5.0 +; new eq + bias = sum(abs(ug-vg)/(abs(ug)+abs(vg))) + MG = (1. - (bias/dimsizes(ug)))*5. + + M_biomass = sprintf("%.2f", MG) + system("sed s#M_biomass#"+M_biomass+"# table.html > table.html.new") + system("mv -f table.html.new table.html") + print (M_biomass) + +; plot correlation coef + + gRes = True + gRes@txFontHeightF = 0.02 + gRes@txAngleF = 90 + + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")" + + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) +;-------------------------------------------------------------------- + +;(a) ob + + title = ob_name+" "+ model_grid + resg@tiMainString = title + + plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg) + +;(b) model + + title = "Model "+ model_name + resg@tiMainString = title + + plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg) + +;(c) model-ob + + zz = datamod + zz = datamod - dataob + title = "Model_"+model_name+" - Observed" + + resg@cnMinLevelValF = -10. ; Min level + resg@cnMaxLevelValF = 10. ; Max level + resg@cnLevelSpacingF = 2. ; interval + resg@tiMainString = title + + 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) +;------------------------------------------------------------------------ + temp_name = "biomass." + model_name + system("mkdir -p " + temp_name) + system("cp table.html " + temp_name) + system("mv *.png " + temp_name) + system("tar cf "+ temp_name +".tar " + temp_name) +;------------------------------------------------------------------------ +end