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 Axis and Grid objects 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.

  • 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.

[1]:
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.

[2]:
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
[2]:
<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.864 1.704 1.402 1.999 ... 1.163 0.956 0.6621 1.311
    theta    (z) float64 0.6931 1.099 1.386 1.609 ... 2.079 2.197 2.303 2.398

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.

[3]:
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();
_images/transform_5_0.svg

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.

[4]:
# 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
[4]:
<xarray.DataArray 'phi' (theta: 20)>
array([       nan,        nan,        nan,        nan,        nan,
       1.82634284, 1.76385252, 1.69701395, 1.531294  , 1.49504957,
       1.91749189, 1.91003094, 1.77717191, 1.27696796, 0.91885967,
       1.11001084,        nan,        nan,        nan,        nan])
Coordinates:
  * theta    (theta) float64 0.0 0.1579 0.3158 0.4737 ... 2.526 2.684 2.842 3.0

Now let’s see what the result looks like.

[5]:
ax = plot_profile()
phi_transformed.plot(ax=ax, y='theta', marker='.')
[5]:
[<matplotlib.lines.Line2D at 0x11a3ecac0>]
_images/transform_9_1.svg

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

[6]:
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='.')
[6]:
[<matplotlib.lines.Line2D at 0x11a4dad90>]
_images/transform_11_1.svg

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.

[7]:
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='.')
[7]:
[<matplotlib.lines.Line2D at 0x11a605670>]
_images/transform_13_1.svg

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.

[8]:
# 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
[8]:
<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.864 1.704 1.402 1.999 ... 1.163 0.956 0.6621 1.311
    theta    (z) float64 0.6931 1.099 1.386 1.609 ... 2.079 2.197 2.303 2.398
[9]:
# Recreate the grid object with a staggered `center`/`outer` coordinate layout
grid = Grid(ds, coords={'Z':{'center':'z', 'outer':'zc'}},
            periodic=False)
grid
[9]:
<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.

[10]:
ds['theta_outer'] = grid.interp(ds.theta, 'Z', boundary='fill')
ds['theta_outer']
[10]:
<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

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.

[11]:
# 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
/Users/ddeauna/Documents/GitHub/xgcm/xgcm/transform.py:239: 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(
[11]:
<xarray.DataArray 'phi' (theta_outer: 19)>
array([0.        , 0.        , 0.43144147, 0.53592955, 0.53592955,
       0.61431944, 0.77631926, 0.86126781, 1.0464762 , 1.40090257,
       1.78250446, 1.95974   , 1.96822435, 1.5915601 , 1.15854591,
       0.        , 0.        , 0.        , 0.        ])
Coordinates:
  * theta_outer  (theta_outer) float64 0.07895 0.2368 0.3947 ... 2.763 2.921
[12]:
phi_transformed_cons.plot(y='theta_outer', marker='.', yincrease=False)
[12]:
[<matplotlib.lines.Line2D at 0x11a688df0>]
_images/transform_21_1.svg

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.

[13]:
ds.phi.sum().values
[13]:
array(14.66316066)
[14]:
phi_transformed_cons.sum().values
[14]:
array(14.66316066)

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.

[15]:
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]
[16]:
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
[16]:
<xarray.Dataset>
Dimensions:       (axis_nbounds: 2, lev: 75, nvertex: 4, time: 1980, x: 362, y: 294, 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: axis_nbounds, nvertex, x, y, 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...

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

[17]:
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.

[18]:
grid = Grid(ds, coords={'Z': {'center': 'lev'},
                        },
                periodic=False
            )
grid
[18]:
<xgcm.Grid>
Z Axis (not periodic, boundary=None):
  * center   lev
[19]:
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
[19]:
<xarray.DataArray 'so' (time: 1980, y: 294, x: 362, lev: 10)>
dask.array<transpose, shape=(1980, 294, 362, 10), dtype=float64, 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

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

[20]:
salt_on_depth.isel(time=0).sel(lev=50).plot()
[20]:
<matplotlib.collections.QuadMesh at 0x14bf28850>
_images/transform_35_1.svg

Depth to Potential Temperature

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

[21]:
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
[21]:
<xarray.DataArray 'so' (time: 1980, y: 294, x: 362, thetao: 38)>
dask.array<transpose, shape=(1980, 294, 362, 38), dtype=float64, 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
[22]:
salt_on_theta.isel(time=0).sel(thetao=20).plot()
/Users/ddeauna/opt/miniconda3/envs/test_env_xgcm/lib/python3.9/site-packages/numba/np/ufunc/gufunc.py:151: RuntimeWarning: invalid value encountered in _interp_1d_linear
  return self.ufunc(*args, **kwargs)
[22]:
<matplotlib.collections.QuadMesh at 0x14c7109d0>
_images/transform_38_2.svg
[23]:
salt_on_theta.isel(time=0).mean(dim='x').plot(x='y')
/Users/ddeauna/opt/miniconda3/envs/test_env_xgcm/lib/python3.9/site-packages/numba/np/ufunc/gufunc.py:151: RuntimeWarning: invalid value encountered in _interp_1d_linear
  return self.ufunc(*args, **kwargs)
/Users/ddeauna/opt/miniconda3/envs/test_env_xgcm/lib/python3.9/site-packages/dask/array/numpy_compat.py:39: RuntimeWarning: invalid value encountered in true_divide
  x = np.divide(x1, x2, out)
[23]:
<matplotlib.collections.QuadMesh at 0x14cff2970>
_images/transform_39_2.svg

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.

[24]:
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
[24]:
<xgcm.Grid>
Z Axis (not periodic, boundary=None):
  * center   lev --> outer
  * outer    level_outer --> center
X Axis (not periodic, boundary=None):
  * center   x --> right
  * right    x_c --> center
Y Axis (not periodic, boundary=None):
  * center   y --> right
  * right    y_c --> center

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.

[25]:
thickness = grid.diff(ds.level_outer, 'Z')
v_transport =  ds.vo * thickness
v_transport = v_transport.fillna(0.).rename('v_transport')
v_transport
[25]:
<xarray.DataArray 'v_transport' (time: 1980, lev: 75, y_c: 294, x: 362)>
dask.array<where, shape=(1980, 75, 294, 362), dtype=float64, chunksize=(3, 75, 294, 362), chunktype=numpy.ndarray>
Coordinates:
  * lev      (lev) float64 0.5058 1.556 2.668 ... 5.495e+03 5.698e+03 5.902e+03
  * time     (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
    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>
Dimensions without coordinates: y_c, x

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

[26]:
ds['theta'] = grid.interp(ds.thetao, ['Y'], boundary='extend')
ds.theta
[26]:
<xarray.DataArray 'theta' (time: 1980, lev: 75, y_c: 294, x: 362)>
dask.array<mul, shape=(1980, 75, 294, 362), dtype=float32, chunksize=(4, 75, 293, 362), chunktype=numpy.ndarray>
Coordinates:
  * lev      (lev) float64 0.5058 1.556 2.668 ... 5.495e+03 5.698e+03 5.902e+03
  * time     (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
    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>
Dimensions without coordinates: y_c, x

We can transform v_transport to temperature space (target_theta_levels).

[27]:
v_transport_theta = grid.transform(v_transport, 'Z', target_theta_levels,
                                   target_data=ds.theta, method='conservative')
v_transport_theta
/Users/ddeauna/Documents/GitHub/xgcm/xgcm/grid.py:956: UserWarning: The `target data` input is not located on the cell bounds. This method will continue with linear interpolation with repeated boundary values. For most accurate results provide values on cell bounds.
  warnings.warn(
/Users/ddeauna/Documents/GitHub/xgcm/xgcm/transform.py:239: 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(
[27]:
<xarray.DataArray 'v_transport' (time: 1980, y_c: 294, x: 362, theta: 37)>
dask.array<transpose, shape=(1980, 294, 362, 37), dtype=float64, chunksize=(3, 293, 362, 37), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
    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>
  * theta    (theta) float64 -1.5 -0.5 0.5 1.5 2.5 ... 30.5 31.5 32.5 33.5 34.5
Dimensions without coordinates: y_c, x

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:

[28]:
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
[28]:
<xarray.DataArray 'theta' (time: 1980, level_outer: 76, y_c: 294, x: 362)>
dask.array<rechunk-merge, shape=(1980, 76, 294, 362), dtype=float32, chunksize=(4, 76, 293, 362), chunktype=numpy.ndarray>
Coordinates:
  * time         (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
  * level_outer  (level_outer) float64 0.0 1.024 2.103 ... 5.8e+03 6.004e+03
Dimensions without coordinates: y_c, x

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

[29]:
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))

/Users/ddeauna/Documents/GitHub/xgcm/xgcm/transform.py:239: 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(

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

[30]:
v_transport.isel(time=0).sum(dim='lev').plot(robust=True)
[30]:
<matplotlib.collections.QuadMesh at 0x14cce4b50>
_images/transform_54_1.svg
[31]:
v_transport_theta.isel(time=0).sum(dim='theta').plot(robust=True)
[31]:
<matplotlib.collections.QuadMesh at 0x15091cd00>
_images/transform_55_1.svg

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

[32]:
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)
[32]:
<matplotlib.contour.QuadContourSet at 0x14c9fa760>
_images/transform_57_1.svg

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)