Matplotlib三角剖分插值等高面：曲线下方数据屏蔽问题

我想用Matplotlib的三角剖分生成插值等高面，但数据呈曲线形态，无法去除曲线下方的无效数据，希望将外围数据点设为边界。

参考Matplotlib官方教程编写了以下代码：

import matplotlib.tri as tri
import numpy as np
import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(nrows=2)
xi = np.linspace(-10,150,2000)
yi = np.linspace(-10,60,2000)
triang = tri.Triangulation(x_after, y_after)

interpolator = tri.LinearTriInterpolator(triang, strain_after)
Xi, Yi = np.meshgrid(xi, yi)
zi = interpolator(Xi, Yi)

ax1.triplot(triang, 'ro-', lw=5)
ax1.contour(xi, yi, zi, levels=30, linewidths=0.5, colors='k')
cntr1 = ax1.contourf(xi, yi, zi, levels=30, cmap="jet")
fig.colorbar(cntr1, ax=ax1)
ax1.plot(x_after, y_after, 'ko', ms=3)

ax2.tricontour(x_after, y_after, strain_after, levels=30, linewidths=0.5, colors='k')
cntr2 = ax2.tricontourf(x_after, y_after, strain_after, levels=30, cmap="jet")
fig.colorbar(cntr2, ax=ax2)
ax2.plot(x_after, y_after, 'ko', ms=3)

plt.subplots_adjust(hspace=0.5)
plt.show()

当前生成的图像效果如下：

我发现可以用triang.set_mask()方法掩码数据，但不知道如何定义掩码以实现需求。

以下是数据的具体数值：

x_after y_after z_after strain_after
39  117.2757    8.7586  0.1904  7.164
40  119.9474    7.152   0.1862  6.6456
37  111.8319    12.0568 0.1671  6.273
38  114.5314    10.4186 0.1651  5.7309
41  122.7482    5.4811  0.1617  9.1563
36  108.8823    13.4417 0.1421  8.8683
42  125.5035    3.8309  0.141   9.7385
33  99.8064 17.6315 0.1357  9.8613
32  96.8869 18.6449 0.1197  4.4147
35  105.8846    14.6086 0.1079  7.7055
28  84.2221 22.0191 0.1076  6.2098
26  77.8689 23.158  0.1067  7.5833
29  87.354  21.2974 0.1044  11.4365
27  81.0778 22.6443 0.1019  8.3794
24  71.4004 23.7749 0.0968  8.6207
34  102.8772    15.9558 0.0959  18.2025
23  68.2124 23.962  0.0939  7.9201
25  74.6905 23.4465 0.0901  9.0361
30  90.5282 20.398  0.0864  14.1051
31  93.802  19.335  0.0794  10.4563
43  128.3489    2.1002  0.0689  9.0292
22  65.0282 24.1107 0.0654  7.99
21  61.7857 24.0129 0.0543  8.2589
20  58.5831 23.9527 0.0407  9.0087
0   -0.0498 -0.5159 0.0308  7.1073
19  55.3115 23.7794 0.0251  9.6441
5   12.5674 9.3369  0.0203  7.2051
2   4.8147  3.6074  0.0191  8.0103
1   2.363   1.5329  0.0184  7.8285
18  52.0701 23.526  0.016   8.0149
3   7.4067  5.5988  0.0111  8.9994
7   18.2495 12.5836 0.0098  9.771
9   23.9992 15.4145 0.0098  6.7995
16  45.5954 22.5274 0.0098  12.9428
4   9.9776  7.5563  0.0093  6.9804
17  48.9177 23.0669 0.0084  9.3782
13  35.9812 20.0588 0.0066  9.6005
6   15.3578 11.0027 0.0062  9.7801
15  42.2909 21.8663 0.0052  12.0288
11  29.816  17.8723 0.0049  8.9085
8   21.1241 14.0893 0.0032  6.5716
10  26.8691 16.7093 0.0014  6.9672
44  131.1371    0.4155  0.0 11.9578
14  39.0687 20.991  -0.0008 9.9907
12  32.9645 18.9796 -0.0102 9.3389
45  134.083 -1.3928 -0.0616 15.29

解决方法 #

要实现将外围数据点设为边界、掩码掉无效区域的三角形，可按以下步骤操作：

1. 提取并排序外围边界点：从数据中挑选构成曲线外围的点，按曲线顺序排列（确保能形成闭合多边形边界）。
2. 创建多边形路径：用matplotlib.path.Path将边界点转化为可判断内外的路径。
3. 生成掩码：计算每个三角形的重心，判断重心是否在多边形内部，不在内部的三角形标记为需要掩码。
4. 设置掩码：将生成的掩码应用到三角剖分对象上。

修改后的代码如下：

import matplotlib.tri as tri
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.path import Path

# 假设x_after、y_after、strain_after是已加载的数据数组
fig, (ax1, ax2) = plt.subplots(nrows=2)
xi = np.linspace(-10,150,2000)
yi = np.linspace(-10,60,2000)
triang = tri.Triangulation(x_after, y_after)

# --- 添加掩码处理部分 ---
# 提取外围边界点的索引（按曲线顺序排列，确保形成闭合区域）
boundary_indices = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,36,37,38,39,40,41,42,43,44,45]
boundary_x = x_after[boundary_indices]
boundary_y = y_after[boundary_indices]

# 创建多边形路径
path = Path(np.column_stack((boundary_x, boundary_y)))

# 计算每个三角形的重心
tri_centers_x = np.array([triang.x[tri].mean() for tri in triang.triangles])
tri_centers_y = np.array([triang.y[tri].mean() for tri in triang.triangles])
centers = np.column_stack((tri_centers_x, tri_centers_y))

# 判断重心是否在多边形内部，不在的三角形需要掩码（~取反表示外部为True）
mask = ~path.contains_points(centers)

# 应用掩码
triang.set_mask(mask)
# --- 掩码处理结束 ---

interpolator = tri.LinearTriInterpolator(triang, strain_after)
Xi, Yi = np.meshgrid(xi, yi)
zi = interpolator(Xi, Yi)

ax1.triplot(triang, 'ro-', lw=5)
ax1.contour(xi, yi, zi, levels=30, linewidths=0.5, colors='k')
cntr1 = ax1.contourf(xi, yi, zi, levels=30, cmap="jet")
fig.colorbar(cntr1, ax=ax1)
ax1.plot(x_after, y_after, 'ko', ms=3)

ax2.tricontour(triang, strain_after, levels=30, linewidths=0.5, colors='k')
cntr2 = ax2.tricontourf(triang, strain_after, levels=30, cmap="jet")
fig.colorbar(cntr2, ax=ax2)
ax2.plot(x_after, y_after, 'ko', ms=3)

plt.subplots_adjust(hspace=0.5)
plt.show()

说明 #

• 边界点的索引需要根据数据实际形态调整，确保排列后形成的多边形能准确包围有效数据区域。
• 通过判断三角形重心是否在多边形内部，可以精准掩码掉曲线下方的无效三角形，只保留有效区域的插值结果。

内容的提问来源于stack exchange，提问作者Alex