clamp: co2/99.all.ncl@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"

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

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

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

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

    10 begin

    11 ; add line to ascci/html file

    13   nnewlines = dimsizes(newlines)

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

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

    16     return

    17   end if

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

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

    20   nline = nline + nnewlines

    21   return

    22 end

    23 ;****************************************************************************

    25 begin

    27   plot_type = "ps"

    28   plot_type_new = "png"

    30 ;************************************************

    31 ; read model data

    32 ;************************************************

    34 ; from command line inputs

    36 ;--------------------------------------------------

    37 ; edit table.html of current model for movel1_vs_model2

    39   if (isvar("compare")) then

    40      html_name2 = compare+"/table.html"

    41      html_new2  = html_name2 +".new"

    43 ; the following only needs to execute once

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

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

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

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

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

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

    50   end if

    52 ;-------------------------------------

    53 ; edit table.html for current model

    55   html_name = model_name+"/table.html"

    56   html_new  = html_name +".new"

    58 ; the following only needs to execute once

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

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

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

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

    64   x     = fm->CO2

    65   xi    = fm->lon

    66   yi    = fm->lat

    68   xdim  = dimsizes(x)

    69   nlev  = xdim(1)

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

    72 ; get co2 at the lowest level

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

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

    76 ; Model_units [=] kgCO2 / kgDryAir

    77 ; 28.966 = molecular weight of dry air

    78 ; 44.       = molecular weight of CO2

    79 ; u mol = 1e-6 mol

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

    82   y      = y * factor

    84   y@_FillValue = 1.e36

    85   y@units      = "u mol/mol"

    86 ;************************************************

    87 ; read data: observed

    88 ;************************************************

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

    90   fili  = "co2_globalView_98.nc"

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

    92   val   = g->CO2_SEAS

    93   lon   = g->LON

    94   lat   = g->LAT

    95   sta   = chartostring(g->STATION)

    96   delete (g)

    98   ncase = dimsizes(lat)

    99 ;**************************************************************

   100 ; get only the lowest level at each station

   101 ;**************************************************************

   102   lat_tmp = lat

   103   lat_tmp@_FillValue = 1.e+36

   105   do n = 0,ncase-1

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

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

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

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

   110         end if

   111         delete (indexes)

   112      end if

   113   end do

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

   117   lat_ob = lat(indexes)

   118   lon_ob = lon(indexes)

   119   val_ob = val(indexes,:)

   120 ;************************************************************

   121 ; interpolate model data into observed station

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

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

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

   125   yi(0) = -90.

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

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

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

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

   133   val_model_0 = val_model

   134 ;************************************************************

   135 ; remove annual mean

   136 ;************************************************************

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

   139 ;*******************************************************************

   140 ; res for station line plot

   141 ;*******************************************************************

   142 ; for x-axis in xyplot

   143   mon = ispan(1,12,1)

   144   mon@long_name = "month"

   146   res                   = True                      ; plot mods desired

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

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

   150 ; Add a boxed legend using the more simple method

   151   res@pmLegendDisplayMode    = "Always"

   152 ; res@pmLegendWidthF         = 0.1

   153   res@pmLegendWidthF         = 0.08

   154   res@pmLegendHeightF        = 0.06

   155 ; res@pmLegendOrthogonalPosF = -1.17

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

   157   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   159 ; res@pmLegendParallelPosF   =  0.18

   160   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   162 ; res@lgPerimOn             = False

   163   res@lgLabelFontHeightF     = 0.015

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

   165 ;************************************************************

   166 ; number of latitude zone

   167 ;************************************************************

   168   nzone = 4

   170 ; saving data for zone

   171 ; number of rows for zone table (with data)

   172   nrow_zone = nzone

   174 ; number of columns for zone table

   175   ncol_zone = 6

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

   179 do z = 0,nzone-1

   181   if (z .eq. 0) then

   182      zone = "60N-90N"

   183      score_max = 5.0

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

   185   end if

   187   if (z .eq. 1) then

   188      zone = "30N-60N"

   189      score_max = 5.0

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

   191   end if

   193   if (z .eq. 2) then

   194      zone = "EQ-30N"

   195      score_max = 5.0

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

   197   end if

   199   if (z .eq. 3) then

   200      zone = "90S-EQ"

   201      score_max = 5.0

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

   203   end if

   205   npts = dimsizes(ind_z)

   207 ;------------------------------------------------------

   208 ; for metric table computation

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

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

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

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

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

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

   216 ;-----------------------------------------------------

   217 ; for station line plot

   218   npts_str = ""

   219   npts_str = npts

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

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

   224   plot_data!0 = "case"

   225   plot_data!1 = "month"

   226   plot_data!2 = "pts"

   227   plot_data@long_name   = "CO2 Seasonal"

   229   plot_data_0!0 = "month"

   230   plot_data_0!1 = "pts"

   231   plot_data_0@long_name = "CO2"

   232 ;--------------------------------------------------------------------------

   233   do n=0,npts-1

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

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

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

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

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

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

   243 ;----------------------------------------------------------------------

   244 ; for station line plot

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

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

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

   251      plot_name = sta(ind_z(n))

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

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

   255 ;------------------------------------------

   256 ;    for panel plot

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

   259      res@gsnFrame     = False                   ; Do not draw plot

   260      res@gsnDraw      = False                   ; Do not advance frame

   262      pres                            = True     ; panel plot mods desired

   263      pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   264                                                ; indiv. plots in panel

   265      pres@gsnMaximize                = True     ; fill the page

   266 ;------------------------------------------

   267      res@tiMainString = title                           ; add title

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

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

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

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

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

   278      clear (wks)

   279 ;---------------------------------------------------------------------------

   280   end do

   281 ;-------------------------------------------------------------------------

   282 ; for line plot in a zone

   284   plot_name = "All_"+npts_str

   285   title = plot_name + " in "+ zone

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

   288 ;-----------------------------------------

   289 ; for panel plot

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

   292   res@gsnFrame     = False                   ; Do not draw plot

   293   res@gsnDraw      = False                   ; Do not advance frame

   295   pres                            = True     ; panel plot mods desired

   296   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   297                                                ; indiv. plots in panel

   298   pres@gsnMaximize                = True     ; fill the page

   299 ;-----------------------------------------

   300   res@tiMainString = title                                     ; add title

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

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

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

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

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

   311   clear (wks)

   312 ; delete (ind_z)

   313   delete (plot_data)

   314   delete (plot_data_0)

   315 ;---------------------------------------------------------------------------

   316 ; values saved for zone table

   318   amp_ratio_zone = avg(amp_ratio_sta)

   319   ccr_zone       = avg(ccr_sta)

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

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

   323   text4(z,0) = zone

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

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

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

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

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

   330 ;*******************************************************************

   331 ; html table -- station

   332 ;*******************************************************************

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

   335   header = (/"<HTML>" \

   336             ,"<HEAD>" \

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

   338             ,"</HEAD>" \

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

   340             /)

   341   footer = "</HTML>"

   343   table_header = (/ \

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

   345        ,"<tr>" \

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

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

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

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

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

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

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

   353        ,"</tr>" \

   354        /)

   355   table_footer = "</table>"

   356   row_header = "<tr>"

   357   row_footer = "</tr>"

   359   lines = new(50000,string)

   360   nline = 0

   362   set_line(lines,nline,header)

   363   set_line(lines,nline,table_header)

   364 ;-----------------------------------------------

   365 ; row of table

   367   do n = 0,npts-1

   368      set_line(lines,nline,row_header)

   370      txt0 = sta(ind_z(n))

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

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

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

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

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

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

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

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

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

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

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

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

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

   386      set_line(lines,nline,row_footer)

   387   end do

   389 ; last row, summary

   390   set_line(lines,nline,row_header)

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

   393   txt1 = "-"

   394   txt2 = "-"

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

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

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

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

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

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

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

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

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

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

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

   408   set_line(lines,nline,row_footer)

   409 ;-----------------------------------------------

   410   set_line(lines,nline,table_footer)

   411   set_line(lines,nline,footer)

   413 ; Now write to an HTML file.

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

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

   416     asciiwrite(output_html,lines(idx))

   417   else

   418    print ("error?")

   419   end if

   420   delete (idx)

   421 ;-----------------------------------------------------------------

   423   delete (ind_z)

   424   delete (amp_model)

   425   delete (amp_ob)

   426   delete (amp_ratio_sta)

   427   delete (ccr_sta)

   428   delete (M_sta)

   429   delete (score_sta)

   430   clear (wks)

   431 end do

   433 ;*******************************************************************

   434 ; html table -- zone

   435 ;*******************************************************************

   436   output_html = "score+line_vs_ob.html"

   438   header = (/"<HTML>" \

   439             ,"<HEAD>" \

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

   441             ,"</HEAD>" \

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

   443             /)

   444   footer = "</HTML>"

   446   delete (table_header)

   447   table_header = (/ \

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

   449        ,"<tr>" \

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

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

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

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

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

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

   456        ,"</tr>" \

   457        /)

   458   table_footer = "</table>"

   459   row_header = "<tr>"

   460   row_footer = "</tr>"

   462   lines = new(50000,string)

   463   nline = 0

   465   set_line(lines,nline,header)

   466   set_line(lines,nline,table_header)

   467 ;-----------------------------------------------

   468 ;row of table

   470   do n = 0,nrow_zone-1

   471      set_line(lines,nline,row_header)

   473      set_line(lines,nline,"<th><a href=score+line_"+text4(n,0)+".html>"+text4(n,0)+"</th>")

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

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

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

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

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

   480      set_line(lines,nline,row_footer)

   481   end do

   483 ; for the last row

   485      txt0 = "All"

   486      txt1 = sum(stringtofloat(text4(0:3,1)))

   487      txt2 = "-"

   488      txt3 = "-"

   489      txt4 = "-"

   490      txt5 = sum(stringtofloat(text4(0:3,5)))

   492      set_line(lines,nline,row_header)

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

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

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

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

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

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

   501      set_line(lines,nline,row_footer)

   502 ;-----------------------------------------------

   503   set_line(lines,nline,table_footer)

   504   set_line(lines,nline,footer)

   506 ; Now write to an HTML file.

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

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

   509     asciiwrite(output_html,lines(idx))

   510   else

   511    print ("error?")

   512   end if

   513   delete (idx)

   514 ;--------------------------------------------------------------------------

   516   M_co2 = txt5

   518   if (isvar("compare")) then

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

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

   521   end if

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

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

   525 ;***************************************************************************

   526 ; output plots

   527 ;***************************************************************************

   528   output_dir = model_name+"/co2"

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

   531 ;***************************************************************************

   532 end

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