绘制散点图最密集区域的等高线答案

【问题标题】：Plot contours for the densest region of a scatter plot绘制散点图最密集区域的等高线
【发布时间】：2013-10-19 04:03:09
【问题描述】：

我正在生成约 30 万个数据点的散点图，但我遇到的问题是在某些地方过于拥挤，以至于看不到任何结构 - 所以我有了一个想法！

我想让绘图为最密集的部分生成等高线图，并用 scatter() 数据点留下不太密集的区域。

所以我尝试单独计算每个数据点的最近邻距离，然后当该距离达到特定值时，绘制轮廓并填充它，然后当它达到更大的值（密度较小）只是做分散......

我已经尝试并失败了几天，我不确定传统的等高线图是否适用于这种情况。

我会提供代码，但它太乱了，可能只会混淆问题。而且它的计算量如此之大，如果它确实工作，它可能只会让我的电脑崩溃！

提前谢谢大家！

附言我一直在寻找和寻找答案！我相信它甚至不可能出现所有结果！

编辑：所以这样做的想法是查看某些特定点在 300k 样本的结构中的位置。这是一个示例图，我的观点分散在三个差异中。颜色。

我将尝试从我的数据中随机抽取 1000 个数据点并将其作为文本文件上传。干杯堆垛机。 :)

编辑：嘿，这是一些 1000 行的示例数据 - 只有两列 [X,Y]（或上图中的 [g-i,i]）以空格分隔。谢谢你们！ the data

【问题讨论】：

根据这些值的拥挤程度，您可以通过 scatter(x, y, alpha=0.1) 或一些合适的小值来梳理一些结构。按照你的建议，我会建立一个内核密度估计（参见scipy.stats.kde）。
为什么不使用二维直方图来显示数据？
@FriskyGrub 您可以只提供与您的真实数据具有相同类型/形状/等的随机数据 - 您并不总是需要首先发布生成真实数据的复杂步骤地方。让我们更容易提供对您有用的答案。
@RutgerKassies - 这并没有真正以有意义的方式显示该数据，并且受到分箱问题的影响。此外，很难在打印输出中正确表示它。
@FriskyGrub “这并没有真正以有意义的方式显示该数据，并且存在分箱问题。此外，很难在打印输出中正确表示它。” - 你是什么意思它没有意义？直方图是传达分布质量所在位置的一种完全有效的方式。通过在 300,000 个云中绘制每个 x,y 点的确切位置绝对没有任何好处，无论如何它们都相互重叠。找到在打印输出上看起来不错的颜色图并不难。

标签： python numpy matplotlib scipy contour

【解决方案1】：

您可以使用各种 numpy/scipy/matplotlib 工具来实现这一点：

创建一个原始点的scipy.spatial.KDTree 以便快速查找。
使用np.meshgrid 以您想要的轮廓的分辨率创建点网格
使用KDTree.query 创建目标密度内所有位置的掩码
使用矩形 bin 或 plt.hexbin 对数据进行分类。
从分箱数据中绘制轮廓，但使用第 3 步中的掩码过滤出密度较低的区域。
使用掩码的逆向plt.scatter 剩余点。

【讨论】：

我实际上并没有直接尝试过，但这基本上是我最终做的。我求助于使用 hexbin 的“热图”，因为我无法从命令 n^n -_- 中减少轮廓内容的计算时间……可能值得回头看看，这是一个有趣的问题。跨度>

【解决方案2】：

4年后，我终于可以回答这个问题了！这可以使用来自matplotlib.path 的contains_points 来完成。

我使用了来自astropy 的高斯平滑，可以根据需要省略或替换。

import matplotlib.colors as colors
from matplotlib import path
import numpy as np
from matplotlib import pyplot as plt
try:
    from astropy.convolution import Gaussian2DKernel, convolve
    astro_smooth = True
except ImportError as IE:
    astro_smooth = False

np.random.seed(123)
t = np.linspace(-1,1.2,2000)
x = (t**2)+(0.3*np.random.randn(2000))
y = (t**5)+(0.5*np.random.randn(2000))

H, xedges, yedges = np.histogram2d(x,y, bins=(50,40))
xmesh, ymesh = np.meshgrid(xedges[:-1], yedges[:-1])

# Smooth the contours (if astropy is installed)
if astro_smooth:
    kernel = Gaussian2DKernel(stddev=1.)
    H=convolve(H,kernel)

fig,ax = plt.subplots(1, figsize=(7,6)) 
clevels = ax.contour(xmesh,ymesh,H.T,lw=.9,cmap='winter')#,zorder=90)

# Identify points within contours
p = clevels.collections[0].get_paths()
inside = np.full_like(x,False,dtype=bool)
for level in p:
    inside |= level.contains_points(zip(*(x,y)))

ax.plot(x[~inside],y[~inside],'kx')
plt.show(block=False)

【讨论】：

【解决方案3】：

也许有人（比如我）会偶然发现互联网寻找答案。 @FriskyGrub，我喜欢你的平滑方法。 AstroML 库中有一个解决方案，例如 https://www.astroml.org/book_figures/chapter1/fig_S82_scatter_contour.html#book-fig-chapter1-fig-s82-scatter-contour 。我不确定您如何在代码中设置阈值（在其之上包含轮廓中的点而不是散布点），但我设法重现了与您类似的结果：

import matplotlib.pyplot as plt
from astroML.plotting import scatter_contour
np.random.seed(123)
t = np.linspace(-1,1.2,2000)
x = (t**2)+(0.3*np.random.randn(2000))
y = (t**5)+(0.5*np.random.randn(2000))
fig,ax = plt.subplots(1,1,figsize=(6,6))
scatter_contour(x,y, threshold=15, log_counts=True, ax=ax,
            histogram2d_args=dict(bins=15),
            plot_args=dict(marker='+', linestyle='none', color='black',
                          markersize=5),
            contour_args=dict(cmap='winter',),
           filled_contour=False)

( scatter_contour?? 在帮助下提出了很多文档，但基本上正如 kwargs 所建议的那样，histogram2d_args 是 numpy.histogram2d 采用的那些 args，plot_args 是 scatter plt.plot 和 @987654330 采用的那些 args @plt.contour（或plt.contourf）的那些人

最好的祝福

克里斯

【讨论】：