将 Z 维度中的一个数组与另一个数组相加或减去，其中纬度（Y 维度）和经度（X）不同

Question

最终，我想计算模拟空气质量和测量空气质量之间的差异。 我有两个 netcdf 文件。第一个 (A) 包含来自具有纬度 (y1) 索引长度 100 和经度 (x1) 索引长度 200 的模型的空气粒子数据。 从这个A中，我想减去纬度（y2）索引长度1300和经度（x2）索引长度1300的观测数据（B）。B的实际纬度值（以北和东为单位）存在于A中，虽然不完全正确，即 A 中的值是均匀分布的（例如 55.95°、55.85°、55.75° 等），但 B 中的值有 3 位小数，并且通过改变大约 0.001 到 0.003 的增量来分隔。

感觉这应该是直截了当的：获取纬度/经度范围内的 obs 数据（例如 50.5 到 51°N 和 8.1 到 8.2°E），然后从相同纬度/经度范围内的模型数据中减去它。

起初我尝试使用 numpy 改编自 this example 计算“与全球温度的偏差”。但我总是陷入死胡同。

然后，我尝试了无数种类似的东西（这显然是错误的，但我不是编码专家）：

anomaly=[]
for j in range(len(100)):
    for k in range(len(200)):
        for i in range(len(1300)):
            if latitude_model[j] == latitude_observation[i] and longitude_model[k] == longitude_observation[i]:
                departure = model_data[0,0,j,k] - observation_data[i,i] #the first two dimensions of the model data are 'time' and 'level'
                anomaly = np.append(departure)

我的第三种方法是使用 xarray 改编自 this example。 Xarray 将允许使用 method='nearest' 和 tolerance = 0.1 函数，这将有助于处理不匹配的纬度/经度数据（据我所知）。 但是在加载了两个 netcdf 文件后，我什至找不到如何继续代码的入口点。另外，我可能不得不重塑（但如何？）模型数据以匹配观察结果。或者从同一个模型网格中减去观测数据，如果多个观测点落在同一个网格内。

PS：这个问题最终和my other question有关，大概是同一个数据和问题。

Answer 1

在计算模型和观测值之间的偏差时，在进行分析之前匹配网格非常重要。换句话说，预处理始终是必须的。因此，在减去两个文件之前，您要么必须将模型的网格与观察到的网格相匹配，反之亦然，否则，由于存在差异，您的输出将毫无意义。最简单的方法是使用特殊运算符，如 CDO、NCO、NCL 等。

在您的命令行中（尽管 CDO 在 Python 中也可用，但所需的语法与下面的不同）

### Match the grids ###
cdo remapbil,obs.nc model.nc model1.nc

### Subtract the files ###
cdo sub model1.nc obs.nc bias.nc

然后，您可以轻松地在 Python 中映射差异。我更喜欢这种方法，因为它比在 Python 中预处理数据更容易、更轻量。

（从智能手机发送）

Answer 2

如果您想使用 Python（以 CDO 作为后端，需要安装）来执行此操作，您可以使用我的包 nctoolkit (https://nctoolkit.readthedocs.io/en/latest/installing.html)。

所以，如果你的两个文件分别命名为 file1 和 file2。您将首先将它们作为数据集读入。

将 nctoolkit 导入为 nc

data1 = nc.open_data(file1)

data2 = nc.open_data(file2)

然后，您可以将第一个数据集重新网格化，使其与第一个数据集具有相同的网格。这是必要的，以便单元格匹配。

data1.regrid(data2)

您可以从第一个数据集中减去第二个数据集。

data1.sub(data2)

如果您想将其转换为 xarray 对象，您可以这样做：

d1_xr = data1.to_xarray()

或者如果您想要一个 pandas 数据框，请执行以下操作：

d1_df = data1.to_dataframe()

还有一种自动绘图方法，使用全息视图：

df1_df.plot()