Update to DuckDB V1.2.2 + Add Spatial Join Operator

[Maxxen]

This PR updates the workflow and pinned DuckDB version in the v1.2.2 branch to DuckDB v1.2.2. It also adds a new physical operator, the SPATIAL_JOIN.

Spatial Joins

Executing spatial joins in DuckDB did not scale well, as the only method available for the query engine to execute them was through the "blockwise nested-loop join" operator. The blockwise-nl-join operator uses the simplest possible strategy to join two datasets, namely by comparing each item on the left side with the each item on the right side, evaluating the join condition for each pair. This means that joining a dataset of size n with another of size m requires m*n comparisons. In spatial, where the join conditions tend to be pretty costly to evaluate (e.g. polygon-polygon intersection) the effect is that a lot of spatial join queries just hit the wall as the dataset sizes grow, despite DuckDB generally being pretty good at brute-forcing through its parallelized vectorized execution.

Enter, the SPATIAL_JOIN operator!
Compared to the BLOCKWISE_NL_JOIN, the spatial join operator creates a temporary spatial index on-the-fly on the "build side" of the join which it then uses to quickly filter out only those rows whose bounding boxes intersect with the probe-side value and thus "might" be a match. Bounding box intersection checks are comparatively cheap, and so this is a lot more efficient than having to evaluate the actual join predicate on every pairing.

The spatial join operator is instantiated any time there is a INNER/LEFT/RIGHT/OUTER join using one of the following predicates as the join condition:

• ST_Equals
• ST_Intersects
• ST_Touches
• ST_Crosses
• ST_Within
• ST_Contains
• ST_Overlaps
• ST_Covers
• ST_CoveredBy
• ST_ContainsPropery

Limitations:

There are currently some limitations, but Im very much interested in working on resolving these in the future.

• The spatial join operator only supports a single join condition for now
• SEMI/ANTI joins are not yet supported either
• The "build side" of the join needs to be able to fit in memory. However, all memory used to construct the temporary index is tracked by DuckDB, and should respect any set memory limit.
  • This might sound limiting, but as it stands executing spatial joins before this PR tend to hit the combinatorial explosion performance wall way before the build side grows large enough to cause memory issues, so this initial spatial join operator will still allow you to execute larger joins that you probably couldn't before.

Example

Here's a quick example/benchmark:

CREATE TABLE lhs AS
SELECT
    ST_Point(x, y) as geom,
    (y * 50) + x // 10 as id
FROM
    generate_series(0, 1000, 5) r1(x),
    generate_series(0, 1000, 5) r2(y);

-- 40401 rows

CREATE TABLE rhs AS
SELECT
    ST_Buffer(ST_Point(x, y), 5) as geom,
    (y * 50) + x // 10 as id
FROM
    generate_series(0, 500, 10) r1(x),
    generate_series(0, 500, 10) r2(y);

-- 2601 rows

EXPLAIN SELECT * FROM rhs JOIN lhs ON st_intersects(lhs.geom, rhs.geom);
----
┌─────────────────────────────┐
│┌───────────────────────────┐│
││       Physical Plan       ││
│└───────────────────────────┘│
└─────────────────────────────┘
┌───────────────────────────┐
│         PROJECTION        │
│    ────────────────────   │
│            geom           │
│             id            │
│            geom           │
│             id            │
│                           │
│        ~40401 Rows        │
└─────────────┬─────────────┘
┌─────────────┴─────────────┐
│        SPATIAL_JOIN       │
│    ────────────────────   │
│      Join Type: INNER     │
│                           │
│        Conditions:        ├──────────────┐
│ ST_Intersects(geom, geom) │              │
│                           │              │
│        ~40401 Rows        │              │
└─────────────┬─────────────┘              │
┌─────────────┴─────────────┐┌─────────────┴─────────────┐
│         SEQ_SCAN          ││         SEQ_SCAN          │
│    ────────────────────   ││    ────────────────────   │
│         Table: lhs        ││         Table: rhs        │
│   Type: Sequential Scan   ││   Type: Sequential Scan   │
│                           ││                           │
│        Projections:       ││        Projections:       │
│            geom           ││            geom           │
│             id            ││             id            │
│                           ││                           │
│        ~40401 Rows        ││         ~2601 Rows        │
└───────────────────────────┘└───────────────────────────┘

Executing this in the CLI with .timer on gives me :

Run Time (s): real 0.150 user 0.030077 sys 0.001769 on my M3 MacBook Pro.

If we now disable the optimizer with pragma disabled_optimizers = 'extension'; and rerun, we get:

Run Time (s): real 10.470 user 20.675175 sys 0.011965

Now, this is a very contrived example because the geometries are very uniform and half of them trivially match, but it still illustrates the performance improvements.

Other changes

• ST_Extent_Agg is now precise instead of relying on cached bounds
• Fixes a bug causing ST_Translate not to work properly for xy geometries.

[Copilot]

Copilot reviewed 25 out of 29 changed files in this pull request and generated 1 comment.

Files not reviewed (4)

• Makefile: Language not supported
• duckdb: Language not supported
• src/spatial/CMakeLists.txt: Language not supported
• src/spatial/operators/CMakeLists.txt: Language not supported

Comments suppressed due to low confidence (11)

src/spatial/spatial_optimizers.cpp:1

• The removal of spatial_optimizers.cpp is significant; ensure that the optimization logic it contained is fully integrated into the new SpatialJoinOptimizer.

-#include "spatial/spatial_types.hpp"

src/spatial/spatial_extension.cpp:16

• [nitpick] The include directive was updated to reference the new spatial join optimizer; please verify that all dependent modules have been updated accordingly.

#include "spatial/operators/spatial_join_optimizer.hpp"

src/spatial/modules/main/spatial_functions_scalar.cpp:2664

• Changing the bounding box type from double to float in ST_Extent_Approx may reduce precision; please verify that this trade-off is acceptable.

Box2D<float> bbox;

src/spatial/modules/main/spatial_functions_scalar.cpp:5196

• Ensure that sgl::util::hilbert_encode provides equivalent results and performance compared to the previous local HilbertEncode implementation.

const auto h = sgl::util::hilbert_encode(16, hilbert_x, hilbert_y);

src/spatial/modules/main/spatial_functions_aggregate.cpp:145

• Switching the aggregate function input type from geometry_t to string_t should be validated to ensure compatibility with downstream spatial functions.

const auto agg = AggregateFunction::UnaryAggregate<ExtentAggState, string_t, string_t, ExtentAggFunction>(

src/spatial/modules/geos/geos_module.cpp:286

• [nitpick] The stray semicolon was removed from ST_Buffer; this clean-up improves readability and reduces potential confusion.

return lstate.Serialize(result, buffer);

src/spatial/index/rtree/rtree_index_plan_scan.cpp:84

• Verify that extracting the bounding box directly as a float without explicit MathUtil conversion maintains the expected precision for index scans.

static bool TryGetBoundingBox(const Value &value, Box2D<float> &bbox) {

src/spatial/index/rtree/rtree_index.cpp:170

• Transitioning from double to float in bounding box extraction should be reviewed; ensure any precision loss does not affect the spatial index accuracy.

if (!geom_data[i].TryGetCachedBounds(bbox)) {

src/spatial/geometry/geometry_type.hpp:102

• The conversion from double to float in TryGetCachedBounds may introduce precision loss; please confirm that this trade-off is acceptable for spatial operations.

bool TryGetCachedBounds(Box2D<float> &bbox) const {

src/sgl/sgl.hpp:1040

• Changing the default values for 'old_vertex' from {0, 0, 0, 0} to {0, 0, 1, 1} in affine_transform may affect transformation results; please verify that these new defaults are intended.

vertex_xyzm old_vertex = {0, 0, 1, 1};

.github/workflows/StableDistributionPipeline.yml:14

• CI workflow configuration was updated to use DuckDB v1.2.2; ensure all pipeline steps and branch triggers are compatible with the new version.

- v1.2.2

[aborruso]

Hi @Maxxen I have duckdb 1.2.2 and duckdb spatial 3bb37f8.

How to use this great new feature?

Thank you

[brpy]

Great feature @Maxxen, much awaited!

This significantly improves our workloads. Some of our workloads which were not possible earlier are now getting executed in just a few minutes.

Thanks a lot!