在 ClickHouse，我们经常需要下载和探索越来越大的数据集，这总是需要快速的互联网连接。当考虑搬家和出国，并需要确保我有快速的互联网时，我发现了规模庞大的 Ookla 数据集，它提供了全球所有位置的互联网性能数据。该数据集是世界上最大的众包网络测试数据来源，涵盖了固定和移动互联网性能，覆盖了世界上大部分陆地表面。
虽然 Ookla 已经提供了 出色的工具，可以使用此数据集按位置探索互联网性能，但这为探索 ClickHouse 的 Apache Iceberg 支持以及不断扩展的地理功能提供了机会。
进一步缩小范围，我很好奇仅使用 SQL 通过将地理多边形转换为 SVG 来可视化地理数据的可行性。这将带领我们游览 h3 索引、多边形字典、墨卡托投影、颜色插值以及使用 UDF 的质心计算。如果您好奇我们是如何仅使用 SQL 就得到以下结果的，请继续阅读……
注意：ClickHouse 支持多种适用于此数据集的可视化工具，包括但不限于 Superset。我们仅限于使用 SQL 只是为了好玩！
在深入研究数据集之前，我们将快速回顾一下 Apache Iceberg 的价值及其对 ClickHouse 的支持和相关性。对于熟悉此主题的用户，请随意跳转到我们开始查询的部分！
Apache Iceberg
Apache Iceberg 近年来备受关注，它对于将“数据湖”概念演变为“湖仓一体”至关重要。通过为数据提供高性能的表格式，Iceberg 将类似表的语义引入到存储在数据湖中的结构化数据中。更具体地说，这种开放表格式是供应商中立的，并允许数据文件（通常是结构化的 Parquet 格式）作为表公开，并支持模式演变、删除、更新、事务，甚至 ACID 合规性。这代表了从基于 Hadoop 和 Hive 等技术的更传统的旧式数据湖的逻辑演变。
也许最重要的是，这种数据格式允许 Spark 和 Trino 等查询引擎安全地同时使用相同的表，同时还提供分区等功能，这些功能可能会被查询规划器利用来加速查询。
ClickHouse 支持
虽然 ClickHouse 不仅仅是一个查询引擎，但我们承认这些开放标准的重要性，它可以让用户避免供应商锁定，并提供与其他工具进行数据交换的合理格式。正如最近一篇博文 《云数据仓库的解绑》 中所讨论的那样，我们认为湖仓一体是现代数据架构中的一个基本组成部分，与实时分析数据库完美互补。在这种架构中，湖仓一体提供了经济高效的冷存储和通用的数据源，并支持即席查询。
但是，对于更苛刻的用例，例如需要低延迟查询和高并发支持的面向用户的数据密集型应用程序，可以将子集移动到 ClickHouse。这允许用户从 Iceberg 的开放交换格式中受益，用于其冷存储和长期保留，并将 ClickHouse 用作实时数据仓库，以在查询性能至关重要时使用。
使用 ClickHouse 查询 Iceberg
实现此愿景需要能够从 ClickHouse 查询 Iceberg 文件。作为一种湖仓一体格式，目前通过专用的表 iceberg 函数支持此功能。这假设数据托管在与 S3 兼容的服务上，例如 AWS S3 或 Minio。
对于我们的示例，我们已将 Ookla 数据集放在以下公共 S3 存储桶中，供用户实验。
s3://datasets-documentation/ookla/iceberg/
我们可以使用 DESCRIBE 查询来描述 Ookla 数据的模式，其中 iceberg 函数替换了我们的表名。与依赖于对实际数据进行抽样的 Parquet 不同，模式是从 Iceberg 元数据文件中读取的。
DESCRIBE TABLE iceberg('https://datasets-documentation.s3.eu-west-3.amazonaws.com/ookla/iceberg/')
SETTINGS describe_compact_output = 1
┌─name────────────┬─type─────────────┐
│ quadkey │ Nullable(String) │
│ tile │ Nullable(String) │
│ avg_d_kbps │ Nullable(Int32) │
│ avg_u_kbps │ Nullable(Int32) │
│ avg_lat_ms │ Nullable(Int32) │
│ avg_lat_down_ms │ Nullable(Int32) │
│ avg_lat_up_ms │ Nullable(Int32) │
│ tests │ Nullable(Int32) │
│ devices │ Nullable(Int32) │
│ year_month │ Nullable(Date) │
└─────────────────┴──────────────────┘
10 rows in set. Elapsed: 0.216 sec.
同样，可以使用标准 SQL 执行行查询，其中表函数替换了表名。在以下示例中，我们对一些行进行抽样，并计算数据集的总大小。
SELECT *
FROM iceberg('https://datasets-documentation.s3.eu-west-3.amazonaws.com/ookla/iceberg/')
LIMIT 1
FORMAT Vertical
Row 1:
──────
quadkey: 1202021303331311
tile: POLYGON((4.9163818359375 51.2206474303833, 4.921875 51.2206474303833, 4.921875 51.2172068072334, 4.9163818359375 51.2172068072334, 4.9163818359375 51.2206474303833))
avg_d_kbps: 109291
avg_u_kbps: 15426
avg_lat_ms: 12
avg_lat_down_ms: ᴺᵁᴸᴸ
avg_lat_up_ms: ᴺᵁᴸᴸ
tests: 6
devices: 4
year_month: 2021-06-01
1 row in set. Elapsed: 2.100 sec. Processed 8.19 thousand rows, 232.12 KB (3.90 thousand rows/s., 110.52 KB/s.)
Peak memory usage: 4.09 MiB.
SELECT count()
FROM iceberg('https://datasets-documentation.s3.eu-west-3.amazonaws.com/ookla/iceberg/')
┌───count()─┐
│ 128006990 │
└───────────┘
1 row in set. Elapsed: 0.701 sec. Processed 128.01 million rows, 21.55 KB (182.68 million rows/s., 30.75 KB/s.)
Peak memory usage: 4.09 MiB
我们可以看到，Ookla 使用 WKT 格式将地球表面划分为多边形。我们将在博客的后半部分深入探讨这一点；现在，我们将完成我们的简单示例，其中一个聚合计算每个多边形的平均设备数量，以及下载速度的中位数、第 90、99 和 99.9 百分位数。
SELECT
round(avg(devices), 2) AS avg_devices,
arrayMap(m -> round(m / 1000), quantiles(0.5, 0.9, 0.99, 0.999)(avg_d_kbps)) AS download_mbps
FROM iceberg('https://datasets-documentation.s3.eu-west-3.amazonaws.com/ookla/iceberg/')
┌─avg_devices─┬─download_mbps────┐
│ 5.64 │ [51,249,502,762] │
└─────────────┴──────────────────┘
1 row in set. Elapsed: 4.704 sec. Processed 128.01 million rows, 3.75 GB (27.21 million rows/s., 797.63 MB/s.)
Peak memory usage: 22.89 MiB.
以上内容提供了一些关于其实时数据仓库如何运行的简单示例：湖仓一体中的冷 Iceberg 存储和 ClickHouse 中用于实时分析的热查询。
一些限制
以上是一些非常简单的查询。在我们深入探讨如何可视化此数据的更深层次主题之前，值得强调 ClickHouse 的 Iceberg 支持在哪些方面需要进一步改进（并且已计划）。
- 目前，ClickHouse 不会读取具有基于行的删除 (v2) 或其模式已演变的 Iceberg 文件。除非用户显式请求忽略，否则在这些情况下，用户将收到读取异常。
- 在查询评估期间，ClickHouse 当前不会利用分区来修剪扫描的数据。这代表了提高查询执行效率的巨大可能性，并且与 最近对 Parquet 所做的改进保持一致。
- 当前的实现绑定到 S3。未来，我们希望消除这种耦合，并添加对 ClickHouse 支持的任何对象存储的支持。
- 与 ClickHouse 的排序键类似，Iceberg 支持对数据进行排序。通过确保 Iceberg 表中的数据文件被拆分和排序，查询可能会在执行期间跳过文件。这有可能进一步减少扫描的数据，从而提高性能。
探索数据
在介绍了基础知识之后，让我们稍微改变一下本博客的方向，更详细地探索 Ookla 数据，目标是仅使用 SQL 可视化互联网速度。
读取和丰富
如前所述，Ookla 数据将地球划分为大小不同的多边形，从中收集互联网性能统计数据。这些多边形在 WKT 格式中定义，WKT 格式是一种用于表示空间对象的标记语言。这种格式对于 Iceberg 和数据分发目的来说是有意义的，因为底层的 Parquet 文件没有原生的多边形表示。
POLYGON((-160.02685546875 70.6435894914449, -160.021362304688 70.6435894914449, -160.021362304688 70.6417687358462, -160.02685546875 70.6417687358462, -160.02685546875 70.6435894914449))
但是，查询此数据需要转换为数值表示。对于 ClickHouse，可以使用 readWKTPolygon 函数执行此转换。如下所示，这会将文本表示转换为 Polygon 类型，我们稍后可以在此类型上使用 ClickHouse 的大量几何函数。请注意，这些多边形包含经度和纬度值。
SELECT readWKTPolygon(assumeNotNull(tile)) AS poly, toTypeName(poly) as type
FROM iceberg('https://datasets-documentation.s3.eu-west-3.amazonaws.com/ookla/iceberg/')
LIMIT 1
FORMAT Vertical
Row 1:
──────
poly: [[(4.9163818359375,51.2206474303833),(4.921875,51.2206474303833),(4.921875,51.2172068072334),(4.9163818359375,51.2172068072334),(4.9163818359375,51.2206474303833)]]
type: Polygon
1 row in set. Elapsed: 1.920 sec. Processed 8.19 thousand rows, 232.12 KB (4.27 thousand rows/s., 120.91 KB/s.)
Peak memory usage: 4.10 MiB.
要温和地介绍多边形和解析 WKT 数据，我们推荐我们自己的 “Data with Mark” 系列视频。
请注意，以上视频使用了 Ookla 数据集的最新文件，其中提供了质心。
虽然 Ookla 数据集的更新版本包含每个多边形的预计算质心（从 2023 年第二季度开始），但该数据集可以追溯到 2019 年。正如我们稍后将展示的那样，此质心（纬度/经度对）在渲染数据时可能很有用。幸运的是，由于我们的多边形是非自相交且闭合的，因此计算起来相对简单。对于好奇的读者，下面显示了一些示例 Python 代码（来源）。
def centroid(vertices):
x, y = 0, 0
n = len(vertices)
signed_area = 0
# computes the signed area and y and x accumulators - signed_area UDF function
for i in range(len(vertices)):
x0, y0 = vertices[i]
x1, y1 = vertices[(i + 1) % n]
# shoelace formula - maps to area function below
area = (x0 * y1) - (x1 * y0)
signed_area += area
x += (x0 + x1) * area
y += (y0 + y1) * area
# final step of centroid function
signed_area *= 0.5
x /= 6 * signed_area
y /= 6 * signed_area
return x, y
x, y = centroid([(-160.037841796875, 70.6363054807905), (-160.032348632812, 70.6363054807905),(-160.032348632812, 70.6344840663086), (-160.037841796875, 70.6344840663086), (-160.037841796875, 70.6363054807905)])
print(x,y)
python centroid.py
-160.03509089358298 70.63539286610323
如果我们将问题分解为几个函数，则将其转换为可在查询时使用的 SQL UDF 最为简单。
CREATE FUNCTION area AS (v1,v2) -> ((v1.1 * v2.2) - (v2.1 * v1.2))
CREATE FUNCTION signed_area AS polygon ->
(
arrayFold((acc, x) -> (
(acc.1) + ((((x.1).1) + ((polygon[((x.2) + 1) % (length(polygon) + 1)]).1)) * area(x.1, polygon[((x.2) + 1) % (length(polygon) + 1)])),
(acc.2) + ((((x.1).2) + ((polygon[((x.2) + 1) % (length(polygon) + 1)]).2)) * area(x.1, polygon[((x.2) + 1) % (length(polygon) + 1)])),
(acc.3) + area(x.1, polygon[((x.2) + 1) % (length(polygon) + 1)])), arrayZip(polygon, range(1, length(polygon) + 1)), (0.0, 0.0, 0.0)
)
)
CREATE FUNCTION centroid as polygon -> ((signed_area(polygon).1) / ((6 * (signed_area(polygon).3)) * 0.5), (signed_area(polygon).2) / ((6 * (signed_area(polygon).3)) * 0.5))
快速测试证实我们的 UDF 给出的结果与 Python 实现相同。
SELECT centroid([(-160.037841796875, 70.6363054807905), (-160.032348632812, 70.6363054807905), (-160.032348632812, 70.6344840663086), (-160.037841796875, 70.6344840663086), (-160.037841796875, 70.6363054807905)]) AS centroid
┌─centroid────────────────────────────────┐
│ (-160.03509089358298,70.63539286610323) │
└─────────────────────────────────────────┘
1 row in set. Elapsed: 0.017 sec.
使用上述函数，我们现在可以读取我们的多边形数据并计算由经度和纬度组成的质心。正如我们将展示的那样，我们可以在可视化数据时利用这些质心。
加载数据
虽然以下所有查询都可以与 iceberg 表函数一起使用，但这需要为每个查询下载数据，通常会将性能限制为我们的网络带宽。由于我们将运行多个查询，因此将数据加载到 ClickHouse 中以享受前面强调的性能优势是有意义的。我们修改了我们的表模式，以包含转换后的多边形（Ring 类型 - 见下文）和质心（Point 类型），并使用上面的函数提取这些数据。
CREATE TABLE default.ookla
(
`quadkey` Nullable(String),
`polygon` Ring,
`centroid` Point,
`avg_d_kbps` Nullable(Int32),
`avg_u_kbps` Nullable(Int32),
`avg_lat_ms` Nullable(Int32),
`avg_lat_down_ms` Nullable(Int32),
`avg_lat_up_ms` Nullable(Int32),
`tests` Nullable(Int32),
`devices` Nullable(Int32),
`year_month` Nullable(Date)
)
ENGINE = MergeTree ORDER BY tuple()
INSERT INTO ookla SELECT
quadkey,
readWKTPolygon(assumeNotNull(tile))[1] AS polygon,
centroid(polygon) AS centroid,
avg_d_kbps,
avg_u_kbps,
avg_lat_ms,
avg_lat_down_ms,
avg_lat_up_ms,
tests,
devices,
year_month
FROM iceberg('https://datasets-documentation.s3.eu-west-3.amazonaws.com/ookla/iceberg/')
0 rows in set. Elapsed: 549.829 sec. Processed 128.01 million rows, 3.75 GB (232.81 thousand rows/s., 6.82 MB/s.)
Peak memory usage: 2.16 GiB.
读者可能会注意到，我们使用了 readWKTPolygon 函数返回的数组中的第一个元素。多边形可以由多个环组成，每个环都由一个封闭空间组成。第一个条目描述了外部空间，随后的数组条目指定了任何孔洞。由于我们的多边形是闭合的，没有内部孔洞，因此我们只有一个条目 - 允许我们将我们的多边形映射到更简单的 Ring 类型。
我们的示例仅包含固定设备的数据。对于好奇的用户，还存在移动设备的数据。
多边形到 SVG
如前所述，我们现在可以利用多种可视化工具来渲染此数据。但是，这些工具确实会产生开销 - 安装，更重要的是，获得熟悉度和专业知识的时间。鉴于 SVG 格式本质上是一种能够渲染多边形的基于矢量的格式，也许我们可以简单地将这些多边形转换为 SVG，以提供数据的 2D 可视化效果，作为地图？
这需要一种将每个多边形转换为 SVG 多边形条目 的方法。理想情况下，我们还应该能够设置每个多边形的样式以传达指标，例如，下载速度。与以往一样，社区已经考虑过这个问题，ClickHouse 提供了 SVG 函数！
这使我们能够将 ClickHouse 中的多边形转换为 SVG 多边形，以及传递样式
SELECT SVG([(-160.037841796875, 70.6363054807905), (-160.032348632812, 70.6363054807905), (-160.032348632812, 70.6344840663086), (-160.037841796875, 70.6344840663086), (-160.037841796875, 70.6363054807905)], 'fill: rgb(255, 55,38);') AS svg FORMAT Vertical
Row 1:
──────
svg: <polygon points="-160.038,70.6363 -160.032,70.6363 -160.032,70.6345 -160.038,70.6345 -160.038,70.6363" style="fill: rgb(255, 55,38);"/>
1 row in set. Elapsed: 0.002 sec.
点到像素
虽然以上内容很有希望，但如果以 2d 地图形式可视化，我们不能简单地使用坐标作为我们的纬度和经度。这些坐标表示三维表面（即地球表面）上的位置，而不是二维表面。存在许多称为投影的技术，用于将纬度和经度转换为二维表面（例如地图或计算机屏幕）上的点。此过程不可避免地涉及对地球表面表示形式的一些折衷或失真，因为不可能在不拉伸或压缩其表面的情况下展平球体。
地图上常用的一种流行的投影是 墨卡托投影。这有很多优点；最主要的是，它在局部范围内保留了角度和形状，使其成为导航的绝佳选择。它还具有恒向线（行进方向）是直线的优点，从而简化了导航。最后，虽然墨卡托投影准确地表示了赤道附近的区域，但当向两极移动时，它会显着扭曲大小和形状。与赤道地区相比，格陵兰岛和南极洲等陆地显得比实际大得多。
虽然 ClickHouse 没有内置墨卡托函数，但我们可以使用另一个 SQL UDF 解决此问题。在这里，我们将我们的点投影到尺寸为 1024x1024 的像素空间中。
-- coord (lon, lat) format
CREATE OR REPLACE FUNCTION mercator AS (coord, width, height) -> (
((coord.1) + 180) * (width / 360),
(height / 2) - ((width * ln(tan((pi() / 4) + ((((coord.2) * pi()) / 180) / 2)))) / (2 * pi()))
)
SELECT mercator((-160.037841796875, 70.6363054807905), 1024, 1024) AS pixel_coords
┌─pixel_coords──────────────────┐
│ (56.78125,223.79687511938704) │
└───────────────────────────────┘
1 row in set. Elapsed: 0.001 sec.
要将投影应用于整个多边形，我们可以利用 arrayMap 函数。
SELECT polygon, arrayMap(p -> mercator(p, 1024, 1024), polygon) AS pixels
FROM ookla
LIMIT 1
FORMAT Vertical
Row 1:
──────
polygon: [(-51.30615234375,-30.0358110426678),(-51.3006591796875,-30.0358110426678),(-51.3006591796875,-30.0405664305846),(-51.30615234375,-30.0405664305846),(-51.30615234375,-30.0358110426678)]
pixels: [(366.0625,601.6406251005949),(366.078125,601.6406251005949),(366.078125,601.6562501003996),(366.0625,601.6562501003996),(366.0625,601.6406251005949)]
1 row in set. Elapsed: 0.012 sec.
大小限制
借助我们的 SVG 和墨卡托函数，我们完全有能力生成我们的第一个 SVG。以下查询迭代所有 1.28 亿个多边形，并将它们转换为 svg 元素。我们还输出一个包含背景颜色和 viewBox 坐标的开头和结尾 <svg> 标记。
对于此初始查询，我们想确认此渲染技术有效，因此，使用固定样式，并且目前不考虑任何指标。我们将此查询的输出重定向到一个文件。请注意，我们使用输出格式 CustomSeparated 来生成我们的 svg 开头和结尾标记。
clickhouse client --query "SELECT SVG(arrayMap(p -> mercator(p, 1024, 1024), polygon), 'fill: rgb(255, 55,38);')
FROM ookla
FORMAT CustomSeparated
SETTINGS format_custom_result_before_delimiter = '<svg style=\"background-color: rgb(17,17,17);\" viewBox=\"0 0 1024 1024\" xmlns:svg=\"http://www.w3.org/2000/svg\" xmlns=\"http://www.w3.org/2000/svg\">', format_custom_result_after_delimiter = '</svg>', format_custom_escaping_rule = 'Raw'" > world.svg
对于许多人来说，可以预测的是，1.28 亿个多边形会生成一个相当大的 svg 文件！
(base) clickhouse@ClickHouse-MacBook-Pro ookla % ls -lh world.svg
-rw-r--r-- 1 dalemcdiarmid wheel 31G 2 Feb 17:04 world.svg
31GB 的大小，大多数机器上的浏览器或工具都无法渲染它，甚至使用 ImageMagick 等工具将其转换为压缩格式也是不现实的。如果您有关于如何实现此目标的想法，我们欢迎您的建议！
使用字典进行汇总
我们的第一个努力可能是将数据汇总为代表地理区域的更大的多边形。为此，我们需要一种方法来确定我们的小多边形位于哪些更大的多边形中。这可以使用 polygonsWithinCartesian 等函数来实现，但对于 1.28 亿个多边形来说，计算成本会很高。
幸运的是，ClickHouse 具有 多边形字典 来解决这个问题。这些字典允许用户有效地搜索包含特定点的多边形。这需要我们的小多边形中的每一个都由一个单点表示，例如质心 :)
例如，假设我们想将我们的多边形汇总到国家/地区。我们将我们的目标多边形存储在 countries 表中
CREATE TABLE countries
(
`name` String,
`coordinates` MultiPolygon
)
ENGINE = MergeTree
ORDER BY name
INSERT INTO countries
SELECT name, coordinates
FROM (
SELECT JSONExtractString(JSONExtractString(json, 'properties'), 'ADMIN') as name,
JSONExtractString(JSONExtractRaw(json, 'geometry'), 'type') as type,
if(type == 'Polygon', [JSONExtract(JSONExtractRaw(JSONExtractRaw(json, 'geometry'), 'coordinates'), 'Polygon')], JSONExtract(JSONExtractRaw(JSONExtractRaw(json, 'geometry'), 'coordinates'), 'MultiPolygon')) as coordinates
FROM (SELECT arrayJoin(JSONExtractArrayRaw(json, 'features')) as json
FROM url('https://datasets-documentation.s3.eu-west-3.amazonaws.com/countries/countries.geojson', JSONAsString)))
我们的国家/地区多边形以 geoJSON 形式表示，并且可在此处获得。上面的 INSERT SELECT 使用 JSONExtract* 函数从 JSON 中提取这些多边形，并将它们存储在 countries 表中。
我们再次有一个关于此主题的视频，涵盖了查询和后续字典
CREATE DICTIONARY country_polygons
(
`coordinates` MultiPolygon,
`name` String
)
PRIMARY KEY coordinates
SOURCE(CLICKHOUSE(TABLE 'countries'))
LIFETIME(MIN 0 MAX 0)
LAYOUT(POLYGON(STORE_POLYGON_KEY_COLUMN 1))
此字典按我们的多边形的坐标索引，并允许我们使用 dictGet 函数检索与特定纬度和经度关联的多边形。例如，
SELECT dictGet(country_polygons, 'name', (-2.72647, 55.599621)) AS location
┌─location───────┐
│ United Kingdom │
└────────────────┘
1 row in set. Elapsed: 0.012 sec.
请注意此处的检索速度。除了在我们稍后的可视化中很有用之外，此字典还允许我们快速了解不同国家/地区的互联网速度。例如，我们可以快速确定南乔治亚岛可能是居住以获得快速互联网的最佳地点！
SELECT dictGet(country_polygons, 'name', centroid) AS country,
avg(avg_d_kbps) AS download_speed,
sum(devices) AS total_devices
FROM ookla
GROUP BY country
ORDER BY download_speed DESC
LIMIT 10
┌─country──────────────────────────────────┬─────download_speed─┬─total_devices─┐
│ South Georgia and South Sandwich Islands │ 198885 │ 4 │
│ Singapore │ 197264.31481017102 │ 3591446 │
│ Hong Kong S.A.R. │ 166375.523948868 │ 3891226 │
│ China │ 158978.10242450528 │ 53841880 │
│ Thailand │ 155309.49656953348 │ 27506456 │
│ Iceland │ 152290.62935414084 │ 281514 │
│ United States of America │ 149532.95964270137 │ 246966470 │
│ South Korea │ 149115.5647068738 │ 2129690 │
│ United Arab Emirates │ 148936.4285346194 │ 8380900 │
│ Jersey │ 146030.67963711882 │ 100902 │
└──────────────────────────────────────────┴────────────────────┴───────────────┘
10 rows in set. Elapsed: 5.467 sec. Processed 256.01 million rows, 6.66 GB (46.83 million rows/s., 1.22 GB/s.)
Peak memory usage: 129.74 MiB.
最后，正如我们将在下面看到的那样，能够根据国家/地区的名称有效地检索该国家/地区的多边形也可能很有帮助。这可以使用使用相同 countries 表的简单字典来实现。
CREATE DICTIONARY default.country_names
(
`coordinates` MultiPolygon,
`name` String
)
PRIMARY KEY name
SOURCE(CLICKHOUSE(TABLE 'countries'))
LIFETIME(MIN 0 MAX 0)
LAYOUT(COMPLEX_KEY_HASHED())
SELECT dictGet(country_names, 'coordinates', 'Yemen') AS coords
[[[(53.308239457000155,12.11839453700017),(53.310267307000146,12.11144367200015)..
1 row in set. Elapsed: 0.007 sec.
有关字典的更多详细信息，我们推荐以下博文。
全球 SVG
借助我们的字典、质心、墨卡托 UDF 和 SVG 函数，我们可以生成我们的第一个（可用）SVG。
上面的 SVG 为 9.2MB，由以下查询生成，该查询在 11 秒内完成。如键所示，红色表示慢速互联网速度（下载速率），而绿色表示更快的连接速度。
WITH
-- settings determining height, width and color range
1024 AS height, 1024 AS width, [234, 36, 19] AS lc, [0, 156, 31] AS uc,
country_rates AS
(
-- compute the average download speed per country using country_polygons dictionary
SELECT
dictGet(country_polygons, 'name', centroid) AS country,
avg(avg_d_kbps) AS download_speed
FROM ookla
-- exclude 'Antarctica' as polygon isnt fully closed
WHERE country != 'Antarctica'
GROUP BY country
),
(
-- compute min speed over all countries
SELECT min(download_speed)
FROM country_rates
) AS min_speed,
(
-- compute max speed over all countries
SELECT max(download_speed)
FROM country_rates
) AS max_speed,
country_colors AS
(
SELECT
country,
download_speed,
(download_speed - min_speed) / (max_speed - min_speed) AS percent,
-- compute a rgb value based on linear interolation between yellow and red
format('rgb({},{},{})', toUInt8((lc[1]) + (percent * ((uc[1]) - (lc[1])))), toUInt8((lc[2]) + (percent * ((uc[2]) - (lc[2])))), toUInt8((lc[3]) + (percent * ((uc[3]) - (lc[3]))))) AS rgb,
-- get polygon for country name and mercator project every point
arrayMap(p -> arrayMap(r -> arrayMap(x -> mercator(x, width, height), r), p), dictGet(country_names, 'coordinates', country)) AS pixel_poly
FROM country_rates
)
-- convert country polygon to svg
SELECT SVG(pixel_poly, format('fill: {};', rgb))
FROM country_colors
-- send output to file
INTO OUTFILE 'countries.svg'
FORMAT CustomSeparated
SETTINGS format_custom_result_before_delimiter = '<svg style="background-color: rgb(17,17,17);" viewBox="0 0 1024 1024" xmlns:svg="http://www.w3.org/2000/svg" xmlns="http://www.w3.org/2000/svg">', format_custom_result_after_delimiter = '</svg>', format_custom_escaping_rule = 'Raw'
246 rows in set. Elapsed: 1.396 sec. Processed 384.02 million rows, 8.06 GB (275.15 million rows/s., 5.78 GB/s.)
Peak memory usage: 39.54 MiB.
尽管只是组合了上述概念，但这里有很多内容。我们对 SQL 进行了注释，希望这有所帮助，但总结如下
- 我们在第一个 CTE 中声明常量，例如颜色渐变（lc = 黄色，uc = 红色）和图像高度。
- 对于每个国家/地区，我们使用
country_polygons字典来计算平均速度，类似于我们之前的查询。 - 对于上述值，我们确定最大和最小互联网速度。这些用于为每个国家/地区生成颜色，使用我们颜色范围之间的简单 线性插值。
- 对于每个国家/地区，使用其名称，我们使用反向字典
country_names来获取多边形。这些多边形很复杂，有许多内部孔洞以 Rings 形式表示。因此，我们需要使用嵌套的arrayMap函数来使用墨卡托投影函数投影它们。 - 最后，我们使用 SVG 函数和 RGB 值将投影的国家/地区多边形转换为 SVG。这使用 INTO OUTFILE 子句输出到本地文件。
使用 h3 索引进行汇总
虽然上述查询证明了从我们的地理数据渲染 SVG 的概念，但它的汇总过于粗略。按国家/地区可视化互联网速度很有趣，但只是在之前的查询（列出了按国家/地区划分的互联网速度）的基础上略有构建。
我们可以尝试使用来自 公共来源的更大的多边形来填充我们的 country_polygons 字典。但是，这似乎有点限制性，理想情况下，我们应该能够通过指定的分辨率生成我们的矢量，该分辨率控制我们渲染的多边形的大小。
H3 是 Uber 开发的一种地理索引系统，其中地球表面被划分为均匀的六边形单元格网格。此系统是分层的，即，顶层（“父级”）上的每个六边形都可以拆分为七个更小的均匀六边形（“子级”），依此类推。使用索引系统时，我们定义介于 0 和 15 之间的分辨率。分辨率越大，每个六边形越小，因此覆盖地球表面所需的数量也越多。
可以使用一个简单的 ClickHouse 查询和函数 h3NumHexagons来说明分辨率在多边形数量上的差异，该函数返回每个分辨率的六边形数量。
SELECT h3NumHexagons(CAST(number, 'UInt8')) AS num_hexagons
FROM numbers(0, 16)
┌─res─┬────num_hexagons─┐
│ 0 │ 122 │
│ 1 │ 842 │
│ 2 │ 5882 │
│ 3 │ 41162 │
│ 4 │ 88122 │
│ 5 │ 2016842 │
│ 6 │ 14117882 │
│ 7 │ 98825162 │
│ 8 │ 691776122 │
│ 9 │ 4842432842 │
│ 10 │ 33897029882 │
│ 11 │ 237279209162 │
│ 12 │ 1660954464122 │
│ 13 │ 11626681248842 │
│ 14 │ 81386768741882 │
│ 15 │ 569707381193162 │
└─────┴─────────────────┘
16 rows in set. Elapsed: 0.001 sec.
ClickHouse 通过多个函数支持 h3，包括通过 geoToH3 函数将经度和纬度转换为具有指定分辨率的 h3 索引的功能。使用此函数，我们可以将我们的质心转换为 h3 索引，我们可以在该索引上聚合和计算统计数据，例如平均下载速度。聚合后，我们可以使用函数 h3ToGeoBoundary 反转此过程，以将 h3 索引及其六边形转换为多边形。最后，我们使用此墨卡托投影和 SVG 函数来生成我们的 SVG。我们在下面可视化此过程
上面，我们显示了分辨率为 6 的结果，它生成了 14,117,882 个多边形，文件大小为 165MB。这需要 Chrome 在 Macbook M2 上大约 1 分钟才能渲染。通过修改分辨率，我们可以控制我们的多边形的大小和我们可视化中的细节级别。如果我们将分辨率降低到 4，我们的文件大小也会减小到 9.2MB，因为我们的可视化变得更加粗糙，Chrome 能够几乎立即渲染它。
我们的查询在此处
Our query here:
WITH
1024 AS height, 1024 AS width, 6 AS resolution, 30 AS lc, 144 AS uc,
h3_rates AS
(
SELECT
geoToH3(centroid.1, centroid.2, resolution) AS h3,
avg(avg_d_kbps) AS download_speed
FROM ookla
– filter out as these are sporadic and don’t cover the continent
WHERE dictGet(country_polygons, 'name', centroid) != 'Antarctica'
GROUP BY h3
),
(
SELECT min(download_speed) FROM h3_rates
) AS min_speed,
(
SELECT quantile(0.95)(download_speed) FROM h3_rates
) AS max_speed,
h3_colors AS
(
SELECT
-- sqrt gradient
arrayMap(x -> (x.2), arrayConcat(h3ToGeoBoundary(h3), [h3ToGeoBoundary(h3)[1]])) AS longs,
sqrt(download_speed - min_speed) / sqrt(max_speed - min_speed) AS percent,
-- oklch color with gradient on hue
format('oklch(60% 0.23 {})', toUInt8(lc + (percent * (uc - lc)))) AS oklch,
arrayMap(p -> mercator((p.2, p.1), width, height), h3ToGeoBoundary(h3)) AS pixel_poly
FROM h3_rates
-- filter out points crossing 180 meridian
WHERE arrayExists(i -> (abs((longs[i]) - (longs[i + 1])) > 180), range(1, length(longs) + 1)) = 0
)
SELECT SVG(pixel_poly, format('fill: {};', oklch))
FROM h3_colors
INTO OUTFILE 'h3_countries.svg' TRUNCATE
FORMAT CustomSeparated
SETTINGS format_custom_result_before_delimiter = '<svg style="background-color: rgb(17,17,17);" viewBox="0 0 1024 1024" xmlns:svg="http://www.w3.org/2000/svg" xmlns="http://www.w3.org/2000/svg">', format_custom_result_after_delimiter = '</svg>', format_custom_escaping_rule = 'Raw'
55889 rows in set. Elapsed: 10.101 sec. Processed 384.02 million rows, 8.06 GB (38.02 million rows/s., 798.38 MB/s.)
Peak memory usage: 23.73 MiB.
除了将质心转换为 h3 索引及其聚合之外，此查询还有一些其他有趣的部分
- 我们将经度提取到一个单独的数组中，并且第一个条目也附加到末尾。即
arrayMap(x -> (x.2), arrayConcat(h3ToGeoBoundary(h3), [h3ToGeoBoundary(h3)[1]]))。这使我们能够检查连续的经度点对，并过滤掉那些跨越 180 子午线的点（这些点渲染为线条，破坏了我们的视觉效果），即WHERE arrayExists(i -> (abs((longs[i]) - (longs[i + 1])) > 180), range(1, length(longs) + 1)) = 0 - 我们没有使用 RGB 颜色渐变，而是选择使用 Oklch。这提供了 RGB（和 HSL）的几个优势，在提供更符合预期的结果和更平滑的渐变的同时，也更易于生成，因为我们的调色板只需要更改一个维度（色调）。在上面，我们的调色板从红色（低带宽）到绿色（高带宽），色调在 30 到 144 之间呈平方根渐变。
sqrt在这里为我们提供了更均匀的值分布，如果线性分布，这些值自然会聚集在一起，从而使我们能够区分区域。我们在此处使用第 95 个百分位数的 144 也很有帮助。
光栅化作为替代方案
上面的 h3 索引方法灵活且实际，可用于渲染分辨率高达 7 左右的图像。高于此值的值会生成大多数浏览器无法渲染的文件，例如，7 和 8 会生成文件大小分别为 600MB（勉强可以渲染！）和 1.85GB 的文件。在这些更高的分辨率下，某些地理区域（例如非洲）缺少数据也变得更加明显，因为每个质心的多边形变得更小，即更多空间只是空的。
这种方法的替代方法是光栅化图像，为每个质心生成一个像素而不是一个多边形。当我们指定所需的分辨率（通过墨卡托）函数时，这种方法可能也更灵活，这将直接映射到用于表示我们的质心的像素数。对于相同等效数量的多边形，数据表示在磁盘上也应该更小。
此方法的原理很简单。我们使用我们的墨卡托函数投影每个质心，指定分辨率并对 x 和 y 值进行舍入。这些 x 和 y 对中的每一个都代表一个像素，我们可以使用简单的 GROUP BY 计算该像素的统计信息，例如平均下载速率。以下内容使用此方法计算每个 x,y 像素（分辨率为 2048 x 2048）的平均下载速率，并按行和列顺序返回像素。
SELECT
(y * 2048) + x AS pos,
avg(avg_d_kbps) AS download_speed
FROM ookla
GROUP BY
round(mercator(centroid, 2048, 2048).1) AS x,
round(mercator(centroid, 2048, 2048).2) AS y
ORDER BY pos ASC WITH FILL FROM 0 TO 2048 * 2048
使用上面描述的相同方法，以下查询使用下载速度的最小和最大范围计算颜色渐变，并为每个像素输出 RGBA 值。在下面，我们将 alpha 通道固定为 255，并对我们的颜色渐变使用线性函数，其中第 95 个百分位数为上限（以帮助压缩我们的渐变）。在这种情况下，我们使用查询创建一个视图 ookla_as_pixels，以简化可视化此数据的下一步。请注意，我们使用背景颜色 (17, 17, 17, 255) 来表示那些下载速度的空值的像素，即没有数据点，并确保我们的 RGBA 值是 UInt8。
CREATE OR REPLACE VIEW ookla_as_pixels AS
WITH [234, 36, 19] AS lc, [0, 156, 31] AS uc, [17,17,17] As bg,
pixels AS
(
SELECT (y * 2048) + x AS pos,
avg(avg_d_kbps) AS download_speed
FROM ookla
GROUP BY
round(mercator(centroid, 2048, 2048).1) AS x,
round(mercator(centroid, 2048, 2048).2) AS y
ORDER BY
pos WITH FILL FROM 0 TO 2048 * 2048
),
(
SELECT min(download_speed)
FROM pixels
) AS min_speed,
(
SELECT quantile(0.95)(download_speed)
FROM pixels
) AS max_speed,
pixel_colors AS (
SELECT
least(if(isNull(download_speed),0, (download_speed - min_speed) / (max_speed - min_speed)), 1.0) AS percent,
if(isNull(download_speed),bg[1], toUInt8(lc[1] + (percent * (uc[1] - lc[1])))) AS red,
if(isNull(download_speed),bg[2], toUInt8(lc[2] + (percent * (uc[2] - lc[2])))) AS green,
if(isNull(download_speed),bg[3], toUInt8(lc[3] + (percent * (uc[3] - lc[3])))) AS blue
FROM pixels
) SELECT red::UInt8, green::UInt8, blue::UInt8, 255::UInt8 as alpha FROM pixel_colors
将此数据转换为图像最容易使用 Canvas 元素 和一些简单的 js 来完成。我们有效地为每个像素提供 4 个字节，总共 2048 x 2048 x 4 字节，或大约 16MB。使用 putImageData，我们将此数据直接渲染到画布中。
<!doctype html>
<html>
<head>
<meta charset="utf-8">
<link rel="icon" href="favicon.png">
<title>Simple Ookla Visual</title>
</head>
<body>
<div id="error"></div>
<canvas id="canvas" width="2048" height="2048"></canvas>
<script>
async function render(tile) {
const url = `https://127.0.0.1:8123/?user=default&default_format=RowBinary`;
const response = await fetch(url, { method: 'POST', body: 'SELECT * FROM ookla_as_pixels' });
if (!response.ok) {
const text = await response.text();
let err = document.getElementById('error');
err.textContent = text;
err.style.display = 'block';
return;
}
buf = await response.arrayBuffer();
let ctx = tile.getContext('2d');
let image = ctx.createImageData(2048, 2048);
let arr = new Uint8ClampedArray(buf);
for (let i = 0; i < 2048 * 2048 * 4; ++i) {
image.data[i] = arr[i];
}
ctx.putImageData(image, 0, 0, 0, 0, 2048, 2048);
let err = document.getElementById('error');
err.style.display = 'none';
}
const canvas = document.getElementById("canvas");
render(canvas).then((err) => {
if (err) {
let err = document.getElementById('error');
err.textContent = text;
} else {
err.style.display = 'none'
}
})
</script>
</body>
</html>
此 代码由 Alexey Milovidov 提供，他通过他的 ADS-B 数据交互式可视化和分析 将此技术提升到了一个全新的水平。
以上内容利用了 ClickHouse 的 HTTP 接口 和以 RowBinary 格式 返回结果的能力。视图的创建大大简化了我们的 SQL 查询，只需 SELECT * FROM ookla_as_pixels。我们的结果图像。
即使我们不得不使用比仅仅 SQL 稍微多一点的东西，这也代表了我们 SVG 方法的一个不错的替代方案 :)
结论
这篇文章介绍了 Iceberg 数据格式，并讨论了 ClickHouse 中当前的支持状态，其中包含有关如何查询和导入数据的示例。我们探索了以 Iceberg 格式分发的 Ookla 互联网速度数据集，并借此机会探索了 h3 索引、多边形字典、墨卡托投影、颜色插值和使用 UDF 的质心计算 - 所有这些都旨在仅使用 SQL 可视化数据。
虽然我们主要关注下载速度和固定设备，但我们很乐意看到其他人使用类似的方法来探索其他指标。或者，我们欢迎您就如何改进上述内容提出建议！
