1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/all/04.biomass.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,536 @@
1.4 +; ****************************************************
1.5 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.6 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.7 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.8 +;************************************************
1.9 +begin
1.10 +
1.11 + plot_type = "ps"
1.12 + plot_type_new = "png"
1.13 +
1.14 +;------------------------------------------------------
1.15 +; edit table.html of current model for movel1_vs_model2
1.16 +
1.17 + if (isvar("compare")) then
1.18 + html_name2 = compare+"/table.html"
1.19 + html_new2 = html_name2 +".new"
1.20 + end if
1.21 +
1.22 +;------------------------------------------------------
1.23 +; edit table.html for current model
1.24 +
1.25 + html_name = model_name+"/table.html"
1.26 + html_new = html_name +".new"
1.27 +
1.28 +;---------------------------------
1.29 +; get model data: landfrac and area
1.30 +
1.31 + film_l = "lnd_"+ model_grid + ".nc"
1.32 + fm_l = addfile (dirs+film_l,"r")
1.33 +
1.34 + landfrac = fm_l->landfrac
1.35 + area0 = fm_l->area
1.36 +
1.37 + delete (fm_l)
1.38 +
1.39 +; change area from km**2 to m**2
1.40 + area0 = area0 * 1.e6
1.41 +
1.42 +;-----------------------------
1.43 +; take into account landfrac
1.44 +
1.45 + area0 = area0 * landfrac
1.46 +
1.47 +;--------------------------------------------
1.48 +; read model data
1.49 +
1.50 + fm = addfile (dirm+film1,"r")
1.51 +
1.52 + if (BGC .eq. "cn") then
1.53 + data1 = fm->LIVESTEMC
1.54 + data2 = fm->DEADSTEMC
1.55 + data3 = fm->LEAFC
1.56 + datamod0 = data1(0,:,:)
1.57 + datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)
1.58 + end if
1.59 +
1.60 + if (BGC .eq. "casa") then
1.61 + factor_WOODC = 0.7
1.62 + data1 = fm->WOODC
1.63 + data2 = fm->LEAFC
1.64 + datamod0 = data1(0,:,:)
1.65 + datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:)
1.66 + end if
1.67 +
1.68 +; unit: gC/m2
1.69 +
1.70 + xm = fm->lon
1.71 + ym = fm->lat
1.72 +
1.73 + delete (fm)
1.74 +
1.75 +;------------------------------------------------
1.76 +; read amazon mask data
1.77 +
1.78 + dir_m = diro + "biomass/"
1.79 + fil_m = "amazon_mask_"+ model_grid + ".nc"
1.80 + fm = addfile (dir_m+fil_m,"r")
1.81 +
1.82 + mask_amazon0 = fm->mask_amazon
1.83 +
1.84 + delete (fm)
1.85 +
1.86 +;------------------------------------------------
1.87 +; read ob data
1.88 +
1.89 + ob_name = "LC15_Amazon_Biomass"
1.90 +
1.91 + dir_b = diro + "biomass/"
1.92 + fil_b = "amazon_biomass_"+model_grid+".nc"
1.93 + fo = addfile (dir_b+fil_b,"r")
1.94 +
1.95 + dataob = fo->BIOMASS
1.96 + xo = fo->lon
1.97 + yo = fo->lat
1.98 +
1.99 + delete (fo)
1.100 +
1.101 +;************************************************
1.102 +; Units for these variables are:
1.103 +; dataob : MgC/ha
1.104 +; datamod0 : gC/m2
1.105 +; We want to convert these to KgC/m2
1.106 +; ha = 100m*100m = 10,000 m2
1.107 +; MgC/ha*1000/10,000 = KgC/m2
1.108 +
1.109 + factor_aboveground = 0.5
1.110 + factor_unit_ob = 0.1
1.111 + factor_unit_mod = 0.001
1.112 +
1.113 + dataob = dataob * factor_aboveground * factor_unit_ob
1.114 + datamod0 = datamod0 * factor_unit_mod
1.115 +
1.116 + dataob@units = "KgC/m2"
1.117 + datamod0@units = "KgC/m2"
1.118 +
1.119 + dataob@long_name = "Amazon Biomass"
1.120 + datamod0@long_name = "Amazon Biomass"
1.121 +;********************************************************
1.122 +; get subset of datamod0 that match dataob
1.123 +
1.124 + nlon = dimsizes(xo)
1.125 + nlat = dimsizes(yo)
1.126 +
1.127 + ind_lonL = ind(xm .eq. xo(0))
1.128 + ind_lonR = ind(xm .eq. xo(nlon-1))
1.129 + ind_latS = ind(ym .eq. yo(0))
1.130 + ind_latN = ind(ym .eq. yo(nlat-1))
1.131 +
1.132 + datamod = dataob
1.133 + datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)
1.134 +
1.135 + area = dataob
1.136 + area(:,:) = area0(ind_latS:ind_latN,ind_lonL:ind_lonR)
1.137 +
1.138 + mask_amazon = dataob
1.139 + mask_amazon(:,:) = mask_amazon0(ind_latS:ind_latN,ind_lonL:ind_lonR)
1.140 +
1.141 + mask_amazon@units = ""
1.142 +;********************************************************
1.143 +; sum over amazom_mask area:
1.144 +
1.145 +; Peta g = 1.e15 g = 1.e12 Kg
1.146 + factor_unit = 1.e-12
1.147 +
1.148 +; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.)
1.149 +
1.150 + Sum_area = sum(area*mask_amazon)*factor_unit
1.151 +
1.152 + Sum_ob = sum(dataob*area*mask_amazon)*factor_unit
1.153 + Sum_mod = sum(datamod*area*mask_amazon)*factor_unit
1.154 +
1.155 + avg_ob = Sum_ob /Sum_area
1.156 + avg_mod = Sum_mod/Sum_area
1.157 +
1.158 + Sum_biomass_ob = sprintf("%.2f",Sum_ob )
1.159 + Sum_biomass_mod = sprintf("%.2f",Sum_mod)
1.160 +
1.161 +;----------------------------------------------------------------------
1.162 +; contour plot res
1.163 +
1.164 + resg = True ; Use plot options
1.165 + resg@cnFillOn = True ; Turn on color fill
1.166 + resg@gsnSpreadColors = True ; use full colormap
1.167 + resg@cnLinesOn = False ; Turn off contourn lines
1.168 + resg@mpFillOn = False ; Turn off map fill
1.169 + resg@gsnAddCyclic = False
1.170 +
1.171 + resg@gsnSpreadColors = True ; use full colormap
1.172 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.173 +
1.174 + resg@mpMinLatF = -21.1 ; range to zoom in on
1.175 + resg@mpMaxLatF = 13.8
1.176 + resg@mpMinLonF = 277.28
1.177 + resg@mpMaxLonF = 326.43
1.178 +;------------------------------------------------------------------------
1.179 +; mask plot
1.180 +
1.181 + plot_name = "mask_ob"
1.182 +
1.183 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.184 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.185 +
1.186 +;-----------------------------------------
1.187 +; plot area sum
1.188 +
1.189 + gRes = True
1.190 + gRes@txFontHeightF = 0.02
1.191 +; gRes@txAngleF = 90
1.192 +
1.193 + area_sum_text = "(mask area = "+sprintf("%.2f", Sum_area)+"e12 m2)"
1.194 +
1.195 + gsn_text_ndc(wks,area_sum_text,0.50,0.80,gRes)
1.196 +;-----------------------------------------
1.197 +
1.198 + resg@cnMinLevelValF = 0. ; Min level
1.199 + resg@cnMaxLevelValF = 1. ; Max level
1.200 + resg@cnLevelSpacingF = 0.1 ; interval
1.201 +
1.202 + resg@tiMainString = "Amazon Mask: grid = "+ model_grid
1.203 +
1.204 + plot = gsn_csm_contour_map_ce(wks,mask_amazon,resg)
1.205 +
1.206 + delete (wks)
1.207 +
1.208 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.209 + "rm "+plot_name+"."+plot_type)
1.210 +
1.211 +;------------------------------------------------------------------------
1.212 +; contour ob
1.213 +
1.214 + resg@cnMinLevelValF = 0. ; Min level
1.215 + resg@cnMaxLevelValF = 30. ; Max level
1.216 + resg@cnLevelSpacingF = 2. ; interval
1.217 +
1.218 + plot_name = "global_ob"
1.219 + title = ob_name+" "+ model_grid
1.220 + resg@tiMainString = title
1.221 +
1.222 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.223 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.224 +
1.225 + plot = gsn_csm_contour_map_ce(wks,dataob,resg)
1.226 +
1.227 + delete (wks)
1.228 +
1.229 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.230 + "rm "+plot_name+"."+plot_type)
1.231 +
1.232 +;------------------------------------------------------------------------
1.233 +; contour model
1.234 +
1.235 + resg@cnMinLevelValF = 0. ; Min level
1.236 + resg@cnMaxLevelValF = 30. ; Max level
1.237 + resg@cnLevelSpacingF = 2. ; interval
1.238 +
1.239 + plot_name = "global_model"
1.240 + title = "Model "+ model_name
1.241 + resg@tiMainString = title
1.242 +
1.243 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.244 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.245 +
1.246 + plot = gsn_csm_contour_map_ce(wks,datamod,resg)
1.247 +
1.248 + delete (wks)
1.249 +
1.250 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.251 + "rm "+plot_name+"."+plot_type)
1.252 +
1.253 +;------------------------------------------------------------------------
1.254 +; contour model vs ob
1.255 +
1.256 + plot_name = "global_model_vs_ob"
1.257 +
1.258 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.259 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.260 +
1.261 + delete (plot)
1.262 + plot=new(3,graphic) ; create graphic array
1.263 +
1.264 + resg@gsnFrame = False ; Do not draw plot
1.265 + resg@gsnDraw = False ; Do not advance frame
1.266 +
1.267 +;(d) compute correlation coef and M score
1.268 +
1.269 + uu = ndtooned(datamod)
1.270 + vv = ndtooned(dataob)
1.271 +
1.272 + good = ind(.not.ismissing(uu) .and. .not.ismissing(vv))
1.273 +
1.274 + ug = uu(good)
1.275 + vg = vv(good)
1.276 +
1.277 + ccrG = esccr(ug,vg,0)
1.278 +
1.279 + score_max = 5.0
1.280 +
1.281 +; Miomass = (ccrG*ccrG)* score_max
1.282 +; new eq
1.283 + bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
1.284 + Mbiomass = (1. - (bias/dimsizes(ug)))*score_max
1.285 + M_biomass = sprintf("%.2f", Mbiomass)
1.286 +
1.287 + if (isvar("compare")) then
1.288 + system("sed -e '1,/M_biomass/s/M_biomass/"+M_biomass+"/' "+html_name2+" > "+html_new2+";"+ \
1.289 + "mv -f "+html_new2+" "+html_name2)
1.290 + end if
1.291 +
1.292 + system("sed s#M_biomass#"+M_biomass+"# "+html_name+" > "+html_new+";"+ \
1.293 + "mv -f "+html_new+" "+html_name)
1.294 +
1.295 +; plot correlation coef
1.296 +
1.297 + gRes = True
1.298 + gRes@txFontHeightF = 0.02
1.299 + gRes@txAngleF = 90
1.300 +
1.301 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
1.302 +
1.303 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.304 +;--------------------------------------------------------------------
1.305 +
1.306 +;(a) ob
1.307 +
1.308 + title = ob_name+" "+ model_grid
1.309 + resg@tiMainString = title
1.310 +
1.311 + plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg)
1.312 +
1.313 +;(b) model
1.314 +
1.315 + title = "Model "+ model_name
1.316 + resg@tiMainString = title
1.317 +
1.318 + plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg)
1.319 +
1.320 +;(c) model-ob
1.321 +
1.322 + zz = datamod
1.323 + zz = datamod - dataob
1.324 + title = "Model_"+model_name+" - Observed"
1.325 +
1.326 + resg@cnMinLevelValF = -10. ; Min level
1.327 + resg@cnMaxLevelValF = 10. ; Max level
1.328 + resg@cnLevelSpacingF = 2. ; interval
1.329 + resg@tiMainString = title
1.330 +
1.331 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.332 +
1.333 + pres = True ; panel plot mods desired
1.334 +; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.335 + ; indiv. plots in panel
1.336 + pres@gsnMaximize = True ; fill the page
1.337 +
1.338 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.339 +
1.340 + delete (wks)
1.341 + delete (plot)
1.342 + delete (zz)
1.343 +
1.344 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.345 + "rm "+plot_name+"."+plot_type)
1.346 +
1.347 + resg@gsnFrame = True ; draw plot
1.348 + resg@gsnDraw = True ; advance frame
1.349 +;------------------------------------------------------------------------
1.350 +; contour ob : masked
1.351 +
1.352 + resg@cnMinLevelValF = 0. ; Min level
1.353 + resg@cnMaxLevelValF = 30. ; Max level
1.354 + resg@cnLevelSpacingF = 2. ; interval
1.355 +
1.356 + plot_name = "global_mask_ob"
1.357 + title = ob_name+" "+ model_grid
1.358 + resg@tiMainString = title
1.359 +
1.360 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.361 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.362 +;-----------------------------------------
1.363 +; plot average over mask region
1.364 +
1.365 + gRes = True
1.366 + gRes@txFontHeightF = 0.02
1.367 + gRes@txAngleF = 0
1.368 +
1.369 + area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_ob)+" Kg C/m2)"
1.370 +
1.371 + gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)
1.372 +;-----------------------------------------
1.373 + zo = dataob
1.374 + zo = dataob*mask_amazon
1.375 + zo = where((mask_amazon .le. 0.01), zo@_FillValue, zo)
1.376 + plot = gsn_csm_contour_map_ce(wks,zo,resg)
1.377 +
1.378 + delete (wks)
1.379 + delete (plot)
1.380 +
1.381 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.382 + "rm "+plot_name+"."+plot_type)
1.383 +
1.384 +;------------------------------------------------------------------------
1.385 +; contour model: masked
1.386 +
1.387 + resg@cnMinLevelValF = 0. ; Min level
1.388 + resg@cnMaxLevelValF = 30. ; Max level
1.389 + resg@cnLevelSpacingF = 2. ; interval
1.390 +
1.391 + plot_name = "global_mask_model"
1.392 + title = "Model "+ model_name
1.393 + resg@tiMainString = title
1.394 +
1.395 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.396 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.397 +;-----------------------------------------
1.398 +; plot average over mask region
1.399 +
1.400 + gRes = True
1.401 + gRes@txFontHeightF = 0.02
1.402 + gRes@txAngleF = 0
1.403 +
1.404 + area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_mod)+" Kg C/m2)"
1.405 +
1.406 + gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)
1.407 +;-----------------------------------------
1.408 + zm = datamod
1.409 + zm = datamod*mask_amazon
1.410 + zm = where((mask_amazon .le. 0.01), zm@_FillValue, zm)
1.411 + plot = gsn_csm_contour_map_ce(wks,zm,resg)
1.412 +
1.413 + delete (wks)
1.414 + delete (plot)
1.415 +
1.416 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.417 + "rm "+plot_name+"."+plot_type)
1.418 +
1.419 +;------------------------------------------------------------------------
1.420 +; contour model vs ob: masked
1.421 +
1.422 + plot_name = "global_mask_vs_ob"
1.423 +
1.424 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.425 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.426 +
1.427 + plot=new(3,graphic) ; create graphic array
1.428 +
1.429 + resg@gsnFrame = False ; Do not draw plot
1.430 + resg@gsnDraw = False ; Do not advance frame
1.431 +
1.432 +;(d) compute correlation coef and M score
1.433 +
1.434 + delete (good)
1.435 + delete (uu)
1.436 + delete (vv)
1.437 + delete (ug)
1.438 + delete (vg)
1.439 +
1.440 + score_max = 5.
1.441 +
1.442 + uu = ndtooned(zm)
1.443 + vv = ndtooned(zo)
1.444 +
1.445 + good = ind((uu .gt. 0.) .and. (vv .gt. 0.))
1.446 +
1.447 + ug = uu(good)
1.448 + vg = vv(good)
1.449 +
1.450 + ccrG = esccr(ug,vg,0)
1.451 +
1.452 +; Miomass = (ccrG*ccrG)*score_max
1.453 +; new eq
1.454 + bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))
1.455 + Mbiomass2 = (1. - (bias/dimsizes(ug)))*score_max
1.456 + M_biomask = sprintf("%.2f", Mbiomass2)
1.457 +
1.458 + if (isvar("compare")) then
1.459 + system("sed -e '1,/M_biomask/s/M_biomask/"+M_biomask+"/' "+html_name2+" > "+html_new2+";"+ \
1.460 + "mv -f "+html_new2+" "+html_name2)
1.461 + system("sed -e '1,/Sum_biomass_ob/s/Sum_biomass_ob/"+Sum_biomass_ob+"/' "+html_name2+" > "+html_new2+";"+ \
1.462 + "mv -f "+html_new2+" "+html_name2)
1.463 + system("sed -e '1,/Sum_biomass_mod/s/Sum_biomass_mod/"+Sum_biomass_mod+"/' "+html_name2+" > "+html_new2+";"+ \
1.464 + "mv -f "+html_new2+" "+html_name2)
1.465 + end if
1.466 +
1.467 + system("sed s#M_biomask#"+M_biomask+"# "+html_name+" > "+html_new+";"+ \
1.468 + "mv -f "+html_new+" "+html_name)
1.469 + system("sed s#Sum_biomass_ob#"+Sum_biomass_ob+"# "+html_name+" > "+html_new+";"+ \
1.470 + "mv -f "+html_new+" "+html_name)
1.471 + system("sed s#Sum_biomass_mod#"+Sum_biomass_mod+"# "+html_name+" > "+html_new+";"+ \
1.472 + "mv -f "+html_new+" "+html_name)
1.473 +;--------------------------------------------------------------------
1.474 +; plot correlation coef
1.475 +
1.476 + gRes = True
1.477 + gRes@txFontHeightF = 0.02
1.478 + gRes@txAngleF = 90
1.479 +
1.480 + correlation_text = "(correlation coef = "+sprintf("%.2f", ccrG)+")"
1.481 +
1.482 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.483 +;--------------------------------------------------------------------
1.484 +
1.485 +;(a) ob
1.486 +
1.487 + title = ob_name+" "+ model_grid
1.488 + resg@tiMainString = title
1.489 +
1.490 + plot(0) = gsn_csm_contour_map_ce(wks,zo,resg)
1.491 +
1.492 +;(b) model
1.493 +
1.494 + title = "Model "+ model_name
1.495 + resg@tiMainString = title
1.496 +
1.497 + plot(1) = gsn_csm_contour_map_ce(wks,zm,resg)
1.498 +
1.499 +;(c) model-ob
1.500 +
1.501 + zz = zo
1.502 + zz = zm - zo
1.503 + title = "Model_"+model_name+" - Observed"
1.504 +
1.505 + resg@cnMinLevelValF = -10. ; Min level
1.506 + resg@cnMaxLevelValF = 10. ; Max level
1.507 + resg@cnLevelSpacingF = 2. ; interval
1.508 + resg@tiMainString = title
1.509 +
1.510 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.511 +
1.512 + pres = True ; panel plot mods desired
1.513 +; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.514 + ; indiv. plots in panel
1.515 + pres@gsnMaximize = True ; fill the page
1.516 +
1.517 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.518 +
1.519 + delete (wks)
1.520 + delete (plot)
1.521 +
1.522 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.523 + "rm "+plot_name+"."+plot_type)
1.524 +
1.525 +;***************************************************************************
1.526 +; add total score and write to file
1.527 +;***************************************************************************
1.528 + M_total = Mbiomass + Mbiomass2
1.529 +
1.530 + asciiwrite("M_save.biomass", M_total)
1.531 +
1.532 +;***************************************************************************
1.533 +; output plots
1.534 +;***************************************************************************
1.535 + output_dir = model_name+"/biomass"
1.536 +
1.537 + system("mv *.png " + output_dir)
1.538 +;***************************************************************************
1.539 +end