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

    11 ;************************************************

    12 ; read in model data

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

    14   diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

    15 ; fili2 = "b30.061m_401_425_MONS_climo_atm.nc"

    16   fili2 = "b30.061n_1995-2004_MONS_climo_atm.nc"

    18   g     = addfile(diri2+fili2,"r")

    19   x     = g->CO2

    20   xi    = g->lon

    21   yi    = g->lat

    22   xdim  = dimsizes(x)

    23   nlev  = xdim(1)

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

    25 ; printVarSummary (y)

    27 ; get the co2 at the lowest level

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

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

    31 ; Model_units [=] kgCO2 / kgDryAir

    32 ; 28.966 = molecular weight of dry air

    33 ; 44.       = molecular weight of CO2

    34 ; u mol = 1e-6 mol

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

    37   y      = y * factor

    39   y@_FillValue = 1.e36

    40   y@units      = "u mol/mol"

    41 ; y = where(y0 .lt. 287.,y@_FillValue,y)

    42 ; printVarSummary (y)

    43 ; print (min(y)+"/"+max(y))

    45 ;************************************************

    46 ; read in data: observed

    47 ;************************************************

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

    49  fili  = "co2_globalView_98.nc"

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

    51  val   = g->CO2_SEAS

    52  lon   = g->LON

    53  lat   = g->LAT

    54  sta   = chartostring(g->STATION)

    55  delete (g)

    57 ;print (sta(0))

    59  ncase = dimsizes(lat)

    60 ;print (ncase)

    62 ;**************************************************************

    63 ; get only the lowest level at each station

    64 ;**************************************************************

    65  lat_tmp = lat

    66  lat_tmp@_FillValue = 1.e+36

    68  do n = 0,ncase-1

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

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

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

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

    73        end if

    74        delete (indexes)

    75     end if

    76  end do

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

    79 ;print (dimsizes(indexes))

    80 ;print (indexes)

    82  lat_ob = lat(indexes)

    83  lon_ob = lon(indexes)

    84  val_ob = val(indexes,:)

    85 ;printVarSummary (val_ob)

    86 ;print (lat_ob +"/"+lon_ob)

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

    89 ; interpolate model data into observed station

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

    91 ;************************************************************

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

    93   print (yi(0))

    94   yi(0) = -90.

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

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

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

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

   102   val_model_0 = val_model

   103 ; printVarSummary (val_model)

   104 ; print (min(val_model)+"/"+max(val_model))

   106 ; remove annual mean

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

   108 ; print (min(val_model)+"/"+max(val_model))

   110   nzone = 4

   112 ;*******************************************************************

   113 ; for table -- zone

   114 ;*******************************************************************

   115 ; table header name

   116   table_header_name = "Zone"

   118 ; column (not including header column)

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

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

   121   ncol_zone       = dimsizes(col_header_zone )

   123 ; row (not including header row)

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

   125   nrow_zone       = dimsizes(row_header_zone)

   127 ; arrays to be passed to table.

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

   129 ;--------------------------------------------------------------------

   131   table_length_zone = 0.4

   133 ; Table header

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

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

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

   137   text1 = table_header_name

   138   res1               = True

   139   res1@txFontHeightF = 0.03

   140   res1@gsFillColor   = "CornFlowerBlue"

   142 ; Column header (equally space in x2)

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

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

   145   y2    = y1                      ; same as y1

   146   text2 = col_header_zone

   147   res2               = True

   148   res2@txFontHeightF = 0.015

   149   res2@gsFillColor   = "Gray"

   151 ; Row header (equally space in y2)

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

   153   x3    = x1                      ; same as x1

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

   155   text3 = row_header_zone

   156   res3               = True

   157   res3@txFontHeightF = 0.02

   158   res3@gsFillColor   = "Gray"

   160 ; Main table body

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

   162   x4    = x2                      ; Start and end x

   163   y4    = y3                      ; Start and end Y

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

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

   167   color_fill4      = "white"

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

   170   res4               = True       ; Set up resource list

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

   172   res4@txFontHeightF = 0.02

   173   res4@gsFillColor   = color_fill4

   175   delete (color_fill4)

   177 ;*******************************************************************

   178 ; for table -- station

   179 ;*******************************************************************

   180 ; column (not including header column)

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

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

   183   ncol_sta       = dimsizes(col_header_sta )

   184 ;-------------------------------------------------------------------

   186 ; Table header

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

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

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

   190   res5               = True

   191   res5@txFontHeightF = 0.02

   192   res5@gsFillColor   = "CornFlowerBlue"

   194 ; Column header (equally space in x2)

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

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

   197   y6    = y5                      ; same as y1

   198   text6 = col_header_sta

   199   res6               = True

   200   res6@txFontHeightF = 0.012

   201   res6@gsFillColor   = "Gray"

   203 ; Row header (equally space in y2)

   205   res7               = True

   206   res7@txFontHeightF = 0.015

   207   res7@gsFillColor   = "Gray"

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

   209 ; for station line plot

   211 ; for x-axis in xyplot

   212   mon = ispan(1,12,1)

   213   mon@long_name = "month"

   215   plot_type = "ps"

   216   plot_type_new = "png"

   218   res                   = True                      ; plot mods desired

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

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

   222 ; Add a boxed legend using the more simple method

   224   res@pmLegendDisplayMode    = "Always"

   225 ; res@pmLegendWidthF         = 0.1

   226   res@pmLegendWidthF         = 0.08

   227   res@pmLegendHeightF        = 0.06

   228 ; res@pmLegendOrthogonalPosF = -1.17

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

   230   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   232 ; res@pmLegendParallelPosF   =  0.18

   233   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   235 ; res@lgPerimOn             = False

   236   res@lgLabelFontHeightF     = 0.015

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

   238 ;-------------------------------------------------------------------

   240   do z = 0,nzone-1

   242   if (z .eq. 0) then

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

   244      zone = "60N-90N"

   245      score_max = 5.0

   246 ;    index of stations in this zone

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

   248 ;    print (ind_z)

   249 ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

   250 ;    print (val_ob(ind_z,:))

   251 ;    print (val_model(ind_z,:))

   252   end if

   254   if (z .eq. 1) then

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

   256      zone = "30N-60N"

   257      score_max = 5.0

   258 ;    index of stations in this zone

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

   260 ;    print (ind_z)

   261 ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

   262 ;    print (val_ob(ind_z,:))

   263 ;    print (val_model(ind_z,:))

   264   end if

   266   if (z .eq. 2) then

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

   268      zone = "EQ-30N"

   269      score_max = 5.0

   270 ;    index of stations in this zone

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

   272 ;    print (ind_z)

   273 ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

   274 ;    print (val_ob(ind_z,:))

   275 ;    print (val_model(ind_z,:))

   276   end if

   278   if (z .eq. 3) then

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

   280      zone = "90S-EQ"

   281      score_max = 5.0

   282 ;    index of stations in this zone

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

   284 ;    print (ind_z)

   285 ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

   286 ;    print (val_ob(ind_z,:))

   287 ;    print (val_model(ind_z,:))

   288   end if

   290  npts = dimsizes(ind_z)

   291  print (npts)

   292 ;-------------------------------------------------------------------------

   293 ; for metric table computation

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

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

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

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

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

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

   303 ;----------------------

   304 ; for table -- station

   305 ;----------------------

   307 ; row (not including header row)

   308   nrow_sta       = npts

   310 ; Table header

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

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

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

   314   res5               = True

   315   res5@txFontHeightF = 0.02

   316   res5@gsFillColor   = "CornFlowerBlue"

   318 ; Column header (equally space in x2)

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

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

   321   y6    = y1                      ; same as y1

   322   text6 = col_header_sta

   323   res6               = True

   324   res6@txFontHeightF = 0.012

   325   res6@gsFillColor   = "Gray"

   327 ; Row header (equally space in y2)

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

   329   x7    = x1                      ; same as x1

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

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

   332   res7               = True

   333   res7@txFontHeightF = 0.015

   334   res7@gsFillColor   = "Gray"

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

   336 ; for station line plot

   338   npts_str = ""

   339   npts_str = npts

   340 ; print (npts_str)

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

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

   345   plot_data!0 = "case"

   346   plot_data!1 = "month"

   347   plot_data!2 = "pts"

   348   plot_data@long_name   = "CO2 Seasonal"

   350   plot_data_0!0 = "month"

   351   plot_data_0!1 = "pts"

   352   plot_data_0@long_name = "CO2"

   353 ;--------------------------------------------------------------------------

   354  do n=0,npts-1

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

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

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

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

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

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

   363     print (sta(ind_z(n))+"/"+lat(ind_z(n))+"/"+lon(ind_z(n))+"/"+amp_ratio_sta(n)+"/"+ccr_sta(n)+"/"+M_sta(n)+"/"+score_sta(n))

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

   365 ; for station line plot

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

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

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

   372     plot_name = sta(ind_z(n))

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

   374 ;   print (title)

   375 ;   print (plot_name)

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

   378 ;------------------------------------------

   379 ;   for panel plot

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

   382     res@gsnFrame     = False                   ; Do not draw plot

   383     res@gsnDraw      = False                   ; Do not advance frame

   385     pres                            = True     ; panel plot mods desired

   386     pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   387                                                ; indiv. plots in panel

   388     pres@gsnMaximize                = True     ; fill the page

   389 ;------------------------------------------

   390     res@tiMainString = title                           ; add title

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

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

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

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

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

   400     clear (wks)

   401 ;---------------------------------------------------------------------------

   402  end do

   403 ;-------------------------------------------------------------------------

   404 ;   for line plot in a zone

   406     plot_name = "All_"+npts_str

   407     title = plot_name + " in "+ zone

   408 ;   print (title)

   409 ;   print (plot_name)

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

   412 ;-----------------------------------------

   413 ;   for panel plot

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

   416     res@gsnFrame     = False                   ; Do not draw plot

   417     res@gsnDraw      = False                   ; Do not advance frame

   419     pres                            = True     ; panel plot mods desired

   420     pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   421                                                ; indiv. plots in panel

   422     pres@gsnMaximize                = True     ; fill the page

   423 ;-----------------------------------------

   424     res@tiMainString = title                                     ; add title

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

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

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

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

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

   435     clear (wks)

   436 ;   delete (ind_z)

   437     delete (plot_data)

   438     delete (plot_data_0)

   439 ;---------------------------------------------------------------------------

   440 ; for zone table value

   442  amp_ratio_zone = avg(amp_ratio_sta)

   443  ccr_zone       = avg(ccr_sta)

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

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

   447  print (npts+"/"+amp_ratio_zone+"/"+ccr_zone+"/"+M_zone+"/"+score_zone)

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

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

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

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

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

   454 ;---------------------------------------------------------------------------

   455 ; plot station table

   456 ;----------------------------------

   457 ; header value for station table

   458   text5  = zone

   460 ; row value for station table

   461   text7  = sta(ind_z)

   463   if (z .eq. 0) then

   464      table_length_sta = 0.5

   465   end if

   466   if (z .eq. 1) then

   467      table_length_sta = 0.995

   468   end if

   469   if (z .eq. 2) then

   470      table_length_sta = 0.8

   471   end if

   472   if (z .eq. 3) then

   473      table_length_sta = 0.8

   474   end if

   476   x7 = x5

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

   479 ; Main table body

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

   481   x8    = x6                      ; Start and end x

   482   y8    = y7                      ; Start and end Y

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

   485   color_fill8      = text8

   486   color_fill8      = "white"

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

   489   res8               = True       ; Set up resource list

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

   491   res8@txFontHeightF = 0.02

   492   res8@gsFillColor   = color_fill8

   494   delete (color_fill8)

   496 ; table value for station table

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

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

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

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

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

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

   504   plot_name = "table_sta." + zone

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

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

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

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

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

   513   frame(wks)

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

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

   518  delete (ind_z)

   519  delete (amp_model)

   520  delete (amp_ob)

   521  delete (amp_ratio_sta)

   522  delete (ccr_sta)

   523  delete (M_sta)

   524  delete (score_sta)

   525  delete (text7)

   526  delete (text8)

   527  delete (res7)

   528  delete (res8)

   529  clear (wks)

   530  end do

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

   532 ; plot zone table

   533 ;**************************************************

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

   536   text4(4,1) = 0.

   537   text4(4,2) = 0.

   538   text4(4,3) = 0.

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

   541   plot_name = "table_zone"

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

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

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

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

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

   549   frame(wks)

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

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

   553 ;-------------------------------------------------------------------

   554   system("mkdir temp")

   555   system("mv *.png temp")

   556   system("tar cf temp.tar temp")

   557 ;-------------------------------------------------------------------

   558 end

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