如何在不进行三角测量、匹配pcolormesh的情况下绘制粗糙数据集的轮廓？

Question

我有一个非常粗略的数据集，用于使用 matplotlib 的 pcolormesh。 x,y 是 2D numpy 数组，代表当前的统一网格。 z 包含范围从 1-9 的整数，每个数字匹配一个相位。选项 shading='nearest' 将根据 z 选择的颜色置于 (x,y) 的中心。我的颜色图是根据可能的 z 值进行分段的。

vmin, vmax = 1, 9
colors = ['blue', 'orange', 'black', 'gray', 'cyan', 'lime', 'yellow', 'green', 'red']
cmap = ListedColormap(colors)
axes[0].pcolormesh(x, y, z, shading = 'nearest', vmin = vmin, vmax = vmax, cmap = cmap)

这样我就得到了我可以接受的上部子图。

但是，某些相具有共同的属性，这就是我想添加轮廓的原因。例如，我想绘制一个将彩色部分和黑色/灰色部分分开的轮廓。我这里有两个问题：

1. 如果我可以使用 contour< /a> 但我无法做到这一点，请参阅用 contourf，没有它对我的数据进行三角测量（？）。如果我有更多的数据点，这不会成为问题，但我不太可能大幅提高分辨率。即使我可以接受三角测量：不应绘制黄色区域。但由于 z 从 z=8（绿色）跳转到 z = 6（石灰），轮廓会插入一个中间黄色区域。
2. 取决于我们如何解决这个问题：我真的希望能够为连接和断开区域绘制轮廓。

我的一个想法是定义一个涵盖石灰、绿色和青色的新阶段，然后勾勒出该区域的轮廓。数据操作很简单，但是，我不知道之后如何继续使用 matplotlib。此外，我不知道如何识别连接和断开的单元格。

Answer 1

我设法完成了以下我几乎满意的设置：

关键是所谓的 <强>阿尔法形状。它本质上是三角测量，粗略地说是为了确定一组点的边界多边形。 这里是对应的python模块。实施起来非常简单。我之前没有任何关于 Shapely 的经验。此外，我还必须深入研究 matplotlib 的 pcolor源代码。最后我想出了以下脚本，主要代码在底部。

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import numpy.ma as ma
from descartes import PolygonPatch
import alphashape

#various parameters
nrows, ncols = 1, 1

# create segmented colormap with 9 colors (I will need 9 colors in my actual application)
colors = ['blue', 'orange', 'black', 'gray', 'cyan', 'lime', 'yellow', 'green', 'red']
cmap = mpl.colors.ListedColormap(colors)
# set vmin and vmax from 1 to 9
vmin = 1
vmax = 9
# name of the saved imagefile
savename = 'plot_alpha_shape'

# alpha value for the alpha shape, dont confuse it with the opacity from the
# standard mpl kwargs
alpha = 0.5
# grid discretization
dx = 2.0
dy = 1
# create 2d rectangular mesh
x, y = np.meshgrid(dx*np.arange(3), dy*np.arange(4))
# create homogeneous demo data
z = np.ones_like(x)*2
# change some values to make z heterogeneous
z[0] = np.ones(3)
z[1,-1] = 1
# define mask for the contour
mask = z>1
z_masked = ma.masked_array(z, mask = mask)

pcolor_kwargs = dict(shading = 'nearest', vmin = vmin, vmax = vmax, cmap = cmap)
contour_kwargs = dict(fc = 'none', ec = 'k', linewidth = 3)

def get_quadrilateral_vertices(x,y):
  X = interp_grid(x)
  Y = interp_grid(y)
  X = interp_grid(X.T).T
  Y = interp_grid(Y.T).T
  return X, Y

def interp_grid(X):
  dX = np.diff(X, axis=1)/2.
  X = np.hstack((X[:, [0]] - dX[:, [0]],
                  X[:, :-1] + dX,
                  X[:, [-1]] + dX[:, [-1]]))
  return X


def get_xymask(x,y):
  # merge x and y masks in case they are different
  mask = ma.getmaskarray(x) + ma.getmaskarray(y)
  # map mask to the cells in order to merge it with z mask
  xymask = (mask[0:-1, 0:-1] + mask[1:, 1:] +
            mask[0:-1, 1:] + mask[1:, 0:-1])
  return xymask

def execute_masking(x,y,z):
  # get dimensions
  Ny, Nx = x.shape
  xymask = get_xymask(x,y)
  # merge all masks
  # don't plot if C or any of the surrounding vertices are masked.
  mask = ma.getmaskarray(z) + xymask

  unmask = ~mask
  X1 = ma.filled(x[:-1, :-1])[unmask]
  Y1 = ma.filled(y[:-1, :-1])[unmask]
  X2 = ma.filled(x[1:, :-1])[unmask]
  Y2 = ma.filled(y[1:, :-1])[unmask]
  X3 = ma.filled(x[1:, 1:])[unmask]
  Y3 = ma.filled(y[1:, 1:])[unmask]
  X4 = ma.filled(x[:-1, 1:])[unmask]
  Y4 = ma.filled(y[:-1, 1:])[unmask]
  # npoly = len(X1)

  xy = np.stack([X1, Y1, X2, Y2, X3, Y3, X4, Y4], axis=-1)
  # one vertex is duplicate in the original code
  # xy = np.stack([X1, Y1, X2, Y2, X3, Y3, X4, Y4, X1, Y1], axis=-1)
  # transform to array of xy pairs
  verts = xy.reshape((-1, 2))
  z = ma.filled(z[:Ny - 1, :Nx - 1])[unmask]
  return verts, z

def get_masked_data(x,y,z):
  X, Y = get_quadrilateral_vertices(x,y)
  # convert to MA, if necessary.
  z = ma.asarray(z)
  X = ma.asarray(X)
  Y = ma.asarray(Y)
  return execute_masking(X,Y,z)

def plot_vertices(ax, verts):
  verts = verts.T
  ax.plot(*verts, linestyle = '', marker = 'x', color = 'r', ms = 10)


# main code
# x,y,z are cellcentered data
# use get_masked_data and its inner functions to get the vertices of the
# cells used in pcolor
# we are not using zdummy here
verts, zdummy = get_masked_data(x,y,z_masked)
# map vertices to a list of (x,y) tuples, each representing one vertex
contour_data = list(zip(verts[:,0], verts[:,1]))
# create an alpha shape from the vertices
contour_alphashape = alphashape.alphashape(verts, alpha)

# create figure with one subplot
fig = plt.figure(figsize=(15/2.54,nrows*4/2.54), constrained_layout=True)
gs = GridSpec(nrows, ncols, figure=fig)
axes = [fig.add_subplot(gs[j,i]) for j in range(nrows) for i in range(ncols)]

# plot vertices
plot_vertices(axes[0], verts)
# plot pcolor
pmesh1 = axes[0].pcolor(x,y,z, **pcolor_kwargs)
# plot the contour using alphashape
contour = PolygonPatch(contour_alphashape,**contour_kwargs)
axes[0].add_patch(contour)


# save the plot
plt.savefig(savename + '.png')

几乎所有定义的函数都取自matplotlib的 pcolor 和 _pcolorargs。

选择的 Alpha 越大，Alpha 形状就会变得越详细。对于非常小的 alpha，你会得到一个凸包。我附上脚本的结果。

，轮廓与蓝色区域不完全匹配。如果 alpha 变得太大，如果我理解正确的话，alpha 形状将不会返回正确的多边形，这就是为什么我无法使轮廓对齐得更紧密。我认为这也与我的数据的规则间距有关。