Python：如何绘制二维不连续的节点中心数据？

Question

我有一组二维数据和一个由四边形构成的二维网格，描述了一个被分割为补丁的区域。该数据在每个网格节点处定义。数据在补丁边界上存在不连续性，即相同位置处的数据被多次定义。 如何使用Python进行绘图，在节点之间进行线性插值，并以正确的方式呈现每个补丁面上的不连续值？ 以下是三个示例元素或补丁，每个元素都有六个节点值。

节点位置和值数据可以存储在 [Kx3x2] 的数组中，其中 K 是元素数量。例如：

x = np.array( [
[ [0.0, 1.0], [0.0, 1.0], [0.0, 1.0]  ],  #element 0
[ [1.0, 2.0], [1.0, 2.0], [1.0, 2.0]  ],  #element 1
[ [2.0, 3.0], [2.0, 3.0], [2.0, 3.0]  ],  #element 2
] )

y = np.array( [
[ [0.0, 0.0], [0.5, 0.5], [1.0, 1.0]  ],  #element 0
[ [0.0, 1.0], [0.5, 1.5], [1.0, 2.0]  ],  #element 1
[ [1.0, 1.0], [1.5, 1.5], [2.0, 2.0]  ],  #element 2
] )

z = np.array( [
[ [0.0, 0.5], [0.0, 0.8], [0.0, 1.0]  ],  #element 0
[ [0.3, 1.0], [0.6, 1.2], [0.8, 1.3]  ],  #element 1
[ [1.2, 1.5], [1.3, 1.4], [1.5, 1.7]  ],  #element 2
] )

我考虑了pyplot.imshow()。但它无法同时考虑整个域并仍能表示多值不连续节点。对于每个补丁，单独调用imshow()可能有效。但是，我该如何在同一坐标轴上绘制每个补丁图像？ 对于非矩形补丁，这也是imshow()存在问题的情况，而这是我所面临的普通情况。

我考虑了pyplot.pcolormesh()，但它似乎仅适用于以单元为中心的数据。

- NoahR

---

你能否尝试像 imshow 一样显示你的数据，将 alpha 通道设置为 0，以表示你的数据未定义的地方？ - ali_m

Answer 1

一个选项是通过所有元素的三角测量，并使用 matplotlib 的 tripcolor() 函数进行绘制。我现在已经发现了两个有用的演示，分别是这里和这里。

我的全局域的自动三角测量可能会出现问题，但单个四边形的 Delaunay 三角测量效果非常好： 我通过附加每个元素的三角测量来创建全局三角测量。这意味着共享节点实际上在位置数组和值数组中被复制。这允许在元素面上出现不连续的数据。 使用线性插值和所需的间断，可以使用 tripcolor() 函数进行绘制，提供每个节点的节点位置和值。 我有点担心等高线绘图可能会出现问题，因为元素面不再逻辑上连接。但是tricontour()仍然按预期工作。(用三角测量叠加显示在此处) 使用以下代码进行复制：

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

x = np.array( [
[ [0.0, 1.0], [0.0, 1.0], [0.0, 1.0]  ],  #element 0
[ [1.0, 2.0], [1.0, 2.0], [1.0, 2.0]  ],  #element 1
[ [2.0, 3.0], [2.0, 3.0], [2.0, 3.0]  ],  #element 2
] )

y = np.array( [
[ [0.0, 0.0], [0.5, 0.5], [1.0, 1.0]  ],  #element 0
[ [0.0, 1.0], [0.5, 1.5], [1.0, 2.0]  ],  #element 1
[ [1.0, 1.0], [1.5, 1.5], [2.0, 2.0]  ],  #element 2
] )

z = np.array( [
[ [0.0, 0.5], [0.0, 0.8], [0.0, 1.0]  ],  #element 0
[ [0.3, 1.0], [0.6, 1.2], [0.8, 1.3]  ],  #element 1
[ [1.2, 1.5], [1.3, 1.4], [1.5, 1.7]  ],  #element 2
] )



global_num_pts =  z.size
global_x = np.zeros( global_num_pts )
global_y = np.zeros( global_num_pts )
global_z = np.zeros( global_num_pts )
global_triang_list = list()

offset = 0;
num_triangles = 0;

#process triangulation element-by-element
for k in range(z.shape[0]):
    points_x = x[k,...].flatten()
    points_y = y[k,...].flatten()
    z_element = z[k,...].flatten()
    num_points_this_element = points_x.size

    #auto-generate Delauny triangulation for the element, which should be flawless due to quadrilateral element shape
    triang = tri.Triangulation(points_x, points_y)
    global_triang_list.append( triang.triangles + offset ) #offseting triangle indices by start index of this element

    #store results for this element in global triangulation arrays
    global_x[offset:(offset+num_points_this_element)] = points_x
    global_y[offset:(offset+num_points_this_element)] = points_y
    global_z[offset:(offset+num_points_this_element)] = z_element

    num_triangles += triang.triangles.shape[0]
    offset += num_points_this_element


#go back and turn all of the triangle indices into one global triangle array
offset = 0
global_triang = np.zeros( (num_triangles, 3) )
for t in global_triang_list:
    global_triang[ offset:(offset+t.shape[0] )] = t
    offset += t.shape[0]

plt.figure()
plt.gca().set_aspect('equal')

plt.tripcolor(global_x, global_y, global_triang, global_z, shading='gouraud' )
#plt.tricontour(global_x, global_y, global_triang, global_z )
#plt.triplot(global_x, global_y, global_triang, 'go-') #plot just the triangle mesh

plt.xlim((-0.25, 3.25))
plt.ylim((-0.25, 2.25))
plt.show()

- NoahR

---

显然，我必须等待两天才能允许这样做。 - NoahR

没关系，只是提醒一下。两天后它会在列表中向下移动，也许我就看不到它来提醒你了 ;) - tacaswell

我会说，沿着内部元素三角形面的渲染伪影让人不满意。相同两个节点之间的线性插值应该是相同的，因此不应该看到任何内部三角剖分的痕迹。 - NoahR