对数坐标系中的椭圆与多边形补丁

Question

我想做一些自动化的程序来生成阿什比图。这种图其实就是散点图，主要用于材料选择。通常会在材料的分组上加上一些包络线形状。因为数据的变化范围很大，所以在图表中常常使用对数坐标轴。

为了生成这个包络线，我做了两件事：

• 找出一组点的凸包；
• 找出一个椭球，表示每个数据集的两个标准差。

在使用线性坐标的图表时，这些方法都很好用，有时在对数坐标的图表中也能正常工作。这里有两个例子：

但在其他情况下，我的方法就不太管用了，也就是说，阴影区域的计算是正确的，但在对数坐标下的绘图却出错了（线性坐标下的绘图总是没问题）。这里是一些出现问题的例子：

问题：

• 在使用对数坐标轴时，matplotlib的Paths/Patches绘图有没有问题？

• 有没有办法绘制出我想要的图？

编辑

我按照@farenorth的建议添加了一个可运行的代码。我发现负值会造成问题。对于多边形（凸包）来说，这不是问题，因为材料属性通常是正值，但也有一些例外。对于这些情况，我认为使用symlog缩放就足够了。

至于椭圆，我需要更多考虑我的实现方式。我的方法是基于以算术平均值为中心，使用2*inc标准差作为半轴的椭圆。这可能会导致某些区域出现负值（因为数据不一定是围绕均值对称的）。数据文件可以在这里下载。

import numpy as np
from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt
from matplotlib.path import Path
import matplotlib.patches as patches
from matplotlib import rcParams

rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16


def poly_enclose(points, color, inc=1.2, rad=0.3, lw=2):
    """
    Plot the convex hull around a set of points as a 
    shaded polygon.
    """
    hull = ConvexHull(points)


    cent = np.mean(points, 0)
    pts = []
    for pt in points[hull.simplices]:
        pts.append(pt[0].tolist())
        pts.append(pt[1].tolist())

    pts.sort(key=lambda p: np.arctan2(p[1] - cent[1],
                                    p[0] - cent[0]))
    pts = pts[0::2]  # Deleting duplicates
    pts.insert(len(pts), pts[0])


    verts = inc*(np.array(pts)- cent) + cent
    verts2 = np.zeros((3*verts.shape[0]-2,2))
    verts2[0::3] = verts
    verts2[1::3,:] = (1-rad)*verts[0:-1,:] + rad*verts[1:,:]
    verts2[2::3,:] = rad*verts[0:-1,:] + (1-rad)*verts[1:,:]
    verts2[0:-1] = verts2[1:]
    verts2[-1] = verts2[0]

    codes = [Path.MOVETO, Path.LINETO, Path.CURVE3,]
    for j in range(len(pts)-2):
        codes.extend([Path.CURVE3, Path.LINETO, Path.CURVE3,])
    codes.append(Path.CURVE3)


    path = Path(verts2, codes)
    patch = patches.PathPatch(path, facecolor=color, lw=0, alpha=0.2)
    edge = patches.PathPatch(path, edgecolor=color, facecolor='none', lw=lw)
    plt.gca().add_patch(patch)
    plt.gca().add_patch(edge)


def ellip_enclose(points, color, inc=1.2, lw=2, nst=2):
    """
    Plot the minimum ellipse around a set of points.

    Based on: 
    https://github.com/joferkington/oost_paper_code/blob/master/error_ellipse.py
    """

    def eigsorted(cov):
        vals, vecs = np.linalg.eigh(cov)
        order = vals.argsort()[::-1]
        return vals[order], vecs[:,order]

    x = points[:,0]
    y = points[:,1]
    cov = np.cov(x, y)
    vals, vecs = eigsorted(cov)
    theta = np.degrees(np.arctan2(*vecs[:,0][::-1]))
    w, h = 2 * nst * np.sqrt(vals)        
    center = np.mean(points, 0)
    ell = patches.Ellipse(center, width=inc*w, height=inc*h, angle=theta,
                          facecolor=color, alpha=0.2, lw=0)
    edge = patches.Ellipse(center, width=inc*w, height=inc*h, angle=theta,
                          facecolor='none', edgecolor=color, lw=lw)
    plt.gca().add_artist(ell)
    plt.gca().add_artist(edge)


inc = 1.2
rad = 0.3
lw = 2
colors = ['blue', 'green', 'red', 'orange']

## Example 1
plt.figure()
for k in range(4):
    points = 1.5*(np.random.rand(20, 2) - 0.5) + k + 3
    plt.plot(points[:,0], points[:,1], 'o', ms=8, color=colors[k],
             mfc="white", mec=colors[k])
    poly_enclose(points, colors[k], inc=inc, rad=rad, lw=lw)
##    ellip_enclose(points, colors[k], inc=inc, lw=lw)

#plt.xscale('symlog')
#plt.yscale('symlog')
plt.grid(True)

##  Example 2
E = {}
E["poly"] = np.loadtxt('young_poly.txt')
E["metals"] = np.loadtxt('young_metals.txt')
E["comp"] = np.loadtxt('young_comp.txt')
E["ceramic"] = np.loadtxt('young_ceramic.txt')

rho = {}
rho["poly"] = np.loadtxt('dens_poly.txt')
rho["metals"] = np.loadtxt('dens_metals.txt')
rho["comp"] = np.loadtxt('dens_comp.txt')
rho["ceramic"] = np.loadtxt('dens_ceramic.txt')

plt.figure()
for k, key  in enumerate(E.keys()):
    x = rho[key][:,0] * 1000
    y = E[key][:,0] * 1e9
    points = np.vstack([x,y]).T
    poly_enclose(points, colors[k], inc=inc, rad=0.3, lw=lw)
##    ellip_enclose(points, colors[k], inc=1, lw=lw)
    plt.plot(x, y, 'o', ms=8, color=colors[k], mfc="white", mec=colors[k])

##plt.xscale('symlog')
##plt.yscale('symlog')
plt.grid(True)

plt.show()

Answer 1

这样导入ConvexHull，应该在大多数情况下都能正常工作：

from scipy.spatial.qhull import ConvexHull

Answer 2

在你的代码中，有两个地方假设了数据是线性的：一个是用线性平均值（numpy.mean）来计算点云的中心，另一个是将一个常数（inc）乘到每个点和计算出的中心之间的差向量上。

如果你想要直观的证明，可以试试在对数坐标系中绘制点云的平均中心，这个中心应该会偏离中心，往上和往右偏得很远。同样，如果你把点云乘以一个常数，然后把两个版本都画出来，你会发现并没有看到预期中的“变大”的线性扩展，而是更像是点云的点在原始云的右上角像下水道一样流走。

如果你在计算平均值和缩放差向量之前，把你的点转换成对数形式，然后在绘图之前再把点转换回线性形式，那么你会得到和线性图上看到的结果相同。这可以通过在你的函数中添加一个额外的标志来实现，或者让你的原始函数返回形状坐标，而不是直接绘图。我觉得后者看起来更简洁。

例如：

def poly_enclose(...):
    ....
    return patch, edge

log_points = np.log(points)
log_patch, log_edge = poly_enclose(log_points, colors[k], inc=inc, rad=0.3, lw=lw)
lin_patch = np.exp(log_patch)
lin_edge = np.exp(log_edge)
plt.gca().add_patch(lin_patch)
plt.gca().add_patch(lin_edge)