clamp: all/04.biomass.ncl@0c6405ab2ff4

     1 ; ****************************************************

     2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

     3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

     4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

     5 ;************************************************

     6 begin

     8  plot_type     = "ps"

     9  plot_type_new = "png"

    11 ;------------------------------------------------------

    12 ; edit table.html of current model for movel1_vs_model2

    14  if (isvar("compare")) then

    15     html_name2 = compare+"/table.html"

    16     html_new2  = html_name2 +".new"

    17  end if

    19 ;------------------------------------------------------

    20 ; edit table.html for current model

    22  html_name = model_name+"/table.html"

    23  html_new  = html_name +".new"

    25 ;---------------------------------

    26 ; get model data: landfrac and area

    28   film_l   = "lnd_"+ model_grid + ".nc"

    29   fm_l     = addfile (dirs+film_l,"r")

    31   landfrac = fm_l->landfrac

    32   area0    = fm_l->area

    34   delete (fm_l)

    36 ; change area from km**2 to m**2

    37   area0 = area0 * 1.e6

    39 ;-----------------------------

    40 ; take into account landfrac

    42   area0     = area0 * landfrac

    44 ;--------------------------------------------

    45 ; read model data

    47  fm   = addfile (dirm+film1,"r")

    49  if (BGC .eq. "cn") then

    50     data1  = fm->LIVESTEMC

    51     data2  = fm->DEADSTEMC

    52     data3  = fm->LEAFC

    53     datamod0 = data1(0,:,:)

    54     datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)

    55  end if

    57  if (BGC .eq. "casa") then

    58     factor_WOODC = 0.7

    59     data1  = fm->WOODC

    60     data2  = fm->LEAFC

    61     datamod0 = data1(0,:,:)

    62     datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:)

    63  end if

    65 ; unit: gC/m2

    67  xm       = fm->lon

    68  ym       = fm->lat

    70  delete (fm)

    72 ;------------------------------------------------

    73 ; read amazon mask data

    75   dir_m = diro + "biomass/"

    76   fil_m = "amazon_mask_"+ model_grid + ".nc"

    77   fm    = addfile (dir_m+fil_m,"r")

    79   mask_amazon0 = fm->mask_amazon

    81   delete (fm)

    83 ;------------------------------------------------

    84 ; read ob data

    86  ob_name = "LC15_Amazon_Biomass"

    88  dir_b = diro + "biomass/"

    89  fil_b = "amazon_biomass_"+model_grid+".nc"

    90  fo   = addfile (dir_b+fil_b,"r")

    92  dataob   = fo->BIOMASS

    93  xo       = fo->lon

    94  yo       = fo->lat

    96  delete (fo)

    98 ;************************************************

    99 ; Units for these variables are:

   100 ; dataob   : MgC/ha

   101 ; datamod0 : gC/m2

   102 ; We want to convert these to KgC/m2

   103 ; ha = 100m*100m = 10,000 m2

   104 ; MgC/ha*1000/10,000 = KgC/m2

   106   factor_aboveground = 0.5

   107   factor_unit_ob     = 0.1

   108   factor_unit_mod    = 0.001

   110   dataob   = dataob * factor_aboveground * factor_unit_ob

   111   datamod0 = datamod0 * factor_unit_mod

   113   dataob@units      = "KgC/m2"

   114   datamod0@units    = "KgC/m2"

   116   dataob@long_name      = "Amazon Biomass"

   117   datamod0@long_name    = "Amazon Biomass"

   118 ;********************************************************

   119 ; get subset of datamod0 that match dataob

   121   nlon = dimsizes(xo)

   122   nlat = dimsizes(yo)

   124   ind_lonL = ind(xm .eq. xo(0))

   125   ind_lonR = ind(xm .eq. xo(nlon-1))

   126   ind_latS = ind(ym .eq. yo(0))

   127   ind_latN = ind(ym .eq. yo(nlat-1))

   129   datamod  = dataob

   130   datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)

   132   area  = dataob

   133   area(:,:) = area0(ind_latS:ind_latN,ind_lonL:ind_lonR)

   135   mask_amazon  = dataob

   136   mask_amazon(:,:) = mask_amazon0(ind_latS:ind_latN,ind_lonL:ind_lonR)

   138   mask_amazon@units = ""

   139 ;********************************************************

   140 ; sum over amazom_mask area:

   142 ; Peta g = 1.e15 g = 1.e12 Kg

   143   factor_unit = 1.e-12

   145 ; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.)

   147   Sum_area = sum(area*mask_amazon)*factor_unit

   149   Sum_ob  = sum(dataob*area*mask_amazon)*factor_unit

   150   Sum_mod = sum(datamod*area*mask_amazon)*factor_unit

   152   avg_ob  = Sum_ob /Sum_area

   153   avg_mod = Sum_mod/Sum_area

   155   Sum_biomass_ob  = sprintf("%.2f",Sum_ob )

   156   Sum_biomass_mod = sprintf("%.2f",Sum_mod)

   158 ;----------------------------------------------------------------------

   159 ; contour plot res

   161   resg                     = True             ; Use plot options

   162   resg@cnFillOn            = True             ; Turn on color fill

   163   resg@gsnSpreadColors     = True             ; use full colormap

   164   resg@cnLinesOn           = False            ; Turn off contourn lines

   165   resg@mpFillOn            = False            ; Turn off map fill

   166   resg@gsnAddCyclic        = False

   168   resg@gsnSpreadColors      = True            ; use full colormap

   169   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

   171   resg@mpMinLatF            = -21.1      ; range to zoom in on

   172   resg@mpMaxLatF            =  13.8

   173   resg@mpMinLonF            =  277.28

   174   resg@mpMaxLonF            =  326.43

   175 ;------------------------------------------------------------------------

   176 ; mask plot

   178   plot_name = "mask_ob"

   180   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   181   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   183 ;-----------------------------------------

   184 ; plot area sum

   186   gRes  = True

   187   gRes@txFontHeightF = 0.02

   188 ; gRes@txAngleF = 90

   190   area_sum_text = "(mask area = "+sprintf("%.2f", Sum_area)+"e12 m2)"

   192   gsn_text_ndc(wks,area_sum_text,0.50,0.80,gRes)

   193 ;-----------------------------------------

   195   resg@cnMinLevelValF      = 0.              ; Min level

   196   resg@cnMaxLevelValF      = 1.              ; Max level

   197   resg@cnLevelSpacingF     = 0.1             ; interval

   199   resg@tiMainString        = "Amazon Mask: grid = "+ model_grid

   201   plot = gsn_csm_contour_map_ce(wks,mask_amazon,resg)

   203   delete (wks)

   205   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   206          "rm "+plot_name+"."+plot_type)

   208 ;------------------------------------------------------------------------

   209 ; contour ob

   211   resg@cnMinLevelValF       = 0.              ; Min level

   212   resg@cnMaxLevelValF       = 30.             ; Max level

   213   resg@cnLevelSpacingF      = 2.              ; interval

   215   plot_name = "global_ob"

   216   title     = ob_name+" "+ model_grid

   217   resg@tiMainString  = title

   219   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   220   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   222   plot = gsn_csm_contour_map_ce(wks,dataob,resg)

   224   delete (wks)

   226   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   227          "rm "+plot_name+"."+plot_type)

   229 ;------------------------------------------------------------------------

   230 ; contour model

   232   resg@cnMinLevelValF       = 0.              ; Min level

   233   resg@cnMaxLevelValF       = 30.             ; Max level

   234   resg@cnLevelSpacingF      = 2.              ; interval

   236   plot_name = "global_model"

   237   title     = "Model "+ model_name

   238   resg@tiMainString  = title

   240   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   241   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   243   plot = gsn_csm_contour_map_ce(wks,datamod,resg)

   245   delete (wks)

   247   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   248          "rm "+plot_name+"."+plot_type)

   250 ;------------------------------------------------------------------------

   251 ; contour model vs ob

   253   plot_name = "global_model_vs_ob"

   255   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   256   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   258   delete (plot)

   259   plot=new(3,graphic)                        ; create graphic array

   261   resg@gsnFrame             = False          ; Do not draw plot

   262   resg@gsnDraw              = False          ; Do not advance frame

   264 ;(d) compute correlation coef and M score

   266   uu = ndtooned(datamod)

   267   vv = ndtooned(dataob)

   269   good = ind(.not.ismissing(uu) .and. .not.ismissing(vv))

   271   ug = uu(good)

   272   vg = vv(good)

   274   ccrG = esccr(ug,vg,0)

   276   score_max = 5.0

   278 ; Miomass = (ccrG*ccrG)* score_max

   279 ; new eq

   280   bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))

   281   Mbiomass  = (1. - (bias/dimsizes(ug)))*score_max

   282   M_biomass = sprintf("%.2f", Mbiomass)

   284   if (isvar("compare")) then

   285      system("sed -e '1,/M_biomass/s/M_biomass/"+M_biomass+"/' "+html_name2+" > "+html_new2+";"+ \

   286             "mv -f "+html_new2+" "+html_name2)

   287   end if

   289   system("sed s#M_biomass#"+M_biomass+"# "+html_name+" > "+html_new+";"+ \

   290          "mv -f "+html_new+" "+html_name)

   292 ; plot correlation coef

   294   gRes  = True

   295   gRes@txFontHeightF = 0.02

   296   gRes@txAngleF = 90

   298   correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"

   300   gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)

   301 ;--------------------------------------------------------------------

   303 ;(a) ob

   305   title     = ob_name+" "+ model_grid

   306   resg@tiMainString  = title

   308   plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg)

   310 ;(b) model

   312   title     = "Model "+ model_name

   313   resg@tiMainString  = title

   315   plot(1) = gsn_csm_contour_map_ce(wks,datamod,resg)

   317 ;(c) model-ob

   319   zz = datamod

   320   zz = datamod - dataob

   321   title = "Model_"+model_name+" - Observed"

   323   resg@cnMinLevelValF  = -10.          ; Min level

   324   resg@cnMaxLevelValF  =  10.          ; Max level

   325   resg@cnLevelSpacingF =  2.           ; interval

   326   resg@tiMainString    = title

   328   plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)

   330   pres                            = True        ; panel plot mods desired

   331 ; pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

   332                                                 ; indiv. plots in panel

   333   pres@gsnMaximize                = True        ; fill the page

   335   gsn_panel(wks,plot,(/3,1/),pres)              ; create panel plot

   337   delete (wks)

   338   delete (plot)

   339   delete (zz)

   341   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   342          "rm "+plot_name+"."+plot_type)

   344   resg@gsnFrame             = True            ; draw plot

   345   resg@gsnDraw              = True            ; advance frame

   346 ;------------------------------------------------------------------------

   347 ; contour ob : masked

   349   resg@cnMinLevelValF       = 0.              ; Min level

   350   resg@cnMaxLevelValF       = 30.             ; Max level

   351   resg@cnLevelSpacingF      = 2.              ; interval

   353   plot_name = "global_mask_ob"

   354   title     = ob_name+" "+ model_grid

   355   resg@tiMainString  = title

   357   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   358   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   359 ;-----------------------------------------

   360 ; plot average over mask region

   362   gRes  = True

   363   gRes@txFontHeightF = 0.02

   364   gRes@txAngleF = 0

   366   area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_ob)+" Kg C/m2)"

   368   gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)

   369 ;-----------------------------------------

   370   zo = dataob

   371   zo = dataob*mask_amazon

   372   zo = where((mask_amazon .le. 0.01), zo@_FillValue, zo)

   373   plot = gsn_csm_contour_map_ce(wks,zo,resg)

   375   delete (wks)

   376   delete (plot)

   378   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   379          "rm "+plot_name+"."+plot_type)

   381 ;------------------------------------------------------------------------

   382 ; contour model: masked

   384   resg@cnMinLevelValF       = 0.              ; Min level

   385   resg@cnMaxLevelValF       = 30.             ; Max level

   386   resg@cnLevelSpacingF      = 2.              ; interval

   388   plot_name = "global_mask_model"

   389   title     = "Model "+ model_name

   390   resg@tiMainString  = title

   392   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   393   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   394 ;-----------------------------------------

   395 ; plot average over mask region

   397   gRes  = True

   398   gRes@txFontHeightF = 0.02

   399   gRes@txAngleF = 0

   401   area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_mod)+" Kg C/m2)"

   403   gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes)

   404 ;-----------------------------------------

   405   zm = datamod

   406   zm = datamod*mask_amazon

   407   zm = where((mask_amazon .le. 0.01), zm@_FillValue, zm)

   408   plot = gsn_csm_contour_map_ce(wks,zm,resg)

   410   delete (wks)

   411   delete (plot)

   413   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   414          "rm "+plot_name+"."+plot_type)

   416 ;------------------------------------------------------------------------

   417 ; contour model vs ob: masked

   419   plot_name = "global_mask_vs_ob"

   421   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   422   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   424   plot=new(3,graphic)                        ; create graphic array

   426   resg@gsnFrame             = False          ; Do not draw plot

   427   resg@gsnDraw              = False          ; Do not advance frame

   429 ;(d) compute correlation coef and M score

   431   delete (good)

   432   delete (uu)

   433   delete (vv)

   434   delete (ug)

   435   delete (vg)

   437   score_max = 5.

   439   uu = ndtooned(zm)

   440   vv = ndtooned(zo)

   442   good = ind((uu .gt. 0.) .and. (vv .gt. 0.))

   444   ug = uu(good)

   445   vg = vv(good)

   447   ccrG = esccr(ug,vg,0)

   449 ; Miomass = (ccrG*ccrG)*score_max

   450 ; new eq

   451   bias = sum(abs(ug-vg)/(abs(ug)+abs(vg)))

   452   Mbiomass2  = (1. - (bias/dimsizes(ug)))*score_max

   453   M_biomask = sprintf("%.2f", Mbiomass2)

   455   if (isvar("compare")) then

   456      system("sed -e '1,/M_biomask/s/M_biomask/"+M_biomask+"/' "+html_name2+" > "+html_new2+";"+ \

   457             "mv -f "+html_new2+" "+html_name2)

   458      system("sed -e '1,/Sum_biomass_ob/s/Sum_biomass_ob/"+Sum_biomass_ob+"/' "+html_name2+" > "+html_new2+";"+ \

   459             "mv -f "+html_new2+" "+html_name2)

   460      system("sed -e '1,/Sum_biomass_mod/s/Sum_biomass_mod/"+Sum_biomass_mod+"/' "+html_name2+" > "+html_new2+";"+ \

   461             "mv -f "+html_new2+" "+html_name2)

   462   end if

   464   system("sed s#M_biomask#"+M_biomask+"# "+html_name+" > "+html_new+";"+ \

   465          "mv -f "+html_new+" "+html_name)

   466   system("sed s#Sum_biomass_ob#"+Sum_biomass_ob+"# "+html_name+" > "+html_new+";"+ \

   467          "mv -f "+html_new+" "+html_name)

   468   system("sed s#Sum_biomass_mod#"+Sum_biomass_mod+"# "+html_name+" > "+html_new+";"+ \

   469          "mv -f "+html_new+" "+html_name)

   470 ;--------------------------------------------------------------------

   471 ; plot correlation coef

   473   gRes  = True

   474   gRes@txFontHeightF = 0.02

   475   gRes@txAngleF = 90

   477   correlation_text = "(correlation coef = "+sprintf("%.2f", ccrG)+")"

   479   gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)

   480 ;--------------------------------------------------------------------

   482 ;(a) ob

   484   title     = ob_name+" "+ model_grid

   485   resg@tiMainString  = title

   487   plot(0) = gsn_csm_contour_map_ce(wks,zo,resg)

   489 ;(b) model

   491   title     = "Model "+ model_name

   492   resg@tiMainString  = title

   494   plot(1) = gsn_csm_contour_map_ce(wks,zm,resg)

   496 ;(c) model-ob

   498   zz = zo

   499   zz = zm - zo

   500   title = "Model_"+model_name+" - Observed"

   502   resg@cnMinLevelValF  = -10.          ; Min level

   503   resg@cnMaxLevelValF  =  10.          ; Max level

   504   resg@cnLevelSpacingF =  2.           ; interval

   505   resg@tiMainString    = title

   507   plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)

   509   pres                            = True        ; panel plot mods desired

   510 ; pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

   511                                                 ; indiv. plots in panel

   512   pres@gsnMaximize                = True        ; fill the page

   514   gsn_panel(wks,plot,(/3,1/),pres)              ; create panel plot

   516   delete (wks)

   517   delete (plot)

   519   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \

   520          "rm "+plot_name+"."+plot_type)

   522 ;***************************************************************************

   523 ; add total score and write to file

   524 ;***************************************************************************

   525   M_total = Mbiomass + Mbiomass2

   527   asciiwrite("M_save.biomass", M_total)

   529 ;***************************************************************************

   530 ; output plots

   531 ;***************************************************************************

   532   output_dir = model_name+"/biomass"

   534   system("mv *.png " + output_dir)

   535 ;***************************************************************************

   536 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Mon, 26 Jan 2009 22:08:20 -0500
changeset 0	0c6405ab2ff4
child 1	4be95183fbcd
permissions	-rw-r--r--