clamp: co2/99.all.ncl.1@4be95183fbcd

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

     2 ; using gsn_table for all

     3 ; combine 19.metric_plot.ncl and 24.lines.ncl

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

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

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

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

     8 ;load "/fis/cgd/cseg/people/jeff/clamp/co2/metrics_table.ncl"

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

    10 begin

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

    13 ; read model data

    14 ;************************************************

    16 ; from command line inputs

    18 ;--------------------------------------------------

    19 ; edit table.html of current model for movel1_vs_model2

    21  if (isvar("compare")) then

    22     html_name2 = compare+"/table.html"

    23     html_new2  = html_name2 +".new"

    25 ; the following only needs to execute once

    26 ;   system("sed 1,/model_nameA/s//"+model_name+"/ "+html_name2+" > "+html_new2+";"+ \

    27 ;          "mv -f "+html_new2+" "+html_name2+";"+ \

    28 ;          "sed 1,/model_nameB/s//"+model_name+"/ "+html_name2+" > "+html_new2+";"+ \

    29 ;          "mv -f "+html_new2+" "+html_name2+";"+ \

    30 ;          "sed s#"+modeln+"#"+model_name+"# "+html_name2+" > "+html_new2+";"+ \

    31 ;          "mv -f "+html_new2+" "+html_name2)

    32  end if

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

    35 ; edit table.html for current model

    37  html_name = model_name+"/table.html"

    38  html_new  = html_name +".new"

    40 ; the following only needs to execute once

    41 ; system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \

    42 ;        "mv -f "+html_new+" "+html_name)

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

    45   fm    = addfile(dirm+film3,"r")

    47   x     = fm->CO2

    48   xi    = fm->lon

    49   yi    = fm->lat

    51   xdim  = dimsizes(x)

    52   nlev  = xdim(1)

    53   y     = x(:,0,:,:)

    55 ; get co2 at the lowest level

    56   y     = x(:,nlev-1,:,:)

    58 ; change to unit of observed (u mol/mol)

    59 ; Model_units [=] kgCO2 / kgDryAir

    60 ; 28.966 = molecular weight of dry air

    61 ; 44.       = molecular weight of CO2

    62 ; u mol = 1e-6 mol

    64   factor = (28.966/44.) * 1e6

    65   y      = y * factor

    67   y@_FillValue = 1.e36

    68   y@units      = "u mol/mol"

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

    70 ; read data: observed

    71 ;************************************************

    72  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/co2/ob/"

    73  fili  = "co2_globalView_98.nc"

    74  g     = addfile (diri+fili,"r")

    75  val   = g->CO2_SEAS

    76  lon   = g->LON

    77  lat   = g->LAT

    78  sta   = chartostring(g->STATION)

    79  delete (g)

    81  ncase = dimsizes(lat)

    82 ;**************************************************************

    83 ; get only the lowest level at each station

    84 ;**************************************************************

    85  lat_tmp = lat

    86  lat_tmp@_FillValue = 1.e+36

    88  do n = 0,ncase-1

    89     if (.not. ismissing(lat_tmp(n))) then

    90        indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)

    91        if (dimsizes(indexes) .gt. 1) then

    92           lat_tmp(indexes(1:)) = lat_tmp@_FillValue

    93        end if

    94        delete (indexes)

    95     end if

    96  end do

    98  indexes = ind(.not. ismissing(lat_tmp))

   100  lat_ob = lat(indexes)

   101  lon_ob = lon(indexes)

   102  val_ob = val(indexes,:)

   103 ;************************************************************

   104 ; interpolate model data into observed station

   105 ; note: model is 0-360E, 90S-90N

   106 ;************************************************************

   107 ; to be able to handle observation at (-89.98,-24.80)

   108   yi(0) = -90.

   110   i = ind(lon_ob .lt. 0.)

   111   lon_ob(i) = lon_ob(i) + 360.

   113   yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)

   115   val_model = yo(pts|:,time|:)

   116   val_model_0 = val_model

   118 ; remove annual mean

   119   val_model = val_model - conform(val_model,dim_avg(val_model),0)

   121   nzone = 4

   123 ;*******************************************************************

   124 ; for table -- zone

   125 ;*******************************************************************

   126 ; table header name

   127   table_header_name = "Zone"

   129 ; column (not including header column)

   130   col_header_zone = (/"Stations","Amplitude Ratio", \

   131                       "Correlation Coef","M score","Combined Score"/)

   132   ncol_zone       = dimsizes(col_header_zone )

   134 ; row (not including header row)

   135   row_header_zone = (/"60N-90N","30N-60N","EQ-30N","90S-EQ","Total"/)

   136   nrow_zone       = dimsizes(row_header_zone)

   138 ; arrays to be passed to table.

   139   value_zone = new ((/nrow_zone, ncol_zone/),string )

   140 ;--------------------------------------------------------------------

   142   table_length_zone = 0.4

   144 ; Table header

   145   ncr1  = (/1,1/)               ; 1 row, 1 column

   146   x1    = (/0.005,0.15/)        ; Start and end X

   147   y1    = (/0.900,0.995/)       ; Start and end Y

   148   text1 = table_header_name

   149   res1               = True

   150   res1@txFontHeightF = 0.03

   151   res1@gsFillColor   = "CornFlowerBlue"

   153 ; Column header (equally space in x2)

   154   ncr2  = (/1,ncol_zone/)         ; 1 rows, 5 columns

   155   x2    = (/x1(1),0.995/)          ; start from end of x1

   156   y2    = y1                      ; same as y1

   157   text2 = col_header_zone

   158   res2               = True

   159   res2@txFontHeightF = 0.015

   160   res2@gsFillColor   = "Gray"

   162 ; Row header (equally space in y2)

   163   ncr3  = (/nrow_zone,1/)         ; 5 rows, 1 columns

   164   x3    = x1                      ; same as x1

   165   y3    = (/1.0-table_length_zone,0.900/) ; end at start of y1

   166   text3 = row_header_zone

   167   res3               = True

   168   res3@txFontHeightF = 0.02

   169   res3@gsFillColor   = "Gray"

   171 ; Main table body

   172   ncr4  = (/nrow_zone,ncol_zone/) ; 5 rows, 5 columns

   173   x4    = x2                      ; Start and end x

   174   y4    = y3                      ; Start and end Y

   175   text4 = new((/nrow_zone,ncol_zone/),string)

   177   color_fill4      = new((/nrow_zone,ncol_zone/),string)

   178   color_fill4      = "white"

   179   color_fill4(:,ncol_zone-1) = "grey"

   181   res4               = True       ; Set up resource list

   182 ; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.

   183   res4@txFontHeightF = 0.02

   184   res4@gsFillColor   = color_fill4

   186   delete (color_fill4)

   188 ;*******************************************************************

   189 ; for table -- station

   190 ;*******************************************************************

   191 ; column (not including header column)

   192   col_header_sta = (/"Latitude","Longitude","Amplitude Ratio", \

   193                       "Correlation Coef","M score","Combined Score"/)

   194   ncol_sta       = dimsizes(col_header_sta )

   195 ;-------------------------------------------------------------------

   197 ; Table header

   198   ncr5  = (/1,1/)               ; 1 row, 1 column

   199   x5    = (/0.005,0.15/)        ; Start and end X

   200   y5    = (/0.900,0.995/)       ; Start and end Y

   201   res5               = True

   202   res5@txFontHeightF = 0.02

   203   res5@gsFillColor   = "CornFlowerBlue"

   205 ; Column header (equally space in x2)

   206   ncr6  = (/1,ncol_sta/)         ; 1 rows, 5 columns

   207   x6    = (/x5(1),0.995/)          ; start from end of x1

   208   y6    = y5                      ; same as y1

   209   text6 = col_header_sta

   210   res6               = True

   211   res6@txFontHeightF = 0.012

   212   res6@gsFillColor   = "Gray"

   214 ; Row header (equally space in y2)

   216   res7               = True

   217   res7@txFontHeightF = 0.015

   218   res7@gsFillColor   = "Gray"

   219 ;--------------------------------------------------------------

   220 ; for station line plot

   222 ; for x-axis in xyplot

   223   mon = ispan(1,12,1)

   224   mon@long_name = "month"

   226   plot_type = "ps"

   227   plot_type_new = "png"

   229   res                   = True                      ; plot mods desired

   230   res@xyLineThicknesses = (/2.0,2.0,2.0/)           ; make 2nd lines thicker

   231   res@xyLineColors      = (/"red","black"/)  ; change line color

   233 ; Add a boxed legend using the more simple method

   235   res@pmLegendDisplayMode    = "Always"

   236 ; res@pmLegendWidthF         = 0.1

   237   res@pmLegendWidthF         = 0.08

   238   res@pmLegendHeightF        = 0.06

   239 ; res@pmLegendOrthogonalPosF = -1.17

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

   241   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   243 ; res@pmLegendParallelPosF   =  0.18

   244   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   246 ; res@lgPerimOn             = False

   247   res@lgLabelFontHeightF     = 0.015

   248   res@xyExplicitLegendLabels = (/"b30.061n","observed"/)

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

   251   do z = 0,nzone-1

   253   if (z .eq. 0) then

   254 ;    maximum score for the zone, 60N-90N

   255      zone = "60N-90N"

   256      score_max = 5.0

   257 ;    index of stations in this zone

   258      ind_z = ind(lat_ob .ge. 60.)

   259   end if

   261   if (z .eq. 1) then

   262 ;    maximum score for the zone, 30N-60N

   263      zone = "30N-60N"

   264      score_max = 5.0

   265 ;    index of stations in this zone

   266      ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)

   267   end if

   269   if (z .eq. 2) then

   270 ;    maximum score for the zone, EQ-30N

   271      zone = "EQ-30N"

   272      score_max = 5.0

   273 ;    index of stations in this zone

   274      ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)

   275   end if

   277   if (z .eq. 3) then

   278 ;    maximum score for the zone, 90S-EQ

   279      zone = "90S-EQ"

   280      score_max = 5.0

   281 ;    index of stations in this zone

   282      ind_z = ind(lat_ob .lt. 0. )

   283   end if

   285  npts = dimsizes(ind_z)

   287 ;-------------------------------------------------------------------------

   288 ; for metric table computation

   290  amp_ob        = new((/npts/),float)

   291  amp_model     = new((/npts/),float)

   293  amp_ratio_sta = new((/npts/),float)

   294  ccr_sta       = new((/npts/),float)

   295  M_sta         = new((/npts/),float)

   296  score_sta     = new((/npts/),float)

   298 ;----------------------

   299 ; for table -- station

   300 ;----------------------

   302 ; row (not including header row)

   303   nrow_sta       = npts

   305 ; Table header

   306   ncr5  = (/1,1/)               ; 1 row, 1 column

   307   x5    = (/0.005,0.15/)        ; Start and end X

   308   y5    = (/0.900,0.995/)       ; Start and end Y

   309   res5               = True

   310   res5@txFontHeightF = 0.02

   311   res5@gsFillColor   = "CornFlowerBlue"

   313 ; Column header (equally space in x2)

   314   ncr6  = (/1,ncol_sta/)         ; 1 rows, 5 columns

   315   x6    = (/x1(1),0.995/)          ; start from end of x1

   316   y6    = y1                      ; same as y1

   317   text6 = col_header_sta

   318   res6               = True

   319   res6@txFontHeightF = 0.012

   320   res6@gsFillColor   = "Gray"

   322 ; Row header (equally space in y2)

   323   ncr7  = (/nrow_sta,1/)         ; 5 rows, 1 columns

   324   x7    = x1                      ; same as x1

   325 ; y7    = (/1.0-table_length_sta,0.900/) ; end at start of y1

   326   text7 = new((/nrow_sta/),string)

   327   res7               = True

   328   res7@txFontHeightF = 0.015

   329   res7@gsFillColor   = "Gray"

   330 ;-------------------------------------------------------------------------

   331 ; for station line plot

   333   npts_str = ""

   334   npts_str = npts

   336   plot_data   = new((/2,12,npts/),float)

   337   plot_data_0 = new((/12,npts/),float)

   339   plot_data!0 = "case"

   340   plot_data!1 = "month"

   341   plot_data!2 = "pts"

   342   plot_data@long_name   = "CO2 Seasonal"

   344   plot_data_0!0 = "month"

   345   plot_data_0!1 = "pts"

   346   plot_data_0@long_name = "CO2"

   347 ;--------------------------------------------------------------------------

   348  do n=0,npts-1

   349     amp_ob(n)    = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:))

   350     amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))

   352     amp_ratio_sta(n) = amp_model(n)/amp_ob(n)

   353     ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)

   354     M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)

   355     score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max

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

   358 ; for station line plot

   360     plot_data(0,:,n) = (/val_model(ind_z(n),:)/)

   361     plot_data(1,:,n) = (/val_ob(ind_z(n),:)/)

   363     plot_data_0(:,n) = (/val_model_0(ind_z(n),:)/)

   365     plot_name = sta(ind_z(n))

   366     title = plot_name+"("+lat(ind_z(n))+","+lon(ind_z(n))+")"

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

   369 ;------------------------------------------

   370 ;   for panel plot

   372     plot=new(2,graphic)                        ; create graphic array

   373     res@gsnFrame     = False                   ; Do not draw plot

   374     res@gsnDraw      = False                   ; Do not advance frame

   376     pres                            = True     ; panel plot mods desired

   377     pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   378                                                ; indiv. plots in panel

   379     pres@gsnMaximize                = True     ; fill the page

   380 ;------------------------------------------

   381     res@tiMainString = title                           ; add title

   383     plot(0)=gsn_csm_xy(wks,mon,plot_data(:,:,n),res)   ; create plot 1

   385     plot(1)=gsn_csm_xy(wks,mon,plot_data_0(:,n),res) ; create plot 2

   387     gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot

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

   390     system("rm "+plot_name+"."+plot_type)

   392     clear (wks)

   393 ;---------------------------------------------------------------------------

   394  end do

   395 ;-------------------------------------------------------------------------

   396 ;   for line plot in a zone

   398     plot_name = "All_"+npts_str

   399     title = plot_name + " in "+ zone

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

   402 ;-----------------------------------------

   403 ;   for panel plot

   405     plot=new(2,graphic)                        ; create graphic array

   406     res@gsnFrame     = False                   ; Do not draw plot

   407     res@gsnDraw      = False                   ; Do not advance frame

   409     pres                            = True     ; panel plot mods desired

   410     pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   411                                                ; indiv. plots in panel

   412     pres@gsnMaximize                = True     ; fill the page

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

   414     res@tiMainString = title                                     ; add title

   416     plot(0) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data),res)   ; create plot 1

   418     plot(1) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data_0),res) ; create plot 2

   420     gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot

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

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

   425     clear (wks)

   426 ;   delete (ind_z)

   427     delete (plot_data)

   428     delete (plot_data_0)

   429 ;---------------------------------------------------------------------------

   430 ; for zone table value

   432  amp_ratio_zone = avg(amp_ratio_sta)

   433  ccr_zone       = avg(ccr_sta)

   434  M_zone   = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts)

   435  score_zone   = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max

   437   text4(z,0) = sprintf("%5.2f", npts)

   438   text4(z,1) = sprintf("%5.2f", amp_ratio_zone)

   439   text4(z,2) = sprintf("%5.2f", ccr_zone)

   440   text4(z,3) = sprintf("%5.2f", M_zone)

   441   text4(z,4) = sprintf("%5.2f", score_zone)

   442 ;---------------------------------------------------------------------------

   443 ; plot station table

   444 ;----------------------------------

   445 ; header value for station table

   446   text5  = zone

   448 ; row value for station table

   449   text7  = sta(ind_z)

   451   if (z .eq. 0) then

   452      table_length_sta = 0.5

   453   end if

   454   if (z .eq. 1) then

   455      table_length_sta = 0.995

   456   end if

   457   if (z .eq. 2) then

   458      table_length_sta = 0.8

   459   end if

   460   if (z .eq. 3) then

   461      table_length_sta = 0.8

   462   end if

   464   x7 = x5

   465   y7 = (/1.0-table_length_sta,0.900/) ; end at start of y1

   467 ; Main table body

   468   ncr8  = (/nrow_sta,ncol_sta/) ; 5 rows, 5 columns

   469   x8    = x6                      ; Start and end x

   470   y8    = y7                      ; Start and end Y

   471   text8 = new((/nrow_sta,ncol_sta/),string)

   473   color_fill8      = text8

   474   color_fill8      = "white"

   475   color_fill8(:,ncol_sta-1) = "grey"

   477   res8               = True       ; Set up resource list

   478 ; res8@gsnDebug      = True       ; Useful to print NDC row,col values used.

   479   res8@txFontHeightF = 0.02

   480   res8@gsFillColor   = color_fill8

   482   delete (color_fill8)

   484 ; table value for station table

   485   text8(:,0) = sprintf("%5.2f", (/lat(ind_z)/))

   486   text8(:,1) = sprintf("%5.2f", (/lon(ind_z)/))

   487   text8(:,2) = sprintf("%5.2f", (/amp_ratio_sta/))

   488   text8(:,3) = sprintf("%5.2f", (/ccr_sta/))

   489   text8(:,4) = sprintf("%5.2f", (/M_sta/))

   490   text8(:,5) = sprintf("%5.2f", (/score_sta/))

   492   plot_name = "table_sta." + zone

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

   496   gsn_table(wks,ncr5,x5,y5,text5,res5)

   497   gsn_table(wks,ncr6,x6,y6,text6,res6)

   498   gsn_table(wks,ncr7,x7,y7,text7,res7)

   499   gsn_table(wks,ncr8,x8,y8,text8,res8)

   501   frame(wks)

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

   504   system("rm "+plot_name+"."+plot_type)

   506  delete (ind_z)

   507  delete (amp_model)

   508  delete (amp_ob)

   509  delete (amp_ratio_sta)

   510  delete (ccr_sta)

   511  delete (M_sta)

   512  delete (score_sta)

   513  delete (text7)

   514  delete (text8)

   515  delete (res7)

   516  delete (res8)

   517  clear (wks)

   518  end do

   519 ;**************************************************

   520 ; plot zone table

   521 ;**************************************************

   523   text4(4,0) = sum(stringtofloat(text4(0:3,0)))

   524   text4(4,1) = 0.

   525   text4(4,2) = 0.

   526   text4(4,3) = 0.

   527   text4(4,4) = sum(stringtofloat(text4(0:3,4)))

   529   plot_name = "table_zone"

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

   532   gsn_table(wks,ncr1,x1,y1,text1,res1)

   533   gsn_table(wks,ncr2,x2,y2,text2,res2)

   534   gsn_table(wks,ncr3,x3,y3,text3,res3)

   535   gsn_table(wks,ncr4,x4,y4,text4,res4)

   537   frame(wks)

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

   540   system("rm "+plot_name+"."+plot_type)

   542   M_co2 = text4(4,4)

   544   if (isvar("compare")) then

   545      system("sed 1,/M_co2/s//"+M_co2+"/ "+html_name2+" > "+html_new2+";"+ \

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

   547   end if

   549   system("sed s#M_co2#"+M_co2+"# "+html_name+" > "+html_new+";"+ \

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

   551 ;***************************************************************************

   552 ; output plots

   553 ;***************************************************************************

   554   output_dir = model_name+"/co2"

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

   557 ;***************************************************************************

   558 end

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