ST_Intersection_Agg not returning intersects

[bartzwemmer]

When testing a function, I ran into the situation where ST_Intersection_Agg returns "POLYGON EMPTY". When I view the data in e.g. QGIS, there is clearly an overlap.

Full code:

import duckdb

DB = duckdb.connect()
DB.install_extension("spatial")
DB.load_extension("spatial")

def main() -> None:
    # Prepare some test data
    """Creates and prints a Polars DataFrame with ID and WKT geometry."""
    DB.execute("""
        CREATE TABLE geometries AS 
        SELECT * FROM ST_Read('data/test_overlap.gml')
        """)
    DB.sql("SELECT * FROM geometries").show()

    # Test for intersections
    intersections = DB.execute("""
                               SELECT 
                                ST_AsText(ST_Intersection_Agg(geometryProperty)) AS overlap
                               FROM geometries
                               WHERE geometryProperty IS NOT NULL
                               """).fetchall()
    
    if intersections[0][0]  == "POLYGON EMPTY":
        print("No intersections found")
        return

    # generate_map(intersections) # do something with the returned overlap here


if __name__ == "__main__":
    main()

And the GML file i'm reading:

<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
     gml:id="aFeatureCollection"
     xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
     xsi:schemaLocation="http://ogr.maptools.org/ test_overlap.xsd"
     xmlns:ogr="http://ogr.maptools.org/"
     xmlns:gml="http://www.opengis.net/gml/3.2">
  <gml:boundedBy><gml:Envelope srsName="urn:ogc:def:crs:EPSG::3857"><gml:lowerCorner>568455.285329327 6817178.14118306</gml:lowerCorner><gml:upperCorner>568665.336942971 6817321.76081046</gml:upperCorner></gml:Envelope></gml:boundedBy>
                                                                                                                 
  <ogr:featureMember>
    <ogr:test_overlap gml:id="test_overlap.0">
      <gml:boundedBy><gml:Envelope srsName="urn:ogc:def:crs:EPSG::3857"><gml:lowerCorner>568455.285329327 6817178.14118306</gml:lowerCorner><gml:upperCorner>568561.576507316 6817278.10550521</gml:upperCorner></gml:Envelope></gml:boundedBy>
      <ogr:geometryProperty><gml:Polygon srsName="urn:ogc:def:crs:EPSG::3857" gml:id="test_overlap.geom.0"><gml:exterior><gml:LinearRing><gml:posList>568455.285329327 6817247.73659721 568521.71731557 6817278.10550521 568561.576507316 6817204.71397755 568500.206005739 6817178.14118306 568455.285329327 6817247.73659721</gml:posList></gml:LinearRing></gml:exterior></gml:Polygon></ogr:geometryProperty>
      <ogr:id>1</ogr:id>
    </ogr:test_overlap>
  </ogr:featureMember>
  <ogr:featureMember>
    <ogr:test_overlap gml:id="test_overlap.1">
      <gml:boundedBy><gml:Envelope srsName="urn:ogc:def:crs:EPSG::3857"><gml:lowerCorner>568551.453537983 6817203.44860639</gml:lowerCorner><gml:upperCorner>568665.336942971 6817316.06664021</gml:upperCorner></gml:Envelope></gml:boundedBy>
      <ogr:geometryProperty><gml:Polygon srsName="urn:ogc:def:crs:EPSG::3857" gml:id="test_overlap.geom.1"><gml:exterior><gml:LinearRing><gml:posList>568551.453537983 6817286.33041779 568624.212380059 6817316.06664021 568665.336942971 6817238.87899905 568601.435699061 6817203.44860639 568551.453537983 6817286.33041779</gml:posList></gml:LinearRing></gml:exterior></gml:Polygon></ogr:geometryProperty>
      <ogr:id>2</ogr:id>
    </ogr:test_overlap>
  </ogr:featureMember>
  <ogr:featureMember>
    <ogr:test_overlap gml:id="test_overlap.2">
      <gml:boundedBy><gml:Envelope srsName="urn:ogc:def:crs:EPSG::3857"><gml:lowerCorner>568490.083036407 6817238.87899905</gml:lowerCorner><gml:upperCorner>568555.249651483 6817321.76081046</gml:upperCorner></gml:Envelope></gml:boundedBy>
      <ogr:geometryProperty><gml:Polygon srsName="urn:ogc:def:crs:EPSG::3857" gml:id="test_overlap.geom.2"><gml:exterior><gml:LinearRing><gml:posList>568490.083036407 6817314.16858346 568504.634804822 6817238.87899905 568555.249651483 6817244.5731693 568528.676856986 6817321.76081046 568490.083036407 6817314.16858346</gml:posList></gml:LinearRing></gml:exterior></gml:Polygon></ogr:geometryProperty>
      <ogr:id>3</ogr:id>
    </ogr:test_overlap>
  </ogr:featureMember>
</ogr:FeatureCollection>

Dependencies:

requires-python = ">=3.14"
dependencies = [
    "duckdb>=1.4.2",
    "folium>=0.20.0",
    "pyarrow>=22.0.0",
]

I tried several data formats (GeoJSON, GML, GeoParquet) and projections, but all with the same result. From the description in the docs, I was expecting a polygon in the shape of the parts of the two features that overlap each other:

Computes the intersection of a set of geometries

[Maxxen]

Hello!
Thanks for opening this issue.
This function is implemented by taking the intersection of each input "one by one", so ordering within the aggregate (which DuckDB allows you to control) matters as intersection is not associative. This could probably be better explained in the docs.

[bartzwemmer]

Thanks for the prompt reply. Regarding the one by one comparison; I have 3 features in my test. It would take feature 1, compare it against 2 and 3. Then take feature 2, compare it against 1 and 3. Finally, it would take feature 3 and compare it against 1 and 2. This can happen in any particular order, without an ORDER BY.

With an updated query like so:

SELECT 
ST_AsText(ST_Intersection_Agg(geometryProperty ORDER BY id)) AS overlap
FROM geometries;

it would go in order of the id column. I tried this with the gml_id (varchar), id (int32) and geometryProperty (geometry) columns, but it still doesn't find the overlap. I always get [('POLYGON EMPTY',)] as a result from the .fetchall(). What am I missing?