在Python中为所选网格区域的外边缘绘制轮廓的方法

问题描述

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-np.pi/2,np.pi/2,30)
y = np.linspace(-np.pi/2,30)
x,y = np.meshgrid(x,y)

z = np.sin(x**2+y**2)[:-1,:-1]

fig,ax = plt.subplots()
ax.pcolormesh(x,y,z)

哪个给出此图像：

现在可以说我要突出显示某些网格框的边缘：

highlight = (z > 0.9)

我可以使用轮廓函数，但这会导致轮廓“平滑”。我只想突出显示网格框边缘之后的区域边缘。

我最近来的就是添加这样的内容：

highlight = np.ma.masked_less(highlight,1)

ax.pcolormesh(x,highlight,facecolor = 'None',edgecolors = 'w')

给出以下情节：

哪一个很近，但是我真正想要的是仅突出显示“甜甜圈”的外边缘和内边缘。

所以从本质上讲，我正在寻找轮廓和pcolormesh函数的某种混合形式-遵循某些值的轮廓，但遵循“步骤”而不是点对点连接的网格箱。这有道理吗？

侧面说明：在pcolormesh参数中，我有edgecolors = 'w'，但边缘仍然显示为蓝色。那里发生了什么事？

编辑： JohanC使用add_iso_line（）的初始答案适用于提出的问题。但是，我正在使用的实际数据是非常不规则的x，y网格，无法转换为1D（add_iso_line()所必需。

我正在使用从极坐标（rho，phi）转换为笛卡尔（x，y）的数据。 JohanC提出的2D解决方案似乎不适用于以下情况：

import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage

def pol2cart(rho,phi):
    x = rho * np.cos(phi)
    y = rho * np.sin(phi)
    return(x,y)

phi = np.linspace(0,2*np.pi,30)
rho = np.linspace(0,2,30)

pp,rr = np.meshgrid(phi,rho)

xx,yy = pol2cart(rr,pp)

z = np.sin(xx**2 + yy**2)

scale = 5
zz = ndimage.zoom(z,scale,order=0)

fig,ax = plt.subplots()
ax.pcolormesh(xx,yy,z[:-1,:-1])

xlim = ax.get_xlim()
ylim = ax.get_ylim()
xmin,xmax = xx.min(),xx.max()
ymin,ymax = yy.min(),yy.max()
ax.contour(np.linspace(xmin,xmax,zz.shape[1]) + (xmax-xmin)/z.shape[1]/2,np.linspace(ymin,ymax,zz.shape[0]) + (ymax-ymin)/z.shape[0]/2,np.where(zz < 0.9,1),levels=[0.5],colors='red')
ax.set_xlim(*xlim)
ax.set_ylim(*ylim)

解决方法

This post显示了绘制此类线条的方法。由于适应当前的pcolormesh并不容易，因此以下代码演示了可能的适应方法。请注意，x和y的2d版本已重命名，因为线段需要1d版本。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection

x = np.linspace(-np.pi / 2,np.pi / 2,30)
y = np.linspace(-np.pi / 2,30)
xx,yy = np.meshgrid(x,y)

z = np.sin(xx ** 2 + yy ** 2)[:-1,:-1]

fig,ax = plt.subplots()
ax.pcolormesh(x,y,z)

def add_iso_line(ax,value,color):
    v = np.diff(z > value,axis=1)
    h = np.diff(z > value,axis=0)

    l = np.argwhere(v.T)
    vlines = np.array(list(zip(np.stack((x[l[:,0] + 1],y[l[:,1]])).T,np.stack((x[l[:,1] + 1])).T)))
    l = np.argwhere(h.T)
    hlines = np.array(list(zip(np.stack((x[l[:,0]],1] + 1])).T,1] + 1])).T)))
    lines = np.vstack((vlines,hlines))
    ax.add_collection(LineCollection(lines,lw=1,colors=color))

add_iso_line(ax,0.9,'r')
plt.show()

这里是第二个答案的改编，它只能用于2d数组：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from scipy import ndimage

x = np.linspace(-np.pi / 2,30)
x,y = np.meshgrid(x,y)

z = np.sin(x ** 2 + y ** 2)

scale = 5
zz = ndimage.zoom(z,scale,order=0)

fig,z[:-1,:-1] )
xlim = ax.get_xlim()
ylim = ax.get_ylim()
xmin,xmax = x.min(),x.max()
ymin,ymax = y.min(),y.max()
ax.contour(np.linspace(xmin,xmax,zz.shape[1]) + (xmax-xmin)/z.shape[1]/2,np.linspace(ymin,ymax,zz.shape[0]) + (ymax-ymin)/z.shape[0]/2,np.where(zz < 0.9,1),levels=[0.5],colors='red')
ax.set_xlim(*xlim)
ax.set_ylim(*ylim)
plt.show()

我将尝试重构add_iso_line方法，以使其更加清楚地进行优化。因此，首先，必须要做的部分：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection

x = np.linspace(-np.pi/2,np.pi/2,30)
y = np.linspace(-np.pi/2,y)
z = np.sin(x**2+y**2)[:-1,z)
xlim,ylim = ax.get_xlim(),ax.get_ylim()
highlight = (z > 0.9)

现在highlight是一个二进制数组，如下所示：之后，我们可以提取True单元的索引，查找False邻域并标识“红”线的位置。我对使用矢量化的方式（例如在add_iso_line中的方式）感到不满意，因此仅使用简单的循环即可：

lines = []
cells = zip(*np.where(highlight))
for x,y in cells:
    if x == 0 or highlight[x - 1,y] == 0: lines.append(([x,y],[x,y + 1]))
    if x == highlight.shape[0] or highlight[x + 1,y] == 0: lines.append(([x + 1,[x + 1,y + 1]))
    if y == 0 or highlight[x,y - 1] == 0: lines.append(([x,y]))
    if y == highlight.shape[1] or highlight[x,y + 1] == 0: lines.append(([x,y + 1],y + 1]))

最后，我调整线的大小并居中，以适合pcolormesh：

lines = (np.array(lines) / highlight.shape - [0.5,0.5]) * [xlim[1] - xlim[0],ylim[1] - ylim[0]]
ax.add_collection(LineCollection(lines,colors='r'))
plt.show()

总而言之，这与JohanC解决方案非常相似，并且通常较慢。幸运的是，我们可以大大减少cells的数量，仅使用python-opencv包提取轮廓：

import cv2
highlight = highlight.astype(np.uint8)
contours,hierarchy = cv2.findContours(highlight,cv2.RETR_TREE,cv2.CHAIN_APPROX_NONE)
cells = np.vstack(contours).squeeze()

这是正在检查的单元格的图示：

contour matplotlib python