clamp: co2/42.metric+lines.ncl@4be95183fbcd

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

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

     3 ; ***********************************************

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

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

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

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

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

     9 begin

    10 ;************************************************

    11 ; read in data: observed

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

    13  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/co2/"

    14  fili  = "co2_globalView_98.nc"

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

    16  val   = g->CO2_SEAS

    17  lon   = g->LON

    18  lat   = g->LAT

    19  sta   = chartostring(g->STATION)

    20  delete (g)

    22 ;print (sta(0))

    24  ncase = dimsizes(lat)

    25 ;print (ncase)

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

    28 ; get only the lowest level at each station

    29 ;**************************************************************

    30  lat_tmp = lat

    31  lat_tmp@_FillValue = 1.e+36

    33  do n = 0,ncase-1

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

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

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

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

    38        end if

    39        delete (indexes)

    40     end if

    41  end do

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

    44 ;print (dimsizes(indexes))

    45 ;print (indexes)

    47  lat_ob = lat(indexes)

    48  lon_ob = lon(indexes)

    49  val_ob = val(indexes,:)

    50 ;printVarSummary (val_ob)

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

    53 ;************************************************

    54 ; read in model data

    55 ;************************************************

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

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

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

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

    61   x     = g->CO2

    62   xi    = g->lon

    63   yi    = g->lat

    64   xdim  = dimsizes(x)

    65   nlev  = xdim(1)

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

    67 ; printVarSummary (y)

    69 ; get the co2 at the lowest level

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

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

    73 ; Model_units [=] kgCO2 / kgDryAir

    74 ; 28.966 = molecular weight of dry air

    75 ; 44.       = molecular weight of CO2

    76 ; u mol = 1e-6 mol

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

    79   y      = y * factor

    81   y@_FillValue = 1.e36

    82   y@units      = "u mol/mol"

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

    84 ; printVarSummary (y)

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

    87 ; interpolate into observed station

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

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

    91   print (yi(0))

    92   yi(0) = -90.

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

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

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

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

   100   val_model_0 = val_model

   101 ; printVarSummary (val_model)

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

   104 ; remove annual mean

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

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

   108   nzone = 4

   110 ;*******************************************************************

   111 ; for table -- zone

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

   113 ; table header name

   114   table_header_name = "Zone"

   116 ; column (not including header column)

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

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

   119   ncol_zone       = dimsizes(col_header_zone )

   121 ; row (not including header row)

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

   123   nrow_zone       = dimsizes(row_header_zone)

   125 ; arrays to be passed to table.

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

   127 ;--------------------------------------------------------------------

   129   table_length_zone = 0.4

   131 ; Table header

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

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

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

   135   text1 = table_header_name

   136   res1               = True

   137   res1@txFontHeightF = 0.03

   138   res1@gsFillColor   = "CornFlowerBlue"

   140 ; Column header (equally space in x2)

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

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

   143   y2    = y1                      ; same as y1

   144   text2 = col_header_zone

   145   res2               = True

   146   res2@txFontHeightF = 0.015

   147   res2@gsFillColor   = "Gray"

   149 ; Row header (equally space in y2)

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

   151   x3    = x1                      ; same as x1

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

   153   text3 = row_header_zone

   154   res3               = True

   155   res3@txFontHeightF = 0.02

   156   res3@gsFillColor   = "Gray"

   158 ; Main table body

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

   160   x4    = x2                      ; Start and end x

   161   y4    = y3                      ; Start and end Y

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

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

   165   color_fill4      = "white"

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

   168   res4               = True       ; Set up resource list

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

   170   res4@txFontHeightF = 0.02

   171   res4@gsFillColor   = color_fill4

   173   delete (color_fill4)

   175 ;*******************************************************************

   176 ; for table -- station

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

   178 ; table header name

   179   table_header_name = "Zone"

   181 ; column (not including header column)

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

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

   184   ncol_zone       = dimsizes(col_header_zone )

   186 ; row (not including header row)

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

   188   nrow_zone       = dimsizes(row_header_zone)

   190 ; arrays to be passed to table.

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

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

   194   table_length_zone = 0.4

   196 ; Table header

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

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

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

   200   text1 = table_header_name

   201   res1               = True

   202   res1@txFontHeightF = 0.03

   203   res1@gsFillColor   = "CornFlowerBlue"

   205 ; Column header (equally space in x2)

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

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

   208   y2    = y1                      ; same as y1

   209   text2 = col_header_zone

   210   res2               = True

   211   res2@txFontHeightF = 0.015

   212   res2@gsFillColor   = "Gray"

   214 ; Row header (equally space in y2)

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

   216   x3    = x1                      ; same as x1

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

   218   text3 = row_header_zone

   219   res3               = True

   220   res3@txFontHeightF = 0.02

   221   res3@gsFillColor   = "Gray"

   223 ; Main table body

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

   225   x4    = x2                      ; Start and end x

   226   y4    = y3                      ; Start and end Y

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

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

   230   color_fill4      = "white"

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

   233   res4               = True       ; Set up resource list

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

   235   res4@txFontHeightF = 0.02

   236   res4@gsFillColor   = color_fill4

   238   delete (color_fill4)

   239 ;--------------------------------------------------------------

   240 ; for station line plot

   242 ; for x-axis in xyplot

   243   mon = ispan(1,12,1)

   244   mon@long_name = "month"

   246   plot_type = "ps"

   247   plot_type_new = "png"

   249   res                   = True                      ; plot mods desired

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

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

   253 ; Add a boxed legend using the more simple method

   255   res@pmLegendDisplayMode    = "Always"

   256 ; res@pmLegendWidthF         = 0.1

   257   res@pmLegendWidthF         = 0.08

   258   res@pmLegendHeightF        = 0.06

   259 ; res@pmLegendOrthogonalPosF = -1.17

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

   261   res@pmLegendOrthogonalPosF = -0.30  ;(downward)

   263 ; res@pmLegendParallelPosF   =  0.18

   264   res@pmLegendParallelPosF   =  0.23  ;(rightward)

   266 ; res@lgPerimOn             = False

   267   res@lgLabelFontHeightF     = 0.015

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

   269 ;-------------------------------------------------------------------

   271   do z = 0,nzone-1

   273   if (z .eq. 0) then

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

   275      zone = "60N-90N"

   276      score_max = 5.0

   277 ;    index of stations in this zone

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

   279 ;    print (ind_z)

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

   281 ;    print (val_ob(ind_z,:))

   282 ;    print (val_model(ind_z,:))

   283   end if

   285   if (z .eq. 1) then

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

   287      zone = "30N-60N"

   288      score_max = 5.0

   289 ;    index of stations in this zone

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

   291 ;    print (ind_z)

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

   293 ;    print (val_ob(ind_z,:))

   294 ;    print (val_model(ind_z,:))

   295   end if

   297   if (z .eq. 2) then

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

   299      zone = "EQ-30N"

   300      score_max = 5.0

   301 ;    index of stations in this zone

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

   303 ;    print (ind_z)

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

   305 ;    print (val_ob(ind_z,:))

   306 ;    print (val_model(ind_z,:))

   307   end if

   309   if (z .eq. 3) then

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

   311      zone = "90S-EQ"

   312      score_max = 5.0

   313 ;    index of stations in this zone

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

   315 ;    print (ind_z)

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

   317 ;    print (val_ob(ind_z,:))

   318 ;    print (val_model(ind_z,:))

   319   end if

   321  npts = dimsizes(ind_z)

   322  print (npts)

   323 ;-------------------------------------------------------------------------

   324 ; for metric table

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

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

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

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

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

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

   333 ;-------------------------------------------------------------------------

   334 ; for station line plot

   336   npts_str = ""

   337   npts_str = npts

   338 ; print (npts_str)

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

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

   343   plot_data!0 = "case"

   344   plot_data!1 = "month"

   345   plot_data!2 = "pts"

   346   plot_data@long_name   = "CO2 Seasonal"

   348   plot_data_0!0 = "month"

   349   plot_data_0!1 = "pts"

   350   plot_data_0@long_name = "CO2"

   351 ;--------------------------------------------------------------------------

   352  do n=0,npts-1

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

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

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

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

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

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

   361     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))

   362 ;----------------------------------------------------------------------

   363 ; for station line plot

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

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

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

   370     plot_name = sta(ind_z(n))

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

   372 ;   print (title)

   373 ;   print (plot_name)

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

   376 ;------------------------------------------

   377 ;   for panel plot

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

   380     res@gsnFrame     = False                   ; Do not draw plot

   381     res@gsnDraw      = False                   ; Do not advance frame

   383     pres                            = True     ; panel plot mods desired

   384     pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   385                                                ; indiv. plots in panel

   386     pres@gsnMaximize                = True     ; fill the page

   387 ;------------------------------------------

   388     res@tiMainString = title                           ; add title

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

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

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

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

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

   398     clear (wks)

   399 ;---------------------------------------------------------------------------

   400  end do

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

   402 ;   for line plot in a zone

   404     plot_name = "All_"+npts_str

   405     title = plot_name + " in "+ zone

   406 ;   print (title)

   407 ;   print (plot_name)

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

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

   411 ;   for panel plot

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

   414     res@gsnFrame     = False                   ; Do not draw plot

   415     res@gsnDraw      = False                   ; Do not advance frame

   417     pres                            = True     ; panel plot mods desired

   418     pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

   419                                                ; indiv. plots in panel

   420     pres@gsnMaximize                = True     ; fill the page

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

   422     res@tiMainString = title                                     ; add title

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

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

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

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

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

   433     clear (wks)

   434 ;   delete (ind_z)

   435     delete (plot_data)

   436     delete (plot_data_0)

   437 ;---------------------------------------------------------------------------

   438 ; for zone table

   440  amp_ratio_zone = avg(amp_ratio_sta)

   441  ccr_zone       = avg(ccr_sta)

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

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

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

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

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

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

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

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

   452 ;---------------------------------------------------------------------------

   453 ; column for station table

   454   case_sta   = (/"Latitude","Longitude","Amplitude Ratio", \

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

   456   nCase_sta  = dimsizes(case_sta )

   458 ; row for station table

   459   var_sta  = sta(ind_z)

   460   nVar_sta = dimsizes(var_sta)

   462 ; arrays to be passed to diagram.

   463   case_value_sta = new ((/nCase_sta, nVar_sta/),float )

   465   case_value_sta(0,:) = (/lat(ind_z)/)

   466   case_value_sta(1,:) = (/lon(ind_z)/)

   467   case_value_sta(2,:) = (/amp_ratio_sta/)

   468   case_value_sta(3,:) = (/ccr_sta/)

   469   case_value_sta(4,:) = (/M_sta/)

   470   case_value_sta(5,:) = (/score_sta/)

   472 ;**************************************************

   473 ; plot station table

   474 ;**************************************************

   475   tt_opt        = True

   476   tt_opt@pltType= "png" ; "eps" [default], "pdf", "ps"

   477                         ; "png", "gif" [if you have ImageMajik 'convert']

   479   tt_opt@tableTitle = "Zone "+ zone

   480   plot_name         = "table_sta." + zone

   482 ; metrics_table(plot_name, var_sta, case_sta , case_value_sta, tt_opt)

   484  delete (ind_z)

   485  delete (amp_model)

   486  delete (amp_ob)

   487  delete (amp_ratio_sta)

   488  delete (ccr_sta)

   489  delete (M_sta)

   490  delete (score_sta)

   491  delete (var_sta)

   492  delete (case_value_sta)

   493  end do

   494 ;**************************************************

   495 ; plot zone table

   496 ;**************************************************

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

   499   text4(4,1) = 0.

   500   text4(4,2) = 0.

   501   text4(4,3) = 0.

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

   504   plot_name = "table_zone"

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

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

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

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

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

   512   frame(wks)

   514 end

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