clamp: co2/99.all.ncl.2@0c6405ab2ff4

     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 procedure set_line(lines:string,nline:integer,newlines:string)

    11 begin

    12 ; add line to ascci/html file

    14   nnewlines = dimsizes(newlines)

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

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

    17     return

    18   end if

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

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

    21   nline = nline + nnewlines

    22   return

    23 end

    24 ;****************************************************************************

    26 begin

    28   plot_type = "ps"

    29   plot_type_new = "png"

    31 ;************************************************

    32 ; read model data

    33 ;************************************************

    35 ; from command line inputs

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

    38 ; edit table.html of current model for movel1_vs_model2

    40   if (isvar("compare")) then

    41      html_name2 = compare+"/table.html"

    42      html_new2  = html_name2 +".new"

    44 ; the following only needs to execute once

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

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

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

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

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

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

    51   end if

    53 ;-------------------------------------

    54 ; edit table.html for current model

    56   html_name = model_name+"/table.html"

    57   html_new  = html_name +".new"

    59 ; the following only needs to execute once

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

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

    62 ;------------------------------------------------

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

    65   x     = fm->CO2

    66   xi    = fm->lon

    67   yi    = fm->lat

    69   xdim  = dimsizes(x)

    70   nlev  = xdim(1)

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

    73 ; get co2 at the lowest level

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

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

    77 ; Model_units [=] kgCO2 / kgDryAir

    78 ; 28.966 = molecular weight of dry air

    79 ; 44.       = molecular weight of CO2

    80 ; u mol = 1e-6 mol

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

    83   y      = y * factor

    85   y@_FillValue = 1.e36

    86   y@units      = "u mol/mol"

    87 ;************************************************

    88 ; read data: observed

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

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

    91   fili  = "co2_globalView_98.nc"

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

    93   val   = g->CO2_SEAS

    94   lon   = g->LON

    95   lat   = g->LAT

    96   sta   = chartostring(g->STATION)

    97   delete (g)

    99   ncase = dimsizes(lat)

   100 ;**************************************************************

   101 ; get only the lowest level at each station

   102 ;**************************************************************

   103   lat_tmp = lat

   104   lat_tmp@_FillValue = 1.e+36

   106   do n = 0,ncase-1

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

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

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

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

   111         end if

   112         delete (indexes)

   113      end if

   114   end do

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

   118   lat_ob = lat(indexes)

   119   lon_ob = lon(indexes)

   120   val_ob = val(indexes,:)

   121 ;************************************************************

   122 ; interpolate model data into observed station

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

   124 ;************************************************************

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

   126   yi(0) = -90.

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

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

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

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

   134   val_model_0 = val_model

   135 ;************************************************************

   136 ; remove annual mean

   137 ;************************************************************

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

   140 ;*******************************************************************

   141 ; res for station line plot

   142 ;*******************************************************************

   143 ; for x-axis in xyplot

   144   mon = ispan(1,12,1)

   145   mon@long_name = "month"

   147   res                   = True                      ; plot mods desired

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

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

   151 ; Add a boxed legend using the more simple method

   152   res@pmLegendDisplayMode    = "Always"

   153 ; res@pmLegendWidthF         = 0.1

   154   res@pmLegendWidthF         = 0.08

   155   res@pmLegendHeightF        = 0.06

   156 ; res@pmLegendOrthogonalPosF = -1.17

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

   158   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   160 ; res@pmLegendParallelPosF   =  0.18

   161   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   163 ; res@lgPerimOn             = False

   164   res@lgLabelFontHeightF     = 0.015

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

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

   167 ; number of latitude zone

   168 ;************************************************************

   169   nzone = 4

   171 ; number of rows for zone table

   172   nrow_zone = nzone + 1

   174 ; number of columns for zone table

   175   ncol_zone = 6

   177 ; saving data for zone

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

   180 do z = 0,nzone-1

   182   if (z .eq. 0) then

   183      zone = "60N-90N"

   184      score_max = 5.0

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

   186   end if

   188   if (z .eq. 1) then

   189      zone = "30N-60N"

   190      score_max = 5.0

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

   192   end if

   194   if (z .eq. 2) then

   195      zone = "EQ-30N"

   196      score_max = 5.0

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

   198   end if

   200   if (z .eq. 3) then

   201      zone = "90S-EQ"

   202      score_max = 5.0

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

   204   end if

   206   npts = dimsizes(ind_z)

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

   209 ; for metric table computation

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

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

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

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

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

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

   217 ;-----------------------------------------------------

   218 ; for station line plot

   219   npts_str = ""

   220   npts_str = npts

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

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

   225   plot_data!0 = "case"

   226   plot_data!1 = "month"

   227   plot_data!2 = "pts"

   228   plot_data@long_name   = "CO2 Seasonal"

   230   plot_data_0!0 = "month"

   231   plot_data_0!1 = "pts"

   232   plot_data_0@long_name = "CO2"

   233 ;--------------------------------------------------------------------------

   234   do n=0,npts-1

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

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

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

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

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

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

   244 ;----------------------------------------------------------------------

   245 ; for station line plot

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

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

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

   252      plot_name = sta(ind_z(n))

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

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

   256 ;------------------------------------------

   257 ;    for panel plot

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

   260      res@gsnFrame     = False                   ; Do not draw plot

   261      res@gsnDraw      = False                   ; Do not advance frame

   263      pres                            = True     ; panel plot mods desired

   264      pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   265                                                ; indiv. plots in panel

   266      pres@gsnMaximize                = True     ; fill the page

   267 ;------------------------------------------

   268      res@tiMainString = title                           ; add title

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

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

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

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

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

   279      clear (wks)

   280 ;---------------------------------------------------------------------------

   281   end do

   282 ;-------------------------------------------------------------------------

   283 ; for line plot in a zone

   285   plot_name = "All_"+npts_str

   286   title = plot_name + " in "+ zone

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

   289 ;-----------------------------------------

   290 ; for panel plot

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

   293   res@gsnFrame     = False                   ; Do not draw plot

   294   res@gsnDraw      = False                   ; Do not advance frame

   296   pres                            = True     ; panel plot mods desired

   297   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   298                                                ; indiv. plots in panel

   299   pres@gsnMaximize                = True     ; fill the page

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

   301   res@tiMainString = title                                     ; add title

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

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

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

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

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

   312   clear (wks)

   313 ; delete (ind_z)

   314   delete (plot_data)

   315   delete (plot_data_0)

   316 ;---------------------------------------------------------------------------

   317 ; values saved for zone table

   319   amp_ratio_zone = avg(amp_ratio_sta)

   320   ccr_zone       = avg(ccr_sta)

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

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

   324   text4(z,0) = zone

   325   text4(z,1) = sprintf("%.0f", npts)

   326   text4(z,2) = sprintf("%.2f", amp_ratio_zone)

   327   text4(z,3) = sprintf("%.2f", ccr_zone)

   328   text4(z,4) = sprintf("%.2f", M_zone)

   329   text4(z,5) = sprintf("%.2f", score_zone)

   331 ;*******************************************************************

   332 ; html table -- station

   333 ;*******************************************************************

   334   output_html = "score+line_"+zone+".html"

   336   header = (/"<HTML>" \

   337             ,"<HEAD>" \

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

   339             ,"</HEAD>" \

   340             ,"<H1>CO2 in Zone "+zone+": Model "+model_name+"</H1>" \

   341             /)

   342   footer = "</HTML>"

   344   table_header = (/ \

   345         "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \

   346        ,"<tr>" \

   347        ,"   <th bgcolor=DDDDDD >Site Name</th>" \

   348        ,"   <th bgcolor=DDDDDD >Latitude</th>" \

   349        ,"   <th bgcolor=DDDDDD >Longitude</th>" \

   350        ,"   <th bgcolor=DDDDDD >Amplitude Ratio</th>" \

   351        ,"   <th bgcolor=DDDDDD >Correleration Coef</th>" \

   352        ,"   <th bgcolor=DDDDDD >M Score</th>" \

   353        ,"   <th bgcolor=DDDDDD >Combined Score</th>" \

   354        ,"</tr>" \

   355        /)

   356   table_footer = "</table>"

   357   row_header = "<tr>"

   358   row_footer = "</tr>"

   360   lines = new(50000,string)

   361   nline = 0

   363   set_line(lines,nline,header)

   364   set_line(lines,nline,table_header)

   365 ;-----------------------------------------------

   366 ; row of table

   368   do n = 0,npts-1

   369      set_line(lines,nline,row_header)

   371      txt0 = sta(ind_z(n))

   372      txt1 = sprintf("%5.2f", (/lat(ind_z(n))/))

   373      txt2 = sprintf("%5.2f", (/lon(ind_z(n))/))

   374      txt3 = sprintf("%5.2f", (/amp_ratio_sta(n)/))

   375      txt4 = sprintf("%5.2f", (/ccr_sta(n)/))

   376      txt5 = sprintf("%5.2f", (/M_sta(n)/))

   377      txt6 = sprintf("%5.2f", (/score_sta(n)/))

   379      set_line(lines,nline,"<th><a href="+txt0+".png>"+txt0+"</a></th>")

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

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

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

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

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

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

   387      set_line(lines,nline,row_footer)

   388   end do

   390 ; last row, summary

   391   set_line(lines,nline,row_header)

   393   txt0 = "All_"+sprintf("%.0f", (/npts/))

   394   txt1 = "-"

   395   txt2 = "-"

   396   txt3 = sprintf("%5.2f", (/amp_ratio_zone/))

   397   txt4 = sprintf("%5.2f", (/ccr_zone/))

   398   txt5 = sprintf("%5.2f", (/M_zone/))

   399   txt6 = sprintf("%5.2f", (/score_zone/))

   401   set_line(lines,nline,"<th><a href="+txt0+".png>"+txt0+"</a></th>")

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

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

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

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

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

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

   409   set_line(lines,nline,row_footer)

   410 ;-----------------------------------------------

   411   set_line(lines,nline,table_footer)

   412   set_line(lines,nline,footer)

   414 ; Now write to an HTML file.

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

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

   417     asciiwrite(output_html,lines(idx))

   418   else

   419    print ("error?")

   420   end if

   421   delete (idx)

   422 ;-----------------------------------------------------------------

   424   delete (ind_z)

   425   delete (amp_model)

   426   delete (amp_ob)

   427   delete (amp_ratio_sta)

   428   delete (ccr_sta)

   429   delete (M_sta)

   430   delete (score_sta)

   431   clear (wks)

   432 end do

   434 ;*******************************************************************

   435 ; html table -- zone

   436 ;*******************************************************************

   437   output_html = "score_zone.html"

   439   header = (/"<HTML>" \

   440             ,"<HEAD>" \

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

   442             ,"</HEAD>" \

   443             ,"<H1>CO2 in Latitude Zone: Model "+model_name+"</H1>" \

   444             /)

   445   footer = "</HTML>"

   447   delete (table_header)

   448   table_header = (/ \

   449         "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \

   450        ,"<tr>" \

   451        ,"   <th bgcolor=DDDDDD >Zone</th>" \

   452        ,"   <th bgcolor=DDDDDD >Number of Site</th>" \

   453        ,"   <th bgcolor=DDDDDD >Amplitide Ratio</th>" \

   454        ,"   <th bgcolor=DDDDDD >Correleration Coef</th>" \

   455        ,"   <th bgcolor=DDDDDD >M Score</th>" \

   456        ,"   <th bgcolor=DDDDDD >Combined Score</th>" \

   457        ,"</tr>" \

   458        /)

   459   table_footer = "</table>"

   460   row_header = "<tr>"

   461   row_footer = "</tr>"

   463   lines = new(50000,string)

   464   nline = 0

   466   set_line(lines,nline,header)

   467   set_line(lines,nline,table_header)

   468 ;-----------------------------------------------

   469 ;row of table

   471 ; sum for the last row

   472   text4(4,0) = "All"

   473   text4(4,1) = sum(stringtofloat(text4(0:3,1)))

   474   text4(4,2) = "-"

   475   text4(4,3) = "-"

   476   text4(4,4) = "-"

   477   text4(4,5) = sum(stringtofloat(text4(0:3,5)))

   479   do n = 0,nrow_zone-1

   480      set_line(lines,nline,row_header)

   482      set_line(lines,nline,"<th>"+text4(n,0)+"</th>")

   483      set_line(lines,nline,"<th>"+text4(n,1)+"</th>")

   484      set_line(lines,nline,"<th>"+text4(n,2)+"</th>")

   485      set_line(lines,nline,"<th>"+text4(n,3)+"</th>")

   486      set_line(lines,nline,"<th>"+text4(n,4)+"</th>")

   487      set_line(lines,nline,"<th>"+text4(n,5)+"</th>")

   489      set_line(lines,nline,row_footer)

   490   end do

   491 ;-----------------------------------------------

   492   set_line(lines,nline,table_footer)

   493   set_line(lines,nline,footer)

   495 ; Now write to an HTML file.

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

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

   498     asciiwrite(output_html,lines(idx))

   499   else

   500    print ("error?")

   501   end if

   502   delete (idx)

   503 ;--------------------------------------------------------------------------

   505   M_co2 = text4(4,5)

   507   if (isvar("compare")) then

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

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

   510   end if

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

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

   514 ;***************************************************************************

   515 ; output plots

   516 ;***************************************************************************

   517   output_dir = model_name+"/co2"

   519   system("mv *.png *.html " + output_dir)

   520 ;***************************************************************************

   521 end

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