Transforming Vertical Coordinates

A common need in the analysis of ocean and atmospheric data is to transform the vertical coordinate from its original coordinate (e.g. depth) to a new coordinate (e.g. density). Xgcm supports this sort of one-dimensional coordinate transform on a single Axis of an xgcm.Grid object using the transform method. Two algorithms are implemented:

• Linear interpolation: Linear interpolation is designed to interpolate intensive quantities (e.g. temperature) from one coordinate to another. This method is suitable when the target coordinate is monotonically increasing or decreasing and the data variable is intensive. For example, you want to visualize oxygen on density surfaces from a z-coordinate ocean model.

  • Logarithmic interpolation: Logarithmic interpolation (which is linear interpolation done after applying a logarithm to the target coordinate) is also available. This method is suitable when variation of the intensive quantity is best related to the logarithm of the target coordinate, rather than the target coordinate itself. For example, you want to analyze data from a sigma-coordinate atmospheric model on isobaric (constant pressure) surfaces.

• Conservative remapping: This algorithm is designed to conserve extensive quantities (e.g. transport, heat content). It requires knowledge of cell bounds in both the source and target coordinate. It also handles non-monotonic target coordinates.

On this page, we explain how to use these coordinate transformation capabilities.

from xgcm import Grid
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt

1D Toy Data Example

First we will create a simple, one-dimensional dataset to illustrate how the transform function works. This dataset contains

• a coordinate called z, representing the original depth coordinate

• a data variable called theta, a function of z, which we want as our new vertical coordinate

• a data variable called phi, also a function of z, which represents the data we want to transform into this new coordinate space

In an oceanic context theta might be density and phi might be oxygen. In an atmospheric context, theta might be potential temperature and phi might be potential vorticity.

z = np.arange(2, 12)
theta = xr.DataArray(np.log(z), dims=['z'], coords={'z': z})
phi = xr.DataArray(np.flip(np.log(z)*0.5+ np.random.rand(len(z))), dims=['z'], coords={'z':z})
ds = xr.Dataset({'phi': phi, 'theta': theta})
ds

<xarray.Dataset>
Dimensions:  (z: 10)
Coordinates:
  * z        (z) int64 2 3 4 5 6 7 8 9 10 11
Data variables:
    phi      (z) float64 1.633 1.165 1.233 1.637 ... 1.657 1.358 1.423 1.133
    theta    (z) float64 0.6931 1.099 1.386 1.609 ... 2.079 2.197 2.303 2.398

xarray.Dataset

• Dimensions:
  • z: 10
• Coordinates: (1)
  • z
    (z)
    int64
    2 3 4 5 6 7 8 9 10 11

    array([ 2,  3,  4,  5,  6,  7,  8,  9, 10, 11])

• Data variables: (2)
  • phi
    (z)
    float64
    1.633 1.165 1.233 ... 1.423 1.133

    array([1.63287665, 1.16507553, 1.23304642, 1.63660282, 1.53322135,
           1.19352428, 1.65746418, 1.35807084, 1.42273616, 1.13307975])

  • theta
    (z)
    float64
    0.6931 1.099 1.386 ... 2.303 2.398

    array([0.69314718, 1.09861229, 1.38629436, 1.60943791, 1.79175947,
           1.94591015, 2.07944154, 2.19722458, 2.30258509, 2.39789527])

• Attributes: (0)

Let’s plot this data. Note that, for a simple 1D profile, we can easily visualize phi in theta space by simply plotting phi vs. theta. This is essentially a form of linear interpolation, performed automatically by matplotlib when it draws lines between the discrete points of our data.

def plot_profile():
    fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=[8,5])
    ds.theta.plot(ax=ax1, y='z', marker='.', yincrease=False)
    ds.phi.plot(ax=ax2, y='z', marker='.', yincrease=False)
    ds.swap_dims({'z': 'theta'}).phi.plot(ax=ax3, y='theta', marker='.', yincrease=False)
    fig.subplots_adjust(wspace=0.5)
    return ax3

plot_profile();

Linear transformation

Ok now lets transform phi to theta coordinates using linear interpolation. A key part of this is to define specific theta levels onto which we want to interpolate the data.

# First create an xgcm grid object
grid = Grid(ds, coords={'Z': {'center':'z'}}, periodic=False)

# define the target values in density, linearly spaced
theta_target = np.linspace(0, 3, 20)

# and transform
phi_transformed = grid.transform(ds.phi, 'Z', theta_target, target_data=ds.theta)
phi_transformed

<xarray.DataArray 'phi' (theta: 20)>
array([       nan,        nan,        nan,        nan,        nan,
       1.52174096, 1.33957156, 1.16664694, 1.20395286, 1.29590695,
       1.5814605 , 1.56436105, 1.30629331, 1.56431596, 1.36623482,
       1.22265461,        nan,        nan,        nan,        nan])
Coordinates:
  * theta    (theta) float64 0.0 0.1579 0.3158 0.4737 ... 2.526 2.684 2.842 3.0

xarray.DataArray

'phi'

• theta: 20

• nan nan nan nan nan 1.522 1.34 ... 1.564 1.366 1.223 nan nan nan nan

  array([       nan,        nan,        nan,        nan,        nan,
         1.52174096, 1.33957156, 1.16664694, 1.20395286, 1.29590695,
         1.5814605 , 1.56436105, 1.30629331, 1.56431596, 1.36623482,
         1.22265461,        nan,        nan,        nan,        nan])

• Coordinates: (1)
  • theta
    (theta)
    float64
    0.0 0.1579 0.3158 ... 2.842 3.0

    array([0.      , 0.157895, 0.315789, 0.473684, 0.631579, 0.789474, 0.947368,
           1.105263, 1.263158, 1.421053, 1.578947, 1.736842, 1.894737, 2.052632,
           2.210526, 2.368421, 2.526316, 2.684211, 2.842105, 3.      ])

• Attributes: (0)

Now let’s see what the result looks like.

ax = plot_profile()
phi_transformed.plot(ax=ax, y='theta', marker='.')

[<matplotlib.lines.Line2D at 0x7fc864c01c90>]

Not too bad. We can increase the number of interpolation levels to capture more of the small scale structure.

target_theta = np.linspace(0,3, 100)
phi_transformed = grid.transform(ds.phi, 'Z', target_theta, target_data=ds.theta)
ax = plot_profile()
phi_transformed.plot(ax=ax, y='theta', marker='.')

[<matplotlib.lines.Line2D at 0x7fc864af5450>]

Note that by default, values of theta_target which lie outside the range of theta have been masked (set to NaN). To disable this behavior, you can pass mask_edges=False; values outside the range of theta will be filled with the nearest valid value.

target_theta = np.linspace(0,3, 60)
phi_transformed = grid.transform(ds.phi, 'Z', target_theta, target_data=ds.theta, mask_edges=False)
ax = plot_profile()
phi_transformed.plot(ax=ax, y='theta', marker='.')

[<matplotlib.lines.Line2D at 0x7fc85b30f3a0>]

Conservative transformation

Conservative transformation is designed to preseve the total sum of phi over the Z axis. It presumes that phi is an extensive quantity, i.e. a quantity that is already volume weighted, with respect to the Z axis: for example, units of Kelvins * meters for heat content, rather than just Kelvins. The conservative method requires more input data at the moment. You have to not only specify the coordinates of the cell centers, but also the cell faces (or bounds/boundaries). In xgcm we achieve this by defining the bounding coordinates as the outer axis position. The target theta values are likewise intepreted as cell boundaries in theta-space. In this way, conservative transformation is similar to calculating a histogram.

# define the cell bounds in depth
zc = np.arange(1,12)+0.5

# add them to the existing dataset
ds = ds.assign_coords({'zc': zc})
ds

<xarray.Dataset>
Dimensions:  (z: 10, zc: 11)
Coordinates:
  * z        (z) int64 2 3 4 5 6 7 8 9 10 11
  * zc       (zc) float64 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5
Data variables:
    phi      (z) float64 1.633 1.165 1.233 1.637 ... 1.657 1.358 1.423 1.133
    theta    (z) float64 0.6931 1.099 1.386 1.609 ... 2.079 2.197 2.303 2.398

xarray.Dataset

• Dimensions:
  • z: 10
  • zc: 11
• Coordinates: (2)
  • z
    (z)
    int64
    2 3 4 5 6 7 8 9 10 11

    array([ 2,  3,  4,  5,  6,  7,  8,  9, 10, 11])

  • zc
    (zc)
    float64
    1.5 2.5 3.5 4.5 ... 9.5 10.5 11.5

    array([ 1.5,  2.5,  3.5,  4.5,  5.5,  6.5,  7.5,  8.5,  9.5, 10.5, 11.5])

• Data variables: (2)
  • phi
    (z)
    float64
    1.633 1.165 1.233 ... 1.423 1.133

    array([1.63287665, 1.16507553, 1.23304642, 1.63660282, 1.53322135,
           1.19352428, 1.65746418, 1.35807084, 1.42273616, 1.13307975])

  • theta
    (z)
    float64
    0.6931 1.099 1.386 ... 2.303 2.398

    array([0.69314718, 1.09861229, 1.38629436, 1.60943791, 1.79175947,
           1.94591015, 2.07944154, 2.19722458, 2.30258509, 2.39789527])

• Attributes: (0)

# Recreate the grid object with a staggered `center`/`outer` coordinate layout
grid = Grid(ds, coords={'Z':{'center':'z', 'outer':'zc'}},
            periodic=False)
grid

<xgcm.Grid>
Z Axis (not periodic, boundary=None):
  * center   z --> outer
  * outer    zc --> center

Currently the target_data(theta in this case) has to be located on the outer coordinate for the conservative method (compared to the center for the linear method).

We can easily interpolate theta on the outer coordinate with the grid object.

ds['theta_outer'] = grid.interp(ds.theta, 'Z', boundary='fill')
ds['theta_outer']

<xarray.DataArray 'theta_outer' (zc: 11)>
array([0.34657359, 0.89587973, 1.24245332, 1.49786614, 1.70059869,
       1.86883481, 2.01267585, 2.13833306, 2.24990484, 2.35024018,
       1.19894764])
Coordinates:
  * zc       (zc) float64 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5

xarray.DataArray

'theta_outer'

• zc: 11

• 0.3466 0.8959 1.242 1.498 1.701 1.869 2.013 2.138 2.25 2.35 1.199

  array([0.34657359, 0.89587973, 1.24245332, 1.49786614, 1.70059869,
         1.86883481, 2.01267585, 2.13833306, 2.24990484, 2.35024018,
         1.19894764])

• Coordinates: (1)
  • zc
    (zc)
    float64
    1.5 2.5 3.5 4.5 ... 9.5 10.5 11.5

    array([ 1.5,  2.5,  3.5,  4.5,  5.5,  6.5,  7.5,  8.5,  9.5, 10.5, 11.5])

• Attributes: (0)

Now lets transform the data using the conservative method. Note that the target values will now be interpreted as cell bounds and not cell centers as before.

# define the target values in density
theta_target = np.linspace(0,3, 20)

# and transform
phi_transformed_cons = grid.transform(ds.phi,
                                      'Z',
                                      theta_target,
                                      method='conservative',
                                      target_data=ds.theta_outer)
phi_transformed_cons

/home/jthielen/develop/xgcm/xgcm/transform.py:247: FutureWarning: ``output_sizes`` should be given in the ``dask_gufunc_kwargs`` parameter. It will be removed as direct parameter in a future version.
  out = xr.apply_ufunc(

<xarray.DataArray 'phi' (theta_outer: 19)>
array([0.        , 0.        , 0.37785116, 0.46936054, 0.46936054,
       0.48939379, 0.53079432, 0.62434095, 0.91765916, 1.18077243,
       1.46775848, 1.57323631, 1.6610291 , 2.16457924, 2.03956194,
       0.        , 0.        , 0.        , 0.        ])
Coordinates:
  * theta_outer  (theta_outer) float64 0.07895 0.2368 0.3947 ... 2.763 2.921

xarray.DataArray

'phi'

• theta_outer: 19

• 0.0 0.0 0.3779 0.4694 0.4694 0.4894 ... 2.165 2.04 0.0 0.0 0.0 0.0

  array([0.        , 0.        , 0.37785116, 0.46936054, 0.46936054,
         0.48939379, 0.53079432, 0.62434095, 0.91765916, 1.18077243,
         1.46775848, 1.57323631, 1.6610291 , 2.16457924, 2.03956194,
         0.        , 0.        , 0.        , 0.        ])

• Coordinates: (1)
  • theta_outer
    (theta_outer)
    float64
    0.07895 0.2368 ... 2.763 2.921

    array([0.078947, 0.236842, 0.394737, 0.552632, 0.710526, 0.868421, 1.026316,
           1.184211, 1.342105, 1.5     , 1.657895, 1.815789, 1.973684, 2.131579,
           2.289474, 2.447368, 2.605263, 2.763158, 2.921053])

• Attributes: (0)

phi_transformed_cons.plot(y='theta_outer', marker='.', yincrease=False)

[<matplotlib.lines.Line2D at 0x7fc85b1c52a0>]

There is no point in comparing phi_transformed_cons directly to phi or the results of linear interoplation, since here we have reinterpreted phi as an extensive quantity. However, we can verify that the sum of the two quantities over the Z axis is exactly the same.

ds.phi.sum().values

array(13.96569797)

phi_transformed_cons.sum().values

array(13.96569797)

Logarithmic Interpolation

Since logarithmic interpolation is most often used when tranforming to pressure (isobaric) coordinates, let’s use a new set of example data:

ds = xr.Dataset(
    {
        'temperature': ('sigma', [271.75452, 272.79956, 274.8517, 279.22043, 296.48782]),
        'pressure': ('sigma', [100180.625, 96250.0, 87369.14, 72186.66, 53718.586]),
    },
    {
        'sigma': [0.9969, 0.9558, 0.8631, 0.7046, 0.5117]
    }
)

Let’s now transform temperature to pressure coordinates. As before, a key part of this is to define specific pressure levels onto which we want to interpolate the data.

# Create an xgcm grid object
grid = Grid(ds, coords={'Z': {'center':'sigma'}}, periodic=False)

# Define isobaric levels of interest (a few standard levels)
isobaric_target_levels = np.array([1.0e5, 8.5e4, 7.0e4])

# and transform
temperature_isobaric = grid.transform(ds.temperature, 'Z', isobaric_target_levels, target_data=ds.pressure, method='log')
temperature_isobaric

<xarray.DataArray 'temperature' (pressure: 3)>
array([271.80163709, 275.48086957, 281.01789905])
Coordinates:
  * pressure  (pressure) float64 1e+05 8.5e+04 7e+04

xarray.DataArray

'temperature'

• pressure: 3

• 271.8 275.5 281.0

  array([271.80163709, 275.48086957, 281.01789905])

• Coordinates: (1)
  • pressure
    (pressure)
    float64
    1e+05 8.5e+04 7e+04

    array([100000.,  85000.,  70000.])

• Attributes: (0)

And finally, let’s see what the result looks like.

fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=[8,5])
ds.pressure.plot(ax=ax1, y='sigma', marker='.', yincrease=False)
ds.temperature.plot(ax=ax2, y='sigma', marker='.', yincrease=False)
ds.swap_dims({'sigma': 'pressure'}).temperature.plot(ax=ax3, y='pressure', marker='.', yincrease=False)
temperature_isobaric.plot(ax=ax3, y='pressure', marker='.')

fig.subplots_adjust(wspace=0.7)

plt.show();

Realistic Data Example

To illustrate these features in a more realistic example, we use data from the CNRM CMIP6 model. These data are available from the Pangeo Cloud Data Library. We can see that this is a full, global, 4D ocean dataset.

import intake
col = intake.open_esm_datastore("https://storage.googleapis.com/cmip6/pangeo-cmip6.json")
cat = col.search(
    source_id = 'CNRM-ESM2-1',
    member_id = 'r1i1p1f2',
    experiment_id = 'historical',
    variable_id= ['thetao','so','vo','areacello'],
    grid_label = 'gn'
)
ddict = cat.to_dataset_dict(zarr_kwargs={'consolidated':True, 'use_cftime':True}, aggregate=False)

--> The keys in the returned dictionary of datasets are constructed as follows:
        'activity_id.institution_id.source_id.experiment_id.member_id.table_id.variable_id.grid_label.zstore.dcpp_init_year.version'

100.00% [4/4 00:00<00:00]

thetao = ddict['CMIP.CNRM-CERFACS.CNRM-ESM2-1.historical.r1i1p1f2.Omon.thetao.gn.gs://cmip6/CMIP6/CMIP/CNRM-CERFACS/CNRM-ESM2-1/historical/r1i1p1f2/Omon/thetao/gn/v20181206/.nan.20181206']
so = ddict['CMIP.CNRM-CERFACS.CNRM-ESM2-1.historical.r1i1p1f2.Omon.so.gn.gs://cmip6/CMIP6/CMIP/CNRM-CERFACS/CNRM-ESM2-1/historical/r1i1p1f2/Omon/so/gn/v20181206/.nan.20181206']
vo = ddict['CMIP.CNRM-CERFACS.CNRM-ESM2-1.historical.r1i1p1f2.Omon.vo.gn.gs://cmip6/CMIP6/CMIP/CNRM-CERFACS/CNRM-ESM2-1/historical/r1i1p1f2/Omon/vo/gn/v20181206/.nan.20181206']
areacello = ddict['CMIP.CNRM-CERFACS.CNRM-ESM2-1.historical.r1i1p1f2.Ofx.areacello.gn.gs://cmip6/CMIP6/CMIP/CNRM-CERFACS/CNRM-ESM2-1/historical/r1i1p1f2/Ofx/areacello/gn/v20181206/.nan.20181206'].areacello

vo = vo.rename({'y':'y_c', 'lon':'lon_v', 'lat':'lat_v', 'bounds_lon':'bounds_lon_v', 'bounds_lat':'bounds_lat_v'})

ds = xr.merge([thetao,so,vo], compat='override')
ds = ds.assign_coords(areacello=areacello.fillna(0))
ds

<xarray.Dataset>
Dimensions:       (y: 294, x: 362, nvertex: 4, lev: 75, axis_nbounds: 2,
                   time: 1980, y_c: 294)
Coordinates: (12/13)
    bounds_lat    (y, x, nvertex) float64 dask.array<chunksize=(294, 362, 4), meta=np.ndarray>
    bounds_lon    (y, x, nvertex) float64 dask.array<chunksize=(294, 362, 4), meta=np.ndarray>
    lat           (y, x) float64 dask.array<chunksize=(294, 362), meta=np.ndarray>
  * lev           (lev) float64 0.5058 1.556 2.668 ... 5.698e+03 5.902e+03
    lev_bounds    (lev, axis_nbounds) float64 dask.array<chunksize=(75, 2), meta=np.ndarray>
    lon           (y, x) float64 dask.array<chunksize=(294, 362), meta=np.ndarray>
    ...            ...
    time_bounds   (time, axis_nbounds) object dask.array<chunksize=(1980, 2), meta=np.ndarray>
    bounds_lat_v  (y_c, x, nvertex) float64 dask.array<chunksize=(294, 362, 4), meta=np.ndarray>
    bounds_lon_v  (y_c, x, nvertex) float64 dask.array<chunksize=(294, 362, 4), meta=np.ndarray>
    lat_v         (y_c, x) float64 dask.array<chunksize=(294, 362), meta=np.ndarray>
    lon_v         (y_c, x) float64 dask.array<chunksize=(294, 362), meta=np.ndarray>
    areacello     (y, x) float32 dask.array<chunksize=(294, 362), meta=np.ndarray>
Dimensions without coordinates: y, x, nvertex, axis_nbounds, y_c
Data variables:
    thetao        (time, lev, y, x) float32 dask.array<chunksize=(4, 75, 294, 362), meta=np.ndarray>
    so            (time, lev, y, x) float32 dask.array<chunksize=(5, 75, 294, 362), meta=np.ndarray>
    vo            (time, lev, y_c, x) float32 dask.array<chunksize=(3, 75, 294, 362), meta=np.ndarray>
Attributes: (12/57)
    CMIP6_CV_version:        cv=6.2.3.0-7-g2019642
    Conventions:             CF-1.7 CMIP-6.2
    EXPID:                   CNRM-ESM2-1_historical_r1i1p1f2
    activity_id:             CMIP
    arpege_minor_version:    6.3.2
    branch_method:           standard
    ...                      ...
    xios_commit:             1442-shuffle
    status:                  2019-11-05;created;by nhn2@columbia.edu
    netcdf_tracking_ids:     hdl:21.14100/9c34b796-c31d-4c1f-be90-21d032267f6...
    version_id:              v20181206
    intake_esm_varname:      None
    intake_esm_dataset_key:  CMIP.CNRM-CERFACS.CNRM-ESM2-1.historical.r1i1p1f...

xarray.Dataset

• Dimensions:
  • y: 294
  • x: 362
  • nvertex: 4
  • lev: 75
  • axis_nbounds: 2
  • time: 1980
  • y_c: 294
• Coordinates: (13)
  • bounds_lat
    (y, x, nvertex)
    float64
    dask.array<chunksize=(294, 362, 4), meta=np.ndarray>

    Array Chunk Bytes 3.25 MiB 3.25 MiB Shape (294, 362, 4) (294, 362, 4) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray

  • bounds_lon
    (y, x, nvertex)
    float64
    dask.array<chunksize=(294, 362, 4), meta=np.ndarray>

    Array Chunk Bytes 3.25 MiB 3.25 MiB Shape (294, 362, 4) (294, 362, 4) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray

  • lat
    (y, x)
    float64
    dask.array<chunksize=(294, 362), meta=np.ndarray>

    bounds :

    bounds_lat

    long_name :

    Latitude

    standard_name :

    latitude

    units :

    degrees_north

    Array Chunk Bytes 831.47 kiB 831.47 kiB Shape (294, 362) (294, 362) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray

  • lev
    (lev)
    float64
    0.5058 1.556 ... 5.902e+03

    bounds :

    lev_bounds

    long_name :

    Vertical T levels

    name :

    lev

    positive :

    down

    units :

    m

    array([5.057600e-01, 1.555855e+00, 2.667682e+00, 3.856280e+00, 5.140361e+00,
           6.543034e+00, 8.092518e+00, 9.822750e+00, 1.177368e+01, 1.399104e+01,
           1.652532e+01, 1.942980e+01, 2.275762e+01, 2.655830e+01, 3.087456e+01,
           3.574020e+01, 4.118002e+01, 4.721189e+01, 5.385064e+01, 6.111284e+01,
           6.902168e+01, 7.761116e+01, 8.692943e+01, 9.704131e+01, 1.080303e+02,
           1.200000e+02, 1.330758e+02, 1.474062e+02, 1.631645e+02, 1.805499e+02,
           1.997900e+02, 2.211412e+02, 2.448906e+02, 2.713564e+02, 3.008875e+02,
           3.338628e+02, 3.706885e+02, 4.117938e+02, 4.576256e+02, 5.086399e+02,
           5.652923e+02, 6.280260e+02, 6.972586e+02, 7.733683e+02, 8.566789e+02,
           9.474479e+02, 1.045854e+03, 1.151991e+03, 1.265861e+03, 1.387377e+03,
           1.516364e+03, 1.652568e+03, 1.795671e+03, 1.945295e+03, 2.101027e+03,
           2.262422e+03, 2.429025e+03, 2.600380e+03, 2.776039e+03, 2.955570e+03,
           3.138565e+03, 3.324641e+03, 3.513446e+03, 3.704657e+03, 3.897982e+03,
           4.093159e+03, 4.289952e+03, 4.488155e+03, 4.687581e+03, 4.888070e+03,
           5.089479e+03, 5.291683e+03, 5.494575e+03, 5.698061e+03, 5.902058e+03])

  • lev_bounds
    (lev, axis_nbounds)
    float64
    dask.array<chunksize=(75, 2), meta=np.ndarray>

    units :

    m

    Array Chunk Bytes 1.17 kiB 1.17 kiB Shape (75, 2) (75, 2) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray

  • lon
    (y, x)
    float64
    dask.array<chunksize=(294, 362), meta=np.ndarray>

    bounds :

    bounds_lon

    long_name :

    Longitude

    standard_name :

    longitude

    units :

    degrees_east

    Array Chunk Bytes 831.47 kiB 831.47 kiB Shape (294, 362) (294, 362) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray

  • time
    (time)
    object
    1850-01-16 12:00:00 ... 2014-12-...

    axis :

    T

    bounds :

    time_bounds

    long_name :

    Time axis

    standard_name :

    time

    time_origin :

    1850-01-01 00:00:00

    array([cftime.DatetimeGregorian(1850, 1, 16, 12, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(1850, 2, 15, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(1850, 3, 16, 12, 0, 0, 0, has_year_zero=False),
           ...,
           cftime.DatetimeGregorian(2014, 10, 16, 12, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2014, 11, 16, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2014, 12, 16, 12, 0, 0, 0, has_year_zero=False)],
          dtype=object)

  • time_bounds
    (time, axis_nbounds)
    object
    dask.array<chunksize=(1980, 2), meta=np.ndarray>

    Array Chunk Bytes 30.94 kiB 30.94 kiB Shape (1980, 2) (1980, 2) Count 2 Tasks 1 Chunks Type object numpy.ndarray

  • bounds_lat_v
    (y_c, x, nvertex)
    float64
    dask.array<chunksize=(294, 362, 4), meta=np.ndarray>

    Array Chunk Bytes 3.25 MiB 3.25 MiB Shape (294, 362, 4) (294, 362, 4) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray

  • bounds_lon_v
    (y_c, x, nvertex)
    float64
    dask.array<chunksize=(294, 362, 4), meta=np.ndarray>

    Array Chunk Bytes 3.25 MiB 3.25 MiB Shape (294, 362, 4) (294, 362, 4) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray

  • lat_v
    (y_c, x)
    float64
    dask.array<chunksize=(294, 362), meta=np.ndarray>

    bounds :

    bounds_lat

    long_name :

    Latitude

    standard_name :

    latitude

    units :

    degrees_north

    Array Chunk Bytes 831.47 kiB 831.47 kiB Shape (294, 362) (294, 362) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray

  • lon_v
    (y_c, x)
    float64
    dask.array<chunksize=(294, 362), meta=np.ndarray>

    bounds :

    bounds_lon

    long_name :

    Longitude

    standard_name :

    longitude

    units :

    degrees_east

    Array Chunk Bytes 831.47 kiB 831.47 kiB Shape (294, 362) (294, 362) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray

  • areacello
    (y, x)
    float32
    dask.array<chunksize=(294, 362), meta=np.ndarray>

    cell_methods :

    area: sum

    description :

    Cell areas for any grid used to report ocean variables and variables which are requested as used on the model ocean grid (e.g. hfsso, which is a downward heat flux from the atmosphere interpolated onto the ocean grid). These cell areas should be defined to enable exact calculation of global integrals (e.g., of vertical fluxes of energy at the surface and top of the atmosphere).

    history :

    none

    long_name :

    Grid-Cell Area

    online_operation :

    once

    standard_name :

    cell_area

    units :

    m2

    Array Chunk Bytes 415.73 kiB 415.73 kiB Shape (294, 362) (294, 362) Count 5 Tasks 1 Chunks Type float32 numpy.ndarray

• Data variables: (3)
  • thetao
    (time, lev, y, x)
    float32
    dask.array<chunksize=(4, 75, 294, 362), meta=np.ndarray>

    cell_measures :

    area: areacello volume: volcello

    cell_methods :

    area: mean where sea time: mean

    description :

    Diagnostic should be contributed even for models using conservative temperature as prognostic field.

    history :

    none

    interval_operation :

    1800 s

    interval_write :

    1 month

    long_name :

    Sea Water Potential Temperature

    online_operation :

    average

    standard_name :

    sea_water_potential_temperature

    units :

    degC

    Array Chunk Bytes 58.88 GiB 121.80 MiB Shape (1980, 75, 294, 362) (4, 75, 294, 362) Count 496 Tasks 495 Chunks Type float32 numpy.ndarray

  • so
    (time, lev, y, x)
    float32
    dask.array<chunksize=(5, 75, 294, 362), meta=np.ndarray>

    cell_measures :

    area: areacello volume: volcello

    cell_methods :

    area: mean where sea time: mean

    description :

    Sea Water Salinity

    history :

    none

    interval_operation :

    1800 s

    interval_write :

    1 month

    long_name :

    Sea Water Salinity

    online_operation :

    average

    standard_name :

    sea_water_salinity

    units :

    0.001

    Array Chunk Bytes 58.88 GiB 152.25 MiB Shape (1980, 75, 294, 362) (5, 75, 294, 362) Count 397 Tasks 396 Chunks Type float32 numpy.ndarray

  • vo
    (time, lev, y_c, x)
    float32
    dask.array<chunksize=(3, 75, 294, 362), meta=np.ndarray>

    cell_methods :

    time: mean

    description :

    Prognostic y-ward velocity component resolved by the model.

    history :

    none

    interval_operation :

    1800 s

    interval_write :

    1 month

    long_name :

    Sea Water Y Velocity

    online_operation :

    average

    standard_name :

    sea_water_y_velocity

    units :

    m s-1

    Array Chunk Bytes 58.88 GiB 91.35 MiB Shape (1980, 75, 294, 362) (3, 75, 294, 362) Count 661 Tasks 660 Chunks Type float32 numpy.ndarray

• Attributes: (57)

  CMIP6_CV_version :

  cv=6.2.3.0-7-g2019642

  Conventions :

  CF-1.7 CMIP-6.2

  EXPID :

  CNRM-ESM2-1_historical_r1i1p1f2

  activity_id :

  CMIP

  arpege_minor_version :

  6.3.2

  branch_method :

  standard

  branch_time_in_child :

  0.0

  branch_time_in_parent :

  0.0

  contact :

  contact.cmip@meteo.fr

  creation_date :

  2018-06-10T22:12:36Z

  data_specs_version :

  01.00.21

  description :

  CMIP6 historical

  dr2xml_md5sum :

  789c0010d555326fb79e73c3e20c5d0e

  dr2xml_version :

  1.9

  experiment :

  all-forcing simulation of the recent past

  experiment_id :

  historical

  external_variables :

  areacello volcello

  forcing_index :

  2

  frequency :

  mon

  further_info_url :

  https://furtherinfo.es-doc.org/CMIP6.CNRM-CERFACS.CNRM-ESM2-1.historical.none.r1i1p1f2

  grid :

  native ocean tri-polar grid with 105 k ocean cells

  grid_label :

  gn

  history :

  none

  initialization_index :

  1

  institution :

  CNRM (Centre National de Recherches Meteorologiques, Toulouse 31057, France), CERFACS (Centre Europeen de Recherche et de Formation Avancee en Calcul Scientifique, Toulouse 31057, France)

  institution_id :

  CNRM-CERFACS

  license :

  CMIP6 model data produced by CNRM-CERFACS is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License (https://creativecommons.org/licenses). Consult https://pcmdi.llnl.gov/CMIP6/TermsOfUse for terms of use governing CMIP6 output, including citation requirements and proper acknowledgment. Further information about this data, including some limitations, can be found via the further_info_url (recorded as a global attribute in this file) and at http://www.umr-cnrm.fr/cmip6/. The data producers and data providers make no warranty, either express or implied, including, but not limited to, warranties of merchantability and fitness for a particular purpose. All liabilities arising from the supply of the information (including any liability arising in negligence) are excluded to the fullest extent permitted by law.

  mip_era :

  CMIP6

  nemo_gelato_commit :

  49095b3accd5d4c_6524fe19b00467a

  nominal_resolution :

  100 km

  parent_activity_id :

  CMIP

  parent_experiment_id :

  piControl

  parent_mip_era :

  CMIP6

  parent_source_id :

  CNRM-ESM2-1

  parent_time_units :

  days since 1850-01-01 00:00:00

  parent_variant_label :

  r1i1p1f2

  physics_index :

  1

  product :

  model-output

  realization_index :

  1

  realm :

  ocean

  references :

  http://www.umr-cnrm.fr/cmip6/references

  source :

  CNRM-ESM2-1 (2017): aerosol: TACTIC_v2 atmos: Arpege 6.3 (T127; Gaussian Reduced with 24572 grid points in total distributed over 128 latitude circles (with 256 grid points per latitude circle between 30degN and 30degS reducing to 20 grid points per latitude circle at 88.9degN and 88.9degS); 91 levels; top level 78.4 km) atmosChem: REPROBUS-C_v2 land: Surfex 8.0c ocean: Nemo 3.6 (eORCA1, tripolar primarily 1deg; 362 x 294 longitude/latitude; 75 levels; top grid cell 0-1 m) ocnBgchem: Pisces 2.s seaIce: Gelato 6.1

  source_id :

  CNRM-ESM2-1

  source_type :

  AOGCM BGC AER CHEM

  sub_experiment :

  none

  sub_experiment_id :

  none

  table_id :

  Omon

  title :

  CNRM-ESM2-1 model output prepared for CMIP6 / CMIP historical

  tracking_id :

  hdl:21.14100/9c34b796-c31d-4c1f-be90-21d032267f6e hdl:21.14100/e6ff74ce-1ace-4151-abf2-3e25386a8449 hdl:21.14100/3e3933f8-8edb-4296-928b-a84a5ba62ed3 hdl:21.14100/a55236c8-6a3e-4a77-97d4-48c7336aa825 hdl:21.14100/7031a63f-4d30-4672-8ab6-f516f604ee41 hdl:21.14100/72029dae-43de-4097-abae-b24a749387b2 hdl:21.14100/f8a5f02d-55c6-4c6a-a37f-9e01c3b68989

  variable_id :

  thetao

  variant_label :

  r1i1p1f2

  xios_commit :

  1442-shuffle

  status :

  2019-11-05;created;by nhn2@columbia.edu

  netcdf_tracking_ids :

  hdl:21.14100/9c34b796-c31d-4c1f-be90-21d032267f6e hdl:21.14100/e6ff74ce-1ace-4151-abf2-3e25386a8449 hdl:21.14100/3e3933f8-8edb-4296-928b-a84a5ba62ed3 hdl:21.14100/a55236c8-6a3e-4a77-97d4-48c7336aa825 hdl:21.14100/7031a63f-4d30-4672-8ab6-f516f604ee41 hdl:21.14100/72029dae-43de-4097-abae-b24a749387b2 hdl:21.14100/f8a5f02d-55c6-4c6a-a37f-9e01c3b68989

  version_id :

  v20181206

  intake_esm_varname :

  None

  intake_esm_dataset_key :

  CMIP.CNRM-CERFACS.CNRM-ESM2-1.historical.r1i1p1f2.Omon.thetao.gn.gs://cmip6/CMIP6/CMIP/CNRM-CERFACS/CNRM-ESM2-1/historical/r1i1p1f2/Omon/thetao/gn/v20181206/.nan.20181206

The grid is missing an outer coordinate for the Z axis, so we will construct one. This will be needed for conservative interpolation.

import cf_xarray
level_outer_data = cf_xarray.bounds_to_vertices(ds.lev_bounds, 'axis_nbounds').load().data

ds = ds.assign_coords({'level_outer': level_outer_data})

Linear Interpolation

Depth to Depth

To illustrate linear interpolation, we will first interpolate salinity onto a uniformly spaced vertical grid.

grid = Grid(ds, coords={'Z': {'center': 'lev'},
                        },
                periodic=False
            )
grid

<xgcm.Grid>
Z Axis (not periodic, boundary=None):
  * center   lev

target_depth_levels = np.arange(0,500,50)
salt_on_depth = grid.transform(ds.so, 'Z', target_depth_levels, target_data=None, method='linear')
salt_on_depth

<xarray.DataArray 'so' (time: 1980, y: 294, x: 362, lev: 10)>
dask.array<transpose, shape=(1980, 294, 362, 10), dtype=float32, chunksize=(5, 294, 362, 10), chunktype=numpy.ndarray>
Coordinates:
    lat        (y, x) float64 dask.array<chunksize=(294, 362), meta=np.ndarray>
    lon        (y, x) float64 dask.array<chunksize=(294, 362), meta=np.ndarray>
  * time       (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
    areacello  (y, x) float32 dask.array<chunksize=(294, 362), meta=np.ndarray>
  * lev        (lev) int64 0 50 100 150 200 250 300 350 400 450
Dimensions without coordinates: y, x

xarray.DataArray

'so'

• time: 1980
• y: 294
• x: 362
• lev: 10

• dask.array<chunksize=(5, 294, 362, 10), meta=np.ndarray>

  Array Chunk Bytes 7.85 GiB 20.30 MiB Shape (1980, 294, 362, 10) (5, 294, 362, 10) Count 1983 Tasks 396 Chunks Type float32 numpy.ndarray

• Coordinates: (5)
  • lat
    (y, x)
    float64
    dask.array<chunksize=(294, 362), meta=np.ndarray>

    Array Chunk Bytes 831.47 kiB 831.47 kiB Shape (294, 362) (294, 362) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray

  • lon
    (y, x)
    float64
    dask.array<chunksize=(294, 362), meta=np.ndarray>

    Array Chunk Bytes 831.47 kiB 831.47 kiB Shape (294, 362) (294, 362) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray

  • time
    (time)
    object
    1850-01-16 12:00:00 ... 2014-12-...

    array([cftime.DatetimeGregorian(1850, 1, 16, 12, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(1850, 2, 15, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(1850, 3, 16, 12, 0, 0, 0, has_year_zero=False),
           ...,
           cftime.DatetimeGregorian(2014, 10, 16, 12, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2014, 11, 16, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2014, 12, 16, 12, 0, 0, 0, has_year_zero=False)],
          dtype=object)

  • areacello
    (y, x)
    float32
    dask.array<chunksize=(294, 362), meta=np.ndarray>

    Array Chunk Bytes 415.73 kiB 415.73 kiB Shape (294, 362) (294, 362) Count 5 Tasks 1 Chunks Type float32 numpy.ndarray

  • lev
    (lev)
    int64
    0 50 100 150 ... 300 350 400 450

    array([  0,  50, 100, 150, 200, 250, 300, 350, 400, 450])

• Attributes: (0)

Note that the computation is lazy. (No data has been downloaded or computed yet.) We can trigger computation by plotting something.

salt_on_depth.isel(time=0).sel(lev=50).plot()

<matplotlib.collections.QuadMesh at 0x7fc8587b8550>

Depth to Potential Temperature

We can also interpolate salinity onto temperature surface through linear interpolation.

target_theta_levels = np.arange(-2, 36)
salt_on_theta = grid.transform(ds.so, 'Z', target_theta_levels, target_data=ds.thetao, method='linear')
salt_on_theta

<xarray.DataArray 'so' (time: 1980, y: 294, x: 362, thetao: 38)>
dask.array<transpose, shape=(1980, 294, 362, 38), dtype=float32, chunksize=(4, 294, 362, 38), chunktype=numpy.ndarray>
Coordinates:
    lat        (y, x) float64 dask.array<chunksize=(294, 362), meta=np.ndarray>
    lon        (y, x) float64 dask.array<chunksize=(294, 362), meta=np.ndarray>
  * time       (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
    areacello  (y, x) float32 dask.array<chunksize=(294, 362), meta=np.ndarray>
  * thetao     (thetao) int64 -2 -1 0 1 2 3 4 5 6 ... 27 28 29 30 31 32 33 34 35
Dimensions without coordinates: y, x

xarray.DataArray

'so'

• time: 1980
• y: 294
• x: 362
• thetao: 38

• dask.array<chunksize=(4, 294, 362, 38), meta=np.ndarray>

  Array Chunk Bytes 29.83 GiB 61.71 MiB Shape (1980, 294, 362, 38) (4, 294, 362, 38) Count 7131 Tasks 792 Chunks Type float32 numpy.ndarray

• Coordinates: (5)
  • lat
    (y, x)
    float64
    dask.array<chunksize=(294, 362), meta=np.ndarray>

    Array Chunk Bytes 831.47 kiB 831.47 kiB Shape (294, 362) (294, 362) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray

  • lon
    (y, x)
    float64
    dask.array<chunksize=(294, 362), meta=np.ndarray>

    Array Chunk Bytes 831.47 kiB 831.47 kiB Shape (294, 362) (294, 362) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray

  • time
    (time)
    object
    1850-01-16 12:00:00 ... 2014-12-...

    array([cftime.DatetimeGregorian(1850, 1, 16, 12, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(1850, 2, 15, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(1850, 3, 16, 12, 0, 0, 0, has_year_zero=False),
           ...,
           cftime.DatetimeGregorian(2014, 10, 16, 12, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2014, 11, 16, 0, 0, 0, 0, has_year_zero=False),
           cftime.DatetimeGregorian(2014, 12, 16, 12, 0, 0, 0, has_year_zero=False)],
          dtype=object)

  • areacello
    (y, x)
    float32
    dask.array<chunksize=(294, 362), meta=np.ndarray>

    Array Chunk Bytes 415.73 kiB 415.73 kiB Shape (294, 362) (294, 362) Count 5 Tasks 1 Chunks Type float32 numpy.ndarray

  • thetao
    (thetao)
    int64
    -2 -1 0 1 2 3 ... 30 31 32 33 34 35

    array([-2, -1,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15,
           16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
           34, 35])

• Attributes: (0)

salt_on_theta.isel(time=0).sel(thetao=20).plot()

/home/jthielen/miniconda3/envs/test_env_xgcm/lib/python3.10/site-packages/numba/np/ufunc/gufunc.py:170: RuntimeWarning: invalid value encountered in _interp_1d_linear
  return self.ufunc(*args, **kwargs)

<matplotlib.collections.QuadMesh at 0x7fc85a7b3760>

salt_on_theta.isel(time=0).mean(dim='x').plot(x='y')

/home/jthielen/miniconda3/envs/test_env_xgcm/lib/python3.10/site-packages/numba/np/ufunc/gufunc.py:170: RuntimeWarning: invalid value encountered in _interp_1d_linear
  return self.ufunc(*args, **kwargs)

<matplotlib.collections.QuadMesh at 0x7fc85a40d1e0>

Conservative Interpolation

To do conservative interpolation, we will attempt to calculate the meridional overturning in temperature space. Note that this is not a perfectly precise calculation. However, it’s sufficient to illustrate the basic principles of the calculation.

Create another grid object for conservative interpolation.

grid = Grid(ds, coords={'Z': {'center': 'lev', 'outer': 'level_outer'},
                        'X': {'center': 'x', 'right': 'x_c'},
                        'Y': {'center': 'y', 'right': 'y_c'}
                        },
            periodic=False,
            )
grid

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [27], in <cell line: 1>()
----> 1 grid = Grid(ds, coords={'Z': {'center': 'lev', 'outer': 'level_outer'},
      2                         'X': {'center': 'x', 'right': 'x_c'},
      3                         'Y': {'center': 'y', 'right': 'y_c'}
      4                         },
      5             periodic=False,
      6             )
      7 grid

File ~/develop/xgcm/xgcm/grid.py:1271, in Grid.__init__(self, ds, check_dims, periodic, default_shifts, face_connections, coords, metrics, boundary, fill_value)
   1269 for pos, dim in positions.items():
   1270     if not (dim in ds.variables or dim in ds.dims):
-> 1271         raise ValueError(
   1272             f"Could not find dimension `{dim}` (for the `{pos}` position on axis `{axis}`) in input dataset."
   1273         )
   1274     if dim not in ds.dims:
   1275         raise ValueError(
   1276             f"Input `{dim}` (for the `{pos}` position on axis `{axis}`) is not a dimension in the input datasets `ds`."
   1277         )

ValueError: Could not find dimension `x_c` (for the `right` position on axis `X`) in input dataset.

To use conservative interpolation, we have to go from an intensive quantity (velocity) to an extensive one (velocity times cell thickness). We fill any missing values with 0, since they don’t contribute to the transport.

thickness = grid.diff(ds.level_outer, 'Z')
v_transport =  ds.vo * thickness
v_transport = v_transport.fillna(0.).rename('v_transport')
v_transport

We also need to interpolate theta or thetao, our target data for interpolation, to the same horizontal position as v_transport. This means moving from cell center to cell corner. This step introduces some considerable errors, particularly near the boundaries of bathymetry. (Xgcm currently has no special treatment for internal boundary conditions–see issue 222.)

ds['theta'] = grid.interp(ds.thetao, ['Y'], boundary='extend')
ds.theta

We can transform v_transport to temperature space (target_theta_levels).

v_transport_theta = grid.transform(v_transport, 'Z', target_theta_levels,
                                   target_data=ds.theta, method='conservative')
v_transport_theta

Notice that this produced a warning. The conservative transformation method natively needs target_data to be provided on the cell bounds (here level_outer). Since transforming onto tracer coordinates is a very common scenario, xgcm uses linear interpolation to infer the values on the outer axis position.

To demonstrate how to provide target_data on the outer grid position, we reproduce the steps xgcm executes internally:

theta_outer = grid.interp(ds.theta,['Z'], boundary='extend')
# the data cannot be chunked along the transformation axis
theta_outer = theta_outer.chunk({'level_outer': -1}).rename('theta')
theta_outer

When we apply the transformation we can see that the results in this case are equivalent:

v_transport_theta_manual = grid.transform(v_transport, 'Z', target_theta_levels,
                                   target_data=theta_outer, method='conservative')

# Warning: this step takes a long time to compute. We will only compare the first time value
xr.testing.assert_allclose(v_transport_theta_manual.isel(time=0), v_transport_theta.isel(time=0))

Now we verify visually that the vertically integrated transport is conserved under this transformation.

v_transport.isel(time=0).sum(dim='lev').plot(robust=True)

v_transport_theta.isel(time=0).sum(dim='theta').plot(robust=True)

Finally, we attempt to plot a crude meridional overturning streamfunction for a single timestep.

dx = 110e3 * np.cos(np.deg2rad(ds.lat_v))
(v_transport_theta.isel(time=0) * dx).sum(dim='x').cumsum(dim='theta').plot.contourf(x='y_c', levels=31)

Logarithmic Interpolation

As noted previously, logarithmic interpolation is most often used to interpolate data from atmospheric models with a non-isobaric vertical coordinate (such as sigma or hybrid sigma) to isobaric levels suitable for analysis. And so, in place of the previous 4D ocean dataset, let’s generate a synthetic 3D atmospheric dataset (loosely based on the Sanders 1971 Analytic Model) to use to explore logarithmic interpolation:

def generate_analytic_model(
    nx=250,
    ny=150,
    nsigma=30,
    L=3e6,
    y_max=1.8e6,
    y_scale_factor=0.333,
    alpha=0.722,
    terrain_rise=1.5e3,
    pressure_variability=1e3,
    pressure_sea_level_mean=1e5,
    pressure_top=5e3,
    scale_height=8e3,
    temperature_variability=5,
    temperature_mean_surface=300,
    temperature_mean_tropopause=228.5,
    a=1.12e-5,
    b=2.3e-3,
    k=20
):
    """Generate sythetic data for an atmospheric trough over a Gaussian terrain.

    Parameters
    ----------
    nx : int
        Count of grid points in x direction
    ny : int
        Count of grid points in y direction
    nsigma : int
        Count of vertical sigma coordinate levels
    L : float
        Zonal wavelength of trough [meters]
    y_max : float
        Meridional extent of domain [meters]
    alpha : float
        Vertical control parameter for tropopause location [dimensionless]
    terrain_rise : float
        Max height of the terrain [meters]
    pressure_variability : float
        Wave amplitude of sea-level pressure field [Pa]
    pressure_sea_level_mean : float
        Domain average of sea-level pressure field [Pa]
    pressure_top : float
        Uniform pressure at top of domain [Pa]
    scale_height : float
        Control parameter for conversion of sea-level pressure to surface pressure [meters]
    temperature_variability : float
        Wave amplitude of temperature perturbation field [K]
    temperature_mean_surface : float
        Horizontal mean of temperature at ground level [K]
    temperature_mean_tropopause : float
        Horizontal mean of temperature at model top [K]
    a : float
        Meridional control parameter for temperature field shape [K / m]
    b : float
        Meridional control parameter for temperature field shape [Pa / m]
    k : float
        Vertical control parameter for mean temperature profile sharpness in vertical

    """
    # Constants
    R = 287.05
    g = 9.81

    # Define coordinates (and broadcasted versions)
    x = np.linspace(-L / 2, L / 2, nx)
    y = np.linspace(0, y_max, ny)
    sigma = np.linspace(1, 0, nsigma)

    x_2d, y_2d = np.meshgrid(x, y)
    x_3d, y_3d, sigma_3d = np.meshgrid(x, y, sigma)

    # Sea-level pressure (2D) as a trough shape
    pressure_sea_level = (
        pressure_sea_level_mean
        - b * y_2d * y_scale_factor
        - pressure_variability * np.cos(2 * np.pi / L * x_2d) * np.cos(2 * np.pi / L * y_2d * y_scale_factor)
    )

    # Terrain height as an offset Gaussian shape (uniform in meridional direction)
    geopotential_height_surface = np.exp(-((x_2d + L / 3) / (L / 3))**2) * terrain_rise

    # Surface pressure (2D) from sea-level pressure and terrain height
    pressure_surface = pressure_sea_level * np.exp(
        -geopotential_height_surface / scale_height
    )

    # Pressure (3D) from definition of sigma coordinate
    pressure = sigma_3d * (pressure_surface - pressure_top)[..., None] + pressure_top

    # Trough component of temperature
    temperature_pertubation = (
        -(1 + alpha * np.log(pressure_sea_level_mean / pressure))
        * (
            a * y_3d * y_scale_factor
            + temperature_variability * np.cos(2 * np.pi / L * x_3d) * np.cos(2 * np.pi / L * y_3d * y_scale_factor)
        )
    )

    # Vertical component of temperature
    temperature_mean = (
        temperature_mean_tropopause
        + (
            np.log(1 + np.exp(k * (sigma + alpha - 1)))
            / np.log(1 + np.exp(k * alpha))
        ) * (temperature_mean_surface - temperature_mean_tropopause)
    )

    # Combine and calcuate temperature and geopotential 3D
    temperature = temperature_mean[None, None] + temperature_pertubation
    geopotential = g * geopotential_height_surface[..., None] - np.concatenate(
        (
            np.zeros_like(x_2d)[..., None],
            np.cumsum(
                (
                    R
                    * (temperature[..., :-1] + temperature[..., 1:])
                    / (sigma_3d[..., :-1] + sigma_3d[..., 1:])
                    * np.diff(sigma_3d, axis=-1)
                ),
                axis=-1
            )
        ),
        axis=-1
    )

    # Return dataset!
    return xr.Dataset(
        {
            'temperature': (('y', 'x', 'sigma'), temperature),
            'geopotential': (('y', 'x', 'sigma'), geopotential),
            'pressure': (('y', 'x', 'sigma'), pressure),
        },
        {
            'y': y,
            'x': x,
            'sigma': sigma
        }
    ), {
        'p_init': pressure_sea_level,
        't_init': temperature_pertubation[..., 0]
    }

ds, stuff = generate_analytic_model()

ds

<xarray.Dataset>
Dimensions:       (y: 150, x: 250, sigma: 30)
Coordinates:
  * y             (y) float64 0.0 1.208e+04 2.416e+04 ... 1.788e+06 1.8e+06
  * x             (x) float64 -1.5e+06 -1.488e+06 ... 1.488e+06 1.5e+06
  * sigma         (sigma) float64 1.0 0.9655 0.931 ... 0.06897 0.03448 0.0
Data variables:
    temperature   (y, x, sigma) float64 305.5 302.2 298.9 ... 215.3 214.1 212.2
    geopotential  (y, x, sigma) float64 1.146e+04 1.452e+04 ... 3.418e+05
    pressure      (y, x, sigma) float64 8.728e+04 8.444e+04 ... 8.238e+03 5e+03

xarray.Dataset

• Dimensions:
  • y: 150
  • x: 250
  • sigma: 30
• Coordinates: (3)
  • y
    (y)
    float64
    0.0 1.208e+04 ... 1.788e+06 1.8e+06

    array([      0.      ,   12080.536913,   24161.073826,   36241.610738,
             48322.147651,   60402.684564,   72483.221477,   84563.758389,
             96644.295302,  108724.832215,  120805.369128,  132885.90604 ,
            144966.442953,  157046.979866,  169127.516779,  181208.053691,
            193288.590604,  205369.127517,  217449.66443 ,  229530.201342,
            241610.738255,  253691.275168,  265771.812081,  277852.348993,
            289932.885906,  302013.422819,  314093.959732,  326174.496644,
            338255.033557,  350335.57047 ,  362416.107383,  374496.644295,
            386577.181208,  398657.718121,  410738.255034,  422818.791946,
            434899.328859,  446979.865772,  459060.402685,  471140.939597,
            483221.47651 ,  495302.013423,  507382.550336,  519463.087248,
            531543.624161,  543624.161074,  555704.697987,  567785.234899,
            579865.771812,  591946.308725,  604026.845638,  616107.38255 ,
            628187.919463,  640268.456376,  652348.993289,  664429.530201,
            676510.067114,  688590.604027,  700671.14094 ,  712751.677852,
            724832.214765,  736912.751678,  748993.288591,  761073.825503,
            773154.362416,  785234.899329,  797315.436242,  809395.973154,
            821476.510067,  833557.04698 ,  845637.583893,  857718.120805,
            869798.657718,  881879.194631,  893959.731544,  906040.268456,
            918120.805369,  930201.342282,  942281.879195,  954362.416107,
            966442.95302 ,  978523.489933,  990604.026846, 1002684.563758,
           1014765.100671, 1026845.637584, 1038926.174497, 1051006.711409,
           1063087.248322, 1075167.785235, 1087248.322148, 1099328.85906 ,
           1111409.395973, 1123489.932886, 1135570.469799, 1147651.006711,
           1159731.543624, 1171812.080537, 1183892.61745 , 1195973.154362,
           1208053.691275, 1220134.228188, 1232214.765101, 1244295.302013,
           1256375.838926, 1268456.375839, 1280536.912752, 1292617.449664,
           1304697.986577, 1316778.52349 , 1328859.060403, 1340939.597315,
           1353020.134228, 1365100.671141, 1377181.208054, 1389261.744966,
           1401342.281879, 1413422.818792, 1425503.355705, 1437583.892617,
           1449664.42953 , 1461744.966443, 1473825.503356, 1485906.040268,
           1497986.577181, 1510067.114094, 1522147.651007, 1534228.187919,
           1546308.724832, 1558389.261745, 1570469.798658, 1582550.33557 ,
           1594630.872483, 1606711.409396, 1618791.946309, 1630872.483221,
           1642953.020134, 1655033.557047, 1667114.09396 , 1679194.630872,
           1691275.167785, 1703355.704698, 1715436.241611, 1727516.778523,
           1739597.315436, 1751677.852349, 1763758.389262, 1775838.926174,
           1787919.463087, 1800000.      ])

  • x
    (x)
    float64
    -1.5e+06 -1.488e+06 ... 1.5e+06

    array([-1500000.      , -1487951.807229, -1475903.614458, ...,  1475903.614458,
            1487951.807229,  1500000.      ])

  • sigma
    (sigma)
    float64
    1.0 0.9655 0.931 ... 0.03448 0.0

    array([1.      , 0.965517, 0.931034, 0.896552, 0.862069, 0.827586, 0.793103,
           0.758621, 0.724138, 0.689655, 0.655172, 0.62069 , 0.586207, 0.551724,
           0.517241, 0.482759, 0.448276, 0.413793, 0.37931 , 0.344828, 0.310345,
           0.275862, 0.241379, 0.206897, 0.172414, 0.137931, 0.103448, 0.068966,
           0.034483, 0.      ])

• Data variables: (3)
  • temperature
    (y, x, sigma)
    float64
    305.5 302.2 298.9 ... 214.1 212.2

    array([[[305.49123009, 302.19569012, 298.90423051, ..., 241.65188259,
             242.72995545, 244.33361423],
            [305.49580349, 302.20021263, 298.90869992, ..., 241.65128288,
             242.72786361, 244.32857962],
            [305.49681135, 302.20109377, 298.90944956, ..., 241.64231743,
             242.71670887, 244.31347899],
            ...,
            [304.95927571, 301.66458149, 298.37403211, ..., 241.32987889,
             242.50115771, 244.31347899],
            [304.96389752, 301.66931839, 298.37888806, ..., 241.34200166,
             242.51445753, 244.32857962],
            [304.96538598, 301.67084531, 298.38045477, ..., 241.346012  ,
             242.51886976, 244.33361423]],
    
           [[305.44188207, 302.14526217, 298.85268567, ..., 241.53376782,
             242.60169126, 244.19054582],
            [305.44639823, 302.14972754, 298.85709808, ..., 241.5331358 ,
             242.59957745, 244.18551138],
            [305.44734927, 302.15055199, 298.85779116, ..., 241.52413856,
             242.5884012 , 244.17041129],
    ...
            [294.88053902, 291.34362661, 287.80251302, ..., 215.55693741,
             214.30393625, 212.45393151],
            [294.88216839, 291.34529231, 287.80421631, ..., 215.56088955,
             214.30823432, 212.45873669],
            [294.88276485, 291.34590079, 287.80483721, ..., 215.56223885,
             214.30968935, 212.46033875]],
    
           [[294.25349055, 290.71563687, 287.17358684, ..., 215.03081378,
             213.85584025, 212.1913367 ],
            [294.24641956, 290.70856772, 287.16652006, ..., 215.02587554,
             213.85196477, 212.1897749 ],
            [294.23832217, 290.70044851, 287.15837877, ..., 215.01839033,
             213.84531158, 212.1850905 ],
            ...,
            [294.79421775, 291.25526371, 287.71203823, ..., 215.33869856,
             214.06547028, 212.1850905 ],
            [294.79581018, 291.25689155, 287.71370269, ..., 215.34255374,
             214.06966197, 212.1897749 ],
            [294.7963952 , 291.25748829, 287.71431154, ..., 215.34387123,
             214.07108189, 212.1913367 ]]])

  • geopotential
    (y, x, sigma)
    float64
    1.146e+04 1.452e+04 ... 3.418e+05

    array([[[1.14600535e+04, 1.45203435e+04, 1.76575387e+04, ...,
             2.04421342e+05, 2.50768611e+05, 3.90580209e+05],
            [1.15972780e+04, 1.46576138e+04, 1.77948559e+04, ...,
             2.04560259e+05, 2.50907271e+05, 3.90716823e+05],
            [1.17327389e+04, 1.47930842e+04, 1.79303348e+04, ...,
             2.04692692e+05, 2.51037778e+05, 3.90839794e+05],
            ...,
            [3.20252275e+01, 3.08696169e+03, 6.21861776e+03, ...,
             1.92645735e+05, 2.38940301e+05, 3.78680443e+05],
            [3.01660906e+01, 3.08514968e+03, 6.21685582e+03, ...,
             1.92650487e+05, 2.38947486e+05, 3.78695780e+05],
            [2.84066326e+01, 3.08340541e+03, 6.21512769e+03, ...,
             1.92650898e+05, 2.38948703e+05, 3.78699708e+05]],
    
           [[1.14600535e+04, 1.45198411e+04, 1.76565040e+04, ...,
             2.04355113e+05, 2.50678808e+05, 3.90412519e+05],
            [1.15972780e+04, 1.46571108e+04, 1.77938200e+04, ...,
             2.04493994e+05, 2.50817425e+05, 3.90549085e+05],
            [1.17327389e+04, 1.47925806e+04, 1.79292978e+04, ...,
             2.04626391e+05, 2.50947892e+05, 3.90672009e+05],
    ...
            [3.20252275e+01, 2.98422956e+03, 6.00684591e+03, ...,
             1.78484011e+05, 2.19614532e+05, 3.42115378e+05],
            [3.01660906e+01, 2.98238701e+03, 6.00502095e+03, ...,
             1.78484358e+05, 2.19615669e+05, 3.42119128e+05],
            [2.84066326e+01, 2.98063362e+03, 6.00327398e+03, ...,
             1.78483369e+05, 2.19614948e+05, 3.42119285e+05]],
    
           [[1.14600535e+04, 1.44059375e+04, 1.74219939e+04, ...,
             1.89445841e+05, 2.30483146e+05, 3.52779988e+05],
            [1.15972780e+04, 1.45430908e+04, 1.75590734e+04, ...,
             1.89578190e+05, 2.30614652e+05, 3.52909933e+05],
            [1.17327389e+04, 1.46784700e+04, 1.76943678e+04, ...,
             1.89707364e+05, 2.30742472e+05, 3.53034499e+05],
            ...,
            [3.20252275e+01, 2.98334985e+03, 6.00503284e+03, ...,
             1.78363682e+05, 2.19450505e+05, 3.41805728e+05],
            [3.01660906e+01, 2.98150693e+03, 6.00320710e+03, ...,
             1.78363977e+05, 2.19451569e+05, 3.41809341e+05],
            [2.84066326e+01, 2.97975343e+03, 6.00145988e+03, ...,
             1.78362971e+05, 2.19450825e+05, 3.41809452e+05]]])

  • pressure
    (y, x, sigma)
    float64
    8.728e+04 8.444e+04 ... 5e+03

    array([[[ 87277.73323609,  84440.57002105,  81603.40680602, ...,
              10674.32643008,   7837.16321504,   5000.        ],
            [ 87124.98439728,  84293.08838358,  81461.19236988, ...,
              10663.7920274 ,   7831.8960137 ,   5000.        ],
            [ 86973.90920732,  84147.22268293,  81320.53615854, ...,
              10653.37304878,   7826.68652439,   5000.        ],
            ...,
            [100957.52078754,  97648.64076038,  94339.76073322, ...,
              11617.76005431,   8308.88002716,   5000.        ],
            [100960.866917  ,  97651.87150607,  94342.87609514, ...,
              11617.99082186,   8308.99541093,   5000.        ],
            [100963.44864026,  97654.36420439,  94345.27976852, ...,
              11618.16887174,   8309.08443587,   5000.        ]],
    
           [[ 87269.70716152,  84432.82070768,  81595.93425383, ...,
              10673.77290769,   7836.88645385,   5000.        ],
            [ 87116.97235401,  84285.35261766,  81453.73288132, ...,
              10663.23947269,   7831.61973635,   5000.        ],
            [ 86965.91101054,  84139.50028604,  81313.08956154, ...,
              10652.821449  ,   7826.4107245 ,   5000.        ],
    ...
            [ 98908.06832214,  95669.85906965,  92431.64981717, ...,
              11476.41850498,   8238.20925249,   5000.        ],
            [ 98910.71514152,  95672.4146194 ,  92434.11409728, ...,
              11476.60104424,   8238.30052212,   5000.        ],
            [ 98913.03393082,  95674.65345044,  92436.27297007, ...,
              11476.76096075,   8238.38048037,   5000.        ]],
    
           [[ 85490.34737605,  82714.81815618,  79939.28893632, ...,
              10551.05843973,   7775.52921986,   5000.        ],
            [ 85340.91052972,  82570.53430455,  79800.15807939, ...,
              10540.75245033,   7770.37622516,   5000.        ],
            [ 85193.47928229,  82428.18689324,  79662.8945042 , ...,
              10530.58477809,   7765.29238904,   5000.        ],
            ...,
            [ 98890.83443521,  95653.21945469,  92415.60447416, ...,
              11475.22996105,   8237.61498052,   5000.        ],
            [ 98893.47321175,  95655.76723893,  92418.06126611, ...,
              11475.41194564,   8237.70597282,   5000.        ],
            [ 98895.78906903,  95658.00323907,  92420.2174091 , ...,
              11475.57165993,   8237.78582997,   5000.        ]]])

• Attributes: (0)

This synthetic model gives a low-pressure trough centered in the domain, with a sloping terrain in the zonal direction. Here is what that terrain (surface geopotential height) looks like:

(ds['geopotential'] / 9.81).rename('geopotential_height').sel(sigma=1, y=0).plot()

[<matplotlib.lines.Line2D at 0x7fc85a69d4b0>]

If we were to inspect a given vertical level of these data, it would be difficult to interpret due to the terrain-following nature of the sigma coordinate:

ds['geopotential'].isel(sigma=15).plot.contourf()

<matplotlib.contour.QuadContourSet at 0x7fc854a50df0>

And so, let’s interpolate to some common meteorological upper-air levels (750 hPa, 500 hPa, and 300 hPa) for plotting and analysis:

grid = Grid(
    ds,
    coords={
        'X': {'center': 'x'},
        'Y': {'center': 'y'},
        'Z': {'center': 'sigma'}
    },
    periodic=False
)

isobaric_levels = np.array([7.5e4, 5.0e4, 3.0e4])

geopotential_isobaric = grid.transform(
    ds['geopotential'],
    'Z',
    isobaric_levels,
    target_data=ds['pressure'],
    method='log'
)
temperature_isobaric = grid.transform(
    ds['temperature'],
    'Z',
    isobaric_levels,
    target_data=ds['pressure'],
    method='log'
)

ds_isobaric = xr.merge([geopotential_isobaric, temperature_isobaric])

ds_isobaric

<xarray.Dataset>
Dimensions:       (y: 150, x: 250, pressure: 3)
Coordinates:
  * y             (y) float64 0.0 1.208e+04 2.416e+04 ... 1.788e+06 1.8e+06
  * x             (x) float64 -1.5e+06 -1.488e+06 ... 1.488e+06 1.5e+06
  * pressure      (pressure) float64 7.5e+04 5e+04 3e+04
Data variables:
    geopotential  (y, x, pressure) float64 2.529e+04 6.027e+04 ... 9.459e+04
    temperature   (y, x, pressure) float64 291.3 262.7 242.8 ... 240.8 221.8

xarray.Dataset

• Dimensions:
  • y: 150
  • x: 250
  • pressure: 3
• Coordinates: (3)
  • y
    (y)
    float64
    0.0 1.208e+04 ... 1.788e+06 1.8e+06

    array([      0.      ,   12080.536913,   24161.073826,   36241.610738,
             48322.147651,   60402.684564,   72483.221477,   84563.758389,
             96644.295302,  108724.832215,  120805.369128,  132885.90604 ,
            144966.442953,  157046.979866,  169127.516779,  181208.053691,
            193288.590604,  205369.127517,  217449.66443 ,  229530.201342,
            241610.738255,  253691.275168,  265771.812081,  277852.348993,
            289932.885906,  302013.422819,  314093.959732,  326174.496644,
            338255.033557,  350335.57047 ,  362416.107383,  374496.644295,
            386577.181208,  398657.718121,  410738.255034,  422818.791946,
            434899.328859,  446979.865772,  459060.402685,  471140.939597,
            483221.47651 ,  495302.013423,  507382.550336,  519463.087248,
            531543.624161,  543624.161074,  555704.697987,  567785.234899,
            579865.771812,  591946.308725,  604026.845638,  616107.38255 ,
            628187.919463,  640268.456376,  652348.993289,  664429.530201,
            676510.067114,  688590.604027,  700671.14094 ,  712751.677852,
            724832.214765,  736912.751678,  748993.288591,  761073.825503,
            773154.362416,  785234.899329,  797315.436242,  809395.973154,
            821476.510067,  833557.04698 ,  845637.583893,  857718.120805,
            869798.657718,  881879.194631,  893959.731544,  906040.268456,
            918120.805369,  930201.342282,  942281.879195,  954362.416107,
            966442.95302 ,  978523.489933,  990604.026846, 1002684.563758,
           1014765.100671, 1026845.637584, 1038926.174497, 1051006.711409,
           1063087.248322, 1075167.785235, 1087248.322148, 1099328.85906 ,
           1111409.395973, 1123489.932886, 1135570.469799, 1147651.006711,
           1159731.543624, 1171812.080537, 1183892.61745 , 1195973.154362,
           1208053.691275, 1220134.228188, 1232214.765101, 1244295.302013,
           1256375.838926, 1268456.375839, 1280536.912752, 1292617.449664,
           1304697.986577, 1316778.52349 , 1328859.060403, 1340939.597315,
           1353020.134228, 1365100.671141, 1377181.208054, 1389261.744966,
           1401342.281879, 1413422.818792, 1425503.355705, 1437583.892617,
           1449664.42953 , 1461744.966443, 1473825.503356, 1485906.040268,
           1497986.577181, 1510067.114094, 1522147.651007, 1534228.187919,
           1546308.724832, 1558389.261745, 1570469.798658, 1582550.33557 ,
           1594630.872483, 1606711.409396, 1618791.946309, 1630872.483221,
           1642953.020134, 1655033.557047, 1667114.09396 , 1679194.630872,
           1691275.167785, 1703355.704698, 1715436.241611, 1727516.778523,
           1739597.315436, 1751677.852349, 1763758.389262, 1775838.926174,
           1787919.463087, 1800000.      ])

  • x
    (x)
    float64
    -1.5e+06 -1.488e+06 ... 1.5e+06

    array([-1500000.      , -1487951.807229, -1475903.614458, ...,  1475903.614458,
            1487951.807229,  1500000.      ])

  • pressure
    (pressure)
    float64
    7.5e+04 5e+04 3e+04

    array([75000., 50000., 30000.])

• Data variables: (2)
  • geopotential
    (y, x, pressure)
    float64
    2.529e+04 6.027e+04 ... 9.459e+04

    array([[[ 25287.83676273,  60269.40608511, 102606.3143656 ],
            [ 25270.14491781,  60267.55577153, 102614.86483933],
            [ 25251.91182001,  60264.46296755, 102621.15250807],
            ...,
            [ 26410.20506226,  60157.23776212, 101657.13023413],
            [ 26411.61399329,  60159.20026952, 101660.29089329],
            [ 26412.16747767,  60159.82125901, 101661.24049674]],
    
           [[ 25277.2894549 ,  60251.9328587 , 102576.53375009],
            [ 25259.62206473,  60250.10702578, 102585.10831535],
            [ 25241.41333033,  60247.03863595, 102591.42005169],
            ...,
            [ 26397.95006478,  60137.94806826, 101625.40547668],
            [ 26399.35848713,  60139.91002913, 101628.56555178],
            [ 26399.91162376,  60140.53065735, 101629.5147873 ]],
    
           [[ 25266.66299891,  60234.3260638 , 102546.52292979],
            [ 25249.02028812,  60232.52495387, 102555.12185708],
            [ 25230.83619202,  60229.48130557, 102561.45807817],
            ...,
    ...
            ...,
            [ 23760.11049494,  55954.73766707,  94708.12010479],
            [ 23760.5827102 ,  55955.260656  ,  94708.95794593],
            [ 23760.78318779,  55955.36076844,  94709.09338721]],
    
           [[ 22995.81352946,  56434.38393864,  96037.52032452],
            [ 22982.48001333,  56437.4772376 ,  96053.15699085],
            [ 22969.0674971 ,  56440.12072135,  96067.89902699],
            ...,
            [ 23738.0025745 ,  55919.45215389,  94649.49172875],
            [ 23738.46426402,  55919.95868659,  94650.30265163],
            [ 23738.66087487,  55920.05295276,  94650.42876817]],
    
           [[ 22976.85430342,  56402.45485677,  95982.66571965],
            [ 22963.57101275,  56405.5996659 ,  95998.35491335],
            [ 22950.21388009,  56408.30371983,  96013.16543872],
            ...,
            [ 23715.87188693,  55884.12669425,  94590.79241342],
            [ 23716.32302722,  55884.61673075,  94591.57634846],
            [ 23716.51576283,  55884.70513665,  94591.69311637]]])

  • temperature
    (y, x, pressure)
    float64
    291.3 262.7 242.8 ... 240.8 221.8

    array([[[291.27521451, 262.66950882, 242.77958028],
            [291.42881863, 262.76509415, 242.80774249],
            [291.57717925, 262.85502899, 242.82968173],
            ...,
            [279.2524916 , 255.043627  , 240.56558187],
            [279.25570078, 255.04922016, 240.57414443],
            [279.25565056, 255.0503969 , 240.57684069]],
    
           [[291.22875063, 262.60678176, 242.69666201],
            [291.38235673, 262.70236873, 242.72482504],
            [291.53071948, 262.79230539, 242.74676531],
            ...,
            [279.20495804, 254.98011489, 240.48203364],
            [279.20816695, 254.98570774, 240.49059587],
            [279.20811662, 254.98688438, 240.49329202]],
    
           [[291.18192136, 262.54356203, 242.61309297],
            [291.33552957, 262.63915081, 242.64125704],
            [291.48389483, 262.72908977, 242.66319897],
            ...,
    ...
            ...,
            [268.76060606, 241.02228804, 222.11796104],
            [268.76041146, 241.02334894, 222.12057768],
            [268.75922362, 241.02301235, 222.12129079]],
    
           [[280.95631955, 248.74066244, 224.31232276],
            [281.11739619, 248.84408999, 224.35483975],
            [281.27539884, 248.94481137, 224.3949731 ],
            ...,
            [268.67154604, 240.90319391, 221.96119698],
            [268.67131132, 240.90420127, 221.96374326],
            [268.67011007, 240.90384681, 221.96443291]],
    
           [[280.86971777, 248.62343285, 224.15767476],
            [281.03081114, 248.72688099, 224.20021661],
            [281.1888565 , 248.82765808, 224.24042113],
            ...,
            [268.58237005, 240.78394305, 221.80422498],
            [268.5820951 , 240.78489672, 221.80670071],
            [268.58088042, 240.78452434, 221.80736683]]])

• Attributes: (0)

An additional benefit of now having our data in isobaric coordinates is that the form of kinematics and dynamics formulas are more straightforward compared to non-isobaric forms. To demonstrate this, let’s plot the 300 hPa Temperature and Geostrophic Wind:

# Add inner coords for derivatives
ds_isobaric.coords['x_inner'] = (ds_isobaric['x'].values[:-1] + ds_isobaric['x'].values[1:]) / 2
ds_isobaric.coords['y_inner'] = (ds_isobaric['y'].values[:-1] + ds_isobaric['y'].values[1:]) / 2

# Create new grid
grid_isobaric = Grid(
    ds_isobaric,
    coords={
        'X': {'center': 'x', 'inner': 'x_inner'},
        'Y': {'center': 'y', 'inner': 'y_inner'},
        'Z': {'center': 'pressure'}
    },
    periodic=False
)

# Calculate geostrophic wind components (without metrics)
f = 1.2e-4  # f-plane approximation of Coriolis force
u_wind = grid_isobaric.interp(
    -grid_isobaric.diff(ds_isobaric['geopotential'], 'Y', to='inner') / grid_isobaric.diff(ds_isobaric['y'], 'Y', to='inner') / f,
    to='inner',
    axis='X'
)
v_wind = grid_isobaric.interp(
    grid_isobaric.diff(ds_isobaric['geopotential'], 'X', to='inner') / grid_isobaric.diff(ds_isobaric['x'], 'X', to='inner') / f,
    to='inner',
    axis='Y'
)

# Plot
level = 3.0e4
fig, ax = plt.subplots(1, 1, figsize=[8,5])
ds_isobaric['temperature'].sel(pressure=level).plot.contourf(ax=ax, alpha=0.3)
ds_isobaric['geopotential'].sel(pressure=level).plot.contour(ax=ax, colors='k')
barb_reduce = slice(10, -10, 25)
ax.barbs(
    u_wind.x_inner.values[barb_reduce],
    u_wind.y_inner.values[barb_reduce],
    u_wind.sel(pressure=level).values[barb_reduce, barb_reduce],
    v_wind.sel(pressure=level).values[barb_reduce, barb_reduce]
)
ax.set_aspect('equal')

plt.show();

Performance

By default xgcm performs some simple checks when using method='linear'. It checks if the last value of the data is larger than the first, and if not, the data is flipped.This ensures that monotonically decreasing variables, like temperature are interpolated correctly. These checks have a performance penalty (~30% in some preliminary tests).

If you have manually flipped your data and ensured that its monotonically increasing, you can switch the checks off to get even better performance.

grid.transform(..., method='linear', bypass_checks=True)