python - 处理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和Y设置为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和X中出现的每个Y值，并在整个网格上重建函数。这是基于较弱的假设，即原始数据位于矩形网格上，但是不需要其他假设。当数据中缺少全波段时，这在您的特定情况下将无济于事，但是如果您的数据中包含nans补丁，它们将对您有帮助。因此，我正在提供这种情况的解决方案。为您

这是一个使用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)

如果没有nan的精简网格的尺寸小于原始尺寸，即数据中包含nan的完整行或列，则此解决方案将失效。但是，如果不是这种情况，那么它将为您提供如下所示的漂亮结果：



这也意味着，如果沿一条直线猜测原始问题中的X和Y值，例如，通过知道两个网格等距，则可以固定X的第一行和第一列Y，并使用上面的最新代码：它将为您生成完整的网格，其结果类似于本文中的第一幅图。