clamp: biomass/99.all.ncl.0@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 = "i01.03cn_1545-1569_ANN_climo.nc"

    17 ;model_name = "i01.03cn"

    18 ;model_grid = "T42"

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

    21 ;model_name = "b30.061n"

    22 ;model_grid = "T31"

    24 ;film = "newcn05_ncep_1i_ANN_climo_lnd.nc"

    25 ;model_name = "newcn"

    26 ;model_grid = "1.9"

    28 ;film = "i01.06cn_1798-2004_ANN_climo.nc"

    29 ;model_name = "06cn"

    30 ;model_grid = "T42"

    32 ;film = "i01.06casa_1798-2004_ANN_climo.nc"

    33 ;model_name = "06casa"

    34 ;model_grid = "T42"

    35 ;note: use 98.all.ncl

    37  film = "i01.10cn_1948-2004_ANN_climo.nc"

    38  model_name = "10cn"

    39  model_grid = "T42"

    41 ;film = "i01.10casa_1948-2004_ANN_climo.nc"

    42 ;model_name = "10casa"

    43 ;model_grid = "T42"

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

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

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

    48 ;unit: gC/m2

    49  data1  = fm->LIVESTEMC

    50  data2  = fm->DEADSTEMC

    51  data3  = fm->LEAFC

    52  datamod0 = data1(0,:,:)

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

    55  xm       = fm->lon

    56  ym       = fm->lat

    57 ;************************************************

    58 ; read in ob data

    59 ;************************************************

    60  ob_name = "LC15_Amazon_Biomass"

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

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

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

    66  dataob   = fo->BIOMASS

    67  xo       = fo->lon

    68  yo       = fo->lat

    69 ;************************************************

    70 ; Units for these variables are:

    71 ; dataob   : MgC/ha

    72 ; datamod0 : gC/m2

    73 ; We want to convert these to KgC/m2

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

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

    77   factor_aboveground = 0.5

    78   factor_unit_ob     = 0.1

    79   factor_unit_mod    = 0.001

    81   dataob   = dataob * factor_aboveground * factor_unit_ob

    82   datamod0 = datamod0 * factor_unit_mod

    84   dataob@units      = "KgC/m2"

    85   datamod0@units    = "KgC/m2"

    87   dataob@long_name      = "Amazon Biomass"

    88   datamod0@long_name    = "Amazon Biomass"

    89 ;********************************************************

    90 ; get subset of datamod0 that match dataob

    92   nlon = dimsizes(xo)

    93   nlat = dimsizes(yo)

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

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

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

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

   100   datamod  = dataob

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

   103 ; print (datamod)

   105 ;----------------------------------------------------------------------

   106 ; global contour

   108 ;res

   109   resg                     = True             ; Use plot options

   110   resg@cnFillOn            = True             ; Turn on color fill

   111   resg@gsnSpreadColors     = True             ; use full colormap

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

   113 ; resg@lbLabelAutoStride   = True

   114   resg@cnLinesOn           = False            ; Turn off contourn lines

   115   resg@mpFillOn            = False            ; Turn off map fill

   116   resg@gsnAddCyclic        = False

   118   resg@gsnSpreadColors      = True            ; use full colormap

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

   120   resg@cnMinLevelValF       = 0.              ; Min level

   121   resg@cnMaxLevelValF       = 30.             ; Max level

   122   resg@cnLevelSpacingF      = 2.              ; interval

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

   125   resg@mpMaxLatF            =  13.8

   126   resg@mpMinLonF            =  277.28

   127   resg@mpMaxLonF            =  326.43

   128 ;------------------------------------------------------------------------

   129 ;global contour ob

   131   plot_name = "global_ob"

   132   title     = ob_name+" "+ model_grid

   133   resg@tiMainString  = title

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

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

   138   plot = gsn_csm_contour_map_ce(wks,dataob,resg)

   139   frame(wks)

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

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

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

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

   146   clear (wks)

   147 ;------------------------------------------------------------------------

   148 ;global contour model

   150   plot_name = "global_model"

   151   title     = "Model "+ model_name

   152   resg@tiMainString  = title

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

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

   157   plot = gsn_csm_contour_map_ce(wks,datamod,resg)

   158   frame(wks)

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

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

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

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

   165   clear (wks)

   166 ;------------------------------------------------------------------------

   167 ;global contour model vs ob

   169   plot_name = "global_model_vs_ob"

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

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

   174   delete (plot)

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

   177   resg@gsnFrame             = False          ; Do not draw plot

   178   resg@gsnDraw              = False          ; Do not advance frame

   180 ;(d) compute correlation coef and M score

   182   uu = ndtooned(datamod)

   183   vv = ndtooned(dataob)

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

   187   ug = uu(good)

   188   vg = vv(good)

   190   ccrG = esccr(ug,vg,0)

   191 ; print (ccrG)

   193   MG   = (ccrG*ccrG)* 5.0

   194 ; new eq

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

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

   197   print (MG)

   199 ; plot correlation coef

   201   gRes  = True

   202   gRes@txFontHeightF = 0.02

   203   gRes@txAngleF = 90

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

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

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

   210 ;(a) ob

   212   title     = ob_name+" "+ model_grid

   213   resg@tiMainString  = title

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

   217 ;(b) model

   219   title     = "Model "+ model_name

   220   resg@tiMainString  = title

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

   224 ;(c) model-ob

   226   zz = datamod

   227   zz = datamod - dataob

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

   230   resg@cnMinLevelValF  = -10.          ; Min level

   231   resg@cnMaxLevelValF  =  10.          ; Max level

   232   resg@cnLevelSpacingF =  2.           ; interval

   233   resg@tiMainString    = title

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

   237   pres                            = True        ; panel plot mods desired

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

   239                                                 ; indiv. plots in panel

   240   pres@gsnMaximize                = True        ; fill the page

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

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

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

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

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

   249   frame (wks)

   250   clear (wks)

   251   delete (plot)

   252 ;------------------------------------------------------------------------

   253 ; temp_name = "temp." + model_name

   254 ; system("mkdir -p " + temp_name)

   255 ; system("mv *.png " + temp_name)

   256 ; system("tar cf "+ temp_name +".tar " + temp_name)

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

   258 end

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