Generate missing C-grid coordinates¶

Climate model output is usually provided with geometric information like the position of the cell center and boundaries. Gridded observational datasets often lack these informations and provide only the position of either the gridcell center or cell boundary.

This makes the calculation of common vector calculus operators like the gradient difficult, and results depend on the particular method used.

xgcm can infer the cell boundary or cell center location depending on the geometry of the gridded source dataset. This enables consistent and easily reproducible calculations across model and observational datasets.

The autogenerate module enables the user to apply the same methods to both model output and observational products, which enables better comparison and a unified workflow using different sources of data.

In this case xgcm can infer the missing coordinates to enable the creation of a grid object. Below we will illustrate how to infer coordinates for several example datasets (nonperiodic and periodic) and show how the resulting dataset can be used to perform common calculations like gradients and distance/area weighted averages on observational data.

Non periodic 1D example: Ocean temperature profile¶

First let’s import xarray and xgcm

[1]:

import numpy as np
import xarray as xr
from xgcm import Grid
from xgcm.autogenerate import generate_grid_ds
import matplotlib.pyplot as plt
%matplotlib inline

Below we will create a synthetic temperature profile with decreasing temperature at depth, with unevenly spaced depth intervals (commonly found in hydrographic data).

[2]:

# create a synthetic ocean temperature profile with uneven spacing in depth
z = np.hstack([np.arange(1,10, 1), np.arange(10,200, 10), np.arange(200,700, 20)])
# Create synthetic temperature profile with maximum temperature gradient at mid depth (e.g. the thermocline)
temp = ((np.cos(np.pi*z/700) + 1) + np.exp(-z/350) / 2) * 10

# Convert to xarray.Dataset
ds = xr.DataArray(temp, dims=['depth'], coords={'depth':z}).to_dataset(name='temperature')

ds.temperature.plot(y='depth', yincrease=False)

[2]:

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

Infer the cell boundaries using xgcm.autogenerate¶

generate_grid_ds can infer the missing cell positions based on the given position (defaults to cell center) and the axis, which is defined by passing a dictionary with the physical axis as key and the dataset dimensions belonging to that axis as values.

[3]:

# Generate 'full' dataset, which includes additional coordinate `depth_left` and appropriate attributes.
ds_full = generate_grid_ds(ds, {'Z':'depth'})
print(ds_full)
print(ds.depth.data)
print(ds_full.depth_left.data)

<xarray.Dataset>
Dimensions:      (depth: 53, depth_left: 53)
Coordinates:
  * depth        (depth) int64 1 2 3 4 5 6 7 8 ... 560 580 600 620 640 660 680
  * depth_left   (depth_left) float64 0.5 1.5 2.5 3.5 ... 630.0 650.0 670.0
Data variables:
    temperature  (depth) float64 24.99 24.97 24.96 24.94 ... 1.164 0.9193 0.7567
[  1   2   3   4   5   6   7   8   9  10  20  30  40  50  60  70  80  90
 100 110 120 130 140 150 160 170 180 190 200 220 240 260 280 300 320 340
 360 380 400 420 440 460 480 500 520 540 560 580 600 620 640 660 680]
[5.00e-01 1.50e+00 2.50e+00 3.50e+00 4.50e+00 5.50e+00 6.50e+00 7.50e+00
 8.50e+00 9.50e+00 1.50e+01 2.50e+01 3.50e+01 4.50e+01 5.50e+01 6.50e+01
 7.50e+01 8.50e+01 9.50e+01 1.05e+02 1.15e+02 1.25e+02 1.35e+02 1.45e+02
 1.55e+02 1.65e+02 1.75e+02 1.85e+02 1.95e+02 2.10e+02 2.30e+02 2.50e+02
 2.70e+02 2.90e+02 3.10e+02 3.30e+02 3.50e+02 3.70e+02 3.90e+02 4.10e+02
 4.30e+02 4.50e+02 4.70e+02 4.90e+02 5.10e+02 5.30e+02 5.50e+02 5.70e+02
 5.90e+02 6.10e+02 6.30e+02 6.50e+02 6.70e+02]

We see now that a new dimension depth_left was created, with cell boundaries shifted towards the surface

The default behaviour of generate_grid_ds is to extrapolate the grid position to the ‘left’ (e.g. towards the surface for a depth profile), assuming that the spacing in the two cells closest to the boundary (here: the first and second cell) is equal. Particular geometries might require adjustments of the boundary treatment, by specifying e.g. pad=0 to ensure the topmost cell boundary is located at the sea surface.

Finally we can create the xgcm.Grid object like we would from model output (see for example here)

[4]:

# Create grid object
grid = Grid(ds_full, periodic=False)
print(grid)

<xgcm.Grid>
Z Axis (not periodic, boundary=None):
  * center   depth --> left
  * left     depth_left --> center

Calculate vertical gradient¶

Now we have all the tools we need to calculate the vertical gradient just like with numerical model output

[5]:

# Calculate vertical distances located on the cellboundary
ds_full.coords['dzc'] = grid.diff(ds_full.depth, 'Z', boundary='extrapolate')
# Calculate vertical distances located on the cellcenter
ds_full.coords['dzt'] = grid.diff(ds_full.depth_left, 'Z', boundary='extrapolate')

[6]:

# note that the temperature gradient is located on the `depth_left` dimension,
# but no temperature information is available, to infer the gradient in the topmost grid cell.
# The following will pad with nan towards the surface. Alternatively the values could be padded with
# with a particular value or linearly extrapolated.

dtemp_dz = grid.diff(ds.temperature, 'Z', boundary='fill', fill_value=np.nan) / ds_full.dzc
print(dtemp_dz)


fig, ax1 = plt.subplots()
ax1.invert_yaxis()
ax2 = ax1.twiny()

ds.temperature.plot(ax=ax1, y='depth', color='C0')
ax1.set_xlabel('temperature [deg C]', color='C0')

dtemp_dz.plot(ax=ax2, y='depth_left', color='C1')
ax2.set_xlabel('vertical temperature gradient [deg C / m]', color='C1');

<xarray.DataArray (depth_left: 53)>
array([        nan, -0.01452675, -0.01468758, -0.01484852, -0.01500955,
       -0.01517068, -0.01533191, -0.01549322, -0.01565462, -0.01581609,
       -0.01670564, -0.01832588, -0.01994686, -0.02156426, -0.02317378,
       -0.02477119, -0.02635229, -0.02791296, -0.02944913, -0.03095681,
       -0.03243211, -0.0338712 , -0.03527037, -0.03662601, -0.0379346 ,
       -0.03919276, -0.04039723, -0.04154487, -0.04263268, -0.04413764,
       -0.04591988, -0.04741622, -0.04860995, -0.04948706, -0.05003634,
       -0.05024945, -0.05012097, -0.04964842, -0.04883232, -0.04767612,
       -0.04618617, -0.04437168, -0.04224465, -0.0398197 , -0.03711403,
       -0.0341472 , -0.030941  , -0.02751928, -0.02390771, -0.02013362,
       -0.01622573, -0.01221393, -0.00812903])
Coordinates:
  * depth_left  (depth_left) float64 0.5 1.5 2.5 3.5 ... 610.0 630.0 650.0 670.0
    dzc         (depth_left) int64 1 1 1 1 1 1 1 1 1 ... 20 20 20 20 20 20 20 20

[6]:

Text(0.5, 0, 'vertical temperature gradient [deg C / m]')

Depth weighted average¶

Another common operation for many climate datasets is a weighted mean along an unevenly spaced dimension. Using the grid spacing for the tracer cells earlier this becomes trivial.

[7]:

mean_temp = ds_full.temperature.mean('depth')
mean_temp_weighted = (ds_full.temperature * ds_full.dzt).sum('depth') / ds_full.dzt.sum('depth')

print(mean_temp.data)
print(mean_temp_weighted.data)

16.271481181388435
12.320585088179214

Periodic 2D example¶

Below we will show how to apply these methods similarly to a global surface wind field, which is periodic in the longitudinal (‘x’) direction. The dataset is stored in a zenodo archive.

[8]:

# download the data
import urllib.request
import shutil

url = 'https://zenodo.org/record/4421428/files/'
file_name = 'uvwnd.10m.gauss.2018.nc'
with urllib.request.urlopen(url + file_name) as response, open(file_name, 'wb') as out_file:
    shutil.copyfileobj(response, out_file)

# open the data
ds = xr.open_dataset(file_name)
ds

[9]:

ds_full = generate_grid_ds(ds, {'X':'lon', 'Y':'lat'})
ds_full

As in the depth profile the longitude and latitude values are extended to the left when the defaults are used. However, since the latitude is not periodic we can specify the position to extend to as outer (more details here). This will extend the latitudinal positions both to the left and right, avoiding missing values later on.

[10]:

grid = Grid(ds_full, periodic=['X'])
grid

[10]:

<xgcm.Grid>
X Axis (periodic, boundary=None):
  * center   lon --> left
  * left     lon_left --> center
T Axis (not periodic, boundary=None):
  * center   time
Y Axis (not periodic, boundary=None):
  * center   lat --> left
  * left     lat_left --> center

Now we compute the difference (in degrees) along the longitude and latitude for both the cell center and the cell face. Since we are not taking the difference of a data variable across the periodic boundary, we need to specify the boundary_discontinutity in order to avoid the introduction of artefacts at the boundary.

[11]:

dlong = grid.diff(ds_full.lon, 'X', boundary_discontinuity=360)
dlonc = grid.diff(ds_full.lon_left, 'X', boundary_discontinuity=360)
dlonc_wo_discontinuity = grid.diff(ds_full.lon_left, 'X')

dlatg = grid.diff(ds_full.lat, 'Y', boundary='fill', fill_value=np.nan)
dlatc = grid.diff(ds_full.lat_left, 'Y', boundary='fill', fill_value=np.nan)

dlonc.plot()
dlonc_wo_discontinuity.plot(linestyle='--')

[11]:

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

../../../_images/repos_xgcm_xgcm-examples_05_autogenerate_24_1.png

Without adding the boundary_discontinuity input, the last cell distance is calculated incorectly.

The values we just calculated are actually not cell distances, but instead differences in longitude and latitude (in degrees). In order to calculate operators like the gradient dlon and dlat have to be converted into approximate cartesian distances on a globe.

[12]:

def dll_dist(dlon, dlat, lon, lat):
        """Converts lat/lon differentials into distances in meters

        PARAMETERS
        ----------
        dlon : xarray.DataArray longitude differentials
        dlat : xarray.DataArray latitude differentials
        lon  : xarray.DataArray longitude values
        lat  : xarray.DataArray latitude values

        RETURNS
        -------
        dx  : xarray.DataArray distance inferred from dlon
        dy  : xarray.DataArray distance inferred from dlat
        """

        distance_1deg_equator = 111000.0
        dx = dlon * xr.ufuncs.cos(xr.ufuncs.deg2rad(lat)) * distance_1deg_equator
        dy = ((lon * 0) + 1) * dlat * distance_1deg_equator
        return dx, dy

ds_full.coords['dxg'], ds_full.coords['dyg'] = dll_dist(dlong, dlatg, ds_full.lon, ds_full.lat)
ds_full.coords['dxc'], ds_full.coords['dyc'] = dll_dist(dlonc, dlatc, ds_full.lon, ds_full.lat)

Based on the distances we can estimate the area of each grid cell and compute the area-weighted meridional average temperature

[13]:

ds_full.coords['area_c'] = ds_full.dxc * ds_full.dyc

Area weighted average¶

By using the approximated cell area, we can easily compute the area weighted average of the zonal wind.

[16]:

u = ds_full.uwnd.mean('time')

u_mean = u.mean('lat')
u_mean_weighted = (u * u.area_c).sum('lat') / (u.area_c).sum('lat')


u_mean.plot(label='unweighted mean')
u_mean_weighted.plot(label='area weighted mean')
plt.legend()
plt.ylabel('zonal wind [m/s]')
plt.gca().autoscale('tight')

../../../_images/repos_xgcm_xgcm-examples_05_autogenerate_36_0.png

Thes two lines show the importance of applying a weighted mean, when the grid spacing is irregular (e.g. datasets gridded on a regular lat-lon grid).