Basic Plotting with MOM6

This notebook will take you through basic plotting of MOM6 data. We use Xarray’s plotting functionality but it is just a thin wrapper around the popular matplotlib library.

Before we get started, you may choose to run this notebook on LEAP-Pangeo hub or Binder!

This notebook is inspired by the MOM6-AnalysisCodebook docs

import xarray as xr
import cartopy.crs as ccrs
import cmocean
import matplotlib.pyplot as plt

ds = xr.open_dataset("./data/thetao_ePBL_data_timemean_1978-2002.nc")
ds

<xarray.Dataset>
Dimensions:  (lat: 180, lon: 360, time: 1, z_l: 75)
Coordinates:
  * lat      (lat) float64 -89.5 -88.5 -87.5 -86.5 -85.5 ... 86.5 87.5 88.5 89.5
  * lon      (lon) float64 0.5 1.5 2.5 3.5 4.5 ... 355.5 356.5 357.5 358.5 359.5
  * time     (time) object 1988-01-01 09:00:00
  * z_l      (z_l) float64 1.0 3.0 5.0 7.0 ... 4.668e+03 5.339e+03 6.098e+03
Data variables:
    thetao   (time, z_l, lat, lon) float32 ...
Attributes:
    title:                 OM4p25_JRA55do1.4_r6_cycle1_MLE1_zgrid
    associated_files:      areacello: 19580101.ocean_static.nc
    grid_type:             regular
    grid_tile:             N/A
    code_release_version:  19.1
    git_hash:              c97c1f6afab0df3545ee82c337d51accc4205bc3
    creationtime:          Wed Dec 28 08:43:57 2022
    hostname:              pp206
    history:               Wed Oct 11 11:21:04 2023: ncks -v thetao ePBL_data...
    comment:               FMS time averaging, version 3.0, precision=double
    NCO:                   netCDF Operators version 5.1.4 (Homepage = http://...

xarray.Dataset

• Dimensions:
  • lat: 180
  • lon: 360
  • time: 1
  • z_l: 75
• Coordinates: (4)
  • lat
    (lat)
    float64
    -89.5 -88.5 -87.5 ... 88.5 89.5

    long_name :

    latitude

    units :

    degrees_N

    axis :

    Y

    bounds :

    lat_bnds

    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

    long_name :

    longitude

    units :

    degrees_E

    axis :

    X

    bounds :

    lon_bnds

    array([  0.5,   1.5,   2.5, ..., 357.5, 358.5, 359.5])

  • time
    (time)
    object
    1988-01-01 09:00:00

    long_name :

    time

    axis :

    T

    calendar_type :

    JULIAN

    bounds :

    time_bnds

    cell_methods :

    time: mean

    array([cftime.DatetimeJulian(1988, 1, 1, 9, 0, 0, 0, has_year_zero=False)],
          dtype=object)

  • z_l
    (z_l)
    float64
    1.0 3.0 5.0 ... 5.339e+03 6.098e+03

    units :

    meters

    long_name :

    Depth at cell center

    axis :

    Z

    positive :

    down

    edges :

    z_i

    array([1.000000e+00, 3.000000e+00, 5.000000e+00, 7.000000e+00, 9.000000e+00,
           1.100000e+01, 1.300000e+01, 1.500000e+01, 1.700000e+01, 1.900000e+01,
           2.100000e+01, 2.300000e+01, 2.500000e+01, 2.700000e+01, 2.900000e+01,
           3.100000e+01, 3.300000e+01, 3.500000e+01, 3.700000e+01, 3.900000e+01,
           4.100500e+01, 4.301500e+01, 4.502500e+01, 4.704000e+01, 4.906500e+01,
           5.110500e+01, 5.316500e+01, 5.524500e+01, 5.735500e+01, 5.951000e+01,
           6.172000e+01, 6.400000e+01, 6.637000e+01, 6.886000e+01, 7.150500e+01,
           7.434500e+01, 7.743500e+01, 8.083500e+01, 8.461500e+01, 8.887000e+01,
           9.371000e+01, 9.926500e+01, 1.056950e+02, 1.131800e+02, 1.219400e+02,
           1.322350e+02, 1.443650e+02, 1.586800e+02, 1.755900e+02, 1.955750e+02,
           2.191800e+02, 2.470400e+02, 2.798900e+02, 3.185650e+02, 3.640200e+02,
           4.173500e+02, 4.798000e+02, 5.527850e+02, 6.379050e+02, 7.369650e+02,
           8.520100e+02, 9.853400e+02, 1.139540e+03, 1.317505e+03, 1.522475e+03,
           1.758080e+03, 2.028360e+03, 2.337815e+03, 2.691455e+03, 3.094835e+03,
           3.554110e+03, 4.076095e+03, 4.668325e+03, 5.339115e+03, 6.097620e+03])

• Data variables: (1)
  • thetao
    (time, z_l, lat, lon)
    float32
    ...

    units :

    degC

    long_name :

    Sea Water Potential Temperature

    cell_measures :

    volume: volcello area: areacello

    standard_name :

    sea_water_potential_temperature

    interp_method :

    conserve_order1

    cell_methods :

    area:mean z_l:mean yh:mean xh:mean time: mean

    time_avg_info :

    average_T1,average_T2,average_DT

    [4860000 values with dtype=float32]

• Indexes: (4)
  • lat
    PandasIndex

    PandasIndex(Index([-89.5, -88.5, -87.5, -86.5, -85.5, -84.5, -83.5, -82.5, -81.5, -80.5,
           ...
            80.5,  81.5,  82.5,  83.5,  84.5,  85.5,  86.5,  87.5,  88.5,  89.5],
          dtype='float64', name='lat', length=180))

  • lon
    PandasIndex

    PandasIndex(Index([               0.5, 1.5000000000000002,                2.5,
                          3.5,                4.5,                5.5,
                          6.5,  7.499999999999999,                8.5,
                          9.5,
           ...
           350.50000000000006,              351.5,              352.5,
                        353.5,              354.5,              355.5,
                        356.5,              357.5,              358.5,
           359.50000000000006],
          dtype='float64', name='lon', length=360))

  • time
    PandasIndex

    PandasIndex(CFTimeIndex([1988-01-01 09:00:00],
                dtype='object', length=1, calendar='julian', freq=None))

  • z_l
    PandasIndex

    PandasIndex(Index([               1.0,                3.0,                5.0,
                          7.0,                9.0,               11.0,
                         13.0,               15.0,               17.0,
                         19.0,               21.0,               23.0,
                         25.0,               27.0,               29.0,
                         31.0,               33.0,               35.0,
                         37.0,               39.0, 41.004999999999995,
                       43.015,  45.02499999999999,  47.03999999999999,
                       49.065,             51.105,  53.16499999999999,
            55.24499999999999,  57.35499999999999,  59.50999999999999,
                        61.72,               64.0,              66.37,
                        68.86,             71.505,             74.345,
                       77.435,  80.83500000000001,  84.61500000000001,
                        88.87,              93.71,  99.26499999999999,
                      105.695, 113.17999999999999,             121.94,
           132.23499999999999, 144.36499999999998, 158.67999999999998,
           175.58999999999997,            195.575, 219.17999999999998,
                       247.04,             279.89, 318.56500000000005,
           364.02000000000004,             417.35, 479.80000000000007,
            552.7850000000001,            637.905,            736.965,
                       852.01,  985.3400000000001,            1139.54,
                     1317.505,           1522.475,            1758.08,
                      2028.36, 2337.8149999999996,           2691.455,
                     3094.835,            3554.11, 4076.0950000000003,
                     4668.325,           5339.115,            6097.62],
          dtype='float64', name='z_l'))

• Attributes: (11)

  title :

  OM4p25_JRA55do1.4_r6_cycle1_MLE1_zgrid

  associated_files :

  areacello: 19580101.ocean_static.nc

  grid_type :

  regular

  grid_tile :

  N/A

  code_release_version :

  19.1

  git_hash :

  c97c1f6afab0df3545ee82c337d51accc4205bc3

  creationtime :

  Wed Dec 28 08:43:57 2022

  hostname :

  pp206

  history :

  Wed Oct 11 11:21:04 2023: ncks -v thetao ePBL_data_timemean_1978_2002.nc thetao_ePBL_data_timemean_1978-2002.nc Fri Aug 11 10:13:48 2023: ncra ocean_monthly_z_1x1deg.1958-1962.ann.nc ocean_monthly_z_1x1deg.1963-1967.ann.nc ocean_monthly_z_1x1deg.1968-1972.ann.nc ocean_monthly_z_1x1deg.1973-1977.ann.nc ocean_monthly_z_1x1deg.1978-1982.ann.nc ocean_monthly_z_1x1deg.1983-1987.ann.nc ocean_monthly_z_1x1deg.1988-1992.ann.nc ocean_monthly_z_1x1deg.1993-1997.ann.nc ocean_monthly_z_1x1deg.1998-2002.ann.nc ocean_monthly_z_1x1deg.2003-2007.ann.nc ocean_monthly_z_1x1deg.2008-2012.ann.nc ocean_monthly_z_1x1deg.2013-2017.ann.nc ePBL_data_timemean_1978_2002.nc fregrid --standard_dimension --input_mosaic ocean_mosaic.nc --input_file annual --interp_method conserve_order1 --remap_file fregrid_remap_file_360_by_180.nc --nlon 360 --nlat 180 --scalar_field volcello,thetao,so,obvfsq,uo,vo,umo,vmo,vhml,N2,time_bnds --output_file out.nc

  comment :

  FMS time averaging, version 3.0, precision=double

  NCO :

  netCDF Operators version 5.1.4 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)

Let us first plot thetao from the file

ds["thetao"].isel(time=0, z_l=0).plot(cmap=cmocean.cm.thermal, figsize=[10, 6])

<matplotlib.collections.QuadMesh at 0x7ff4b9ce3340>

Now, let us use a map projection using cartopy. We will also add our own colorbar to the plot.

plt.figure(figsize=[14, 6])
subplot_kws = dict(projection=ccrs.PlateCarree(), facecolor="white")
p = (
    ds["thetao"]
    .isel(time=0, z_l=0)
    .plot(
        vmin=-2,
        vmax=32,
        cmap=cmocean.cm.thermal,
        subplot_kws=subplot_kws,
        transform=ccrs.PlateCarree(),
        add_labels=False,
        add_colorbar=False,
    )
)

# add separate colorbar
cb = plt.colorbar(p, ticks=[0, 5, 10, 15, 20, 25, 30], shrink=0.95)
cb.ax.tick_params(labelsize=18)

You can pick any projection from the cartopy list but, whichever projection you use, you still have to set transform to ccrs.PlateCarree(). Let us use a different projection and add coastlines and gridlines to the plot.

plt.figure(figsize=[14, 6])
subplot_kws = dict(projection=ccrs.Robinson(), facecolor="white")
p = (
    ds["thetao"]
    .isel(time=0, z_l=0)
    .plot(
        vmin=-2,
        vmax=32,
        cmap=cmocean.cm.thermal,
        subplot_kws=subplot_kws,
        transform=ccrs.PlateCarree(),
        add_labels=False,
        add_colorbar=False,
    )
)

p.axes.coastlines()
p.axes.gridlines(color="black", alpha=0.5, linestyle="--")

# add separate colorbar
cb = plt.colorbar(p, ticks=[0, 5, 10, 15, 20, 25, 30], shrink=0.95)
cb.ax.tick_params(labelsize=18)

/home/runner/miniconda3/envs/DGM2lines/lib/python3.9/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

Another projection!

subplot_kws = dict(
    projection=ccrs.SouthPolarStereo(central_longitude=-120.0), facecolor="white"
)

plt.figure(figsize=[12, 8])
p = (
    ds["thetao"]
    .isel(time=0, z_l=0)
    .plot(
        vmin=-2,
        vmax=24,
        cmap=cmocean.cm.thermal,
        subplot_kws=subplot_kws,
        transform=ccrs.PlateCarree(),
        add_labels=False,
        add_colorbar=False,
    )
)

p.axes.gridlines(color="black", alpha=0.5, linestyle="--")
p.axes.set_extent([-300, 60, -40, -90], ccrs.PlateCarree())

# add separate colorbar
cb = plt.colorbar(
    p, ticks=[-2, 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24], shrink=0.99
)
cb.ax.tick_params(labelsize=18)