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