设计目标#

xgcm的诞生源于对 GCMs 数据分析的迫切需求。在诸如 CMIP6 (Coupled Model Intercomparison Project Phase 6) 这样的大型项目中，不同模式所提供的计算输出数据，无论是在水平方向还是垂直方向上，往往具有不同的网格结构。为了分析不同数值模式之间的计算结果差异，研究人员需要一个能够在模型原生网格上进行计算以及减小插值带来的误差的软件工具。

xgcm的设计目标是提供一个通用的数据处理分析工具，即使用户不了解底层网格的细节，也能像处理分析方程一样，轻松应用诸如梯度和散度等常见算子来进行数据处理。该工具能够自动识别速度变量所在的网格类型（例如 B 或 C 网格），并尽可能接近模型内部代码的方式执行数值计算。

xgcm基于xarray数据结构，这是一种坐标和元数据丰富的多维数组数据表示，非常适合于GCMs数据的分析。xarray提供了便捷的数据索引和分组、支持坐标的数据转换，以及基于dask的并行、核外数组计算能力。在这些功能基础上，xgcm增加了对有限体积Arakawa网格的支持，并提供了适用于此类网格的微分和积分算子。

xgcm的另一个主要动机源于海洋、大气和气候模型分辨率的快速提升。虽然现代超级计算机可以轻松生成TB甚至PB级的数据集，但常用的后处理工作流程在处理这些海量数据时显得力不从心。因此，一个灵活的、开源的、基于Python的GCM分析框架势在必行，xgcm旨在提升该领域的整体生产力，并加速海洋与大气领域的科学发现。

核心功能#

xgcm旨在提供一组核心的、通用的后处理功能，主要包括：

1. 网格拓扑 (Grid Topology)：xgcm能够处理常见的Arakawa网格类型（如 A、B、C 网格），并能自动识别变量所在的网格位置。这些网格通常用于海洋、大气模型中，其特点在于变量（如速度、压力、温度）在网格上的不同位置定义。理解网格拓扑是进行正确数值计算模拟的前提。

   

2. 网格度量 (Grid Metrics)：xgcm可以计算网格单元的距离、面积和体积等度量信息。这些度量信息对于正确计算诸如梯度、散度、通量等物理量至关重要。

3. 网格通用函数 (Grid Ufuncs)：xgcm提供了基于网格的通用函数，如差分和插值。这些函数能够自动处理不同网格类型的变量，从而简化了用户操作。xgcm 提供的 diff 和 interp 函数，用户无需考虑网格的细节，就可以实现变量在网格上的微分和插值运算。

4. 边界条件 (Boundary Conditions)：xgcm允许用户指定模型的边界条件，从而确保计算结果的准确性。

5. 垂直坐标转换 (Transforming Vertical Coordinates)：xgcm正在开发用于垂直网格转换的功能，例如将变量从几何坐标转换为密度坐标。

6. 矢量微积分算子：xgcm计划实现常见的矢量微积分算子，如梯度 (grad)、散度 (div) 和旋度 (curl)。这些算子的实现将进一步提高xgcm的功能，使得用户可以更加轻松地进行矢量场的分析。

xgcm的核心理念在于保持功能的简洁性和通用性。该工具避免添加过于专业化的功能，如Okubo-Weiss参数或剪切和应变场，而是鼓励用户在其他软件包中基于xgcm构建这些功能。这种方法使得xgcm专注于提供基础、通用的网格处理工具，同时保证了其灵活性和可扩展性。

应用案例#

1. MITgcm数值模式计算结果分析处理

1import xarray as xr
2import numpy as np
3import xgcm
4from matplotlib import pyplot as plt
5%matplotlib inline
6plt.rcParams[‘figure.figsize’] = (10,6)
7
8# download the data
9import urllib.request
10import shutil
11
12url = ‘https://zenodo.org/record/4421428/files/‘
13file_name = ‘mitgcm_example_dataset_v2.nc’
14with urllib.request.urlopen(url + file_name) as response, open(file_name, ‘wb’) as out_file:
15 shutil.copyfileobj(response, out_file)
16
17# open the data
18ds = xr.open_dataset(file_name)
19ds

1u_transport = ds.UVEL * ds.dyG * ds.hFacW * ds.drF
2v_transport = ds.VVEL * ds.dxG * ds.hFacS * ds.drF
3div_uv = (grid.diff(u_transport, ‘X’) + grid.diff(v_transport, ‘Y’)) / ds.rA
4div_uv[0, 0].plot()

1zeta = (-grid.diff(ds.UVEL * ds.dxC, ‘Y’) + grid.diff(ds.VVEL * ds.dyC, ‘X’))/ds.rAz
2zeta_bt = (zeta * ds.drF).sum(dim=‘Z’)
3zeta_bt.plot(vmax=2e-4)

1ke = 0.5*(grid.interp((ds.UVELds.hFacW)**2, ‘X’) + grid.interp((ds.VVELds.hFacS)**2, ‘Y’))
2ke[0,0].where(ds.maskC[0]).plot()

2. MITgcm ECCOv4r3数据分析处理

1from matplotlib import pyplot as plt
2import cartopy as cart
3import pyresample
4
5class LLCMapper:
6
7 def init(self, ds, dx=0.25, dy=0.25):
8 # Extract LLC 2D coordinates
9 lons_1d = ds.XC.values.ravel()
10 lats_1d = ds.YC.values.ravel()
11
12 # Define original grid
13 self.orig_grid = pyresample.geometry.SwathDefinition(lons=lons_1d, lats=lats_1d)
14
15 # Longitudes latitudes to which we will we interpolate
16 lon_tmp = np.arange(-180, 180, dx) + dx/2
17 lat_tmp = np.arange(-90, 90, dy) + dy/2
18
19 # Define the lat lon points of the two parts.
20 self.new_grid_lon, self.new_grid_lat = np.meshgrid(lon_tmp, lat_tmp)
21 self.new_grid = pyresample.geometry.GridDefinition(lons=self.new_grid_lon,
22 lats=self.new_grid_lat)
23
24 def call(self, da, ax=None, projection=cart.crs.Robinson(), lon_0=-60, **plt_kwargs):
25
26 assert set(da.dims) == set([‘face’, ‘j’, ‘i’]), “da must have dimensions [‘face’, ‘j’, ‘i’]“
27
28 if ax is None:
29 fig, ax = plt.subplots(figsize=(12, 6), subplot_kw={‘projection’: projection})
30 else:
31 m = plt.axes(projection=projection)
32
33 field = pyresample.kd_tree.resample_nearest(self.orig_grid, da.values,
34 self.new_grid,
35 radius_of_influence=100000,
36 fill_value=None)
37
38 vmax = plt_kwargs.pop(‘vmax’, field.max())
39 vmin = plt_kwargs.pop(‘vmin’, field.min())
40
41
42 x,y = self.new_grid_lon, self.new_grid_lat
43
44 # Find index where data is splitted for mapping
45 split_lon_idx = round(x.shape[1]/(360/(lon_0 if lon_0>0 else lon_0+360)))
46
47
48 p = ax.pcolormesh(x[:,], y[:,], field[:,],
49 vmax=vmax, vmin=vmin, transform=cart.crs.PlateCarree(), zorder=1, **plt_kwargs)
50 p = ax.pcolormesh(x[:,split_lon_idx:], y[:,split_lon_idx:], field[:,split_lon_idx:],
51 vmax=vmax, vmin=vmin, transform=cart.crs.PlateCarree(), zorder=2, **plt_kwargs)
52
53 ax.add_feature(cart.feature.LAND, facecolor=‘0.5’, zorder=3)
54 label = ”
55 if da.name is not None:
56 label = da.name
57 if ‘units’ in da.attrs:
58 label += ’ [%s]’ % da.attrs[‘units’]
59 cb = plt.colorbar(p, shrink=0.4, label=label)
60 return ax

1sst = ds.THETA.sel(time=‘2000-01-15’, k=0)
2mapper(sst, cmap=‘RdBu_r’)

未来开发计划

xgcm的未来开发路线计划主要集中在以下几个方面：

1. 矢量微积分算子: 实现梯度、散度、旋度等矢量微积分算子，并使其能像分析方程一样被轻松使用。
2. 自动网格度量处理: 进一步完善自动网格度量处理，提高用户操作的便捷性，并避免因网格不匹配导致的计算错误。
3. 垂直坐标转换: 实现通用的垂直坐标转换功能，从而使用户可以更灵活地处理和分析数据。
4. 更好的文档和示例: 提供更详尽的文档和示例，从而帮助用户更好地理解和使用 xgcm。
5. 与其他软件包的集成: 加强与其他 Python 科学计算软件包（如 scipy、numpy 等）的集成，从而进一步扩展 xgcm 的功能。

相关链接

• https://xgcm.readthedocs.io/en/latest/index.html