clamp: fire/23.table+tseries.ncl@4be95183fbcd

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

     2 ;using model biome vlass

3 ;

     4 ; required command line input parameters:

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

6 ;

     7 ; histogram normalized by rain and compute correleration

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

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

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

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

    12 ;**************************************************************

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

    14 begin

    15 ; add line to ascci/html file

    17   nnewlines = dimsizes(newlines)

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

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

    20     return

    21   end if

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

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

    24   nline = nline + nnewlines

    25   return

    26 end

    27 ;**************************************************************

    28 ; Main code.

    29 begin

    31   plot_type     = "ps"

    32   plot_type_new = "png"

    34 ;---------------------------------------------------

    35 ; model name and grid

    37   model_grid = "T42"

    39   model_name  = "cn"

    40   model_name1 = "i01.06cn"

    41   model_name2 = "i01.10cn"

    43 ;---------------------------------------------------

    44 ; get biome data: model

    46   biome_name_mod = "Model PFT Class"

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

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

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

    52   classmod = fm->CLASS_PFT

    54   delete (fm)

    56 ; model data has 17 land-type classes

    58   nclass_mod = 17

    60 ;--------------------------------------------------

    61 ; get model data: landmask, landfrac and area

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

    64   film = "lnd_T42.nc"

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

    67   landmask = fm->landmask

    68   landfrac = fm->landfrac

    69   area     = fm->area

    71   delete (fm)

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

    74   area = area * 1.e6

    76 ;----------------------------------------------------

    77 ; read data: time series, model

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

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

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

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

    85  delete (fm)

    87 ; Units for these variables are:

    88 ; g C/m^2/s

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

    92   nsec_per_month = 60*60*24*30

    94   data_mod = data_mod * nsec_per_month

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

    97 ;----------------------------------------------------

    98 ; read data: time series, observed

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

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

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

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

   106  delete (fm)

   108  ob_name = "GFEDv2"

   110 ; Units for these variables are:

   111 ; g C/m^2/month

   113 ;---------------------------------------------------

   114 ; take into account landfrac

   116  area     = area * landfrac

   117  data_mod = data_mod * conform(data_mod, landfrac, (/2,3/))

   118  data_ob  = data_ob  * conform(data_ob,  landfrac, (/2,3/))

   120  delete (landfrac)

   122 ;-------------------------------------------------------------

   123 ; html table1 data

   125 ; column (not including header column)

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

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

   129                ,"Correlation Coefficient" \

   130                ,"Ratio model/observed" \

   131                ,"M_score" \

   132                ,"Timeseries plot" \

   133                /)

   135   ncol = dimsizes(col_head)

   137 ; row (not including header row)

   139 ; using model biome class:

   140   row_head  = (/"Not Vegetated" \

   141                ,"Needleleaf Evergreen Temperate Tree" \

   142                ,"Needleleaf Evergreen Boreal Tree" \

   143 ;              ,"Needleleaf Deciduous Boreal Tree" \

   144                ,"Broadleaf Evergreen Tropical Tree" \

   145                ,"Broadleaf Evergreen Temperate Tree" \

   146                ,"Broadleaf Deciduous Tropical Tree" \

   147                ,"Broadleaf Deciduous Temperate Tree" \

   148 ;              ,"Broadleaf Deciduous Boreal Tree" \

   149 ;              ,"Broadleaf Evergreen Shrub" \

   150                ,"Broadleaf Deciduous Temperate Shrub" \

   151                ,"Broadleaf Deciduous Boreal Shrub" \

   152                ,"C3 Arctic Grass" \

   153                ,"C3 Non-Arctic Grass" \

   154                ,"C4 Grass" \

   155                ,"Corn" \

   156 ;              ,"Wheat" \

   157                ,"All Biome" \

   158                /)

   159   nrow = dimsizes(row_head)

   161 ; arrays to be passed to table.

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

   164 ;*****************************************************************

   165 ; (A) get time-mean

   166 ;*****************************************************************

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

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

   170   delete (x)

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

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

   174   delete (x)

   176 ;----------------------------------------------------

   177 ; compute correlation coef

   179   landmask_1d = ndtooned(landmask)

   180   data_mod_1d = ndtooned(data_mod_m)

   181   data_ob_1d  = ndtooned(data_ob_m)

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

   185   cc = esccr(data_mod_1d(good),data_ob_1d(good),0)

   187   delete (landmask_1d)

   188   delete (data_mod_1d)

   189   delete (data_ob_1d)

   190   delete (good)

   192 ;----------------------------------------------------

   193 ; compute M_global

   195   score_max = 1.

   197   Mscore1 = cc * cc * score_max

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

   201 ;----------------------------------------------------

   202 ; global res

   204   resg                      = True             ; Use plot options

   205   resg@cnFillOn             = True             ; Turn on color fill

   206   resg@gsnSpreadColors      = True             ; use full colormap

   207   resg@cnLinesOn            = False            ; Turn off contourn lines

   208   resg@mpFillOn             = False            ; Turn off map fill

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

   211 ;----------------------------------------------------

   212 ; global contour: model vs ob

   214   plot_name = "global_model_vs_ob"

   216   wks = gsn_open_wks (plot_type,plot_name)

   217   gsn_define_colormap(wks,"gui_default")

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

   221   resg@gsnFrame             = False          ; Do not draw plot

   222   resg@gsnDraw              = False          ; Do not advance frame

   224 ;----------------------

   225 ; plot correlation coef

   227   gRes               = True

   228   gRes@txFontHeightF = 0.02

   229   gRes@txAngleF      = 90

   231   correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")"

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

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

   236 ; plot ob

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

   240   title     = ob_name

   241   resg@tiMainString  = title

   243   resg@cnMinLevelValF       = 1.

   244   resg@cnMaxLevelValF       = 10.

   245   resg@cnLevelSpacingF      = 1.

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

   249 ;-----------------------

   250 ; plot model

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

   254   title     = "Model "+ model_name

   255   resg@tiMainString  = title

   257   resg@cnMinLevelValF       = 1.

   258   resg@cnMaxLevelValF       = 10.

   259   resg@cnLevelSpacingF      = 1.

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

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

   264 ; plot model-ob

   266   resg@cnMinLevelValF  = -8.

   267   resg@cnMaxLevelValF  =  2.

   268   resg@cnLevelSpacingF =  1.

   270   zz = data_ob_m

   271   zz = data_mod_m - data_ob_m

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

   273   resg@tiMainString    = title

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

   277 ; plot panel

   279   pres                            = True        ; panel plot mods desired

   280   pres@gsnMaximize                = True        ; fill the page

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

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

   285         "rm "+plot_name+"."+plot_type)

   287   clear (wks)

   288   delete (plot)

   290   delete (data_ob_m)

   291   delete (data_mod_m)

   292   delete (zz)

   294   resg@gsnFrame             = True          ; Do advance frame

   295   resg@gsnDraw              = True          ; Do draw plot

   297 ;*******************************************************************

   298 ; (B) Time series : per biome

   299 ;*******************************************************************

   301  data_n = 2

   303  dsizes = dimsizes(data_mod)

   304  nyear  = dsizes(0)

   305  nmonth = dsizes(1)

   306  ntime  = nyear * nmonth

   308  year_start = 1997

   309  year_end   = 2004

   311 ;-------------------------------------------

   312 ; Calculate "nice" bins for binning the data

   314 ; using model biome class

   315   nclass = nclass_mod

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

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

   320 ; particular range, and then take an average.

   322   nx = dimsizes(range) - 1

   324 ;-------------------------------------------

   325 ; put data into bins

   327 ; using observed biome class

   328 ; base  = ndtooned(classob)

   329 ; using model biome class

   330   base  = ndtooned(classmod)

   332 ; output

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

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

   337 ; Loop through each range, using base.

   339   do i=0,nx-1

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

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

   343      else

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

   345      end if

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

   347 ;    for area

   349      data = ndtooned(area)

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

   352         area_bin(i) = sum(data(idx))

   353      else

   354         area_bin(i) = area_bin@_FillValue

   355      end if

   357 ;#############################################################

   358 ; using model biome class:

   359 ;     set the following 4 classes to _FillValue:

   360 ;     (3)Needleleaf Deciduous Boreal Tree,

   361 ;     (8)Broadleaf Deciduous Boreal Tree,

   362 ;     (9)Broadleaf Evergreen Shrub,

   363 ;     (16)Wheat

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

   366         area_bin(i) = area_bin@_FillValue

   367      end if

   368 ;#############################################################

   370   delete (data)

   372 ;---------------------

   373 ; for data_mod and data_ob

   375   do n = 0,data_n-1

   377      t = -1

   378      do m = 0,nyear-1

   379      do k = 0,nmonth-1

   381         t = t + 1

   383         if (n.eq.0) then

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

   385         end if

   387         if (n.eq.1) then

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

   389         end if

   391 ;       Calculate average

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

   394            yvalues(t,n,i) = avg(data(idx))

   395         else

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

   397         end if

   399 ;#############################################################

   400 ; using model biome class:

   401 ;     set the following 4 classes to _FillValue:

   402 ;     (3)Needleleaf Deciduous Boreal Tree,

   403 ;     (8)Broadleaf Deciduous Boreal Tree,

   404 ;     (9)Broadleaf Evergreen Shrub,

   405 ;     (16)Wheat

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

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

   409         end if

   410 ;#############################################################

   412      end do

   413      end do

   415      delete(data)

   416   end do

   418     delete(idx)

   419   end do

   421   delete (base)

   422   delete (data_mod)

   423   delete (data_ob)

   425 ;----------------------------------------------------------------

   426 ; get area_good

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

   430   area_g = area_bin(good)

   432   n_biome = dimsizes(good)

   434 ;----------------------------------------------------------------

   435 ; data for tseries plot

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

   439   yvalues_g@units = "TgC/month"

   441 ; change unit to Tg C/month

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

   443   factor_unit = 1.e-12

   445   yvalues_g = yvalues(:,:,good) * conform(yvalues_g,area_g,2) * factor_unit

   447   delete (good)

   449 ;-------------------------------------------------------------------

   450 ; general settings for line plot

   452   res                   = True

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

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

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

   457   res@trXMinF   = year_start

   458   res@trXMaxF   = year_end + 1

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

   461 ; res@vpWidthF  = 0.8

   462   res@vpWidthF  = 0.75

   464   res@tiMainFontHeightF = 0.025       ; size of title

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

   468 ; res@gsnMaximize = True

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

   471 ; Add a boxed legend using the simple method

   473   res@pmLegendDisplayMode    = "Always"

   474 ; res@pmLegendWidthF         = 0.1

   475   res@pmLegendWidthF         = 0.08

   476   res@pmLegendHeightF        = 0.06

   477   res@pmLegendOrthogonalPosF = -1.17

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

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

   481 ; res@pmLegendParallelPosF   =  0.18

   482   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   483   res@pmLegendParallelPosF   =  0.73  ;(rightward)

   484   res@pmLegendParallelPosF   =  0.83  ;(rightward)

   486 ; res@lgPerimOn             = False

   487   res@lgLabelFontHeightF     = 0.015

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

   490 ;*******************************************************************

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

   492 ;*******************************************************************

   494 ; x-axis in time series plot

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

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

   498   timeI = ispan(1,ntime,1)

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

   500   timeF@long_name = "year"

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

   503   plot_data@long_name = "TgC/month"

   505 ;----------------------------------------------

   506 ; time series plot : per biome

   508   do m = 0, n_biome-1

   510      plot_name = "monthly_biome_"+ m

   512      wks = gsn_open_wks (plot_type,plot_name)

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

   515      res@tiMainString = title

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

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

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

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

   523             "rm "+plot_name+"."+plot_type)

   525      clear (wks)

   526      delete (plot)

   528   end do

   530 ;------------------------------------------

   531 ; data for table : per biome

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

   534   unit_factor = 1.e-3

   536   score_max = 1.

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

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

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

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

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

   545   do m = 0, n_biome-1

   547      tmp_ob  = yvalues_g(:,0,m)

   548      tmp_mod = yvalues_g(:,1,m)

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

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

   553      cc = esccr(tmp_mod,tmp_ob,0)

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

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

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

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

   562      delete (good)

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

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

   566      text(m,2) = sprintf("%.2f",cc)

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

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

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

   570   end do

   572   delete (tmp_ob)

   573   delete (tmp_mod)

   575 ;--------------------------------------------

   576 ; time series plot: all biome

   578      plot_name = "monthly_global"

   580      wks = gsn_open_wks (plot_type,plot_name)

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

   583      res@tiMainString = title

   585      do k = 0,ntime-1

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

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

   588      end do

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

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

   593             "rm "+plot_name+"."+plot_type)

   595      clear (wks)

   596      delete (plot)

   598 ;------------------------------------------

   599 ; data for table : global

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

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

   604   cc = esccr(tmp_mod,tmp_ob,0)

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

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

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

   610   text(nrow-1,2) = sprintf("%.2f",cc)

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

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

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

   615 ;**************************************************

   616 ; create html table

   617 ;**************************************************

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

   621   header = (/"<HTML>" \

   622             ,"<HEAD>" \

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

   624             ,"</HEAD>" \

   625             ,header_text \

   626             /)

   627   footer = "</HTML>"

   629   table_header = (/ \

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

   631        ,"<tr>" \

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

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

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

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

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

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

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

   639        ,"</tr>" \

   640        /)

   641   table_footer = "</table>"

   642   row_header = "<tr>"

   643   row_footer = "</tr>"

   645   lines = new(50000,string)

   646   nline = 0

   648   set_line(lines,nline,header)

   649   set_line(lines,nline,table_header)

   650 ;-----------------------------------------------

   651 ;row of table

   653   do n = 0,nrow-1

   654      set_line(lines,nline,row_header)

   656      txt0  = row_head(n)

   657      txt1  = text(n,0)

   658      txt2  = text(n,1)

   659      txt3  = text(n,2)

   660      txt4  = text(n,3)

   661      txt5  = text(n,4)

   662      txt6  = text(n,5)

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

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

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

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

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

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

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

   672      set_line(lines,nline,row_footer)

   673   end do

   674 ;-----------------------------------------------

   675   set_line(lines,nline,table_footer)

   676   set_line(lines,nline,footer)

   678 ; Now write to an HTML file.

   680   output_html = "table_fire.html"

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

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

   684     asciiwrite(output_html,lines(idx))

   685   else

   686    print ("error?")

   687   end if

   689   delete (idx)

   691 end

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