如何使用Matplotlib在极坐标系中绘制带有等高线密度线的散点图？

Question

我正在尝试制作一个极坐标系下的散点图，并将等高线叠加在点云上。我知道如何使用numpy.histogram2d在直角坐标系中完成这个操作：

# Simple case: scatter plot with density contours in cartesian coordinates

import matplotlib.pyplot as pl
import numpy as np

np.random.seed(2015)
N = 1000
shift_value = -6.

x1 = np.random.randn(N) + shift_value
y1 = np.random.randn(N) + shift_value

fig, ax = pl.subplots(nrows=1,ncols=1)

ax.scatter(x1,y1,color='hotpink')

H, xedges, yedges = np.histogram2d(x1,y1)
extent = [xedges[0],xedges[-1],yedges[0],yedges[-1]]
cset1 = ax.contour(H,extent=extent)

# Modify xlim and ylim to be a bit more consistent with what's next
ax.set_xlim(xmin=-10.,xmax=+10.)
ax.set_ylim(ymin=-10.,ymax=+10.)

输出在这里:

然而，当我尝试将我的代码转换为极坐标时，轮廓线变得扭曲了。下面是我的代码和生成的（错误）输出:

# Case with polar coordinates; the contour lines are distorted

np.random.seed(2015)
N = 1000
shift_value = -6.

def CartesianToPolar(x,y):
    r = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y,x)

    return theta, r

x2 = np.random.randn(N) + shift_value
y2 = np.random.randn(N) + shift_value

theta2, r2 = CartesianToPolar(x2,y2)

fig2 = pl.figure()
ax2 = pl.subplot(projection="polar")
ax2.scatter(theta2, r2, color='hotpink')

H, xedges, yedges = np.histogram2d(x2,y2)


theta_edges, r_edges = CartesianToPolar(xedges[:-1],yedges[:-1])
ax2.contour(theta_edges, r_edges,H)

这里是错误输出：

有没有办法让等高线在正确的比例尺下显示？

编辑以回应评论中提出的建议。

编辑2：有人建议这个问题可能是这个问题的重复。虽然我认识到这些问题很相似，但我的问题特别涉及将点的密度等高线绘制在散点图上。另一个问题是关于如何绘制与点坐标一起指定的任何量的等高线级别的。

- BrownianRatchet

---

嗨，保持两个图之间的移位值恒定可能有助于使问题更加清晰。 - camz

当然，你是对的，已经完成了！ - BrownianRatchet

可能是Polar contour plot in Matplotlib的重复问题。 - ali_m

另外：https://dev59.com/q2ox5IYBdhLWcg3wpFz5 - ali_m

Answer 1

问题在于你只转换了数组的边缘。通过仅转换边缘的x和y坐标，实际上相当于转换了跨越2D数组的对角线的坐标。该直线具有非常小的theta值范围，并且您将该范围应用于整个网格。

快速（但不正确）修复

在大多数情况下，您可以将整个网格（即2D数组的x和y，生成2D数组的theta和r）转换为极坐标。

而不是：

H, xedges, yedges = np.histogram2d(x2,y2)
theta_edges, r_edges = CartesianToPolar(xedges[:-1],yedges[:-1])

类似于这样做：

H, xedges, yedges = np.histogram2d(x2,y2)
xedges, yedges = np.meshgrid(xedges[:-1],yedges[:-1]
theta_edges, r_edges = CartesianToPolar(xedges, yedges)

作为一个完整的示例：

import numpy as np
import matplotlib.pyplot as plt

def main():
    x2, y2 = generate_data()
    theta2, r2 = cart2polar(x2,y2)

    fig2 = plt.figure()
    ax2 = fig2.add_subplot(111, projection="polar")
    ax2.scatter(theta2, r2, color='hotpink')

    H, xedges, yedges = np.histogram2d(x2,y2)

    xedges, yedges = np.meshgrid(xedges[:-1], yedges[:-1])
    theta_edges, r_edges = cart2polar(xedges, yedges)
    ax2.contour(theta_edges, r_edges, H)

    plt.show()

def generate_data():
    np.random.seed(2015)
    N = 1000
    shift_value = -6.

    x2 = np.random.randn(N) + shift_value
    y2 = np.random.randn(N) + shift_value
    return x2, y2

def cart2polar(x,y):
    r = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y,x)

    return theta, r

main()

然而，您可能会注意到这看起来有些不正确。 这是因为ax.contour隐含地假定输入数据位于常规网格上。 尽管我们已经提供了笛卡尔坐标系下的常规网格，但是在极坐标系下并没有常规网格。 实际上，它假设我们传递给它了一个在极坐标系下的常规网格。 我们可以重新采样网格，但有一种更简单的方法。

正确的解决方案

要正确绘制2D直方图，请在极坐标空间中计算直方图。

例如，可以执行类似以下操作：

theta2, r2 = cart2polar(x2,y2)
H, theta_edges, r_edges = np.histogram2d(theta2, r2)
ax2.contour(theta_edges[:-1], r_edges[:-1], H)

作为一个完整的示例：

import numpy as np
import matplotlib.pyplot as plt

def main():
    x2, y2 = generate_data()
    theta2, r2 = cart2polar(x2,y2)

    fig2 = plt.figure()
    ax2 = fig2.add_subplot(111, projection="polar")
    ax2.scatter(theta2, r2, color='hotpink')

    H, theta_edges, r_edges = np.histogram2d(theta2, r2)
    ax2.contour(theta_edges[:-1], r_edges[:-1], H)

    plt.show()

def generate_data():
    np.random.seed(2015)
    N = 1000
    shift_value = -6.

    x2 = np.random.randn(N) + shift_value
    y2 = np.random.randn(N) + shift_value
    return x2, y2

def cart2polar(x,y):
    r = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y,x)

    return theta, r

main()

最后，您可能会注意到上述结果略微偏移。这与基于单元格的网格约定（x [0,0]，y [0,0] 给出单元格的中心）与基于边缘的网格约定（x [0,0]，y [0,0] 给出单元格的左下角）。ax.contour 期望事物是以单元格为中心，但您提供的是对齐边缘的 x 和 y 值。

只有半个单元格的位移，但如果您想修复它，请执行以下操作：

def centers(bins):
    return np.vstack([bins[:-1], bins[1:]]).mean(axis=0)

H, theta_edges, r_edges = np.histogram2d(theta2, r2)
theta_centers, r_centers = centers(theta_edges), centers(r_edges)

ax2.contour(theta_centers, r_centers, H)

- Joe Kington

---

非常感谢，这解决了我的初始问题！然而，当我将其应用于实际数据（非高斯分布），密度轮廓线与点的位置略微偏移。你是否知道这个问题？ - BrownianRatchet