clamp: all/04.biomass.ncl@4be95183fbcd

     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+film11,"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      = "kg C m~S~-2~N~"

   114   datamod0@units    = "kg C m~S~-2~N~"

   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)

   207   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   209 ;------------------------------------------------------------------------

   210 ; contour ob

   212   resg@cnMinLevelValF       = 0.              ; Min level

   213   resg@cnMaxLevelValF       = 30.             ; Max level

   214   resg@cnLevelSpacingF      = 2.              ; interval

   216   plot_name = "global_ob"

   217   title     = ob_name+" "+ model_grid

   218   resg@tiMainString  = title

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

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

   223   plot = gsn_csm_contour_map_ce(wks,dataob,resg)

   225   delete (wks)

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

   228          "rm "+plot_name+"."+plot_type)

   229   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   231 ;------------------------------------------------------------------------

   232 ; contour model

   234   resg@cnMinLevelValF       = 0.              ; Min level

   235   resg@cnMaxLevelValF       = 30.             ; Max level

   236   resg@cnLevelSpacingF      = 2.              ; interval

   238   plot_name = "global_model"

   239   title     = "Model "+ model_name + " (2000)"

   240   resg@tiMainString  = title

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

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

   245   plot = gsn_csm_contour_map_ce(wks,datamod,resg)

   247   delete (wks)

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

   250          "rm "+plot_name+"."+plot_type)

   251   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   253 ;------------------------------------------------------------------------

   254 ; contour model vs ob

   256   plot_name = "global_model_vs_ob"

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

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

   261   delete (plot)

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

   264   resg@gsnFrame             = False          ; Do not draw plot

   265   resg@gsnDraw              = False          ; Do not advance frame

   267 ;(d) compute correlation coef and M score

   269   uu = ndtooned(datamod)

   270   vv = ndtooned(dataob)

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

   274   ug = uu(good)

   275   vg = vv(good)

   277   ccrG = esccr(ug,vg,0)

   279   score_max = 10.0

   281 ; Miomass = (ccrG*ccrG)* score_max

   282 ; new eq

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

   284   print ("bias=" + bias + "dimsizes(ug)=" + dimsizes(ug))

   285   print ("other bias=" + sum(abs(ug-vg)/(ug+vg)))

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

   287   print ("Mbiomass=" + Mbiomass)

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

   289   print ("M_biomass=" + M_biomass)

   291   if (isvar("compare")) then

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

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

   294   end if

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

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

   299 ; plot correlation coef

   301   gRes  = True

   302   gRes@txFontHeightF = 0.02

   303   gRes@txAngleF = 90

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

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

   308 ;--------------------------------------------------------------------

   310 ;(a) ob

   312   title     = ob_name+" "+ model_grid

   313   resg@tiMainString  = title

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

   317 ;(b) model

   319   title     = "Model "+ model_name + " (2000)"

   320   resg@tiMainString  = title

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

   324 ;(c) model-ob

   326   zz = datamod

   327   zz = datamod - dataob

   328   title = "Model "+model_name+" (2000) - Observed"

   330   resg@cnMinLevelValF  = -10.          ; Min level

   331   resg@cnMaxLevelValF  =  10.          ; Max level

   332   resg@cnLevelSpacingF =  2.           ; interval

   333   resg@tiMainString    = title

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

   337   pres                            = True        ; panel plot mods desired

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

   339                                                 ; indiv. plots in panel

   340   pres@gsnMaximize                = True        ; fill the page

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

   344   delete (wks)

   345   delete (plot)

   346   delete (zz)

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

   349          "rm "+plot_name+"."+plot_type)

   350   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   352   resg@gsnFrame             = True            ; draw plot

   353   resg@gsnDraw              = True            ; advance frame

   354 ;------------------------------------------------------------------------

   355 ; contour ob : masked

   357   resg@cnMinLevelValF       = 0.              ; Min level

   358   resg@cnMaxLevelValF       = 30.             ; Max level

   359   resg@cnLevelSpacingF      = 2.              ; interval

   361   plot_name = "global_mask_ob"

   362   title     = ob_name+" "+ model_grid

   363   resg@tiMainString  = title

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

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

   367 ;-----------------------------------------

   368 ; plot average over mask region

   370   gRes  = True

   371   gRes@txFontHeightF = 0.02

   372   gRes@txAngleF = 0

   374   area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_ob)+" kg C m~S~-2~N~)"

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

   377 ;-----------------------------------------

   378   zo = dataob

   379   zo = dataob*mask_amazon

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

   381   plot = gsn_csm_contour_map_ce(wks,zo,resg)

   383   delete (wks)

   384   delete (plot)

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

   387          "rm "+plot_name+"."+plot_type)

   388   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   390 ;------------------------------------------------------------------------

   391 ; contour model: masked

   393   resg@cnMinLevelValF       = 0.              ; Min level

   394   resg@cnMaxLevelValF       = 30.             ; Max level

   395   resg@cnLevelSpacingF      = 2.              ; interval

   397   plot_name = "global_mask_model"

   398   title     = "Model "+ model_name + " (2000)"

   399   resg@tiMainString  = title

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

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

   403 ;-----------------------------------------

   404 ; plot average over mask region

   406   gRes  = True

   407   gRes@txFontHeightF = 0.02

   408   gRes@txAngleF = 0

   410   area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_mod)+" kg C m~S~-2~N~)"

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

   413 ;-----------------------------------------

   414   zm = datamod

   415   zm = datamod*mask_amazon

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

   417   plot = gsn_csm_contour_map_ce(wks,zm,resg)

   419   delete (wks)

   420   delete (plot)

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

   423          "rm "+plot_name+"."+plot_type)

   424   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   426 ;------------------------------------------------------------------------

   427 ; contour model vs ob: masked

   429   plot_name = "global_mask_vs_ob"

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

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

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

   436   resg@gsnFrame             = False          ; Do not draw plot

   437   resg@gsnDraw              = False          ; Do not advance frame

   439 ;(d) compute correlation coef and M score

   441   delete (good)

   442   delete (uu)

   443   delete (vv)

   444   delete (ug)

   445   delete (vg)

   447   score_max = 0.

   449   uu = ndtooned(zm)

   450   vv = ndtooned(zo)

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

   454   ug = uu(good)

   455   vg = vv(good)

   457   ccrG = esccr(ug,vg,0)

   459 ; Miomass = (ccrG*ccrG)*score_max

   460 ; new eq

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

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

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

   465   if (isvar("compare")) then

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

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

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

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

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

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

   472   end if

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

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

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

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

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

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

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

   481 ; plot correlation coef

   483   gRes  = True

   484   gRes@txFontHeightF = 0.02

   485   gRes@txAngleF = 90

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

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

   490 ;--------------------------------------------------------------------

   492 ;(a) ob

   494   title     = ob_name+" "+ model_grid

   495   resg@tiMainString  = title

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

   499 ;(b) model

   501   title     = "Model "+ model_name + " (2000)"

   502   resg@tiMainString  = title

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

   506 ;(c) model-ob

   508   zz = zo

   509   zz = zm - zo

   510   title = "Model "+model_name+" (2000) - Observed"

   512   resg@cnMinLevelValF  = -10.          ; Min level

   513   resg@cnMaxLevelValF  =  10.          ; Max level

   514   resg@cnLevelSpacingF =  2.           ; interval

   515   resg@tiMainString    = title

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

   519   pres                            = True        ; panel plot mods desired

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

   521                                                 ; indiv. plots in panel

   522   pres@gsnMaximize                = True        ; fill the page

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

   526   delete (wks)

   527   delete (plot)

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

   530          "rm "+plot_name+"."+plot_type)

   531   ;system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   533 ;***************************************************************************

   534 ; add total score and write to file

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

   536   M_total = Mbiomass + Mbiomass2

   538   asciiwrite("M_save.biomass", M_total)

   540 ;***************************************************************************

   541 ; output plots

   542 ;***************************************************************************

   543   output_dir = model_name+"/biomass"

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

   546 ;***************************************************************************

   547 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Thu, 26 Mar 2009 14:02:21 -0400
changeset 1	4be95183fbcd
parent 0	0c6405ab2ff4
permissions	-rw-r--r--