交互式大图，含约2000万个样本点和数GB数据

Question

我遇到了一个问题（关于我的内存）：它无法保存我想要绘制的数据。不过，我的硬盘空间是足够的。有没有什么办法可以避免我的数据集被“遮蔽”？

具体来说，我在处理数字信号处理，需要使用高采样率。我的框架（GNU Radio）为了节省硬盘空间，会把数值以二进制格式保存。我解压后需要绘图。而且我希望这个图可以放大和互动，这就成了一个问题。

有没有什么优化的办法，或者其他软件/编程语言（比如R之类的）可以处理更大的数据集？其实我希望我的图能显示更多的数据。但我对其他软件没有经验。GNUplot在类似的情况下也失败了。我现在还不太懂R。

import matplotlib.pyplot as plt
import matplotlib.cbook as cbook
import struct

"""
plots a cfile

cfile - IEEE single-precision (4-byte) floats, IQ pairs, binary
txt - index,in-phase,quadrature in plaintext

note: directly plotting with numpy results into shadowed functions
"""

# unpacking the cfile dataset
def unpack_set(input_filename, output_filename):
    index = 0   # index of the samples
    output_filename = open(output_filename, 'wb')

    with open(input_filename, "rb") as f:

        byte = f.read(4)    # read 1. column of the vector

        while byte != "":
        # stored Bit Values
            floati = struct.unpack('f', byte)   # write value of 1. column to a variable
            byte = f.read(4)            # read 2. column of the vector
            floatq = struct.unpack('f', byte)   # write value of 2. column to a variable
            byte = f.read(4)            # next row of the vector and read 1. column
            # delimeter format for matplotlib 
            lines = ["%d," % index, format(floati), ",",  format(floatq), "\n"]
            output_filename.writelines(lines)
            index = index + 1
    output_filename.close
    return output_filename.name

# reformats output (precision configuration here)
def format(value):
    return "%.8f" % value            

# start
def main():

    # specify path
    unpacked_file = unpack_set("test01.cfile", "test01.txt")
    # pass file reference to matplotlib
    fname = str(unpacked_file)
    plt.plotfile(fname, cols=(0,1)) # index vs. in-phase

    # optional
    # plt.axes([0, 0.5, 0, 100000]) # for 100k samples
    plt.grid(True)
    plt.title("Signal-Diagram")
    plt.xlabel("Sample")
    plt.ylabel("In-Phase")

    plt.show();

if __name__ == "__main__":
    main()

像plt.swap_on_disk()这样的功能可以把数据缓存到我的SSD上；)

Answer 1

最近有一个新项目非常适合处理大数据集，那就是Bokeh，这个项目就是为了这个目的而创建的，正是为了处理大数据。

实际上，只有与图表规模相关的数据会被发送到显示的后端。这种方法比Matplotlib的方法要快得多。

Answer 2

#!/usr/bin/env python3

import numpy
import matplotlib.pyplot as plt

x, y, z = numpy.loadtxt('10m.csv', delimiter=',', unpack=True)
plt.figure(figsize=(8, 8), dpi=128)
plt.scatter(x, y, c=z)
# Non-interactive.
#plt.savefig('matplotlib.png')
# Interactive.
plt.show()

from bokeh.io import output_notebook, show
from bokeh.models import HoverTool
from bokeh.transform import linear_cmap
from bokeh.plotting import figure
from bokeh.models import ColumnDataSource
import numpy as np

output_notebook()
N = 1000000
source = ColumnDataSource(data=dict(
    x=np.random.random(size=N) * N,
    y=np.random.random(size=N) * N,
    z=np.random.random(size=N)
))
hover = HoverTool(tooltips=[("z", "@z")])
p = figure()
p.add_tools(hover)
p.circle(
    'x',
    'y',
    source=source,
    color=linear_cmap('z', 'Viridis256', 0, 1.0),
    size=5
)
show(p)

关于在Ubuntu上进行1000万点散点图基准测试的开源交互式绘图软件调查

这个调查是受到一个使用案例的启发，具体内容可以在这里找到：https://stats.stackexchange.com/questions/376361/how-to-find-the-sample-points-that-have-statistically-meaningful-large-outlier-r。我用完全相同的输入文件对几个绘图程序进行了基准测试。

简单来说，我想要做的是：

• 对多维数据进行XY散点图绘制，希望Z值能作为点的颜色
• 用鼠标交互式选择一些看起来有趣的点
• 查看所选点的所有维度（至少包括X、Y和Z），以尝试理解它们为什么在XY散点图中是异常值

这个问题可以用以下简化的测试数据来表示：

python -c 'for i in range(10000000): print(f"{i},{i*2},{i*4}")' > 10m1.csv
echo 5000000,20000000,-1 >> 10m1.csv

文件10m1.csv的前几行大致如下（约239 MB）：

10m1.csv

0,0,0
1,2,4
2,4,8
3,6,12
4,8,16

最后一行，即第10000000行，是异常值，内容如下：

5000000,20000000,-1

所以我们基本上有：

• 一条倾斜度为2的直线，上面有1000万个点
• 加上一个位于图表顶部中心的单个异常点

大致像这样：

Y

^
|
|
|       +       +
|
|             +
|
|           +
|
|         +
|
|       +
|
|     +
|
|   +
|
| +
|
+-------------------> X

这个基准测试的目标是找到图形中点(5000000,20000000)，然后确定它的第三列的值，在我们的测试中是-1。

当我第一次写这个回答时，我使用的是用以下命令生成的10.csv：

python -c 'for i in range(10000000): print(f"{i},{i*2},{i*4}")' > 10m.csv

没有异常值。虽然这测试了性能，但并没有测试选择能力，所以目标是当我找到动力时，将每个测试迁移到10m1.csv。

我还制作了一个包含10个点和一个异常值的示例，以便评估一些无法处理1000万个点的工具的可用性：

i=0;
while [ "$i" -lt 10 ]; do
  echo "$i,$((2 * i)),$((4 * i))"; i=$((i + 1));
done > 11.csv
echo 5,20,-1 >> 11.csv

为了增加乐趣，我还准备了一个更大的10亿点数据集，以防任何程序能够处理1000万个点！CSV文件有点不稳定，所以我转向了HDF5：

#!/usr/bin/env python3

import h5py
import numpy

size = 1000000000

with h5py.File('1b.hdf5', 'w') as f:
    x = numpy.arange(size + 1)
    x[size] =  size / 2
    f.create_dataset('x', data=x, dtype='int64')
    y = numpy.arange(size + 1) * 2
    y[size] =  3 * size / 2
    f.create_dataset('y', data=y, dtype='int64')
    z = numpy.arange(size + 1) * 4
    z[size] = -1
    f.create_dataset('z', data=z, dtype='int64')

这生成了一个约23 GiB的文件，类似于10m1.csv，包含：

• 1亿个点在一条直线上，类似于10m.csv
• 一个位于图表顶部中心的异常点

我还在创建10m1.csv的SQLite版本，因为这可能是实际操作中最合理的格式，因为它允许进行易于理解的SQL查询、明确的索引控制和二进制数值数据：

f=10m.sqlite
rm -f "$f"
n=10000000
time sqlite3 "$f" 'create table t(x integer, y integer, z integer)'
time sqlite3 "$f" "insert into t select value as id, value as x, value * 2 as y, value * 2 as z from generate_series(0, $((n - 1)))"
time sqlite3 "$f" "INSERT INTO t VALUES (?, ?, ?)', ($((n/2)), $((3*n/2)), -1))"
time sqlite3 "$f" 'create index txy on t(x, y)'

我还用n = 10亿运行了那段代码，生成了1b.sqlite。generate_series是我目前找到的最快插入方法：使用Python将大量数据批量插入SQLite

我按(x, y)进行索引，因为这可能会加速查看工具在给定x-y矩形中获取所有点的查询。生成的10m1.sqlite约为367 MB，由于索引的原因，比CSV大。

测试是在Ubuntu 18.10上进行的，除非在子部分中另有说明，使用的是ThinkPad P51笔记本电脑，配备Intel Core i7-7820HQ CPU（4核/8线程）、2个Samsung M471A2K43BB1-CRC RAM（2x 16GiB）和NVIDIA Quadro M1200 4GB GDDR5 GPU。

结果总结

这是我观察到的，考虑到我非常具体的测试用例，以及我是许多被评测软件的首次用户：

它能处理1000万个点吗：

工具	能处理1000万个点吗？	功能多吗？	用户界面感觉好吗？
Vaex	是的，甚至可以处理10亿！	是的。	是的，Jupyter小部件
VisIt	是的，但不能处理1亿	是的，支持2D和3D，专注于交互式。	否
Paraview	否	同上，可能2D功能稍少。	非常好
Mayavi	是的	仅支持3D，良好的交互和脚本支持，但功能较有限。	还可以
gnuplot	在非交互模式下勉强可以。	功能很多，但在交互模式下有限。	还可以
matplotlib	否	同上。	还可以
Bokeh	否，最多支持1百万	是的，易于脚本化。	非常好，Jupyter小部件
PyViz	？	？	？
seaborn	？	？	？
sqlitebrowser	否	可以可视化SQL查询结果	一般

Vaex 2.0.2

https://github.com/vaexio/vaex

安装并运行一个hello world示例，具体可以参考：如何在Vaex中进行交互式2D散点图缩放/点选择？

我用最多1亿个点测试了vaex，它运行得很好，真棒！

它是“以Python脚本为主”，这对可重复性很重要，并且让我可以轻松与其他Python工具对接。

Jupyter的设置有一些复杂，但一旦我用virtualenv运行起来，效果非常好。

要加载我们的CSV文件，在Jupyter中运行：

import vaex
df = vaex.from_csv('10m.csv', names=['x', 'y', 'z'],)
df.plot_widget(df.x, df.y, backend='bqplot')

我们可以立即看到：

现在，我们可以用鼠标缩放、平移和选择点，更新速度非常快，所有操作都在10秒内完成。这里我放大查看了一些单独的点，并选择了其中几个（图像上淡淡的矩形）：

用鼠标选择后，这与使用df.select()方法的效果完全相同。因此，我们可以通过在Jupyter中运行以下命令来提取所选点：

df.to_pandas_df(selection=True)

输出的数据格式为：

        x       y        z   index
0 4525460 9050920 18101840 4525460
1 4525461 9050922 18101844 4525461
2 4525462 9050924 18101848 4525462
3 4525463 9050926 18101852 4525463
4 4525464 9050928 18101856 4525464
5 4525465 9050930 18101860 4525465
6 4525466 9050932 18101864 4525466

由于1000万个点运行得很好，我决定尝试10亿个点……结果也很好！

import vaex
df = vaex.open('1b.hdf5')
df.plot_widget(df.x, df.y, backend='bqplot')

要观察异常值，原始图中看不见的，我们可以参考如何在Vaex交互式Jupyter bqplot绘图小部件中更改点样式以使单个点更大且可见？并使用：

df.plot_widget(df.x, df.y, f='log', shape=128, backend='bqplot')

这将生成：

选择该点后：

我们获得了异常值的完整数据：

   x          y           z
0  500000000  1500000000  -1

这里是创作者的演示，使用了更有趣的数据集和更多功能：https://www.youtube.com/watch?v=2Tt0i823-ec&t=770

不过遗憾的是，它没有内置的sqlite支持：https://github.com/vaexio/vaex/issues/864

在Ubuntu 19.04上测试。

VisIt 2.13.3

网站：https://wci.llnl.gov/simulation/computer-codes/visit

许可证：BSD

由劳伦斯利物浦国家实验室开发，这是一个国家核安全局实验室，所以你可以想象，如果我能让它工作，处理1000万个点对它来说没什么。

安装：没有Debian包，只需从网站下载Linux二进制文件。可以不安装直接运行。另见：https://askubuntu.com/questions/966901/installing-visit

基于VTK，这是许多高性能绘图软件使用的后端库。用C语言编写。

经过3小时的UI调试，我终于让它工作了，确实解决了我的使用案例，详细信息见：https://stats.stackexchange.com/questions/376361/how-to-find-the-sample-points-that-have-statistically-meaningful-large-outlier-r

这是在本帖测试数据上的效果：

放大后的一些选择：

这是选择窗口：

在性能方面，VisIt表现非常好：每个图形操作要么只需很少的时间，要么是立即完成。当我需要等待时，它会显示“处理”消息和剩余工作百分比，GUI不会冻结。

由于1000万个点运行得很好，我还尝试了100000000个点（一个2.7G的CSV文件），但不幸的是它崩溃了/进入了奇怪的状态，我在htop中看到4个VisIt线程占用了我所有的16GiB内存，并且由于malloc失败而死掉了。

最初的入门过程有点痛苦：

• 如果你不是核弹工程师，许多默认设置感觉很糟糕。例如：
  • 默认点大小为1px（与我显示器上的灰尘混淆）
  • 坐标轴范围从0.0到1.0：如何在VisIt绘图程序中显示实际的坐标轴数值，而不是从0.0到1.0的分数？
  • 多窗口设置，选择数据点时会出现讨厌的多个弹出窗口
  • 显示你的用户名和绘图日期（通过“控制”>“注释”>“用户信息”可以去掉）
  • 自动定位的默认设置很糟糕：图例与坐标轴冲突，找不到标题自动化，所以必须手动添加标签并重新定位所有内容
• 功能太多，可能很难找到你想要的东西
• 手册非常有帮助，但它是386页的PDF文档，日期为“2005年10月1.5版”。我想知道他们是否用这个开发了三位一体！而且它是一个漂亮的Sphinx HTML，是在我最初回答这个问题后不久创建的
• 没有Ubuntu包。但预构建的二进制文件确实可以工作。

我将这些问题归因于：

• 它已经存在很长时间，使用了一些过时的GUI理念
• 你不能仅仅通过点击绘图元素来更改它们（例如坐标轴、标题等），而且功能很多，所以很难找到你想要的

我还喜欢一些LLNL的基础设施在这个repo中泄露。例如，看看docs/OfficeHours.txt和该目录中的其他文件！我为“周一早上的人”Brad感到遗憾！哦，回答机的密码是“杀了Ed”，别忘了。

Paraview 5.9.0

网站：https://www.paraview.org/

许可证：BSD

测试于：Ubuntu 20.10。

安装：

sudo apt install paraview

或者从网站下载预构建的最新版本。这是我为这次评测所做的，因为apt中的版本只有5.7.0。我下载了ParaView-5.9.0-MPI-Linux-Python3.8-64bit.tar.gz。

由Kitware和洛斯阿拉莫斯国家实验室开发，后来又由桑迪亚国家实验室（所以是另外两个NNSA实验室），所以我们再次期待它能轻松处理数据。此外，它也是基于VTK并用C++编写，这让人更有信心。

然而我感到失望：由于某种原因，1000万个点让GUI变得非常慢且无响应，无法使用。每当我点击某个东西，比如隐藏线条时，都会花费几十秒。我觉得在某个时刻它就出现了故障，完全不响应。

我对控制良好的、明确的“我现在正在工作，请稍等”的时刻没有意见，但GUI在此期间冻结？这不可接受。

htop显示Paraview使用了8个线程和3GB的内存，所以CPU和内存都没有达到极限。

在GUI方面，Paraview非常漂亮和现代，远比VisIt好，前提是它没有卡顿。

由于10m1.csv让它崩溃，我用11.csv进行了测试，看看除了性能外我是否能解决我的问题，答案是肯定的：

• paraview 11.csv
• 从弹出窗口选择CSV读取器
• 左侧的属性应用
• 在管道浏览器中右键点击CSV
• 添加过滤器 > 字母顺序 > 绘制数据。为什么绘图是一个过滤器？对第一次使用者来说不太直观，相关问题：paraview：从CSV文件绘制数据。我相信这在你理解过滤器可以做的更一般化的事情后会变得有意义，但仍然如此。
• 属性 > 应用
• 取消选择“使用索引作为x轴”
• X数组名称：字段0
• 系列参数去掉字段0和字段2
• 选择字段1并：
  • 线条样式：无
  • 标记样式：十字
  • 标记大小：根据需要增加或减少
• 在绘图上方的“矩形选择（s）”图标
• 选择异常值（点被高亮显示）
• 在绘图过滤器中添加另一个过滤器：“提取选择”
• 应用

最后！！！我得到了一个只包含所选异常值的表格，并显示“字段2”的值为-1：

所以，是的，这并不是轻而易举，但我最终还是成功了。

另一个缺点是，Paraview的功能相比VisIt显得不足，例如：

• 我找不到如何根据第三列的值设置散点的颜色：如何根据第三列的值为Paraview中的散点着色，类似于gnuplot调色板？

Mayavi 4.6.2

网站：https://github.com/enthought/mayavi

开发者：Enthought

安装：

sudo apt-get install libvtk6-dev
python3 -m pip install -u mayavi PyQt5

VTK Python版本。

Mayavi似乎非常专注于3D，我找不到如何在其中进行2D绘图，所以不幸的是这不适合我的使用案例。

不过为了检查性能，我根据https://docs.enthought.com/mayavi/mayavi/auto/example_scatter_plot.html的示例调整了10百万点的示例，它运行得很好，没有延迟：

import numpy as np
from tvtk.api import tvtk
from mayavi.scripts import mayavi2

n = 10000000
pd = tvtk.PolyData()
pd.points = np.linspace((1,1,1),(n,n,n),n)
pd.verts = np.arange(n).reshape((-1, 1))
pd.point_data.scalars = np.arange(n)

@mayavi2.standalone
def main():
   from mayavi.sources.vtk_data_source import VTKDataSource
   from mayavi.modules.outline import Outline
   from mayavi.modules.surface import Surface
   mayavi.new_scene()
   d = VTKDataSource()
   d.data = pd
   mayavi.add_source(d)
   mayavi.add_module(Outline())
   s = Surface()
   mayavi.add_module(s)
   s.actor.property.trait_set(representation='p', point_size=1)
main()

输出：

不过我无法放大到足够看到单个点，近3D平面太远了。也许有办法？

Mayavi的一个酷点是，开发者花了很多精力让你能够通过Python脚本优雅地启动和设置GUI，类似于Matplotlib和gnuplot。这似乎在Paraview中也是可能的，但文档不如Mayavi好。

总体而言，它的功能似乎不如VisIt / Paraview丰富。例如，我无法直接从GUI加载CSV：如何从Mayavi GUI加载CSV文件？

Gnuplot 5.2.2

网站：http://www.gnuplot.info/

gnuplot在我需要快速处理时非常方便，通常是我尝试的第一件事。

安装：

sudo apt-get install gnuplot

对于非交互式使用，它可以合理地处理1000万个点：

#!/usr/bin/env gnuplot
set terminal png size 1024,1024
set output "gnuplot.png"
set key off
set datafile separator ","
plot "10m1.csv" using 1:2:3:3 with labels point

这在7秒内完成：

但如果我尝试进行交互式操作：

#!/usr/bin/env gnuplot
set terminal wxt size 1024,1024
set key off
set datafile separator ","
plot "10m.csv" using 1:2:3 palette

和：

gnuplot -persist main.gnuplot

那么初始渲染和缩放感觉太慢。我甚至看不到矩形选择线！

还要注意的是，对于我的用例，我需要使用超文本标签，如：

plot "10m.csv" using 1:2:3 with labels hypertext

但标签功能存在性能bug，包括非交互式渲染。但我报告了这个问题，Ethan在一天内解决了：https://groups.google.com/forum/#!topic/comp.graphics.apps.gnuplot/qpL8aJIi9ZE

不过我必须说，对于异常值选择，有一个合理的解决方法：只需为所有点添加带有行ID的标签！如果有很多点在一起，你可能无法读取标签。但对于你关心的异常值，你可能会看到！例如，如果我在原始数据中添加一个异常值：

cp 10m.csv 10m1.csv
printf '2500000,10000000,40000000\n' >> 10m1.csv

并将绘图命令修改为：

#!/usr/bin/env gnuplot
set terminal png size 1024,1024
set output "gnuplot.png"
set key off
set datafile separator ","
plot "10.csv" using 1:2:3:3 palette with labels

这显著减慢了绘图速度（在上述修复后40分钟！！！），但产生了合理的输出：

所以通过一些数据过滤，我们最终会达到目标。

Matplotlib 1.5.1, numpy 1.11.1, Python 3.6.7

网站：https://matplotlib.org/

当我的gnuplot脚本变得太复杂时，我通常会尝试Matplotlib。

numpy.loadtxt单独就花了大约10秒，所以我知道这不会顺利：

首先，非交互式尝试给出了不错的输出，但花了3分钟55秒……

然后交互式的在初始渲染和缩放时花了很长时间，无法使用：

注意在这个截图中，缩放选择本应立即缩放并消失，但却在等待缩放计算时长时间停留在屏幕上！

出于某种原因，我不得不注释掉plt.figure(figsize=(8, 8), dpi=128)，否则它会崩溃，显示：

RuntimeError: In set_size: Could not set the fontsize

Bokeh 1.3.1

https://github.com/bokeh/bokeh

Ubuntu 19.04安装：

python3 -m pip install bokeh

然后启动Jupyter：

jupyter notebook

现在如果我绘制1百万个点，一切都运行得很好，界面很棒且快速，包括缩放和悬停信息：

初始视图：

放大后：

但如果我尝试到1000万个点，它就卡住了，htop显示chromium有8个线程占用我所有的内存，处于不可中断的IO状态。

这涉及到引用点的问题：如何引用选定的bokeh数据点

PyViz

https://pyviz.org/

待评估。

集成了Bokeh + datashader + 其他工具。

演示1B数据点的视频：https://www.youtube.com/watch?v=k27MJJLJNT4 “PyViz：在30行Python中可视化10亿数据点的仪表板”由“Anaconda, Inc.”于2018年4月17日发布。

seaborn

https://seaborn.pydata.org/

待评估。

已经有一个关于如何使用seaborn可视化至少5000万行的问答。

sqlitebrowser 3.12.2

https://github.com/sqlitebrowser/sqlitebrowser

我尝试这个工具看看它是否能处理10m1.sqlite，但不幸的是它不能。真可惜！

不过，它可以直接绘制查询结果，这点很酷。

这是它的界面：

在这个图像中，我将10m1.sqlite加载到工具中，然后开始浏览数据。

但它只绘制了加载用于浏览的数据。

你可以点击右下角的按钮“加载所有数据并重新绘制图”，但这会打开一个进度条，每3秒上升1%，所以看起来不太乐观，我放弃了。

在Ubuntu 23.04上测试。

SQL直方图查询

我想知道为什么我找不到一个容易找到的使用这个作为后端的交互式UI工具。基于索引数据库的SQL直方图感觉是处理事情最合理的方式。例如，使用10个步骤并忽略空箱：

div=10
x=0
y=0
x2=10000000
y2=20000000

dx=$(((x2 - x) / div))
dy=$(((y2 - y) / div))
time sqlite3 10m1.sqlite --cmd '.mode csv' <<EOF
select
  floor(x/$dx)*$dx as x,
  floor(y/$dy)*$dy as y,
  count(*) as cnt
from t
where
  x >= $x and x < $x2 and
  y >= $y and y < $y2
group by 1, 2
order by 1, 2
EOF

我们得到：

0,0,1000000
1000000,2000000,1000000
2000000,4000000,1000000
3000000,6000000,1000000
4000000,8000000,1000000
5000000,10000000,1000000
5000000,20000000,1
6000000,12000000,1000000
7000000,14000000,1000000
8000000,16000000,1000000
9000000,18000000,1000000

查询花费6秒，所以它可以勉强处理1000万个点，但无法扩展到10亿。

由于我们已经有了1个点，我们可以在该范围内进行完整列出：

x=5000000
y=20000000
x2=6000000
y2=40000000
time sqlite3 10m1.sqlite --cmd '.mode csv' <<EOF
select *
from t
where
  x >= $x and x < $x2 and
  y >= $y and y < $y2
order by x, y
EOF

这会立即给出最终所需的结果：

5000000,20000000,-1

所以这个GUI可以有一个最大点数限制，如果：

• 如果超过限制，使用热图
• 否则，在该箱中查询所有单独的点并在图中绘制单独的点

如何将SQL扩展到10亿行：R树索引

要扩展到10亿，我们需要r树/空间索引，这允许我们高效地对多个列进行不等式操作。SQLite有它们，但使用起来有点麻烦：

• 所有内容都被转换为32位浮点数，没有64位或整数：https://sqlite-users.sqlite.narkive.com/mom9X2r3/sqlite-is-it-possible-to-support-64-bit-value-in-rtree-module
• 无法重用现有表中的数据
• 只能存储矩形，不能存储点：https://sqlite-users.sqlite.narkive.com/PMLhd3md/r-tree-storage-optimization-for-points

尽管有这些限制，我还是做了一个100百万点的测试，重复的x/y列，创建时间：30分钟，文件大小：5.9 GB

然后进行一次计数限制扫描：

max=100
div=10
x=0
y=0
x2=100000000
y2=200000000
dx=$(((x2 - x) / div))
dy=$(((y2 - y) / div))

cx=0
while [ $cx -lt $x2 ]; do
  cy=0
  while [ $cy -lt $y2 ]; do
    printf "$cx,$cy,"
    sqlite3 100mr.sqlite --cmd '.mode csv' <<EOF
select count(x) from (
  select x from t
  where
    x >= $cx and x < $((cx + dx)) and
    y >= $cy and y < $((cy + dy))
  limit $max
)
EOF
    cy=$((cy+dy))
  done
  cx=$((cx+dx))
done

仅用0.2秒完成，这太惊人了。如果不是因为生成时间太长，它可能会扩展到10亿。

不幸的是，PostgreSQL的索引创建速度也没有更快：如何将使用SQLite R树的简单空间索引移植到Postgres？，不过至少它支持点而不仅仅是矩形。

在Ubuntu 23.04上测试。

Answer 3

你的数据量并不大，但你在绘图时遇到问题，这说明你使用的工具可能有问题。Matplotlib有很多选项，输出效果也不错，但它非常占内存，并且基本上是为了处理小数据量而设计的。不过，还有其他的选择。

举个例子，我用以下代码生成了一个包含2000万数据点的文件'bigdata.bin'：

#!/usr/bin/env python
import numpy
import scipy.io.numpyio

npts=20000000
filename='bigdata.bin'

def main():
    data = (numpy.random.uniform(0,1,(npts,3))).astype(numpy.float32)
    data[:,2] = 0.1*data[:,2]+numpy.exp(-((data[:,1]-0.5)**2.)/(0.25**2))
    fd = open(filename,'wb')
    scipy.io.numpyio.fwrite(fd,data.size,data)
    fd.close()

if __name__ == "__main__":
    main()

这个文件大约有229MB，虽然不算太大，但你提到想处理更大的文件，最终会遇到内存限制。

我们先来关注非交互式的绘图。首先要明白的是，给每个点绘制小图形（比如小十字或小圆圈）会非常糟糕——对于这2000万个点，大部分会重叠，试图渲染这些小图形会导致生成巨大的文件，并且耗费大量时间。我认为这就是Matplotlib默认情况下的一个问题。

而Gnuplot在处理这方面时就没问题：

gnuplot> set term png
gnuplot> set output 'foo.png'
gnuplot> plot 'bigdata.bin' binary format="%3float32" using 2:3 with dots

即使是Matplotlib，只要小心使用（选择光栅后端，用像素标记点），也能表现得不错：

#!/usr/bin/env python
import numpy
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

datatype=[('index',numpy.float32), ('floati',numpy.float32), 
        ('floatq',numpy.float32)]
filename='bigdata.bin'

def main():
    data = numpy.memmap(filename, datatype, 'r') 
    plt.plot(data['floati'],data['floatq'],'r,')
    plt.grid(True)
    plt.title("Signal-Diagram")
    plt.xlabel("Sample")
    plt.ylabel("In-Phase")
    plt.savefig('foo2.png')

if __name__ == "__main__":
    main()  

如果你想要交互式的效果，就需要对数据进行分箱处理，以便动态缩放。我暂时不知道有什么Python工具可以帮助你做到这一点。

另一方面，绘制大数据是一个相当常见的任务，有一些工具可以胜任这个工作。Paraview是我个人的最爱，VisIt也是一个不错的选择。它们主要用于处理3D数据，但特别是Paraview也支持2D，并且非常交互（甚至有Python脚本接口）。唯一需要注意的是，要将数据写入Paraview能轻松读取的文件格式。