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