clamp: fire/24.table+tseries.ncl@0c6405ab2ff4

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

     2 ; landfrac applied to area only.

     3 ; using model biome class

4 ;

     5 ; required command line input parameters:

     6 ;  ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl

7 ;

     8 ; histogram normalized by rain and compute correleration

     9 ;**************************************************************

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

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

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

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

    14 procedure set_line(lines:string,nline:integer,newlines:string)

    15 begin

    16 ; add line to ascci/html file

    18   nnewlines = dimsizes(newlines)

    19   if(nline+nnewlines-1.ge.dimsizes(lines))

    20     print("set_line: bad index, not setting anything.")

    21     return

    22   end if

    23   lines(nline:nline+nnewlines-1) = newlines

    24 ;  print ("lines = " + lines(nline:nline+nnewlines-1))

    25   nline = nline + nnewlines

    26   return

    27 end

    28 ;**************************************************************

    29 ; Main code.

    30 begin

    32   plot_type     = "ps"

    33   plot_type_new = "png"

    35 ;---------------------------------------------------

    36 ; model name and grid

    38   model_grid = "T42"

    40   model_name  = "cn"

    41   model_name1 = "i01.06cn"

    42   model_name2 = "i01.10cn"

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

    45 ; get biome data: model

    47   biome_name_mod = "Model PFT Class"

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

    50   film = "class_pft_"+model_grid+".nc"

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

    53   classmod = fm->CLASS_PFT

    55   delete (fm)

    57 ; model data has 17 land-type classes

    59   nclass_mod = 17

    61 ;--------------------------------------------------

    62 ; get model data: landmask, landfrac and area

    64   dirm = "/fis/cgd/cseg/people/jeff/surface_data/"

    65   film = "lnd_T42.nc"

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

    68   landmask = fm->landmask

    69   landfrac = fm->landfrac

    70   area     = fm->area

    72   delete (fm)

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

    75   area = area * 1.e6

    77 ;---------------------------------------------------

    78 ; take into account landfrac

    80   area     = area * landfrac

    82 ;----------------------------------------------------

    83 ; read data: time series, model

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

    86  film = model_name2 + "_Fire_C_1979-2004_monthly.nc"

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

    89  data_mod = fm->COL_FIRE_CLOSS(18:25,:,:,:)

    91  delete (fm)

    93 ; Units for these variables are:

    94 ; g C/m^2/s

    96 ; change unit to g C/m^2/month

    98   nsec_per_month = 60*60*24*30

   100   data_mod = data_mod * nsec_per_month

   102   data_mod@unit = "gC/m2/month"

   103 ;----------------------------------------------------

   104 ; read data: time series, observed

   106  dirm = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"

   107  film = "Fire_C_1997-2006_monthly_"+ model_grid+".nc"

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

   110  data_ob = fm->FIRE_C(0:7,:,:,:)

   112  delete (fm)

   114  ob_name = "GFEDv2"

   116 ; Units for these variables are:

   117 ; g C/m^2/month

   119 ;-------------------------------------------------------------

   120 ; html table1 data

   122 ; column (not including header column)

   124   col_head  = (/"Observed Fire_Flux (PgC/yr)" \

   125                ,"Model Fire_Flux (PgC/yr)" \

   126                ,"Correlation Coefficient" \

   127                ,"Ratio model/observed" \

   128                ,"M_score" \

   129                ,"Timeseries plot" \

   130                /)

   132   ncol = dimsizes(col_head)

   134 ; row (not including header row)

   136 ; using model biome class:

   137   row_head  = (/"Not Vegetated" \

   138                ,"Needleleaf Evergreen Temperate Tree" \

   139                ,"Needleleaf Evergreen Boreal Tree" \

   140 ;              ,"Needleleaf Deciduous Boreal Tree" \

   141                ,"Broadleaf Evergreen Tropical Tree" \

   142                ,"Broadleaf Evergreen Temperate Tree" \

   143                ,"Broadleaf Deciduous Tropical Tree" \

   144                ,"Broadleaf Deciduous Temperate Tree" \

   145 ;              ,"Broadleaf Deciduous Boreal Tree" \

   146 ;              ,"Broadleaf Evergreen Shrub" \

   147                ,"Broadleaf Deciduous Temperate Shrub" \

   148                ,"Broadleaf Deciduous Boreal Shrub" \

   149                ,"C3 Arctic Grass" \

   150                ,"C3 Non-Arctic Grass" \

   151                ,"C4 Grass" \

   152                ,"Corn" \

   153 ;              ,"Wheat" \

   154                ,"All Biome" \

   155                /)

   156   nrow = dimsizes(row_head)

   158 ; arrays to be passed to table.

   159   text = new ((/nrow, ncol/),string )

   161 ;*****************************************************************

   162 ; (A) get time-mean

   163 ;*****************************************************************

   165   x          = dim_avg_Wrap(data_mod(lat|:,lon|:,month|:,year|:))

   166   data_mod_m = dim_avg_Wrap(       x(lat|:,lon|:,month|:))

   167   delete (x)

   169   x          = dim_avg_Wrap( data_ob(lat|:,lon|:,month|:,year|:))

   170   data_ob_m  = dim_avg_Wrap(       x(lat|:,lon|:,month|:))

   171   delete (x)

   173 ;----------------------------------------------------

   174 ; compute correlation coef: space

   176   landmask_1d = ndtooned(landmask)

   177   data_mod_1d = ndtooned(data_mod_m)

   178   data_ob_1d  = ndtooned(data_ob_m )

   179   area_1d     = ndtooned(area)

   180   landfrac_1d = ndtooned(landfrac)

   182   good = ind(landmask_1d .gt. 0.)

   184   global_mod = sum(data_mod_1d(good)*area_1d(good)) * 1.e-15 * 12.

   185   global_ob  = sum(data_ob_1d(good) *area_1d(good)) * 1.e-15 * 12.

   186   print (global_mod)

   187   print (global_ob)

   189   global_area= sum(area_1d)

   190   global_land= sum(area_1d(good))

   191   print (global_area)

   192   print (global_land)

   194   cc_space = esccr(data_mod_1d(good)*landfrac_1d(good),data_ob_1d(good)*landfrac_1d(good),0)

   196   delete (landmask_1d)

   197   delete (landfrac_1d)

   198 ; delete (area_1d)

   199   delete (data_mod_1d)

   200   delete (data_ob_1d)

   201   delete (good)

   203 ;----------------------------------------------------

   204 ; compute M_global

   206   score_max = 1.

   208   Mscore1 = cc_space * cc_space * score_max

   210   M_global = sprintf("%.2f", Mscore1)

   212 ;----------------------------------------------------

   213 ; global res

   215   resg                      = True             ; Use plot options

   216   resg@cnFillOn             = True             ; Turn on color fill

   217   resg@gsnSpreadColors      = True             ; use full colormap

   218   resg@cnLinesOn            = False            ; Turn off contourn lines

   219   resg@mpFillOn             = False            ; Turn off map fill

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

   222 ;----------------------------------------------------

   223 ; global contour: model vs ob

   225   plot_name = "global_model_vs_ob"

   227   wks = gsn_open_wks (plot_type,plot_name)

   228   gsn_define_colormap(wks,"gui_default")

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

   232   resg@gsnFrame             = False          ; Do not draw plot

   233   resg@gsnDraw              = False          ; Do not advance frame

   235 ;----------------------

   236 ; plot correlation coef

   238   gRes               = True

   239   gRes@txFontHeightF = 0.02

   240   gRes@txAngleF      = 90

   242   correlation_text = "(correlation coef = "+sprintf("%.2f", cc_space)+")"

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

   246 ;-----------------------

   247 ; plot ob

   249   data_ob_m = where(landmask .gt. 0., data_ob_m, data_ob_m@_FillValue)

   251   title     = ob_name

   252   resg@tiMainString  = title

   254   resg@cnMinLevelValF       = 1.

   255   resg@cnMaxLevelValF       = 10.

   256   resg@cnLevelSpacingF      = 1.

   258   plot(0) = gsn_csm_contour_map_ce(wks,data_ob_m,resg)

   260 ;-----------------------

   261 ; plot model

   263   data_mod_m = where(landmask .gt. 0., data_mod_m, data_mod_m@_FillValue)

   265   title     = "Model "+ model_name

   266   resg@tiMainString  = title

   268   resg@cnMinLevelValF       = 1.

   269   resg@cnMaxLevelValF       = 10.

   270   resg@cnLevelSpacingF      = 1.

   272   plot(1) = gsn_csm_contour_map_ce(wks,data_mod_m,resg)

   274 ;-----------------------

   275 ; plot model-ob

   277   resg@cnMinLevelValF  = -8.

   278   resg@cnMaxLevelValF  =  2.

   279   resg@cnLevelSpacingF =  1.

   281   zz = data_ob_m

   282   zz = data_mod_m - data_ob_m

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

   284   resg@tiMainString    = title

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

   288 ; plot panel

   290   pres                            = True        ; panel plot mods desired

   291   pres@gsnMaximize                = True        ; fill the page

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

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

   296         "rm "+plot_name+"."+plot_type)

   298   clear (wks)

   299   delete (plot)

   301   delete (data_ob_m)

   302   delete (data_mod_m)

   303   delete (zz)

   305   resg@gsnFrame             = True          ; Do advance frame

   306   resg@gsnDraw              = True          ; Do draw plot

   308 ;*******************************************************************

   309 ; (B) Time series : per biome

   310 ;*******************************************************************

   312  data_n = 2

   314  dsizes = dimsizes(data_mod)

   315  nyear  = dsizes(0)

   316  nmonth = dsizes(1)

   317  ntime  = nyear * nmonth

   319  year_start = 1997

   320  year_end   = 2004

   322 ;-------------------------------------------

   323 ; Calculate "nice" bins for binning the data

   325 ; using model biome class

   326   nclass = nclass_mod

   328   range  = fspan(0,nclass,nclass+1)

   330 ; Use this range information to grab all the values in a

   331 ; particular range, and then take an average.

   333   nx = dimsizes(range) - 1

   335 ;-------------------------------------------

   336 ; put data into bins

   338 ; using observed biome class

   339 ; base  = ndtooned(classob)

   340 ; using model biome class

   341   base  = ndtooned(classmod)

   343 ; output

   345   area_bin = new((/nx/),float)

   346   yvalues  = new((/ntime,data_n,nx/),float)

   348 ; Loop through each range, using base.

   350   do i=0,nx-1

   352      if (i.ne.(nx-1)) then

   353         idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))

   354      else

   355         idx = ind(base.ge.range(i))

   356      end if

   357 ;---------------------

   358 ;    for area

   360      if (.not.any(ismissing(idx))) then

   361         area_bin(i) = sum(area_1d(idx))

   362      else

   363         area_bin(i) = area_bin@_FillValue

   364      end if

   366 ;#############################################################

   367 ; using model biome class:

   368 ;     set the following 4 classes to _FillValue:

   369 ;     (3)Needleleaf Deciduous Boreal Tree,

   370 ;     (8)Broadleaf Deciduous Boreal Tree,

   371 ;     (9)Broadleaf Evergreen Shrub,

   372 ;     (16)Wheat

   374      if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then

   375         area_bin(i) = area_bin@_FillValue

   376      end if

   377 ;#############################################################

   379 ;---------------------

   380 ; for data_mod and data_ob

   382   do n = 0,data_n-1

   384      t = -1

   385      do m = 0,nyear-1

   386      do k = 0,nmonth-1

   388         t = t + 1

   390         if (n.eq.0) then

   391            data = ndtooned(data_ob(m,k,:,:))

   392         end if

   394         if (n.eq.1) then

   395            data = ndtooned(data_mod(m,k,:,:))

   396         end if

   398 ;       Calculate average

   400         if (.not.any(ismissing(idx))) then

   401            yvalues(t,n,i) = sum(data(idx)*area_1d(idx))

   402         else

   403            yvalues(t,n,i) = yvalues@_FillValue

   404         end if

   406 ;#############################################################

   407 ; using model biome class:

   408 ;     set the following 4 classes to _FillValue:

   409 ;     (3)Needleleaf Deciduous Boreal Tree,

   410 ;     (8)Broadleaf Deciduous Boreal Tree,

   411 ;     (9)Broadleaf Evergreen Shrub,

   412 ;     (16)Wheat

   414         if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then

   415            yvalues(t,n,i) = yvalues@_FillValue

   416         end if

   417 ;#############################################################

   419      end do

   420      end do

   422      delete(data)

   423   end do

   425     delete(idx)

   426   end do

   428   delete (base)

   429   delete (data_mod)

   430   delete (data_ob)

   432   global_bin = sum(area_bin)

   433   print (global_bin)

   435 ;----------------------------------------------------------------

   436 ; get area_good

   438   good = ind(.not.ismissing(area_bin))

   440   area_g = area_bin(good)

   442   n_biome = dimsizes(good)

   444   global_good = sum(area_g)

   445   print (global_good)

   447 ;----------------------------------------------------------------

   448 ; data for tseries plot

   450   yvalues_g = new((/ntime,data_n,n_biome/),float)

   452   yvalues_g@units = "TgC/month"

   454 ; change unit to Tg C/month

   455 ; change unit from g to Tg (Tera gram)

   456   factor_unit = 1.e-12

   458   yvalues_g = yvalues(:,:,good) * factor_unit

   460   delete (good)

   462 ;-------------------------------------------------------------------

   463 ; general settings for line plot

   465   res                   = True

   466   res@xyDashPatterns    = (/0,0/)          ; make lines solid

   467   res@xyLineThicknesses = (/2.0,2.0/)      ; make lines thicker

   468   res@xyLineColors      = (/"blue","red"/) ; line color

   470   res@trXMinF   = year_start

   471   res@trXMaxF   = year_end + 1

   473   res@vpHeightF = 0.4                 ; change aspect ratio of plot

   474 ; res@vpWidthF  = 0.8

   475   res@vpWidthF  = 0.75

   477   res@tiMainFontHeightF = 0.025       ; size of title

   479   res@tmXBFormat  = "f"               ; not to add trailing zeros

   481 ; res@gsnMaximize = True

   483 ;----------------------------------------------

   484 ; Add a boxed legend using the simple method

   486   res@pmLegendDisplayMode    = "Always"

   487 ; res@pmLegendWidthF         = 0.1

   488   res@pmLegendWidthF         = 0.08

   489   res@pmLegendHeightF        = 0.06

   490   res@pmLegendOrthogonalPosF = -1.17

   491 ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)

   492 ; res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   494 ; res@pmLegendParallelPosF   =  0.18

   495   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   496   res@pmLegendParallelPosF   =  0.73  ;(rightward)

   497   res@pmLegendParallelPosF   =  0.83  ;(rightward)

   499 ; res@lgPerimOn             = False

   500   res@lgLabelFontHeightF     = 0.015

   501   res@xyExplicitLegendLabels = (/"observed",model_name/)

   503 ;*******************************************************************

   504 ; (A) time series plot: monthly ( 2 lines per plot)

   505 ;*******************************************************************

   507 ; x-axis in time series plot

   509   timeI = new((/ntime/),integer)

   510   timeF = new((/ntime/),float)

   511   timeI = ispan(1,ntime,1)

   512   timeF = year_start + (timeI-1)/12.

   513   timeF@long_name = "year"

   515   plot_data = new((/2,ntime/),float)

   516   plot_data@long_name = "TgC/month"

   518 ;----------------------------------------------

   519 ; time series plot : per biome

   521   do m = 0, n_biome-1

   523      plot_name = "monthly_biome_"+ m

   525      wks = gsn_open_wks (plot_type,plot_name)

   527      title = "Fire : "+ row_head(m)

   528      res@tiMainString = title

   530      plot_data(0,:) = yvalues_g(:,0,m)

   531      plot_data(1,:) = yvalues_g(:,1,m)

   533      plot = gsn_csm_xy(wks,timeF,plot_data,res)

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

   536             "rm "+plot_name+"."+plot_type)

   538      clear (wks)

   539      delete (plot)

   541   end do

   543 ;------------------------------------------

   544 ; data for table : per biome

   546 ; unit change from TgC/month to PgC/month

   547   unit_factor = 1.e-3

   549   score_max = 1.

   551   tmp_ob    = new((/ntime/),float)

   552   tmp_mod   = new((/ntime/),float)

   554   total_ob  = new((/n_biome/),float)

   555   total_mod = new((/n_biome/),float)

   556   Mscore2   = new((/n_biome/),float)

   558   do m = 0, n_biome-1

   560      tmp_ob  = yvalues_g(:,0,m)

   561      tmp_mod = yvalues_g(:,1,m)

   563      total_ob(m)  = avg(month_to_annual(tmp_ob, 0)) * unit_factor

   564      total_mod(m) = avg(month_to_annual(tmp_mod,0)) * unit_factor

   566      cc_time = esccr(tmp_mod,tmp_ob,0)

   568      ratio = total_mod(m)/total_ob(m)

   570      good = ind(tmp_ob .ne. 0. .and. tmp_mod .ne. 0.)

   572      bias = sum( abs( tmp_mod(good)-tmp_ob(good) )/( abs(tmp_mod(good))+abs(tmp_ob(good)) ) )

   573      Mscore2(m) = (1.- (bias/dimsizes(good)))*score_max

   575      delete (good)

   577      text(m,0) = sprintf("%.2f",total_ob(m))

   578      text(m,1) = sprintf("%.2f",total_mod(m))

   579      text(m,2) = sprintf("%.2f",cc_time)

   580      text(m,3) = sprintf("%.2f",ratio)

   581      text(m,4) = sprintf("%.2f",Mscore2(m))

   582      text(m,5) = "<a href=./monthly_biome_"+m+".png>model_vs_ob</a>"

   583   end do

   585   delete (tmp_ob)

   586   delete (tmp_mod)

   588 ;--------------------------------------------

   589 ; time series plot: all biome

   591      plot_name = "monthly_global"

   593      wks = gsn_open_wks (plot_type,plot_name)

   595      title = "Fire : "+ row_head(n_biome)

   596      res@tiMainString = title

   598      do k = 0,ntime-1

   599         plot_data(0,k) = sum(yvalues_g(k,0,:))

   600         plot_data(1,k) = sum(yvalues_g(k,1,:))

   601      end do

   603      plot = gsn_csm_xy(wks,timeF,plot_data,res)

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

   606             "rm "+plot_name+"."+plot_type)

   608      clear (wks)

   609      delete (plot)

   611 ;------------------------------------------

   612 ; data for table : global

   614   score_max = 1.

   616   tmp_ob  = ndtooned(yvalues_g(:,0,:))

   617   tmp_mod = ndtooned(yvalues_g(:,1,:))

   619   cc_time = esccr(tmp_mod,tmp_ob,0)

   621   ratio = sum(total_mod)/sum(total_ob)

   623   good = ind(tmp_ob .ne. 0. .and. tmp_mod .ne. 0.)

   625   bias = sum( abs( tmp_mod(good)-tmp_ob(good) )/( abs(tmp_mod(good))+abs(tmp_ob(good)) ) )

   626   Mscore3 = (1.- (bias/dimsizes(good)))*score_max

   628   print (Mscore3)

   630   delete (good)

   632   text(nrow-1,0) = sprintf("%.2f",sum(total_ob))

   633   text(nrow-1,1) = sprintf("%.2f",sum(total_mod))

   634   text(nrow-1,2) = sprintf("%.2f",cc_time)

   635   text(nrow-1,3) = sprintf("%.2f",ratio)

   636 ; text(nrow-1,4) = sprintf("%.2f",avg(Mscore2))

   637   text(nrow-1,4) = sprintf("%.2f",    Mscore3)

   638   text(nrow-1,5) = "<a href=./monthly_global.png>model_vs_ob</a>"

   640 ;**************************************************

   641 ; create html table

   642 ;**************************************************

   644   header_text = "<H1>Fire Emission (1997-2004):  Model "+model_name+"</H1>"

   646   header = (/"<HTML>" \

   647             ,"<HEAD>" \

   648             ,"<TITLE>CLAMP metrics</TITLE>" \

   649             ,"</HEAD>" \

   650             ,header_text \

   651             /)

   652   footer = "</HTML>"

   654   table_header = (/ \

   655         "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \

   656        ,"<tr>" \

   657        ,"   <th bgcolor=DDDDDD >Biome Type</th>" \

   658        ,"   <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \

   659        ,"   <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \

   660        ,"   <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \

   661        ,"   <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \

   662        ,"   <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \

   663        ,"   <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \

   664        ,"</tr>" \

   665        /)

   666   table_footer = "</table>"

   667   row_header = "<tr>"

   668   row_footer = "</tr>"

   670   lines = new(50000,string)

   671   nline = 0

   673   set_line(lines,nline,header)

   674   set_line(lines,nline,table_header)

   675 ;-----------------------------------------------

   676 ;row of table

   678   do n = 0,nrow-1

   679      set_line(lines,nline,row_header)

   681      txt0  = row_head(n)

   682      txt1  = text(n,0)

   683      txt2  = text(n,1)

   684      txt3  = text(n,2)

   685      txt4  = text(n,3)

   686      txt5  = text(n,4)

   687      txt6  = text(n,5)

   689      set_line(lines,nline,"<th>"+txt0+"</th>")

   690      set_line(lines,nline,"<th>"+txt1+"</th>")

   691      set_line(lines,nline,"<th>"+txt2+"</th>")

   692      set_line(lines,nline,"<th>"+txt3+"</th>")

   693      set_line(lines,nline,"<th>"+txt4+"</th>")

   694      set_line(lines,nline,"<th>"+txt5+"</th>")

   695      set_line(lines,nline,"<th>"+txt6+"</th>")

   697      set_line(lines,nline,row_footer)

   698   end do

   699 ;-----------------------------------------------

   700   set_line(lines,nline,table_footer)

   701   set_line(lines,nline,footer)

   703 ; Now write to an HTML file.

   705   output_html = "table_fire.html"

   707   idx = ind(.not.ismissing(lines))

   708   if(.not.any(ismissing(idx))) then

   709     asciiwrite(output_html,lines(idx))

   710   else

   711    print ("error?")

   712   end if

   714   delete (idx)

   716 end

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