Known issues: DuckDB spatial joins and CRS tags on geometry columns

docs/geospatial.md

@@ -316,6 +316,14 @@ and SedonaDB. DuckDB does not consume this encoding — pick per
 destination. The CRS tag defaults to `OGC:CRS84`; pass
 `geometry_crs=...` for anything else.
 
+DuckDB compares geometries only when their CRS tags match. Its `ST_Read` tags
+GeoJSON `EPSG:4326`, so to test the raster's points against polygons read that
+way, retag the polygons with `ST_SetCRS(geom, 'OGC:CRS84')`. GeoJSON's
+coordinates are longitude first, which is what `OGC:CRS84` says; `EPSG:4326`
+officially puts latitude first. And for now, join the geometry column to a *table* of polygons
+through `ST_Point(x, y)` or a copied table; see
+[Known issues](limitations.md).
+
 ## Where the array paradigm still earns its keep
 
 The boundary is **weight generation**. Applying a regridding is a join;

docs/limitations.md

@@ -56,6 +56,33 @@ Pick your engine:
       binds as a native `GEOMETRY` with the CRS attached. See
       [Geospatial in SQL](geospatial.md#geoarrow-point-geometry-columns).
 
+    **A spatial join on the geometry column fails.**
+
+    - *Symptom:* joining a registered table's `geometry` column to a
+      table of polygons (`JOIN regions r ON ST_Within(t.geometry,
+      r.geom)`, or the same polygon through a subquery) raises
+      `NotImplementedException: Pushdown Filter Type EXPRESSION_FILTER
+      is not currently supported in PyArrow Scans`. DuckDB's spatial
+      join pushes a filter on the column into its Arrow scan, which
+      can't take it. A constant polygon works.
+    - *What to do:* build the point in the join,
+      `ST_Within(ST_SetCRS(ST_Point(t.x, t.y), 'OGC:CRS84'), r.geom)`,
+      or copy what you need into DuckDB first (`CREATE TABLE ... AS`)
+      and join the copy.
+
+    **Geometries with different CRS tags don't compare.**
+
+    - *Symptom:* `ST_Within` raises "Cannot call function 'ST_Within'
+      with geometries of different coordinate reference systems" when
+      the polygons come from `ST_Read`, which tags GeoJSON
+      `EPSG:4326`, while the registered points carry the default
+      `OGC:CRS84`.
+    - *What to do:* retag the polygons, `ST_SetCRS(geom, 'OGC:CRS84')`.
+      Both tags are WGS84, but `EPSG:4326` officially puts latitude
+      first, while GeoJSON and the registered points put longitude
+      first, as `OGC:CRS84` says. (Tagging the points `EPSG:4326`
+      instead would mislead a later `ST_Transform`.)
+
 === "Polars"
 
     **Float `is_in` literals lose precision.**

tests/test_geometry.py

@@ -112,6 +112,64 @@ def test_duckdb_st_within_on_geometry_column(grid, spatial_con):
     assert got == (expected.size, round(float(expected.mean()), 3))
 
 
+_POLYGON = (
+    "POLYGON ((-58.1 -28.75, -56.9 -28.75, -56.9 -27.9, "
+    "-58.1 -27.9, -58.1 -28.75))"
+)
+
+
+def test_duckdb_spatial_join_needs_a_point_or_a_table(grid, spatial_con):
+    """Joining the geometry column to a table of polygons fails inside
+    DuckDB's Arrow scan, which can't take the filter the spatial join
+    pushes into it. Building the point in the join, or copying the table
+    into DuckDB first, works. When DuckDB fixes this, the first assertion
+    fails: drop the workarounds from the docs then."""
+    duckdb = pytest.importorskip("duckdb")
+    con, _ = spatial_con
+    con.execute(
+        "CREATE TABLE regions AS "
+        f"SELECT ST_SetCRS(ST_GeomFromText('{_POLYGON}'), 'OGC:CRS84') AS geom"
+    )
+    want = int(((grid.y >= -28.75) & (grid.y <= -27.9)).sum()) * grid.sizes["x"]
+
+    with pytest.raises(duckdb.NotImplementedException, match="PyArrow"):
+        con.execute(
+            "SELECT count(*) FROM t JOIN regions r ON ST_Within(t.geometry, r.geom)"
+        ).fetchone()
+    got = con.execute(
+        "SELECT count(*) FROM t JOIN regions r ON "
+        "ST_Within(ST_SetCRS(ST_Point(t.x, t.y), 'OGC:CRS84'), r.geom)"
+    ).fetchone()
+    assert got == (want,)
+    con.execute("CREATE TABLE copied AS FROM t")
+    got = con.execute(
+        "SELECT count(*) FROM copied c JOIN regions r ON ST_Within(c.geometry, r.geom)"
+    ).fetchone()
+    assert got == (want,)
+
+
+def test_duckdb_compares_only_matching_crs_tags(grid):
+    """DuckDB's ST_Read tags GeoJSON as EPSG:4326, and ST_Within refuses
+    to compare it with the points' OGC:CRS84 tag. Retagging the polygon
+    OGC:CRS84, the longitude-first order GeoJSON uses, makes them compare."""
+    duckdb = pytest.importorskip("duckdb")
+    con = duckdb.connect()
+    try:
+        con.execute("INSTALL spatial; LOAD spatial;")
+    except duckdb.Error:
+        pytest.skip("duckdb spatial extension unavailable")
+    xql.register(con, "t", grid, geometry=("x", "y"))
+    polygon = f"ST_SetCRS(ST_GeomFromText('{_POLYGON}'), 'EPSG:4326')"
+    query = "SELECT count(*) FROM t WHERE ST_Within(geometry, {})"
+
+    with pytest.raises(duckdb.BinderException, match="coordinate reference"):
+        con.execute(query.format(polygon)).fetchone()
+    retagged = f"ST_SetCRS({polygon}, 'OGC:CRS84')"
+    inside = (grid.y >= -28.75) & (grid.y <= -27.9)
+    want = int(inside.sum()) * grid.sizes["x"]
+    assert con.execute(query.format(retagged)).fetchone() == (want,)
+
+
 def test_geopandas_consumes_native_points(grid):
     gpd = pytest.importorskip("geopandas")