diff --git a/docs/geospatial.md b/docs/geospatial.md index ca9b4663..86b7a9cb 100644 --- a/docs/geospatial.md +++ b/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; diff --git a/docs/limitations.md b/docs/limitations.md index 50dc16ba..f2607028 100644 --- a/docs/limitations.md +++ b/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.** diff --git a/tests/test_geometry.py b/tests/test_geometry.py index 5e3238d7..6317a42c 100644 --- a/tests/test_geometry.py +++ b/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")