clamp: co2/19.metric_plot.ncl@0c6405ab2ff4

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

     2 ; xy_4.ncl

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

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

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

     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 ; printVarSummary (val_model)

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

   103 ; remove annual mean

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

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

   107   nzone = 4

   109 ; column

   110   case_zone = (/"Stations","Amplitude Ratio", \

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

   112   nCase_zone = dimsizes(case_zone )

   114 ; row

   115   var_zone = (/"60N-90N","30N-60N","EQ-30N","90S-EQ","Total"/)

   116   nVar_zone = dimsizes(var_zone)

   118 ; arrays to be passed to diagram.

   119   case_value_zone = new ((/nCase_zone, nVar_zone/),float )

   121   do z = 0,nzone-1

   123   if (z .eq. 0) then

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

   125      zone = "60N-90N"

   126      score_max = 5.0

   127 ;    index of stations in this zone

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

   129 ;    print (ind_z)

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

   131 ;    print (val_ob(ind_z,:))

   132 ;    print (val_model(ind_z,:))

   133   end if

   135   if (z .eq. 1) then

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

   137      zone = "30N-60N"

   138      score_max = 5.0

   139 ;    index of stations in this zone

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

   141 ;    print (ind_z)

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

   143 ;    print (val_ob(ind_z,:))

   144 ;    print (val_model(ind_z,:))

   145   end if

   147   if (z .eq. 2) then

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

   149      zone = "EQ-30N"

   150      score_max = 5.0

   151 ;    index of stations in this zone

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

   153 ;    print (ind_z)

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

   155 ;    print (val_ob(ind_z,:))

   156 ;    print (val_model(ind_z,:))

   157   end if

   159   if (z .eq. 3) then

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

   161      zone = "90S-EQ"

   162      score_max = 5.0

   163 ;    index of stations in this zone

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

   165 ;    print (ind_z)

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

   167 ;    print (val_ob(ind_z,:))

   168 ;    print (val_model(ind_z,:))

   169   end if

   171  npts = dimsizes(ind_z)

   172  print (npts)

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

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

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

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

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

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

   182  do n=0,npts-1

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

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

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

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

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

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

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

   192  end do

   194  amp_ratio_zone = avg(amp_ratio_sta)

   195  ccr_zone       = avg(ccr_sta)

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

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

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

   201   case_value_zone(0,z) = npts

   202   case_value_zone(1,z) = (/amp_ratio_zone/)

   203   case_value_zone(2,z) = (/ccr_zone/)

   204   case_value_zone(3,z) = (/M_zone/)

   205   case_value_zone(4,z) = (/score_zone/)

   206 ;**************************************************

   207 ; plot station table

   208 ;**************************************************

   209 ; column for station table

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

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

   212   nCase_sta  = dimsizes(case_sta )

   214 ; row for station table

   215   var_sta  = sta(ind_z)

   216   nVar_sta = dimsizes(var_sta)

   218 ; arrays to be passed to diagram.

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

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

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

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

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

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

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

   228 ;**************************************************

   229 ; plot station table

   230 ;**************************************************

   231   tt_opt        = True

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

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

   235   tt_opt@tableTitle = "Zone "+ zone

   236   plot_name         = "table_sta." + zone

   238   metrics_table(plot_name, var_sta, case_sta , case_value_sta, tt_opt)

   240  delete (ind_z)

   241  delete (amp_model)

   242  delete (amp_ob)

   243  delete (amp_ratio_sta)

   244  delete (ccr_sta)

   245  delete (M_sta)

   246  delete (score_sta)

   247  delete (var_sta)

   248  delete (case_value_sta)

   249  end do

   250 ;**************************************************

   251 ; plot zone table

   252 ;**************************************************

   253   case_value_zone(0,4) = case_value_zone(0,0)+case_value_zone(0,1)+case_value_zone(0,2)+case_value_zone(0,3)

   254   case_value_zone(1,4) = 0.

   255   case_value_zone(2,4) = 0.

   256   case_value_zone(3,4) = 0.

   257   case_value_zone(4,4) = case_value_zone(4,0)+case_value_zone(4,1)+case_value_zone(4,2)+case_value_zone(4,3)

   259   tt_opt        = True

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

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

   263   tt_opt@tableTitle = "Zone"

   264   plot_name         = "table_zone"

   266   metrics_table(plot_name, var_zone, case_zone , case_value_zone, tt_opt)

   268 end

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