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

[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>]


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>]


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>]


### 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>]


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

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>


#### 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>

[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>


### 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>

[31]:

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

[31]:

<matplotlib.collections.QuadMesh at 0x15091cd00>


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>


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