Compare sea surface height model output and observations

Comparing the sea-surface height (ssh) from two different resolution runs. Specifically plot the time-mean and standard deviation of ssh and compare to to observations from the AVISO dataset.

import matplotlib.pyplot as plt
import cosima_cookbook as cc
import numpy as np

import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cft
from dask.distributed import Client
import cmocean as cm
client = Client(n_workers=4)
client

Client

Cluster

  • Workers: 4
  • Cores: 48
  • Memory: 202.49 GB

Start a Cosima cookbook database session

session = cc.database.create_session()
#SSH variable in ACCESS-OM2 models
variable = 'sea_level'

# dates to match the AVISO record
start_time = '1993-01-16'
end_time = '2010-12-16'

Load SSH for a 1\(^{\circ}\) run

Here we can specify the rough start and end times using the “start_time” and “end_time” arguments. But, this yields a dataset that has an extra 12 months on either end, so we slice the data to match the AVISO dataset below.

expt = '1deg_jra55v13_iaf_spinup1_B1'  # 1-deg experiment
ssh1 = cc.querying.getvar(expt, variable, session, start_time=start_time, end_time=end_time)
ssh1 = ssh1.sel(time=slice(start_time, end_time)) # slice to the same time range as AVISO
ssh1
Show/Hide data repr Show/Hide attributes
xarray.DataArray
'sea_level'
  • time: 216
  • yt_ocean: 300
  • xt_ocean: 360
  • dask.array<chunksize=(1, 300, 360), meta=np.ndarray>
    Array Chunk
    Bytes 93.31 MB 432.00 kB
    Shape (216, 300, 360) (1, 300, 360)
    Count 700 Tasks 216 Chunks
    Type float32 numpy.ndarray
    360 300 216
    • yt_ocean
      (yt_ocean)
      float64
      -77.88 -77.63 ... 89.32 89.77
      long_name :
      tcell latitude
      units :
      degrees_N
      cartesian_axis :
      Y
      array([-77.876623, -77.629713, -77.381707, ...,  88.872933,  89.324006,
              89.774476])
    • xt_ocean
      (xt_ocean)
      float64
      -279.5 -278.5 -277.5 ... 78.5 79.5
      long_name :
      tcell longitude
      units :
      degrees_E
      cartesian_axis :
      X
      array([-279.5, -278.5, -277.5, ...,   77.5,   78.5,   79.5])
    • time
      (time)
      object
      1993-01-16 12:00:00 ... 2010-12-16 12:00:00
      long_name :
      time
      cartesian_axis :
      T
      calendar_type :
      GREGORIAN
      bounds :
      time_bounds
      array([cftime.DatetimeGregorian(1993, 1, 16, 12, 0, 0, 0, 5, 16),
             cftime.DatetimeGregorian(1993, 2, 15, 0, 0, 0, 0, 0, 46),
             cftime.DatetimeGregorian(1993, 3, 16, 12, 0, 0, 0, 1, 75), ...,
             cftime.DatetimeGregorian(2010, 10, 16, 12, 0, 0, 0, 5, 289),
             cftime.DatetimeGregorian(2010, 11, 16, 0, 0, 0, 0, 1, 320),
             cftime.DatetimeGregorian(2010, 12, 16, 12, 0, 0, 0, 3, 350)],
            dtype=object)
  • long_name :
    effective sea level (eta_t + patm/(rho0*g)) on T cells
    units :
    meter
    valid_range :
    [-1000. 1000.]
    cell_methods :
    time: mean
    time_avg_info :
    average_T1,average_T2,average_DT
    coordinates :
    geolon_t geolat_t
    standard_name :
    sea_surface_height_above_geoid

Load SSH for a 0.25\(^{\circ}\) run

expt = '025deg_jra55v13_ryf8485_gmredi6'  # 0.25-deg experiment
ssh025 = cc.querying.getvar(expt, variable, session, start_time=start_time, end_time=end_time)
ssh025 = ssh025.sel(time=slice(start_time, end_time)) # slice to the same time range as AVISO
ssh025
Show/Hide data repr Show/Hide attributes
xarray.DataArray
'sea_level'
  • time: 216
  • yt_ocean: 1080
  • xt_ocean: 1440
  • dask.array<chunksize=(1, 540, 720), meta=np.ndarray>
    Array Chunk
    Bytes 1.34 GB 1.56 MB
    Shape (216, 1080, 1440) (1, 540, 720)
    Count 2794 Tasks 864 Chunks
    Type float32 numpy.ndarray
    1440 1080 216
    • yt_ocean
      (yt_ocean)
      float64
      -81.08 -80.97 ... 89.84 89.95
      long_name :
      tcell latitude
      units :
      degrees_N
      cartesian_axis :
      Y
      array([-81.077001, -80.971402, -80.865804, ...,  89.736085,  89.841684,
              89.947282])
    • xt_ocean
      (xt_ocean)
      float64
      -279.9 -279.6 ... 79.62 79.88
      long_name :
      tcell longitude
      units :
      degrees_E
      cartesian_axis :
      X
      array([-279.875, -279.625, -279.375, ...,   79.375,   79.625,   79.875])
    • time
      (time)
      object
      1993-01-16 12:00:00 ... 2010-12-16 12:00:00
      long_name :
      time
      cartesian_axis :
      T
      calendar_type :
      NOLEAP
      bounds :
      time_bounds
      array([cftime.DatetimeNoLeap(1993, 1, 16, 12, 0, 0, 0, 6, 16),
             cftime.DatetimeNoLeap(1993, 2, 15, 0, 0, 0, 0, 1, 46),
             cftime.DatetimeNoLeap(1993, 3, 16, 12, 0, 0, 0, 2, 75), ...,
             cftime.DatetimeNoLeap(2010, 10, 16, 12, 0, 0, 0, 2, 289),
             cftime.DatetimeNoLeap(2010, 11, 16, 0, 0, 0, 0, 5, 320),
             cftime.DatetimeNoLeap(2010, 12, 16, 12, 0, 0, 0, 0, 350)], dtype=object)
  • long_name :
    effective sea level (eta_t + patm/(rho0*g)) on T cells
    units :
    meter
    valid_range :
    [-1000. 1000.]
    cell_methods :
    time: mean
    time_avg_info :
    average_T1,average_T2,average_DT
    coordinates :
    geolon_t geolat_t
    standard_name :
    sea_surface_height_above_geoid

Load AVISO observational data

Use xarray to load ‘zos’: the sea surface height variable name. Here we slice off the first three time points so that the data starts in January and spans an even 18 years.

obs = xr.open_dataset('/g/data/hh5/tmp/cosima/observations/original/zos_AVISO_L4_199210-201012.nc')   # AVISO dataset
obs_ssh = obs.zos.sel(time=slice(start_time, end_time))
obs_ssh
Show/Hide data repr Show/Hide attributes
xarray.DataArray
'zos'
  • time: 216
  • lat: 180
  • lon: 360
  • ...
    [13996800 values with dtype=float32]
    • time
      (time)
      datetime64[ns]
      1993-01-16T12:00:00 ... 2010-12-16T12:00:00
      bounds :
      time_bnds
      axis :
      T
      long_name :
      time
      standard_name :
      time
      array(['1993-01-16T12:00:00.000000000', '1993-02-15T00:00:00.000000000',
             '1993-03-16T12:00:00.000000000', ..., '2010-10-16T12:00:00.000000000',
             '2010-11-16T00:00:00.000000000', '2010-12-16T12:00:00.000000000'],
            dtype='datetime64[ns]')
    • lat
      (lat)
      float64
      -89.5 -88.5 -87.5 ... 88.5 89.5
      bounds :
      lat_bnds
      units :
      degrees_north
      axis :
      Y
      long_name :
      latitude
      standard_name :
      latitude
      array([-89.5, -88.5, -87.5, -86.5, -85.5, -84.5, -83.5, -82.5, -81.5, -80.5,
             -79.5, -78.5, -77.5, -76.5, -75.5, -74.5, -73.5, -72.5, -71.5, -70.5,
             -69.5, -68.5, -67.5, -66.5, -65.5, -64.5, -63.5, -62.5, -61.5, -60.5,
             -59.5, -58.5, -57.5, -56.5, -55.5, -54.5, -53.5, -52.5, -51.5, -50.5,
             -49.5, -48.5, -47.5, -46.5, -45.5, -44.5, -43.5, -42.5, -41.5, -40.5,
             -39.5, -38.5, -37.5, -36.5, -35.5, -34.5, -33.5, -32.5, -31.5, -30.5,
             -29.5, -28.5, -27.5, -26.5, -25.5, -24.5, -23.5, -22.5, -21.5, -20.5,
             -19.5, -18.5, -17.5, -16.5, -15.5, -14.5, -13.5, -12.5, -11.5, -10.5,
              -9.5,  -8.5,  -7.5,  -6.5,  -5.5,  -4.5,  -3.5,  -2.5,  -1.5,  -0.5,
               0.5,   1.5,   2.5,   3.5,   4.5,   5.5,   6.5,   7.5,   8.5,   9.5,
              10.5,  11.5,  12.5,  13.5,  14.5,  15.5,  16.5,  17.5,  18.5,  19.5,
              20.5,  21.5,  22.5,  23.5,  24.5,  25.5,  26.5,  27.5,  28.5,  29.5,
              30.5,  31.5,  32.5,  33.5,  34.5,  35.5,  36.5,  37.5,  38.5,  39.5,
              40.5,  41.5,  42.5,  43.5,  44.5,  45.5,  46.5,  47.5,  48.5,  49.5,
              50.5,  51.5,  52.5,  53.5,  54.5,  55.5,  56.5,  57.5,  58.5,  59.5,
              60.5,  61.5,  62.5,  63.5,  64.5,  65.5,  66.5,  67.5,  68.5,  69.5,
              70.5,  71.5,  72.5,  73.5,  74.5,  75.5,  76.5,  77.5,  78.5,  79.5,
              80.5,  81.5,  82.5,  83.5,  84.5,  85.5,  86.5,  87.5,  88.5,  89.5])
    • lon
      (lon)
      float64
      0.5 1.5 2.5 ... 357.5 358.5 359.5
      bounds :
      lon_bnds
      units :
      degrees_east
      axis :
      X
      long_name :
      longitude
      standard_name :
      longitude
      array([  0.5,   1.5,   2.5, ..., 357.5, 358.5, 359.5])
  • standard_name :
    sea_surface_height_above_geoid
    long_name :
    Sea Surface Height Above Geoid
    units :
    m
    original_name :
    maps_of_absolute_dynamic_topography
    history :
    2011-02-11, 12:02:38, AVISO, Aviso2Cmor 2009-01-01 2009-12-31 2011-08-24T20:50:25Z altered by CMOR: Converted units from 'cm' to 'm'. 2011-08-24T20:50:25Z altered by CMOR: Converted type from 'd' to 'f'.
    original_units :
    cm
    cell_methods :
    time: mean
    cell_measures :
    area: areacello
    associated_files :
    baseURL: http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation gridspecFile: gridspec_ocean_fx_Obs-AVISO_obs_r0i0p0.nc areacello: areacello_fx_Obs-AVISO_obs_r0i0p0.nc

Make comparison plot

Use cartopy to plot the time-mean and standard deviation of both of the model outputs and the AVISO data (obs).

plt.figure(figsize=(20,8), dpi=300)

# Define the contour levels and spacing
mlev = np.arange(-1.6, 1.6, 0.02)  # levels for mean ssh
slev = np.arange( 0.0, 0.3, 0.01)  # levels for std ssh


# mean SSH plots

ax = plt.subplot(2, 3, 1, projection=ccrs.Robinson(central_longitude=-100))
p1 = ssh1.mean(dim='time').plot.contourf(cmap=cm.cm.balance, levels=mlev, add_colorbar=False,
                                         extend='both', transform=ccrs.PlateCarree())
ax.set_title('Mean SSH 1$^{\circ}$')

ax = plt.subplot(2, 3, 2, projection=ccrs.Robinson(central_longitude=-100))
p1 = ssh025.mean(dim='time').plot.contourf(cmap=cm.cm.balance, levels=mlev, add_colorbar=False,
                                          extend='both', transform=ccrs.PlateCarree())
ax.set_title('Mean SSH 0.25$^{\circ}$')


ax = plt.subplot(2, 3, 3, projection=ccrs.Robinson(central_longitude=-100))
p1 = obs_ssh.mean(dim='time').plot.contourf(cmap=cm.cm.balance, levels=mlev, add_colorbar=False,
                                            extend='both', transform=ccrs.PlateCarree())
ax.set_title('AVISO Mean SSH')


# std SSH plots
ax2 = plt.subplot(2, 3, 4, projection=ccrs.Robinson(central_longitude=-100))
p2 = ssh1.std(dim='time').plot.contourf(cmap=cm.cm.deep, levels=slev, add_colorbar=False,
                                        extend='max', transform=ccrs.PlateCarree())
ax2.set_title('SSH Standard Deviation 1$^{\circ}$')


ax2 = plt.subplot(2,3,5,projection=ccrs.Robinson(central_longitude=-100))
p2 = ssh025.std(dim='time').plot.contourf(cmap=cm.cm.deep, levels=slev, add_colorbar=False,
                                          extend='max', transform=ccrs.PlateCarree())
ax2.set_title('SSH Standard Deviation 0.25$^{\circ}$')


ax2 = plt.subplot(2,3,6,projection=ccrs.Robinson(central_longitude=-100))
p2 = obs_ssh.std(dim='time').plot.contourf(cmap=cm.cm.deep, levels=slev,add_colorbar=False,
                                           extend='max', transform=ccrs.PlateCarree())
ax2.set_title('AVISO SSH Standard Deviation')

# Colorbars
ax1 = plt.axes([0.92, 0.57, 0.01, 0.3])
cb = plt.colorbar(p1, cax=ax1, orientation='vertical')
cb.ax.set_ylabel('Mean SSH (m)')

ax2 = plt.axes([0.92, 0.14, 0.01, 0.3])
cb = plt.colorbar(p2, cax=ax2, orientation='vertical')
cb.ax.set_ylabel('SSH Standard Deviation (m)');
/g/data3/hh5/public/apps/miniconda3/envs/analysis3-20.01/lib/python3.7/site-packages/xarray/core/nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)
../_images/Compare_SSH_model_obs_0.png

Download python script: Compare_SSH_model_obs.py

Download Jupyter notebook: Compare_SSH_model_obs.ipynb