1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/biomass/99.all.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,262 @@
1.4 +; ****************************************************
1.5 +; combine scatter, histogram, global and zonal plots
1.6 +; compute all correlation coef and M score
1.7 +; ****************************************************
1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.10 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.11 +;************************************************
1.12 +begin
1.13 +
1.14 + plot_type = "ps"
1.15 + plot_type_new = "png"
1.16 +;************************************************
1.17 +; read in model data
1.18 +;************************************************
1.19 +;film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
1.20 +;model_name = "b30.061n"
1.21 +;model_grid = "T31"
1.22 +
1.23 +;film = "newcn05_ncep_1i_ANN_climo_lnd.nc"
1.24 +;model_name = "newcn"
1.25 +;model_grid = "1.9"
1.26 +;--------------------------------------------
1.27 + film = "i01.10cn_1948-2004_ANN_climo.nc"
1.28 + model_name = "10cn"
1.29 + model_grid = "T42"
1.30 +
1.31 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.32 + fm = addfile (dirm+film,"r")
1.33 +
1.34 +;unit: gC/m2
1.35 + data1 = fm->LIVESTEMC
1.36 + data2 = fm->DEADSTEMC
1.37 + data3 = fm->LEAFC
1.38 + datamod0 = data1(0,:,:)
1.39 + datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)
1.40 +;-------------------------------------------
1.41 +;film = "i01.10casa_1948-2004_ANN_climo.nc"
1.42 +;model_name = "10casa"
1.43 +;model_grid = "T42"
1.44 +
1.45 +;dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.46 +;fm = addfile (dirm+film,"r")
1.47 +
1.48 +;unit: gC/m2
1.49 +;factor_WOODC = 0.7
1.50 +;data1 = fm->WOODC
1.51 +;data2 = fm->LEAFC
1.52 +;datamod0 = data1(0,:,:)
1.53 +;datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:)
1.54 +;----------------------------------------------------
1.55 + html_name = "table.html." + model_name
1.56 + html_new = html_name +".new"
1.57 + system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \
1.58 + "mv -f "+html_new+" "+html_name)
1.59 +
1.60 + xm = fm->lon
1.61 + ym = fm->lat
1.62 +;************************************************
1.63 +; read in ob data
1.64 +;************************************************
1.65 + ob_name = "LC15_Amazon_Biomass"
1.66 +
1.67 + diro = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"
1.68 + filo = "amazon_biomass_"+model_grid+".nc"
1.69 + fo = addfile (diro+filo,"r")
1.70 +
1.71 + dataob = fo->BIOMASS
1.72 + xo = fo->lon
1.73 + yo = fo->lat
1.74 +;************************************************
1.75 +; Units for these variables are:
1.76 +; dataob : MgC/ha
1.77 +; datamod0 : gC/m2
1.78 +; We want to convert these to KgC/m2
1.79 +; ha = 100m*100m = 10,000 m2
1.80 +; MgC/ha*1000/10,000 = KgC/m2
1.81 +
1.82 + factor_aboveground = 0.5
1.83 + factor_unit_ob = 0.1
1.84 + factor_unit_mod = 0.001
1.85 +
1.86 + dataob = dataob * factor_aboveground * factor_unit_ob
1.87 + datamod0 = datamod0 * factor_unit_mod
1.88 +
1.89 + dataob@units = "KgC/m2"
1.90 + datamod0@units = "KgC/m2"
1.91 +
1.92 + dataob@long_name = "Amazon Biomass"
1.93 + datamod0@long_name = "Amazon Biomass"
1.94 +;********************************************************
1.95 +; get subset of datamod0 that match dataob
1.96 +
1.97 + nlon = dimsizes(xo)
1.98 + nlat = dimsizes(yo)
1.99 +
1.100 + ind_lonL = ind(xm .eq. xo(0))
1.101 + ind_lonR = ind(xm .eq. xo(nlon-1))
1.102 + ind_latS = ind(ym .eq. yo(0))
1.103 + ind_latN = ind(ym .eq. yo(nlat-1))
1.104 +
1.105 + datamod = dataob
1.106 + datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)
1.107 +
1.108 +; print (datamod)
1.109 +
1.110 +;----------------------------------------------------------------------
1.111 +; global contour
1.112 +
1.113 +;res
1.114 + resg = True ; Use plot options
1.115 + resg@cnFillOn = True ; Turn on color fill
1.116 + resg@gsnSpreadColors = True ; use full colormap
1.117 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.118 +; resg@lbLabelAutoStride = True
1.119 + resg@cnLinesOn = False ; Turn off contourn lines
1.120 + resg@mpFillOn = False ; Turn off map fill
1.121 + resg@gsnAddCyclic = False
1.122 +
1.123 + resg@gsnSpreadColors = True ; use full colormap
1.124 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.125 + resg@cnMinLevelValF = 0. ; Min level
1.126 + resg@cnMaxLevelValF = 30. ; Max level
1.127 + resg@cnLevelSpacingF = 2. ; interval
1.128 +
1.129 + resg@mpMinLatF = -21.1 ; range to zoom in on
1.130 + resg@mpMaxLatF = 13.8
1.131 + resg@mpMinLonF = 277.28
1.132 + resg@mpMaxLonF = 326.43
1.133 +;------------------------------------------------------------------------
1.134 +;global contour ob
1.135 +
1.136 + plot_name = "global_ob"
1.137 + title = ob_name+" "+ model_grid
1.138 + resg@tiMainString = title
1.139 +
1.140 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.141 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.142 +
1.143 + plot = gsn_csm_contour_map_ce(wks,dataob,resg)
1.144 + frame(wks)
1.145 +
1.146 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.147 + "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1.148 + "rm "+plot_name+"-*."+plot_type_new+";"+ \
1.149 + "rm "+plot_name+"."+plot_type)
1.150 +
1.151 + clear (wks)
1.152 +;------------------------------------------------------------------------
1.153 +;global contour model
1.154 +
1.155 + plot_name = "global_model"
1.156 + title = "Model "+ model_name
1.157 + resg@tiMainString = title
1.158 +
1.159 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.160 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.161 +
1.162 + plot = gsn_csm_contour_map_ce(wks,datamod,resg)
1.163 + frame(wks)
1.164 +
1.165 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.166 + "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1.167 + "rm "+plot_name+"-*."+plot_type_new+";"+ \
1.168 + "rm "+plot_name+"."+plot_type)
1.169 +
1.170 + clear (wks)
1.171 +;------------------------------------------------------------------------
1.172 +;global contour model vs ob
1.173 +
1.174 + plot_name = "global_model_vs_ob"
1.175 +
1.176 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.177 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.178 +
1.179 + delete (plot)
1.180 + plot=new(3,graphic) ; create graphic array
1.181 +
1.182 + resg@gsnFrame = False ; Do not draw plot
1.183 + resg@gsnDraw = False ; Do not advance frame
1.184 +
1.185 +;(d) compute correlation coef and M score
1.186 +
1.187 + uu = ndtooned(datamod)
1.188 + vv = ndtooned(dataob)
1.189 +
1.190 + good = ind(.not.ismissing(uu) .and. .not.ismissing(vv))
1.191 +
1.192 + ug = uu(good)
1.193 + vg = vv(good)
1.194 +
1.195 + ccrG = esccr(ug,vg,0)
1.196 +; print (ccrG)
1.197 +
1.198 + MG = (ccrG*ccrG)* 5.0
1.199 +; new eq
1.200 + bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
1.201 + MG = (1. - (bias/dimsizes(ug)))*5.
1.202 +
1.203 + M_biomass = sprintf("%.2f", MG)
1.204 + system("sed s#M_biomass#"+M_biomass+"# "+html_name+" > "+html_new+";"+ \
1.205 + "mv -f "+html_new+" "+html_name)
1.206 + print (M_biomass)
1.207 +
1.208 +; plot correlation coef
1.209 +
1.210 + gRes = True
1.211 + gRes@txFontHeightF = 0.02
1.212 + gRes@txAngleF = 90
1.213 +
1.214 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
1.215 +
1.216 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.217 +;--------------------------------------------------------------------
1.218 +
1.219 +;(a) ob
1.220 +
1.221 + title = ob_name+" "+ model_grid
1.222 + resg@tiMainString = title
1.223 +
1.224 + plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg)
1.225 +
1.226 +;(b) model
1.227 +
1.228 + title = "Model "+ model_name
1.229 + resg@tiMainString = title
1.230 +
1.231 + plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg)
1.232 +
1.233 +;(c) model-ob
1.234 +
1.235 + zz = datamod
1.236 + zz = datamod - dataob
1.237 + title = "Model_"+model_name+" - Observed"
1.238 +
1.239 + resg@cnMinLevelValF = -10. ; Min level
1.240 + resg@cnMaxLevelValF = 10. ; Max level
1.241 + resg@cnLevelSpacingF = 2. ; interval
1.242 + resg@tiMainString = title
1.243 +
1.244 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.245 +
1.246 + pres = True ; panel plot mods desired
1.247 +; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.248 + ; indiv. plots in panel
1.249 + pres@gsnMaximize = True ; fill the page
1.250 +
1.251 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.252 +
1.253 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.254 + "rm "+plot_name+"."+plot_type)
1.255 +
1.256 + clear (wks)
1.257 + delete (plot)
1.258 +;-------------------------------------------------------------------
1.259 + temp_name = "biomass." + model_name
1.260 + system("mkdir -p " + temp_name+";"+ \
1.261 + "cp "+ html_name + " " +temp_name+"/table.html"+";"+ \
1.262 + "mv *.png " + temp_name +";"+ \
1.263 + "tar cf "+ temp_name +".tar " + temp_name)
1.264 +;-------------------------------------------------------------------
1.265 +end