clamp: biomass/99.all.ncl.1@0c6405ab2ff4

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

     2 ; combine scatter, histogram, global and zonal plots

     3 ; compute all correlation coef and M score

     4 ; ****************************************************

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

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

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

     8 ;************************************************

     9 begin

    11  plot_type     = "ps"

    12  plot_type_new = "png"

    13 ;************************************************

    14 ; read in model data

    15 ;************************************************

    16 ;film = "b30.061n_1995-2004_ANN_climo_lnd.nc"

    17 ;model_name = "b30.061n"

    18 ;model_grid = "T31"

    20 ;film = "newcn05_ncep_1i_ANN_climo_lnd.nc"

    21 ;model_name = "newcn"

    22 ;model_grid = "1.9"

    23 ;--------------------------------------------

    24 ;film = "i01.10cn_1948-2004_ANN_climo.nc"

    25 ;model_name = "10cn"

    26 ;model_grid = "T42"

    28 ;dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

    29 ;fm   = addfile (dirm+film,"r")

    31 ;unit: gC/m2

    32 ;data1  = fm->LIVESTEMC

    33 ;data2  = fm->DEADSTEMC

    34 ;data3  = fm->LEAFC

    35 ;datamod0 = data1(0,:,:)

    36 ;datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)

    37 ;-------------------------------------------

    38  film = "i01.10casa_1948-2004_ANN_climo.nc"

    39  model_name = "10casa"

    40  model_grid = "T42"

    42  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

    43  fm   = addfile (dirm+film,"r")

    45 ;unit: gC/m2

    46  factor_WOODC = 0.7

    47  data1  = fm->WOODC

    48  data2  = fm->LEAFC

    49  datamod0 = data1(0,:,:)

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

    51 ;----------------------------------------------------

    53  system("sed s#model_name#"+model_name+"# table.html > table.html.new")

    54  system("mv -f table.html.new table.html")

    56  xm       = fm->lon

    57  ym       = fm->lat

    58 ;************************************************

    59 ; read in ob data

    60 ;************************************************

    61  ob_name = "LC15_Amazon_Biomass"

    63  diro = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"

    64  filo = "amazon_biomass_"+model_grid+".nc"

    65  fo   = addfile (diro+filo,"r")

    67  dataob   = fo->BIOMASS

    68  xo       = fo->lon

    69  yo       = fo->lat

    70 ;************************************************

    71 ; Units for these variables are:

    72 ; dataob   : MgC/ha

    73 ; datamod0 : gC/m2

    74 ; We want to convert these to KgC/m2

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

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

    78   factor_aboveground = 0.5

    79   factor_unit_ob     = 0.1

    80   factor_unit_mod    = 0.001

    82   dataob   = dataob * factor_aboveground * factor_unit_ob

    83   datamod0 = datamod0 * factor_unit_mod

    85   dataob@units      = "KgC/m2"

    86   datamod0@units    = "KgC/m2"

    88   dataob@long_name      = "Amazon Biomass"

    89   datamod0@long_name    = "Amazon Biomass"

    90 ;********************************************************

    91 ; get subset of datamod0 that match dataob

    93   nlon = dimsizes(xo)

    94   nlat = dimsizes(yo)

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

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

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

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

   101   datamod  = dataob

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

   104 ; print (datamod)

   106 ;----------------------------------------------------------------------

   107 ; global contour

   109 ;res

   110   resg                     = True             ; Use plot options

   111   resg@cnFillOn            = True             ; Turn on color fill

   112   resg@gsnSpreadColors     = True             ; use full colormap

   113 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

   114 ; resg@lbLabelAutoStride   = True

   115   resg@cnLinesOn           = False            ; Turn off contourn lines

   116   resg@mpFillOn            = False            ; Turn off map fill

   117   resg@gsnAddCyclic        = False

   119   resg@gsnSpreadColors      = True            ; use full colormap

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

   121   resg@cnMinLevelValF       = 0.              ; Min level

   122   resg@cnMaxLevelValF       = 30.             ; Max level

   123   resg@cnLevelSpacingF      = 2.              ; interval

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

   126   resg@mpMaxLatF            =  13.8

   127   resg@mpMinLonF            =  277.28

   128   resg@mpMaxLonF            =  326.43

   129 ;------------------------------------------------------------------------

   130 ;global contour ob

   132   plot_name = "global_ob"

   133   title     = ob_name+" "+ model_grid

   134   resg@tiMainString  = title

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

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

   139   plot = gsn_csm_contour_map_ce(wks,dataob,resg)

   140   frame(wks)

   142   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   143 ; system("rm "+plot_name+"."+plot_type)

   144   system("rm "+plot_name+"-1."+plot_type_new)

   145   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   147   clear (wks)

   148 ;------------------------------------------------------------------------

   149 ;global contour model

   151   plot_name = "global_model"

   152   title     = "Model "+ model_name

   153   resg@tiMainString  = title

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

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

   158   plot = gsn_csm_contour_map_ce(wks,datamod,resg)

   159   frame(wks)

   161   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   162 ; system("rm "+plot_name+"."+plot_type)

   163   system("rm "+plot_name+"-1."+plot_type_new)

   164   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   166   clear (wks)

   167 ;------------------------------------------------------------------------

   168 ;global contour model vs ob

   170   plot_name = "global_model_vs_ob"

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

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

   175   delete (plot)

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

   178   resg@gsnFrame             = False          ; Do not draw plot

   179   resg@gsnDraw              = False          ; Do not advance frame

   181 ;(d) compute correlation coef and M score

   183   uu = ndtooned(datamod)

   184   vv = ndtooned(dataob)

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

   188   ug = uu(good)

   189   vg = vv(good)

   191   ccrG = esccr(ug,vg,0)

   192 ; print (ccrG)

   194   MG   = (ccrG*ccrG)* 5.0

   195 ; new eq

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

   197   MG   = (1. - (bias/dimsizes(ug)))*5.

   199   M_biomass = sprintf("%.2f", MG)

   200   system("sed s#M_biomass#"+M_biomass+"# table.html > table.html.new")

   201   system("mv -f table.html.new table.html")

   202   print (M_biomass)

   204 ; plot correlation coef

   206   gRes  = True

   207   gRes@txFontHeightF = 0.02

   208   gRes@txAngleF = 90

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

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

   213 ;--------------------------------------------------------------------

   215 ;(a) ob

   217   title     = ob_name+" "+ model_grid

   218   resg@tiMainString  = title

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

   222 ;(b) model

   224   title     = "Model "+ model_name

   225   resg@tiMainString  = title

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

   229 ;(c) model-ob

   231   zz = datamod

   232   zz = datamod - dataob

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

   235   resg@cnMinLevelValF  = -10.          ; Min level

   236   resg@cnMaxLevelValF  =  10.          ; Max level

   237   resg@cnLevelSpacingF =  2.           ; interval

   238   resg@tiMainString    = title

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

   242   pres                            = True        ; panel plot mods desired

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

   244                                                 ; indiv. plots in panel

   245   pres@gsnMaximize                = True        ; fill the page

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

   249   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

   250 ; system("rm "+plot_name+"."+plot_type)

   251 ; system("rm "+plot_name+"-1."+plot_type_new)

   252 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

   254   frame (wks)

   255   clear (wks)

   256   delete (plot)

   257 ;------------------------------------------------------------------------

   258   temp_name = "biomass." + model_name

   259   system("mkdir -p " + temp_name)

   260   system("cp table.html " + temp_name)

   261   system("mv *.png " + temp_name)

   262   system("tar cf "+ temp_name +".tar " + temp_name)

   263 ;------------------------------------------------------------------------

   264 end

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