PostGEESE？ DuckDB 空间扩展简介

Max Gabrielsson

2023-04-28 · 16 分钟

地理空间数据在现代应用程序和数据工程工作流程中变得越来越重要和普遍，其用例范围从基于位置的服务到环境监测。

虽然有很多很棒的专业工具可以处理地理空间数据，但将地理空间功能直接集成到 DuckDB 中有多个优点。首先，您可以利用 DuckDB 丰富的类型系统和扩展（如 JSON 和 ICU）来操作、转换和连接地理空间数据以及常规的非结构化或时间序列数据。其次，涉及几何谓词和关系的 spatial 查询可以很好地转换为 SQL，毕竟 SQL 的全部意义在于表达关系！更不用说 DuckDB 提供的所有其他好处，例如事务语义、高性能多线程矢量化执行和大于内存的数据处理。

因此，我们非常高兴地宣布，DuckDB 现在拥有一个 Spatial extension（空间扩展），其中包含许多功能，可以从 DuckDB CLI 和其他 DuckDB 客户端轻松安装。只需执行

INSTALL spatial;
LOAD spatial;

您就可以开始使用了！

不，我们也不会称它为 GeoDuck，那样太恶心了。

它包含什么？

该扩展的核心是基于“Simple Features”（简单特征）几何模型的 GEOMETRY 类型，以及诸如 ST_Area、ST_Intersects 之类的 accompanying 函数。它还提供用于读取和写入地理空间数据格式以及在坐标参考系统之间转换的方法（详细信息请参见后面的文章！）。虽然我们尚未准备好完全遵守 OGC Simple Feature Access 和 SQL/MM 标准，但如果您在其他数据库系统（例如 PostGIS 或 SpatiaLite）中使用过地理空间功能，那么您应该会感到很熟悉。

大多数已实现的函数都基于基础地理空间库的三大支柱：GEOS、GDAL 和 PROJ，它们分别提供算法、格式转换和坐标参考系统转换。特别是，我们利用 GDAL 提供了一组表和复制函数，这些函数支持从 50 多种不同的地理空间数据格式导入和导出表（到目前为止！），包括最常见的格式，例如 Shapefile、GeoJSON、GeoPackage、KML、GML、WKT、WKB 等。

通过运行以下命令自行检查

最初，我们通过包装现有库来优先提供广泛的功能。我们计划将来以原生方式实现更多的核心功能和算法，以实现更快的性能和更高效的内存管理。

作为朝着这个方向迈出的第一步，我们提供了一组非标准化的专业柱状 DuckDB 原生几何类型，例如 POINT_2D、BOX_2D 等，它们应该可以提供更好的压缩和更快的执行速度，以换取一些灵活性，但是这些工作仍然处于实验阶段。

使用示例

以下演示了如何使用 spatial 扩展来读取和导出多种地理空间数据格式、在不同的坐标参考系统之间转换几何图形以及使用 spatial 属性和谓词函数。虽然此示例可能有些牵强，但我们想展示当前可用功能的强大功能。您可以在 spatial 扩展存储库中找到此示例中使用的数据集。

让我们使用 spatial 扩展提供的 ST_Read 表函数，从 Parquet 格式导入 NYC 出租车乘车数据，以及从 shapefile 导入 accompany 的出租车区域数据。这些出租车区域将 NYC 分成代表区域的多边形，例如纽瓦克机场。然后，我们为乘车数据创建一个表，并为区域数据创建一个表。请注意，ST_Read 生成一个带有 wkb_geometry 列的表，该列包含编码为 WKB（Well-Known Binary）blob 的几何图形数据，然后我们使用 ST_GeomFromWKB 函数将其转换为 GEOMETRY 类型。

如果您不熟悉地理空间生态系统，那么这可能看起来有点多，但请放心，这就是您真正需要开始的一切。简而言之：– Shapefile (.shp, .shx, .dbf) 是一种用于存储几何矢量数据和辅助元数据（例如索引和属性）的常见格式。– WKB (Well Known Binary)，虽然本身不是真正的文件格式，但它是矢量几何图形数据的常见二进制编码，用于例如 GeoParquet 中。它有多种形式，但我们现在只关注“标准”WKB。– GEOMETRY 是 DuckDB 类型，表示 Simple Features 几何对象，该对象基于一组将矢量几何图形数据建模为点、线串、多边形或此类集合的标准。这是 spatial 扩展使用的核心数据类型，也是大多数提供的函数采用和返回的数据类型。

INSTALL spatial;
LOAD spatial;

CREATE TABLE rides AS
    SELECT * 
    FROM 'yellow_tripdata_2010-01-limit1mil.parquet';

-- Load the NYC taxi zone data from a shapefile using the gdal-based ST_Read function
CREATE TABLE zones AS
    SELECT zone, LocationId, borough, geom 
    FROM ST_Read('taxi_zones/taxi_zones.shx');

让我们比较行程距离与上车点和下车点之间的线性距离，以确定出租车司机的工作效率如何（或者数据的质量如何，因为一些差异似乎是负数）。我们将坐标从 "WGS84"（由标识符 EPSG:4326 给出），通常也简称为纬度/经度，转换为 "NAD83 / New York Long Island ftUS"（标识为 ESRI:102718）坐标参考系统，这是一个在纽约附近具有最小失真的投影。然后，我们使用 ST_Distance 函数计算距离。在这种情况下，我们得到以英尺为单位的距离，因为我们已将坐标转换为 NAD83，但是我们可以轻松地将其转换为英里（5280 英尺/英里），这是 rides 数据集中使用的单位，因此我们可以正确地比较它们。

距离短于航空距离的行程可能存在错误，因此我们使用此查询来过滤掉一些坏数据。以下查询利用 DuckDB 引用同一 select 语句中定义的列别名的能力。这是 DuckDB 丰富的 SQL 方言如何简化地理空间分析的一个小例子。

CREATE TABLE cleaned_rides AS
    SELECT 
        ST_Point(pickup_latitude, pickup_longitude) AS pickup_point,
        ST_Point(dropoff_latitude, dropoff_longitude) AS dropoff_point,
        dropoff_datetime::TIMESTAMP - pickup_datetime::TIMESTAMP AS time,
        trip_distance,
        ST_Distance(
            ST_Transform(pickup_point, 'EPSG:4326', 'ESRI:102718'), 
            ST_Transform(dropoff_point, 'EPSG:4326', 'ESRI:102718')) / 5280 
            AS aerial_distance, 
        trip_distance - aerial_distance AS diff 
    FROM rides 
    WHERE diff > 0
    ORDER BY diff DESC;

应该注意的是，这并非完全准确，因为我们使用的 ST_Distance 函数没有考虑地球的曲率。但是，我们将接受它作为我们目的的足够好的近似值。球面和测地线距离计算已在路线图中！

现在，让我们将出租车乘车数据与出租车区域数据连接起来，以获取每次乘车的起点和终点区域。我们使用 ST_Within 函数作为我们的连接条件，以检查上车点或下车点是否在出租车区域多边形内。同样，我们需要将坐标从 WGS84 转换为 NAD83，因为出租车区域数据也使用该投影。像这样的 Spatial join 是地理空间数据处理的核心，但是我们目前没有任何优化措施（例如 spatial 索引）来加快这些查询的速度，这就是为什么我们仅将数据的子集用于以下步骤的原因。

-- Since we don't have spatial indexes yet, use a smaller dataset for the join.
DELETE FROM cleaned_rides WHERE rowid > 5000;

CREATE TABLE joined AS 
    SELECT 
        pickup_point,
        dropoff_point,
        start_zone.zone AS start_zone,
        end_zone.zone AS end_zone, 
        trip_distance,
        time,
    FROM cleaned_rides 
    JOIN zones AS start_zone 
      ON ST_Within(ST_Transform(pickup_point, 'EPSG:4326', 'ESRI:102718'), start_zone.geom) 
    JOIN zones AS end_zone 
      ON ST_Within(ST_Transform(dropoff_point, 'EPSG:4326', 'ESRI:102718'), end_zone.geom);

我们可以使用 GDAL copy 函数将 joined 表导出到 GeoJSONSeq 文件，并传入 GDAL 图层创建选项。由于 GeoJSON 每个记录仅支持单个 GEOMETRY，因此我们使用 ST_MakeLine 函数将上车点和下车点合并为一条线几何图形。 GeoJSON 的默认坐标参考系统是 WGS84，但是坐标对必须是经度/纬度，因此我们需要使用 ST_FlipCoordinates 函数来翻转几何图形。

COPY (
    SELECT 
        ST_MakeLine(pickup_point, dropoff_point)
            .ST_FlipCoordinates()
            .ST_AsWKB()
            AS wkb_geometry,
        start_zone,
        end_zone,
        time::VARCHAR AS trip_time 
    FROM joined) 
TO 'joined.geojsonseq' 
WITH (
    FORMAT gdal,
    DRIVER 'GeoJSONSeq',
    LAYER_CREATION_OPTIONS 'WRITE_BBOX=YES'
);

{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:52:00" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.789923, 40.643515 ], [ -73.97608, 40.680395 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:35:00" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.776445, 40.645422 ], [ -73.98427, 40.670782 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:45:42" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.776878, 40.645065 ], [ -73.992153, 40.662571 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:36:00" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.788028, 40.641508 ], [ -73.97584, 40.670927 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:47:58" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.781855, 40.644749 ], [ -73.980129, 40.663663 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:32:10" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.787494, 40.641559 ], [ -73.974694, 40.673479 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:36:59" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.790138, 40.643342 ], [ -73.982721, 40.662379 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:32:00" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.786952, 40.641248 ], [ -73.97421, 40.676237 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:33:21" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.783892, 40.648514 ], [ -73.979283, 40.669721 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:35:45" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.776643, 40.645272 ], [ -73.978873, 40.66723 ] ] } }

现在我们有了它！我们从 Parquet 中提取了表格数据，将其与 shapefile 中的地理空间数据组合，清理并分析了组合后的数据，并将其输出为人类可读的地理空间格式。当前支持的全部函数及其实现状态可以在 此表中的文档中找到。

下一步是什么？

虽然我们可能需要一段时间才能赶上例如 PostGIS 提供的全部功能集，但是我们相信 DuckDB 的矢量化执行模型和柱状存储格式将为地理空间处理带来一类全新的优化，而我们才刚刚开始探索。因此，提高 Spatial join 和谓词的性能是我们的首要任务。

我们的 GEOMETRY 类型也存在一些限制，我们最终希望解决这些限制，例如，我们不支持额外的 Z 和 M 维度，或者不支持 OGC 标准要求的全部几何图形子类型，例如曲线或多面体表面。

我们还有兴趣在不久的将来支持球形和椭球形计算，可能以专用 GEOGRAPHY 类型的形式。

Wasm 构建也即将推出！

请查看 GitHub 存储库，以获取完整的路线图并查看我们当前正在做的事情。如果您想帮助构建此功能，请在 GitHub 上与我们联系！

结论

DuckDB Spatial 扩展是使 DuckDB 成为数据工程和分析的瑞士军刀的又一步。该扩展提供了一种灵活且熟悉的 GEOMETRY 类型，可以在数千个坐标参考系统之间重新投影，并能够导出和导入 50 多个不同数据源之间的地理空间数据。所有这些都嵌入到一个具有最小运行时依赖项的扩展中。这使 DuckDB 可以无缝地融入您现有的 GIS 工作流程中，而不管您使用的是哪种地理空间数据格式或投影。

我们很高兴听到您对 DuckDB spatial 扩展的看法。现在还处于早期阶段，但是我们希望在未来能够分享更多内容，因为我们将继续取得进展！如果您有任何问题、建议、想法或问题，请随时在 Discord 或 GitHub 上与我们联系！