处理 pcolormesh 中的屏蔽坐标数组

Question

我正在努力可视化一些气候模型的输出。计算是在投影的纬度/经度网格上完成的。由于该模型模拟海冰，因此所有陆地网格单元都被屏蔽。在 Python 中绘制地理信息的标准工具是 Basemap 和 Cartopy，它们都使用 matplotlib 例程。特别是，pcolormesh这是绘图的明显选择。如果没有陆地遮罩，那就很简单了：

X = longitude
Y = latitude
C = variable

fig, ax = plt.subplots()
plt.pcolormesh(X,Y,C)

虽然C允许为掩码数组，但pcolormesh无法处理X和上的掩码数组Y。那么我该如何解决这个问题呢？

举一个简单的例子：

n = 100
X,Y = np.meshgrid(np.linspace(1,5,n),np.linspace(1,5,n))
C = np.sin(X*Y)
fig, ax = plt.subplots()
plt.pcolormesh(X,Y,C)

现在想象我们有一个面具：

X[50:60,:] = np.nan
X[:,50:60] = np.nan
Y[50:60,:] = np.nan
Y[:,50:60] = np.nan
C[50:60,:] = np.nan
C[:,50:60] = np.nan

我必须解决这个问题的第一个想法是只选择有效的条目并重塑X, Y, 和C：

M = np.isnan(X)
X_valid = X[~M]
Y_valid = Y[~M]
C_valid = C[~M]
X_valid.shape = (81,100)
Y_valid.shape = (81,100)
C_valid.shape = (81,100)
plt.pcolormesh(X_valid, Y_valid, C_valid)

像许多幼稚的方法一样，这是行不通的。

理想情况下，生成的图在掩模所在的位置处是空白的。如何才能做到这一点？

Answer 1

我发现您的“天真”方法有两个问题。

首先，通常不应设置坐标数组X和Yto nan，而应仅设置要绘制的函数的值。大多数绘图函数（包括matplotlib其他函数）会自动将这些函数视为缺失值，并用空白代替（将坐标设置为nan，另一方面，可能会干扰涉及插值等的内部例程）。

但是，这仍然不适用于pcolor(mesh). 但这没关系，因为我也不同意你的说法，即它“是阴谋的明显选择”。在我看来，pcolor(mesh)最适合绘制矩阵。对于像你这样的重要情节，类似的东西plt.contourf应该会产生奇迹。它还本质上包括插值，使你的情节更漂亮。它还nan按照我们的预期处理数据点：

n = 100
X, Y = np.meshgrid(np.linspace(1, 5, n), np.linspace(1, 5, n))
C = np.sin(X*Y)
C[50:60,:] = np.nan
C[:,50:60] = np.nan

fig, ax = plt.subplots()
n_levels = 100  # number of contour levels to plot
ax.contourf(X, Y, C, n_levels)

掩蔽之前（左）和之后（右）的结果：

请注意，它contourf代表“填充等高线图”，通过计算输入数据集的水平曲线来工作。这意味着为了获得平滑漂亮的绘图，你需要密集的轮廓线，这就是为什么我选择100条线来绘制。对于您的具体情况，您应该考虑使用关键字参数显式定义级别值levels。

更新：

在评论中，您澄清了您的数据集是给定的，因此您必须处理 和 中的缺失X值Y。这很困难，因为您的输入网格中有漏洞，只有当您非常精确地了解问题的情况时，您才能弥补这一点。

在您的示例中，沿每个维度的坐标都缺少完整区域。这是最好的场景，因为剩余的数据点可能是由meshgrid，仅沿每个维度具有较小的坐标向量。

在这个非常简单的情况下，一个简单的补救措施就是您自己尝试过的：丢弃这些nan值。你几乎是对的，但如果你采用一个形状数组(100,100)并从每个维度切出 10-10 个，你最终会得到一个形状数组(90,90)而不是(81,100)。这就是为什么你的身材看起来如此跳跃。如果你用正确的形状来做，结果会好得多：

n = 100
X, Y = np.meshgrid(np.linspace(1, 5, n), np.linspace(1, 5, n))
C = np.sin(X*Y)
X[50:60, :] = np.nan
X[:, 50:60] = np.nan
Y[50:60, :] = np.nan
Y[:, 50:60] = np.nan
C[50:60, :] = np.nan
C[:, 50:60] = np.nan

endshape = (90, 90)  # needs to be known a priori!

inds = np.logical_not(np.isnan(X) | np.isnan(Y) | np.isnan(C))
X_plot = np.reshape(X[inds], endshape)
Y_plot = np.reshape(Y[inds], endshape)
C_plot = np.reshape(C[inds], endshape)

fig, ax = plt.subplots()
n_levels = 100  # number of contour levels to plot
ax.contourf(X_plot, Y_plot, C_plot, n_levels)

结果显然与丢失的数据接近：执行的插值contourf（或者pcolormesh如果您使用插值）将尝试填补空白，从而扭曲您的数据。您可能会考虑在丢失的数据点上手动绘制白色补丁，但仍然会沿边缘出现一些失真。请注意，我们必须了解缺失点是如何分布的。

为了获得更简单和通用的解决方案，我会尝试猜测底层网格。我的意思是，您应该采用和unique中出现的每个值，并在整个网格上重建您的函数。这是基于一个较弱的假设，即原始数据位于矩形网格上，但不需要其他假设。当数据中缺少完整条带时，这在您的特定情况下没有帮助，但如果您的数据中有 nan 补丁，它们将会有所帮助。所以我也针对这个案例给出了一个解决方案。XY

这是一个用于scipy.interpolate.griddata重建网格的实现（使用插值可能有点矫枉过正，特别是因为我们丢弃了部分结果，但另一个选项是循环整个数据集，我不想这样做）：

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as interp

n = 100
X, Y = np.meshgrid(np.linspace(1, 5, n), np.linspace(1, 5, n))
C = np.sin(X*Y)

# poke a hole into the data
X[40:60, 40:60] = np.nan
Y[40:60, 40:60] = np.nan
C[40:60, 40:60] = np.nan

# indices where nobody is nan
inds = np.logical_not(np.isnan(X) | np.isnan(Y) | np.isnan(C))
X_notnan = X[inds]
Y_notnan = Y[inds]
C_notnan = C[inds]

# construct new mesh
X_vals = np.unique(X[inds])
Y_vals = np.unique(Y[inds])
X_plot, Y_plot = np.meshgrid(X_vals, Y_vals)

# use nearest-neighbour interpolation to match the two meshes
C_plot = interp.griddata(np.array([X_notnan, Y_notnan]).T, C_notnan,
                         (X_plot, Y_plot), method='nearest')
# fill in the nans in C
C_plot[np.logical_not(inds)] = np.nan

fig, ax = plt.subplots()
n_levels = 100  # number of contour levels to plot
ax.contourf(X_plot, Y_plot, C_plot, n_levels)

如果没有 s 的简化网格的nan尺寸小于原始网格的尺寸，即nan数据中存在 s 的完整行或列，则该解决方案将失效。但是，如果情况并非如此，那么它会给你一个漂亮的结果，如下所示：

这也意味着，如果您沿着一条线猜测原始问题中的X和Y值，例如知道两个网格是等距的，那么您可以修复 的第一行和X的第一列Y，并使用上面的最新代码：它应该为您生成完整的网格，产生类似于本文中第一张图的结果。