python - 停止 pyplot.contour 沿不连续点绘制轮廓

Question

我有一个坐标变换的二维地图。每个点的数据是原始坐标系中的轴角，从 0 到 360。我正在尝试使用 pyplot.contour 绘制恒定角度的线，例如 45 度。轮廓沿两个极点之间的 45 度线出现，但轮廓中有一个附加部分沿 0/360 不连续点连接两个极点。这会形成一条非常参差不齐的难看线，因为它基本上只是跟踪一侧数字接近 0 和另一侧接近 360 的像素。

示例：这是使用全彩色地图的图像：

您可以在左侧看到沿蓝色/红色曲线的不连续性。一侧是 360 度，另一侧是 0 度。绘制轮廓时，我得到：

请注意，所有轮廓都连接两个极点，但即使我没有绘制 0 度轮廓，所有其他轮廓都沿着 0 度不连续性（因为 pyplot 认为如果一侧为 0，另一侧为 360，则必须有之间的所有其他角度）。

生成此数据的代码：

import numpy as np
import matplotlib.pyplot as plt

jgal = np.array(
    [
        [-0.054875539726, -0.873437108010, -0.483834985808],
        [0.494109453312, -0.444829589425, 0.746982251810],
        [-0.867666135858, -0.198076386122, 0.455983795705],
    ]
)


def s2v3(rra, rdec, r):
    pos0 = r * np.cos(rra) * np.cos(rdec)
    pos1 = r * np.sin(rra) * np.cos(rdec)
    pos2 = r * np.sin(rdec)
    return np.array([pos0, pos1, pos2])


def v2s3(pos):
    x = pos[0]
    y = pos[1]
    z = pos[2]
    if np.isscalar(x):
        x, y, z = np.array([x]), np.array([y]), np.array([z])
    rra = np.arctan2(y, x)
    low = np.where(rra < 0.0)
    high = np.where(rra > 2.0 * np.pi)
    if len(low[0]):
        rra[low] = rra[low] + (2.0 * np.pi)
    if len(high[0]):
        rra[high] = rra[high] - (2.0 * np.pi)
    rxy = np.sqrt(x ** 2 + y ** 2)
    rdec = np.arctan2(z, rxy)
    r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
    if x.size == 1:
        rra = rra[0]
        rdec = rdec[0]
        r = r[0]
    return rra, rdec, r


def gal2fk5(gl, gb):
    rgl = np.deg2rad(gl)
    rgb = np.deg2rad(gb)
    r = 1.0
    pos = s2v3(rgl, rgb, r)

    pos1 = np.dot(pos.transpose(), jgal).transpose()

    rra, rdec, r = v2s3(pos1)

    dra = np.rad2deg(rra)
    ddec = np.rad2deg(rdec)

    return dra, ddec


def make_coords(resolution=50):
    width = 9
    height = 6
    px = width * resolution
    py = height * resolution
    coords = np.zeros((px, py, 4))
    for ix in range(0, px):
        for iy in range(0, py):
            l = 360.0 / px * ix - 180.0
            b = 180.0 / py * iy - 90.0
            dra, ddec = gal2fk5(l, b)
            coords[ix, iy, 0] = dra
            coords[ix, iy, 1] = ddec
            coords[ix, iy, 2] = l
            coords[ix, iy, 3] = b
    return coords


coords = make_coords()

# now do one of these
# plt.imshow(coords[:,:,0],origin='lower') # color plot
plt.contour(
    coords[:, :, 0], levels=[45, 90, 135, 180, 225, 270, 315]
)  # contour plot with jagged ugliness

plt.show()

我怎么能：

1. 停止 pyplot.contour 沿不连续点绘制轮廓

2. 让 pyplot.contour 认识到 0/360 角度的不连续性根本不是真正的不连续性。

我可以增加基础数据的分辨率，但在我得到一条漂亮的平滑线之前，它开始需要很长时间和大量的内存来绘制。

我还想沿着 0 度绘制轮廓，但如果我能弄清楚如何隐藏不连续性，我可以将它移到不在轮廓附近的其他地方。或者，如果我可以让#2 发生，那将不是问题。

Answer 1

这绝对仍然是一个 hack，但是您可以通过两种方法获得漂亮的平滑轮廓：

1. 绘制相位绝对值的等高线（从 -180° 到 180°），这样就没有不连续性。
2. 在有限区域中绘制两组轮廓，以便靠近极值顶部和底部的数值缺陷不会蔓延。

这是附加到您的示例的完整代码：

Z = np.exp(1j*np.pi*coords[:,:,0]/180.0)
Z *= np.exp(0.25j*np.pi/2.0)   # Shift to get same contours as in your example
X = np.arange(300)
Y = np.arange(450)

N = 2
levels = 90*(0.5 + (np.arange(N) + 0.5)/N)
c1 = plt.contour(X, Y, abs(np.angle(Z)*180/np.pi), levels=levels)
c2 = plt.contour(X, Y, abs(np.angle(Z*np.exp(0.5j*np.pi))*180/np.pi), levels=levels)

可以推广此代码以获得任何“周期性”函数的平滑轮廓。剩下要做的是生成一组具有正确值的新轮廓，以便正确应用颜色图，正确应用标签等。但是，使用 matplotlib 似乎没有一种简单的方法：相关QuadContourSet类做所有事情，我看不到从轮廓构造适当轮廓对象的简单方法c1和c2.

Answer 2

我对完全相同的问题感兴趣。一种解决方案是NaN沿分支切割的轮廓；看这里；另一个是使用matplotx 的 contour()max_jump中的参数。

我将解决方案塑造成一个 Python 包cplot。

import cplot
import numpy as np


def f(z):
    return np.exp(1 / z)


cplot.show(f, (-1.0, +1.0, 400), (-1.0, +1.0, 400))