From b819f1ed5c50abd6c7b16d8e051f749f32457e49 Mon Sep 17 00:00:00 2001 From: Alex Merose Date: Wed, 21 Oct 2026 10:58:25 -0700 Subject: [PATCH 1/3] Add geospatial case 10: rasters, points and polygons, exactly and with H3 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A new case in the geospatial suite joins three shapes of data in DuckDB: the latest NOAA GFS forecast (a raster, read lazily from dynamical.org's public Zarr store), Natural Earth's cities (points) and its 1:10m countries (polygons). It asks how hot each country will be tomorrow and how many people live where it will pass 35 °C, two ways: - exact geometry with DuckDB's spatial extension: ST_Within on the geometry column xarray-sql derives, with bbox_conjuncts pruning a single country's read to its chunks; - H3 hexagon IDs with DuckDB's community h3 extension, which turn each spatial relationship into an equality JOIN. Every step is checked against xarray, shapely or h3-py references. On the 2026-10-21 forecast, the exact country means took 60 s and the H3 ones 13 ms after a 13 s one-time cover, within 0.29 °C for half of the countries. Writing it surfaced two DuckDB behaviors, now in Known issues and pinned by tests: a spatial join on the geometry column of an Arrow scan fails (see #278), and ST_Within won't compare ST_Read's EPSG:4326 tag with the default OGC:CRS84 one. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NyN7idACeDMh94qrnFMu6h --- benchmarks/geospatial/10_h3_join.py | 394 ++++++++++++++++++++++++++++ benchmarks/geospatial/README.md | 11 + docs/geospatial.md | 74 +++++- docs/limitations.md | 24 ++ tests/test_geometry.py | 61 +++++ 5 files changed, 562 insertions(+), 2 deletions(-) create mode 100644 benchmarks/geospatial/10_h3_join.py diff --git a/benchmarks/geospatial/10_h3_join.py b/benchmarks/geospatial/10_h3_join.py new file mode 100644 index 00000000..68402c42 --- /dev/null +++ b/benchmarks/geospatial/10_h3_join.py @@ -0,0 +1,394 @@ +# /// script +# requires-python = ">=3.11" +# dependencies = [ +# "xarray-sql[duckdb]", +# "xarray", +# "numpy", +# "pandas", +# "zarr>=3", +# "fsspec", +# "aiohttp", +# "h3>=4", +# "shapely>=2", +# ] +# +# [tool.uv.sources] +# xarray-sql = { path = "../../", editable = true } +# /// +"""Rasters, points and polygons in one query: exact geometry, then H3. + +Geospatial data comes in different shapes: a forecast is a raster, cities +are points, countries are polygons. Joining them answers everyday +questions, like "how hot will each country be tomorrow?" and "how many +people live where it will pass 35 °C?". This case asks both in DuckDB, an +embedded SQL database, two ways: + +1. **Exact geometry**, with DuckDB's ``spatial`` extension. + ``xql.register(..., geometry=("longitude", "latitude"))`` gives the + raster a point ``geometry`` column, so ``ST_Within(geometry, country)`` + is a point-in-polygon test. A single country's query pairs it with + ``xql.bbox_conjuncts``, so the lazy raster reads only the chunks over + that country. +2. **H3**, Uber's hierarchical grid of hexagons, with DuckDB's community + ``h3`` extension. Every shape becomes hexagon IDs: a grid cell or a city + is the hexagon containing it, and a country is the hexagons that cover + it. Each spatial relationship is then an equality ``JOIN`` on the ID, + which any SQL engine runs as a hash join. + +Exact geometry is exact, and slow on detailed outlines: every grid point +is tested against a country's thousands of vertices. H3 pays that cost +once, to cover the countries with hexagons, and answers every query after +that in milliseconds, within a small error. The case checks each step +against xarray, shapely and h3-py references, and reports how far H3's +answers are from the exact ones. + +Data: + +* the latest NOAA GFS forecast (0.25°, global), from dynamical.org's public + Zarr store; the first 24 hours of 2 m temperature, read lazily; +* Natural Earth's countries (1:10m) and populated places (7,342 cities with + their metropolitan populations), read by DuckDB straight from GeoJSON on + GitHub. + +The forecast changes every six hours, so the numbers do too. +""" + +from __future__ import annotations + +import duckdb +import h3 +import numpy as np +import pandas as pd +import shapely +import xarray as xr + +import xarray_sql as xql +from _harness import ( + CaseSkipped, + assert_grid_close, + measured, + run_case, + show_sql, + timed, +) + +_GFS = "https://data.dynamical.org/noaa/gfs/forecast/latest.zarr" +_NATURAL_EARTH = ( + "/vsicurl/https://raw.githubusercontent.com/nvkelso/natural-earth-vector/" + "master/geojson/" +) +_HOT = 35.0 # °C +# H3 resolution 4 hexagons average 1,770 km², about two GFS grid cells at +# the equator, so every hexagon holds at least one grid cell's center. +_RES = 4 +# Countries whose exact means are checked against shapely: large and small, +# tropical and polar, one-piece and many islands. +_CHECKED = ["India", "Brazil", "Norway", "Indonesia", "Egypt"] + + +def _connect() -> duckdb.DuckDBPyConnection: + con = duckdb.connect() + try: + con.execute("INSTALL spatial; LOAD spatial") + con.execute("INSTALL h3 FROM community; LOAD h3") + con.execute(f""" + CREATE TABLE countries AS + SELECT NAME AS name, geom + FROM ST_Read('{_NATURAL_EARTH}ne_10m_admin_0_countries.geojson'); + CREATE TABLE cities AS + SELECT name, pop_max AS population, latitude, longitude, geom + FROM ST_Read('{_NATURAL_EARTH}ne_10m_populated_places_simple.geojson'); + """) + except duckdb.Error as exc: + raise CaseSkipped( + f"DuckDB extensions or Natural Earth ({exc})" + ) from exc + return con + + +def _outline(con: duckdb.DuckDBPyConnection, name: str): + """A country's outline as WKT, and its bounds (xmin, ymin, xmax, ymax).""" + wkt, *bounds = con.execute( + "SELECT ST_AsText(geom), ST_XMin(geom), ST_YMin(geom)," + " ST_XMax(geom), ST_YMax(geom) FROM countries WHERE name = ?", + [name], + ).fetchone() + return wkt, bounds + + +def _inside(field: xr.DataArray, wkt: str, bounds) -> xr.DataArray: + """The reference: ``field`` inside the outline, NaN elsewhere, cropped + to the rows and columns with any cell inside.""" + xmin, ymin, xmax, ymax = bounds + box = field.sel(latitude=slice(ymax, ymin), longitude=slice(xmin, xmax)) + lon, lat = np.meshgrid(box.longitude, box.latitude) + inside = shapely.contains_xy(shapely.from_wkt(wkt), lon, lat) + return ( + box.where(inside) + .dropna("latitude", how="all") + .dropna("longitude", how="all") + ) + + +def main() -> None: + try: + gfs = xr.open_zarr(_GFS, chunks=None) + except Exception as exc: # noqa: BLE001 — any failure → skip, not crash + raise CaseSkipped(f"GFS unavailable ({exc})") from exc + # The latest forecast's first 24 hours: (lead_time, latitude, longitude). + gfs = gfs[["temperature_2m"]].isel(init_time=-1, lead_time=slice(0, 25)) + print( + f" raster: GFS forecast from {gfs.init_time.values}, " + f"{gfs.sizes['lead_time']} hours × {gfs.sizes['latitude']}" + f"×{gfs.sizes['longitude']}" + ) + + with timed("load DuckDB extensions + Natural Earth"): + con = _connect() + n_countries, n_cities = con.execute( + "SELECT (SELECT count(*) FROM countries), (SELECT count(*) FROM cities)" + ).fetchone() + print(f" vector: {n_countries} countries (polygons), {n_cities:,} cities") + + # DuckDB's ST_Read tags GeoJSON's longitude/latitude as EPSG:4326, and + # ST_Within refuses to compare geometries with different tags, so tag + # the raster's points the same way. (Both name WGS84 longitude and + # latitude; xarray-sql's default tag is OGC:CRS84.) + xql.register( + con, + "gfs", + gfs, + geometry=("longitude", "latitude"), + geometry_crs="EPSG:4326", + ) + + # ── 1. One country, exact geometry, read lazily ────────────────────── + name = "India" + wkt, bounds = _outline(con, name) + # ST_Within alone can't skip chunks: the engine doesn't push functions + # into the scan. The bounding box's BETWEENs can, so the query reads + # only the chunks over the country, and ST_Within refines the box. + sql = f""" + SELECT latitude, longitude, MAX(temperature_2m) AS t_max + FROM gfs + WHERE {xql.bbox_conjuncts(bounds, x="longitude", y="latitude")} + AND ST_Within(geometry, + ST_SetCRS(ST_GeomFromText('{wkt}'), 'EPSG:4326')) + GROUP BY latitude, longitude + ORDER BY latitude DESC, longitude + """ + show_sql(sql.replace(wkt, f"<{name}'s outline>")) + for _ in measured(f"SQL: {name}'s 24-hour high, ST_Within + bbox pruning"): + country = xql.to_dataset(con.sql(sql), template=gfs).t_max + for _ in measured("shapely + xarray reference"): + ref = _inside(gfs.temperature_2m.max("lead_time"), wkt, bounds) + assert_grid_close(f"{name}'s 24-hour high (°C)", country, ref) + print( + f" {name}: {int(country.count()):,} grid cells, mean high " + f"{float(country.mean()):.1f} °C" + ) + + # ── 2. The 24-hour high everywhere, as a table ─────────────────────── + # Every later query uses each grid cell's high, so compute it once into + # a DuckDB table. It keeps xarray-sql's point geometry for the + # point-in-polygon joins below. + with timed("SQL: 24-hour high for every grid cell (reads the globe)"): + con.execute(""" + CREATE TABLE high AS + SELECT latitude, longitude, MAX(temperature_2m) AS t_max, + ANY_VALUE(geometry) AS geometry + FROM gfs GROUP BY latitude, longitude + """) + with timed("xarray reference: gfs.max('lead_time')"): + ref_high = gfs.temperature_2m.max("lead_time").load() + high = xql.to_dataset( + con.sql( + "SELECT latitude, longitude, t_max FROM high" + " ORDER BY latitude DESC, longitude" + ), + template=gfs, + ).t_max + assert_grid_close("24-hour high (°C)", high, ref_high) + + # ── 3. Exact geometry: point-in-polygon joins ──────────────────────── + # A grid cell's area shrinks with the cosine of its latitude, so weight + # each cell by it: an unweighted mean would overweight a country's + # high latitudes, where the grid's points crowd together. + exact_means_sql = """ + SELECT k.name AS country, + SUM(h.t_max * cos(radians(h.latitude))) + / SUM(cos(radians(h.latitude))) AS t_max + FROM high AS h JOIN countries AS k ON ST_Within(h.geometry, k.geom) + GROUP BY k.name + """ + show_sql(exact_means_sql, label="SQL: country mean highs, exact") + for _ in measured( + f"SQL: country mean highs, ST_Within ({high.size:,} points × " + f"{n_countries} countries)" + ): + exact_means = con.sql(exact_means_sql).df().set_index("country").t_max + for name in _CHECKED: + cells = _inside(ref_high, *_outline(con, name)) + weight = np.cos(np.radians(cells.latitude)) + want = float(cells.weighted(weight.fillna(0)).mean()) + np.testing.assert_allclose(exact_means[name], want, rtol=1e-6) + print(f" ✅ PASS exact means match shapely for {', '.join(_CHECKED)}") + + # A city's forecast is its nearest grid cell's. On a regular grid that + # needs no geometry: round the city's coordinates to the grid's 0.25° + # step (and wrap 180° to -180°), then join on equality. + exact_people_sql = f""" + WITH snapped AS ( + SELECT name, population, geom, + ROUND(latitude * 4) / 4 AS lat, + ((ROUND(longitude * 4) / 4 + 180) % 360) - 180 AS lon + FROM cities + ) + SELECT k.name AS country, SUM(c.population) AS people + FROM snapped AS c + JOIN high AS h ON h.latitude = c.lat AND h.longitude = c.lon + JOIN countries AS k ON ST_Within(c.geom, k.geom) + WHERE h.t_max > {_HOT} + GROUP BY k.name + """ + show_sql(exact_people_sql, label="SQL: people above 35 °C, exact") + for _ in measured("SQL: people above 35 °C, nearest cell + ST_Within"): + exact_people = con.sql(exact_people_sql).df().set_index("country") + + # Check the snapping join: each city's high, against indexing the + # reference grid at the same rounded coordinates. + cities = con.sql(""" + SELECT c.latitude, c.longitude, h.t_max + FROM cities AS c JOIN high AS h + ON h.latitude = ROUND(c.latitude * 4) / 4 + AND h.longitude = ((ROUND(c.longitude * 4) / 4 + 180) % 360) - 180 + """).df() + + def round_like_sql(v): # halves away from zero, as SQL's ROUND does + return np.sign(v) * np.floor(np.abs(v) + 0.5) + + lat = round_like_sql(cities.latitude.to_numpy() * 4) / 4 + lon = ( + round_like_sql(cities.longitude.to_numpy() * 4) / 4 + 180 + ) % 360 - 180 + want = ref_high.sel( + latitude=xr.DataArray(lat), longitude=xr.DataArray(lon) + ).values + np.testing.assert_allclose(cities.t_max, want, rtol=1e-6) + print(f" ✅ PASS each of {len(cities):,} cities' high is its grid cell's") + + # ── 4. H3: every shape becomes hexagon IDs ─────────────────────────── + # A hexagon belongs to the country that contains its center. A coastal + # hexagon whose center is at sea goes to a country it overlaps, so + # coastal cities aren't lost; each hexagon still has one country. + # ST_Dump splits a country of several islands into its polygons. + with timed(f"SQL: cover the countries with H3 resolution {_RES} hexagons"): + con.execute(f""" + CREATE TABLE country_cells AS + WITH parts AS ( + SELECT name, UNNEST(ST_Dump(geom)).geom AS geom FROM countries + ), + centered AS ( + SELECT name, + UNNEST(h3_polygon_wkb_to_cells(ST_AsWKB(geom), {_RES})) + AS cell + FROM parts + ), + touching AS ( + SELECT name, + UNNEST(h3_polygon_wkt_to_cells_experimental( + ST_AsText(geom), {_RES}, 'overlap')) AS cell + FROM parts + ) + SELECT name, cell FROM centered + UNION ALL + SELECT MIN(name), cell FROM touching + WHERE cell NOT IN (SELECT cell FROM centered) + GROUP BY cell + """) + with timed("SQL: each grid cell's hexagon, and each hexagon's mean high"): + con.execute(f""" + CREATE TABLE hex_high AS + SELECT h3_latlng_to_cell(latitude, longitude, {_RES}) AS cell, + AVG(t_max) AS t_max + FROM high GROUP BY cell + """) + + # Hexagons have nearly equal areas, so a plain average of a country's + # hexagons is already weighted by area. + h3_means_sql = """ + SELECT k.name AS country, AVG(h.t_max) AS t_max + FROM country_cells AS k JOIN hex_high AS h ON h.cell = k.cell + GROUP BY k.name + """ + show_sql(h3_means_sql, label="SQL: country mean highs, H3") + for _ in measured("SQL: country mean highs, H3 equality JOIN"): + h3_means = con.sql(h3_means_sql).df().set_index("country").t_max + + h3_people_sql = f""" + SELECT k.name AS country, SUM(c.population) AS people + FROM cities AS c + JOIN hex_high AS h + ON h.cell = h3_latlng_to_cell(c.latitude, c.longitude, {_RES}) + JOIN country_cells AS k ON k.cell = h.cell + WHERE h.t_max > {_HOT} + GROUP BY k.name + """ + show_sql(h3_people_sql, label="SQL: people above 35 °C, H3") + for _ in measured("SQL: people above 35 °C, H3 equality JOINs"): + h3_people = con.sql(h3_people_sql).df().set_index("country") + + # Check the hexagons' means against h3-py, H3's Python bindings, + # applied to the reference grid. + with timed("h3-py reference: index every grid cell"): + grid = ref_high.to_dataframe().reset_index() + grid["cell"] = [ + h3.str_to_int(h3.latlng_to_cell(y, x, _RES)) + for y, x in zip(grid.latitude, grid.longitude) + ] + want_hex = grid.groupby("cell").temperature_2m.mean() + got_hex = con.sql("SELECT cell, t_max FROM hex_high").df() + got_hex = got_hex.set_index(got_hex.cell.astype("uint64")).t_max + assert set(got_hex.index) == set(want_hex.index), "different hexagons" + np.testing.assert_allclose(got_hex[want_hex.index], want_hex, rtol=1e-6) + print(f" ✅ PASS {len(want_hex):,} hexagons' mean highs match h3-py") + + # ── How close H3 comes to exact geometry ───────────────────────────── + both = pd.concat({"exact": exact_means, "h3": h3_means}, axis=1).dropna() + diff = (both.h3 - both.exact).abs().sort_values() + print( + f"\n Country mean highs, {len(diff)} countries: H3 is within " + f"{diff.median():.2f} °C of exact for half of them, " + f"{diff.quantile(0.9):.2f} °C for 90%." + ) + print(" Largest differences (small or ragged countries, few hexagons):") + for country in diff.index[-3:][::-1]: + print( + f" {country:<24} exact {both.exact[country]:6.2f} °C," + f" H3 {both.h3[country]:6.2f} °C" + ) + + people = pd.concat( + {"exact": exact_people.people, "h3": h3_people.people}, axis=1 + ).fillna(0) + people = people.sort_values("exact", ascending=False) + print(f"\n People in cities above {_HOT:.0f} °C, top 10 countries:") + print(f"\n {'Country':<24} {'exact':>12} {'H3':>12}") + for country, row in people.head(10).iterrows(): + print(f" {country:<24} {int(row.exact):>12,} {int(row.h3):>12,}") + print( + f" {'Total':<24} {int(people.exact.sum()):>12,}" + f" {int(people.h3.sum()):>12,}" + ) + print( + "\n H3 counts a city as hot when its hexagon's mean high is, and" + " exact geometry\n when its nearest grid cell's is. A coastal" + " city's hexagon averages land with\n cooler sea, so H3 counts" + " fewer hot coastal cities." + ) + + +if __name__ == "__main__": + raise SystemExit( + run_case(main, "Rasters × points × polygons: exact geometry and H3") + ) diff --git a/benchmarks/geospatial/README.md b/benchmarks/geospatial/README.md index c046adc3..80cda53f 100644 --- a/benchmarks/geospatial/README.md +++ b/benchmarks/geospatial/README.md @@ -24,6 +24,7 @@ plain-English definition of the operation, and computes the same numbers. | 07 | `07_reproject_udf.py` | per-pixel CRS transform | scalar **UDF** (`reproject()` from the geo extension), à la PostGIS `ST_Transform` | | 08 | `08_regrid_weights.py` | interpolation to a new grid | sparse-weight table `JOIN` + weighted `GROUP BY` | | 09 | `09_warp.py` | reproject **and** resample (warp) | reproject **UDF** (07) → weight table `JOIN` (08) | +| 10 | `10_h3_join.py` | point-in-polygon per pixel, city and country | `ST_Within` `JOIN`s (DuckDB `spatial`), or H3 cell IDs and equality `JOIN`s (DuckDB `h3`) | Cases 01–06 show operations that are *natively* relational. Cases 07–08 are the "hardest" array operations — reprojection and regridding — and show where a UDF @@ -31,6 +32,9 @@ fits (a per-row coordinate transform) versus where the operation is really a sparse matrix multiply expressed as a `JOIN`. Case 09 composes the two into a full **warp** (GDAL/rasterio `reproject`): the 07 UDF reprojects the target grid, arrays turn the reprojected points into bilinear weights, and the 08 `JOIN` applies them. +Case 10 joins a raster, points and polygons in DuckDB, first with exact +geometry (`ST_Within`) and then through [H3](https://h3geo.org/) hexagon IDs, +which turn each spatial relationship into an equality `JOIN`. See [`docs/geospatial.md`](../../docs/geospatial.md) for the full narrative, including *where the array paradigm still earns its keep* (generating the @@ -60,6 +64,13 @@ interpolation weights — the geometry — which SQL applies but does not comput validates against xarray's `.interp()` at the reprojected points, with Earth Engine's own lon/lat SRTM as a second, cross-CRS check. All three run against Earth Engine using your existing `gcloud` login, and skip cleanly without it. +- **10 H3 join** — the latest **NOAA GFS** forecast (0.25° global) from + [dynamical.org](https://dynamical.org)'s public Zarr store, and + [Natural Earth](https://www.naturalearthdata.com/)'s 1:10m countries and + populated places, which DuckDB reads straight from GeoJSON on GitHub. Needs + DuckDB's `spatial` extension and its community `h3` extension (installed on + first run), and skips cleanly offline. About three minutes, most of it the + exact point-in-polygon join and the h3-py reference. ## Running diff --git a/docs/geospatial.md b/docs/geospatial.md index ca9b4663..04b2e491 100644 --- a/docs/geospatial.md +++ b/docs/geospatial.md @@ -44,7 +44,7 @@ actually runs — and this suite works through nearly all of it: | Transformed Eulerian Mean (circulation diagnostics — zonal means and anomalies) | cases 03, 04 | | Forecast evaluation (scoring forecasts against ground truth) | case 05 | | Regridding and reprojection (resolution and CRS changes) | cases 07, 08, 09 | -| Spatial joins (large polygon-to-polygon joins) | *not covered* — a vector-data problem; the closest analogue here is the raster × vector join in case 06 | +| Spatial joins (large polygon-to-polygon joins) | partly: case 10 joins a raster, points and polygons, exactly and through H3 cells; polygon-to-polygon joins are a vector-data problem and not covered | So the claim isn't that a few cherry-picked operations happen to be relational. It's that an independent survey of the operations geoscience runs at scale, run @@ -63,6 +63,7 @@ through SQL one by one, turns out to be — almost entirely — queries. | Reprojection | per-pixel CRS transform | scalar **UDF** (`ST_Transform`-style) | [`07_reproject_udf.py`](https://github.com/xqlsystems/xarray-sql/blob/main/benchmarks/geospatial/07_reproject_udf.py) | | Regridding | interpolation to a new grid | sparse-weight table `JOIN` | [`08_regrid_weights.py`](https://github.com/xqlsystems/xarray-sql/blob/main/benchmarks/geospatial/08_regrid_weights.py) | | Warp (reproject + resample) | CRS transform *and* interpolation | reproject **UDF** → weight-table `JOIN` | [`09_warp.py`](https://github.com/xqlsystems/xarray-sql/blob/main/benchmarks/geospatial/09_warp.py) | +| Raster × points × polygons | point-in-polygon per pixel | `ST_Within` `JOIN`s, or H3 cell IDs and equality `JOIN`s | [`10_h3_join.py`](https://github.com/xqlsystems/xarray-sql/blob/main/benchmarks/geospatial/10_h3_join.py) | ## 1. A pixel-wise formula is a column expression @@ -206,7 +207,69 @@ the regions are a second SQL table, and the spatial relationship is an ordinary regions keep this simple; arbitrary polygons would follow the same shape, with a point-in-polygon test in the join.) -## 6. The hard cases: where a UDF fits, and where it doesn't +## 6. Rasters, points and polygons: exact geometry, then H3 + +Geospatial data comes in different shapes. A forecast is a raster, cities are +points, countries are polygons, and the interesting questions join them: how +hot will each country be tomorrow, and how many people live where it will pass +35 °C? [`10_h3_join.py`](https://github.com/xqlsystems/xarray-sql/blob/main/benchmarks/geospatial/10_h3_join.py) +asks both in DuckDB, over the latest +[NOAA GFS forecast](https://dynamical.org/catalog/noaa-gfs-forecast/) and +[Natural Earth](https://www.naturalearthdata.com/)'s countries (1:10m) and +cities, two ways. + +**Exact geometry**, with DuckDB's +[`spatial` extension](https://duckdb.org/docs/stable/core_extensions/spatial/overview). +Registering with `geometry=("longitude", "latitude")` gives the raster a point +column, so a country is a point-in-polygon test, and +[`bbox_conjuncts`](#geoarrow-point-geometry-columns) lets the lazy forecast +read only the chunks over it: + +```sql +SELECT latitude, longitude, MAX(temperature_2m) AS t_max +FROM gfs +WHERE "longitude" BETWEEN 68.18 AND 97.40 AND "latitude" BETWEEN 7.97 AND 35.49 + AND ST_Within(geometry, ) +GROUP BY latitude, longitude +``` + +**[H3](https://h3geo.org/)**, Uber's hierarchical grid of hexagons, with +DuckDB's community [`h3` extension](https://duckdb.org/community_extensions/extensions/h3). +Every shape becomes hexagon IDs: a grid cell or a city is the hexagon that +contains it, and a country is the hexagons that cover it. Each spatial +relationship is then an equality `JOIN` on the ID, which any engine runs as a +hash join: + +```sql +SELECT k.name AS country, AVG(h.t_max) AS t_max +FROM country_cells AS k JOIN hex_high AS h ON h.cell = k.cell +GROUP BY k.name +``` + +The two trade exactness for speed. Testing a million grid points against 258 +detailed country outlines takes a minute (60 s); covering the countries with +hexagons takes 13 s once, and the mean high of every country then takes +13 ms. On the forecast of 2026-10-21 the two agreed within 0.29 °C for half of +the countries and 1.24 °C for 90%; the largest differences are small islands +that only a few hexagons cover. Counting people in cities above 35 °C, the two +agree within 10% for eight of the ten countries with the most, but H3 counts far +fewer in two coastal ones (Mozambique, Saudi Arabia): a coastal city's hexagon +averages land with cooler sea. Every step is checked against xarray, shapely and h3-py. + +Three details matter in practice: + +- **Weight by area.** A 0.25° grid's points crowd together toward the poles, so + an unweighted mean over a country's points overweights its high latitudes; + the exact query weights each point by `cos(latitude)`. H3's hexagons have + nearly equal areas, so a plain mean over them is already area-weighted. +- **Cover the coasts.** Covering a country with the hexagons whose *centers* it + contains misses coastal hexagons whose center is at sea, and with them most + coastal cities. The case adds each such hexagon to one country it overlaps. +- **Snap points to a regular grid with arithmetic.** A city's nearest GFS grid + cell is its coordinates rounded to 0.25°, so that join is an equality too, + no geometry needed. + +## 7. The hard cases: where a UDF fits, and where it doesn't Reprojection and regridding are the operations most wedded to the array paradigm. They split cleanly along one line: **is the operation row-independent?** @@ -316,6 +379,13 @@ 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, register with `geometry_crs="EPSG:4326"`: both tags mean WGS84 longitude +and latitude. 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..b70d1263 100644 --- a/docs/limitations.md +++ b/docs/limitations.md @@ -56,6 +56,30 @@ 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:* register with `geometry_crs="EPSG:4326"`. Both tags + mean WGS84 longitude and latitude. + === "Polars" **Float `is_in` literals lose precision.** diff --git a/tests/test_geometry.py b/tests/test_geometry.py index 5e3238d7..e8467cfb 100644 --- a/tests/test_geometry.py +++ b/tests/test_geometry.py @@ -112,6 +112,67 @@ 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 default OGC:CRS84 tag, though both are WGS84 + longitude and latitude. geometry_crs="EPSG:4326" matches it.""" + duckdb = pytest.importorskip("duckdb") + con = duckdb.connect() + try: + con.execute("INSTALL spatial; LOAD spatial;") + except duckdb.Error: + pytest.skip("duckdb spatial extension unavailable") + query = ( + "SELECT count(*) FROM t WHERE ST_Within(geometry, " + f"ST_SetCRS(ST_GeomFromText('{_POLYGON}'), 'EPSG:4326'))" + ) + + xql.register(con, "t", grid, geometry=("x", "y")) + with pytest.raises(duckdb.BinderException, match="coordinate reference"): + con.execute(query).fetchone() + xql.register(con, "t", grid, geometry=("x", "y"), geometry_crs="EPSG:4326") + inside = (grid.y >= -28.75) & (grid.y <= -27.9) + assert con.execute(query).fetchone() == ( + int(inside.sum()) * grid.sizes["x"], + ) + + def test_geopandas_consumes_native_points(grid): gpd = pytest.importorskip("geopandas") From e393efe045d1b7e1199ebb026f2284901119994b Mon Sep 17 00:00:00 2001 From: Alex Merose Date: Wed, 21 Oct 2026 11:07:37 -0700 Subject: [PATCH 2/3] Drop geospatial case 10; keep its DuckDB known issues The geometry and H3 examples now live in the Learn page's tour. Case 10 and its Geospatial section move to the geospatial-case-10 branch, to be evaluated separately. This keeps what writing it turned up: the two DuckDB known issues, the tests that pin them, and a note in the Geospatial page's geometry section. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NyN7idACeDMh94qrnFMu6h --- benchmarks/geospatial/10_h3_join.py | 394 ---------------------------- benchmarks/geospatial/README.md | 11 - docs/geospatial.md | 67 +---- 3 files changed, 2 insertions(+), 470 deletions(-) delete mode 100644 benchmarks/geospatial/10_h3_join.py diff --git a/benchmarks/geospatial/10_h3_join.py b/benchmarks/geospatial/10_h3_join.py deleted file mode 100644 index 68402c42..00000000 --- a/benchmarks/geospatial/10_h3_join.py +++ /dev/null @@ -1,394 +0,0 @@ -# /// script -# requires-python = ">=3.11" -# dependencies = [ -# "xarray-sql[duckdb]", -# "xarray", -# "numpy", -# "pandas", -# "zarr>=3", -# "fsspec", -# "aiohttp", -# "h3>=4", -# "shapely>=2", -# ] -# -# [tool.uv.sources] -# xarray-sql = { path = "../../", editable = true } -# /// -"""Rasters, points and polygons in one query: exact geometry, then H3. - -Geospatial data comes in different shapes: a forecast is a raster, cities -are points, countries are polygons. Joining them answers everyday -questions, like "how hot will each country be tomorrow?" and "how many -people live where it will pass 35 °C?". This case asks both in DuckDB, an -embedded SQL database, two ways: - -1. **Exact geometry**, with DuckDB's ``spatial`` extension. - ``xql.register(..., geometry=("longitude", "latitude"))`` gives the - raster a point ``geometry`` column, so ``ST_Within(geometry, country)`` - is a point-in-polygon test. A single country's query pairs it with - ``xql.bbox_conjuncts``, so the lazy raster reads only the chunks over - that country. -2. **H3**, Uber's hierarchical grid of hexagons, with DuckDB's community - ``h3`` extension. Every shape becomes hexagon IDs: a grid cell or a city - is the hexagon containing it, and a country is the hexagons that cover - it. Each spatial relationship is then an equality ``JOIN`` on the ID, - which any SQL engine runs as a hash join. - -Exact geometry is exact, and slow on detailed outlines: every grid point -is tested against a country's thousands of vertices. H3 pays that cost -once, to cover the countries with hexagons, and answers every query after -that in milliseconds, within a small error. The case checks each step -against xarray, shapely and h3-py references, and reports how far H3's -answers are from the exact ones. - -Data: - -* the latest NOAA GFS forecast (0.25°, global), from dynamical.org's public - Zarr store; the first 24 hours of 2 m temperature, read lazily; -* Natural Earth's countries (1:10m) and populated places (7,342 cities with - their metropolitan populations), read by DuckDB straight from GeoJSON on - GitHub. - -The forecast changes every six hours, so the numbers do too. -""" - -from __future__ import annotations - -import duckdb -import h3 -import numpy as np -import pandas as pd -import shapely -import xarray as xr - -import xarray_sql as xql -from _harness import ( - CaseSkipped, - assert_grid_close, - measured, - run_case, - show_sql, - timed, -) - -_GFS = "https://data.dynamical.org/noaa/gfs/forecast/latest.zarr" -_NATURAL_EARTH = ( - "/vsicurl/https://raw.githubusercontent.com/nvkelso/natural-earth-vector/" - "master/geojson/" -) -_HOT = 35.0 # °C -# H3 resolution 4 hexagons average 1,770 km², about two GFS grid cells at -# the equator, so every hexagon holds at least one grid cell's center. -_RES = 4 -# Countries whose exact means are checked against shapely: large and small, -# tropical and polar, one-piece and many islands. -_CHECKED = ["India", "Brazil", "Norway", "Indonesia", "Egypt"] - - -def _connect() -> duckdb.DuckDBPyConnection: - con = duckdb.connect() - try: - con.execute("INSTALL spatial; LOAD spatial") - con.execute("INSTALL h3 FROM community; LOAD h3") - con.execute(f""" - CREATE TABLE countries AS - SELECT NAME AS name, geom - FROM ST_Read('{_NATURAL_EARTH}ne_10m_admin_0_countries.geojson'); - CREATE TABLE cities AS - SELECT name, pop_max AS population, latitude, longitude, geom - FROM ST_Read('{_NATURAL_EARTH}ne_10m_populated_places_simple.geojson'); - """) - except duckdb.Error as exc: - raise CaseSkipped( - f"DuckDB extensions or Natural Earth ({exc})" - ) from exc - return con - - -def _outline(con: duckdb.DuckDBPyConnection, name: str): - """A country's outline as WKT, and its bounds (xmin, ymin, xmax, ymax).""" - wkt, *bounds = con.execute( - "SELECT ST_AsText(geom), ST_XMin(geom), ST_YMin(geom)," - " ST_XMax(geom), ST_YMax(geom) FROM countries WHERE name = ?", - [name], - ).fetchone() - return wkt, bounds - - -def _inside(field: xr.DataArray, wkt: str, bounds) -> xr.DataArray: - """The reference: ``field`` inside the outline, NaN elsewhere, cropped - to the rows and columns with any cell inside.""" - xmin, ymin, xmax, ymax = bounds - box = field.sel(latitude=slice(ymax, ymin), longitude=slice(xmin, xmax)) - lon, lat = np.meshgrid(box.longitude, box.latitude) - inside = shapely.contains_xy(shapely.from_wkt(wkt), lon, lat) - return ( - box.where(inside) - .dropna("latitude", how="all") - .dropna("longitude", how="all") - ) - - -def main() -> None: - try: - gfs = xr.open_zarr(_GFS, chunks=None) - except Exception as exc: # noqa: BLE001 — any failure → skip, not crash - raise CaseSkipped(f"GFS unavailable ({exc})") from exc - # The latest forecast's first 24 hours: (lead_time, latitude, longitude). - gfs = gfs[["temperature_2m"]].isel(init_time=-1, lead_time=slice(0, 25)) - print( - f" raster: GFS forecast from {gfs.init_time.values}, " - f"{gfs.sizes['lead_time']} hours × {gfs.sizes['latitude']}" - f"×{gfs.sizes['longitude']}" - ) - - with timed("load DuckDB extensions + Natural Earth"): - con = _connect() - n_countries, n_cities = con.execute( - "SELECT (SELECT count(*) FROM countries), (SELECT count(*) FROM cities)" - ).fetchone() - print(f" vector: {n_countries} countries (polygons), {n_cities:,} cities") - - # DuckDB's ST_Read tags GeoJSON's longitude/latitude as EPSG:4326, and - # ST_Within refuses to compare geometries with different tags, so tag - # the raster's points the same way. (Both name WGS84 longitude and - # latitude; xarray-sql's default tag is OGC:CRS84.) - xql.register( - con, - "gfs", - gfs, - geometry=("longitude", "latitude"), - geometry_crs="EPSG:4326", - ) - - # ── 1. One country, exact geometry, read lazily ────────────────────── - name = "India" - wkt, bounds = _outline(con, name) - # ST_Within alone can't skip chunks: the engine doesn't push functions - # into the scan. The bounding box's BETWEENs can, so the query reads - # only the chunks over the country, and ST_Within refines the box. - sql = f""" - SELECT latitude, longitude, MAX(temperature_2m) AS t_max - FROM gfs - WHERE {xql.bbox_conjuncts(bounds, x="longitude", y="latitude")} - AND ST_Within(geometry, - ST_SetCRS(ST_GeomFromText('{wkt}'), 'EPSG:4326')) - GROUP BY latitude, longitude - ORDER BY latitude DESC, longitude - """ - show_sql(sql.replace(wkt, f"<{name}'s outline>")) - for _ in measured(f"SQL: {name}'s 24-hour high, ST_Within + bbox pruning"): - country = xql.to_dataset(con.sql(sql), template=gfs).t_max - for _ in measured("shapely + xarray reference"): - ref = _inside(gfs.temperature_2m.max("lead_time"), wkt, bounds) - assert_grid_close(f"{name}'s 24-hour high (°C)", country, ref) - print( - f" {name}: {int(country.count()):,} grid cells, mean high " - f"{float(country.mean()):.1f} °C" - ) - - # ── 2. The 24-hour high everywhere, as a table ─────────────────────── - # Every later query uses each grid cell's high, so compute it once into - # a DuckDB table. It keeps xarray-sql's point geometry for the - # point-in-polygon joins below. - with timed("SQL: 24-hour high for every grid cell (reads the globe)"): - con.execute(""" - CREATE TABLE high AS - SELECT latitude, longitude, MAX(temperature_2m) AS t_max, - ANY_VALUE(geometry) AS geometry - FROM gfs GROUP BY latitude, longitude - """) - with timed("xarray reference: gfs.max('lead_time')"): - ref_high = gfs.temperature_2m.max("lead_time").load() - high = xql.to_dataset( - con.sql( - "SELECT latitude, longitude, t_max FROM high" - " ORDER BY latitude DESC, longitude" - ), - template=gfs, - ).t_max - assert_grid_close("24-hour high (°C)", high, ref_high) - - # ── 3. Exact geometry: point-in-polygon joins ──────────────────────── - # A grid cell's area shrinks with the cosine of its latitude, so weight - # each cell by it: an unweighted mean would overweight a country's - # high latitudes, where the grid's points crowd together. - exact_means_sql = """ - SELECT k.name AS country, - SUM(h.t_max * cos(radians(h.latitude))) - / SUM(cos(radians(h.latitude))) AS t_max - FROM high AS h JOIN countries AS k ON ST_Within(h.geometry, k.geom) - GROUP BY k.name - """ - show_sql(exact_means_sql, label="SQL: country mean highs, exact") - for _ in measured( - f"SQL: country mean highs, ST_Within ({high.size:,} points × " - f"{n_countries} countries)" - ): - exact_means = con.sql(exact_means_sql).df().set_index("country").t_max - for name in _CHECKED: - cells = _inside(ref_high, *_outline(con, name)) - weight = np.cos(np.radians(cells.latitude)) - want = float(cells.weighted(weight.fillna(0)).mean()) - np.testing.assert_allclose(exact_means[name], want, rtol=1e-6) - print(f" ✅ PASS exact means match shapely for {', '.join(_CHECKED)}") - - # A city's forecast is its nearest grid cell's. On a regular grid that - # needs no geometry: round the city's coordinates to the grid's 0.25° - # step (and wrap 180° to -180°), then join on equality. - exact_people_sql = f""" - WITH snapped AS ( - SELECT name, population, geom, - ROUND(latitude * 4) / 4 AS lat, - ((ROUND(longitude * 4) / 4 + 180) % 360) - 180 AS lon - FROM cities - ) - SELECT k.name AS country, SUM(c.population) AS people - FROM snapped AS c - JOIN high AS h ON h.latitude = c.lat AND h.longitude = c.lon - JOIN countries AS k ON ST_Within(c.geom, k.geom) - WHERE h.t_max > {_HOT} - GROUP BY k.name - """ - show_sql(exact_people_sql, label="SQL: people above 35 °C, exact") - for _ in measured("SQL: people above 35 °C, nearest cell + ST_Within"): - exact_people = con.sql(exact_people_sql).df().set_index("country") - - # Check the snapping join: each city's high, against indexing the - # reference grid at the same rounded coordinates. - cities = con.sql(""" - SELECT c.latitude, c.longitude, h.t_max - FROM cities AS c JOIN high AS h - ON h.latitude = ROUND(c.latitude * 4) / 4 - AND h.longitude = ((ROUND(c.longitude * 4) / 4 + 180) % 360) - 180 - """).df() - - def round_like_sql(v): # halves away from zero, as SQL's ROUND does - return np.sign(v) * np.floor(np.abs(v) + 0.5) - - lat = round_like_sql(cities.latitude.to_numpy() * 4) / 4 - lon = ( - round_like_sql(cities.longitude.to_numpy() * 4) / 4 + 180 - ) % 360 - 180 - want = ref_high.sel( - latitude=xr.DataArray(lat), longitude=xr.DataArray(lon) - ).values - np.testing.assert_allclose(cities.t_max, want, rtol=1e-6) - print(f" ✅ PASS each of {len(cities):,} cities' high is its grid cell's") - - # ── 4. H3: every shape becomes hexagon IDs ─────────────────────────── - # A hexagon belongs to the country that contains its center. A coastal - # hexagon whose center is at sea goes to a country it overlaps, so - # coastal cities aren't lost; each hexagon still has one country. - # ST_Dump splits a country of several islands into its polygons. - with timed(f"SQL: cover the countries with H3 resolution {_RES} hexagons"): - con.execute(f""" - CREATE TABLE country_cells AS - WITH parts AS ( - SELECT name, UNNEST(ST_Dump(geom)).geom AS geom FROM countries - ), - centered AS ( - SELECT name, - UNNEST(h3_polygon_wkb_to_cells(ST_AsWKB(geom), {_RES})) - AS cell - FROM parts - ), - touching AS ( - SELECT name, - UNNEST(h3_polygon_wkt_to_cells_experimental( - ST_AsText(geom), {_RES}, 'overlap')) AS cell - FROM parts - ) - SELECT name, cell FROM centered - UNION ALL - SELECT MIN(name), cell FROM touching - WHERE cell NOT IN (SELECT cell FROM centered) - GROUP BY cell - """) - with timed("SQL: each grid cell's hexagon, and each hexagon's mean high"): - con.execute(f""" - CREATE TABLE hex_high AS - SELECT h3_latlng_to_cell(latitude, longitude, {_RES}) AS cell, - AVG(t_max) AS t_max - FROM high GROUP BY cell - """) - - # Hexagons have nearly equal areas, so a plain average of a country's - # hexagons is already weighted by area. - h3_means_sql = """ - SELECT k.name AS country, AVG(h.t_max) AS t_max - FROM country_cells AS k JOIN hex_high AS h ON h.cell = k.cell - GROUP BY k.name - """ - show_sql(h3_means_sql, label="SQL: country mean highs, H3") - for _ in measured("SQL: country mean highs, H3 equality JOIN"): - h3_means = con.sql(h3_means_sql).df().set_index("country").t_max - - h3_people_sql = f""" - SELECT k.name AS country, SUM(c.population) AS people - FROM cities AS c - JOIN hex_high AS h - ON h.cell = h3_latlng_to_cell(c.latitude, c.longitude, {_RES}) - JOIN country_cells AS k ON k.cell = h.cell - WHERE h.t_max > {_HOT} - GROUP BY k.name - """ - show_sql(h3_people_sql, label="SQL: people above 35 °C, H3") - for _ in measured("SQL: people above 35 °C, H3 equality JOINs"): - h3_people = con.sql(h3_people_sql).df().set_index("country") - - # Check the hexagons' means against h3-py, H3's Python bindings, - # applied to the reference grid. - with timed("h3-py reference: index every grid cell"): - grid = ref_high.to_dataframe().reset_index() - grid["cell"] = [ - h3.str_to_int(h3.latlng_to_cell(y, x, _RES)) - for y, x in zip(grid.latitude, grid.longitude) - ] - want_hex = grid.groupby("cell").temperature_2m.mean() - got_hex = con.sql("SELECT cell, t_max FROM hex_high").df() - got_hex = got_hex.set_index(got_hex.cell.astype("uint64")).t_max - assert set(got_hex.index) == set(want_hex.index), "different hexagons" - np.testing.assert_allclose(got_hex[want_hex.index], want_hex, rtol=1e-6) - print(f" ✅ PASS {len(want_hex):,} hexagons' mean highs match h3-py") - - # ── How close H3 comes to exact geometry ───────────────────────────── - both = pd.concat({"exact": exact_means, "h3": h3_means}, axis=1).dropna() - diff = (both.h3 - both.exact).abs().sort_values() - print( - f"\n Country mean highs, {len(diff)} countries: H3 is within " - f"{diff.median():.2f} °C of exact for half of them, " - f"{diff.quantile(0.9):.2f} °C for 90%." - ) - print(" Largest differences (small or ragged countries, few hexagons):") - for country in diff.index[-3:][::-1]: - print( - f" {country:<24} exact {both.exact[country]:6.2f} °C," - f" H3 {both.h3[country]:6.2f} °C" - ) - - people = pd.concat( - {"exact": exact_people.people, "h3": h3_people.people}, axis=1 - ).fillna(0) - people = people.sort_values("exact", ascending=False) - print(f"\n People in cities above {_HOT:.0f} °C, top 10 countries:") - print(f"\n {'Country':<24} {'exact':>12} {'H3':>12}") - for country, row in people.head(10).iterrows(): - print(f" {country:<24} {int(row.exact):>12,} {int(row.h3):>12,}") - print( - f" {'Total':<24} {int(people.exact.sum()):>12,}" - f" {int(people.h3.sum()):>12,}" - ) - print( - "\n H3 counts a city as hot when its hexagon's mean high is, and" - " exact geometry\n when its nearest grid cell's is. A coastal" - " city's hexagon averages land with\n cooler sea, so H3 counts" - " fewer hot coastal cities." - ) - - -if __name__ == "__main__": - raise SystemExit( - run_case(main, "Rasters × points × polygons: exact geometry and H3") - ) diff --git a/benchmarks/geospatial/README.md b/benchmarks/geospatial/README.md index 80cda53f..c046adc3 100644 --- a/benchmarks/geospatial/README.md +++ b/benchmarks/geospatial/README.md @@ -24,7 +24,6 @@ plain-English definition of the operation, and computes the same numbers. | 07 | `07_reproject_udf.py` | per-pixel CRS transform | scalar **UDF** (`reproject()` from the geo extension), à la PostGIS `ST_Transform` | | 08 | `08_regrid_weights.py` | interpolation to a new grid | sparse-weight table `JOIN` + weighted `GROUP BY` | | 09 | `09_warp.py` | reproject **and** resample (warp) | reproject **UDF** (07) → weight table `JOIN` (08) | -| 10 | `10_h3_join.py` | point-in-polygon per pixel, city and country | `ST_Within` `JOIN`s (DuckDB `spatial`), or H3 cell IDs and equality `JOIN`s (DuckDB `h3`) | Cases 01–06 show operations that are *natively* relational. Cases 07–08 are the "hardest" array operations — reprojection and regridding — and show where a UDF @@ -32,9 +31,6 @@ fits (a per-row coordinate transform) versus where the operation is really a sparse matrix multiply expressed as a `JOIN`. Case 09 composes the two into a full **warp** (GDAL/rasterio `reproject`): the 07 UDF reprojects the target grid, arrays turn the reprojected points into bilinear weights, and the 08 `JOIN` applies them. -Case 10 joins a raster, points and polygons in DuckDB, first with exact -geometry (`ST_Within`) and then through [H3](https://h3geo.org/) hexagon IDs, -which turn each spatial relationship into an equality `JOIN`. See [`docs/geospatial.md`](../../docs/geospatial.md) for the full narrative, including *where the array paradigm still earns its keep* (generating the @@ -64,13 +60,6 @@ interpolation weights — the geometry — which SQL applies but does not comput validates against xarray's `.interp()` at the reprojected points, with Earth Engine's own lon/lat SRTM as a second, cross-CRS check. All three run against Earth Engine using your existing `gcloud` login, and skip cleanly without it. -- **10 H3 join** — the latest **NOAA GFS** forecast (0.25° global) from - [dynamical.org](https://dynamical.org)'s public Zarr store, and - [Natural Earth](https://www.naturalearthdata.com/)'s 1:10m countries and - populated places, which DuckDB reads straight from GeoJSON on GitHub. Needs - DuckDB's `spatial` extension and its community `h3` extension (installed on - first run), and skips cleanly offline. About three minutes, most of it the - exact point-in-polygon join and the h3-py reference. ## Running diff --git a/docs/geospatial.md b/docs/geospatial.md index 04b2e491..4595fdd2 100644 --- a/docs/geospatial.md +++ b/docs/geospatial.md @@ -44,7 +44,7 @@ actually runs — and this suite works through nearly all of it: | Transformed Eulerian Mean (circulation diagnostics — zonal means and anomalies) | cases 03, 04 | | Forecast evaluation (scoring forecasts against ground truth) | case 05 | | Regridding and reprojection (resolution and CRS changes) | cases 07, 08, 09 | -| Spatial joins (large polygon-to-polygon joins) | partly: case 10 joins a raster, points and polygons, exactly and through H3 cells; polygon-to-polygon joins are a vector-data problem and not covered | +| Spatial joins (large polygon-to-polygon joins) | *not covered* — a vector-data problem; the closest analogue here is the raster × vector join in case 06 | So the claim isn't that a few cherry-picked operations happen to be relational. It's that an independent survey of the operations geoscience runs at scale, run @@ -63,7 +63,6 @@ through SQL one by one, turns out to be — almost entirely — queries. | Reprojection | per-pixel CRS transform | scalar **UDF** (`ST_Transform`-style) | [`07_reproject_udf.py`](https://github.com/xqlsystems/xarray-sql/blob/main/benchmarks/geospatial/07_reproject_udf.py) | | Regridding | interpolation to a new grid | sparse-weight table `JOIN` | [`08_regrid_weights.py`](https://github.com/xqlsystems/xarray-sql/blob/main/benchmarks/geospatial/08_regrid_weights.py) | | Warp (reproject + resample) | CRS transform *and* interpolation | reproject **UDF** → weight-table `JOIN` | [`09_warp.py`](https://github.com/xqlsystems/xarray-sql/blob/main/benchmarks/geospatial/09_warp.py) | -| Raster × points × polygons | point-in-polygon per pixel | `ST_Within` `JOIN`s, or H3 cell IDs and equality `JOIN`s | [`10_h3_join.py`](https://github.com/xqlsystems/xarray-sql/blob/main/benchmarks/geospatial/10_h3_join.py) | ## 1. A pixel-wise formula is a column expression @@ -207,69 +206,7 @@ the regions are a second SQL table, and the spatial relationship is an ordinary regions keep this simple; arbitrary polygons would follow the same shape, with a point-in-polygon test in the join.) -## 6. Rasters, points and polygons: exact geometry, then H3 - -Geospatial data comes in different shapes. A forecast is a raster, cities are -points, countries are polygons, and the interesting questions join them: how -hot will each country be tomorrow, and how many people live where it will pass -35 °C? [`10_h3_join.py`](https://github.com/xqlsystems/xarray-sql/blob/main/benchmarks/geospatial/10_h3_join.py) -asks both in DuckDB, over the latest -[NOAA GFS forecast](https://dynamical.org/catalog/noaa-gfs-forecast/) and -[Natural Earth](https://www.naturalearthdata.com/)'s countries (1:10m) and -cities, two ways. - -**Exact geometry**, with DuckDB's -[`spatial` extension](https://duckdb.org/docs/stable/core_extensions/spatial/overview). -Registering with `geometry=("longitude", "latitude")` gives the raster a point -column, so a country is a point-in-polygon test, and -[`bbox_conjuncts`](#geoarrow-point-geometry-columns) lets the lazy forecast -read only the chunks over it: - -```sql -SELECT latitude, longitude, MAX(temperature_2m) AS t_max -FROM gfs -WHERE "longitude" BETWEEN 68.18 AND 97.40 AND "latitude" BETWEEN 7.97 AND 35.49 - AND ST_Within(geometry, ) -GROUP BY latitude, longitude -``` - -**[H3](https://h3geo.org/)**, Uber's hierarchical grid of hexagons, with -DuckDB's community [`h3` extension](https://duckdb.org/community_extensions/extensions/h3). -Every shape becomes hexagon IDs: a grid cell or a city is the hexagon that -contains it, and a country is the hexagons that cover it. Each spatial -relationship is then an equality `JOIN` on the ID, which any engine runs as a -hash join: - -```sql -SELECT k.name AS country, AVG(h.t_max) AS t_max -FROM country_cells AS k JOIN hex_high AS h ON h.cell = k.cell -GROUP BY k.name -``` - -The two trade exactness for speed. Testing a million grid points against 258 -detailed country outlines takes a minute (60 s); covering the countries with -hexagons takes 13 s once, and the mean high of every country then takes -13 ms. On the forecast of 2026-10-21 the two agreed within 0.29 °C for half of -the countries and 1.24 °C for 90%; the largest differences are small islands -that only a few hexagons cover. Counting people in cities above 35 °C, the two -agree within 10% for eight of the ten countries with the most, but H3 counts far -fewer in two coastal ones (Mozambique, Saudi Arabia): a coastal city's hexagon -averages land with cooler sea. Every step is checked against xarray, shapely and h3-py. - -Three details matter in practice: - -- **Weight by area.** A 0.25° grid's points crowd together toward the poles, so - an unweighted mean over a country's points overweights its high latitudes; - the exact query weights each point by `cos(latitude)`. H3's hexagons have - nearly equal areas, so a plain mean over them is already area-weighted. -- **Cover the coasts.** Covering a country with the hexagons whose *centers* it - contains misses coastal hexagons whose center is at sea, and with them most - coastal cities. The case adds each such hexagon to one country it overlaps. -- **Snap points to a regular grid with arithmetic.** A city's nearest GFS grid - cell is its coordinates rounded to 0.25°, so that join is an equality too, - no geometry needed. - -## 7. The hard cases: where a UDF fits, and where it doesn't +## 6. The hard cases: where a UDF fits, and where it doesn't Reprojection and regridding are the operations most wedded to the array paradigm. They split cleanly along one line: **is the operation row-independent?** From 347dcc9d0e8e0069648e17d7fc72f06bfebd1156 Mon Sep 17 00:00:00 2001 From: Alex Merose Date: Wed, 21 Oct 2026 11:15:27 -0700 Subject: [PATCH 3/3] Retag polygons OGC:CRS84 rather than tagging points EPSG:4326 EPSG:4326 officially puts latitude first, and the registered points are longitude first, as OGC:CRS84 says. Retagging the polygons from ST_Read keeps every tag true, so a later ST_Transform needs no always_xy. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NyN7idACeDMh94qrnFMu6h --- docs/geospatial.md | 5 +++-- docs/limitations.md | 7 +++++-- tests/test_geometry.py | 21 +++++++++------------ 3 files changed, 17 insertions(+), 16 deletions(-) diff --git a/docs/geospatial.md b/docs/geospatial.md index 4595fdd2..86b7a9cb 100644 --- a/docs/geospatial.md +++ b/docs/geospatial.md @@ -318,8 +318,9 @@ destination. The CRS tag defaults to `OGC:CRS84`; pass 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, register with `geometry_crs="EPSG:4326"`: both tags mean WGS84 longitude -and latitude. And for now, join the geometry column to a *table* of polygons +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). diff --git a/docs/limitations.md b/docs/limitations.md index b70d1263..f2607028 100644 --- a/docs/limitations.md +++ b/docs/limitations.md @@ -77,8 +77,11 @@ Pick your engine: the polygons come from `ST_Read`, which tags GeoJSON `EPSG:4326`, while the registered points carry the default `OGC:CRS84`. - - *What to do:* register with `geometry_crs="EPSG:4326"`. Both tags - mean WGS84 longitude and latitude. + - *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" diff --git a/tests/test_geometry.py b/tests/test_geometry.py index e8467cfb..6317a42c 100644 --- a/tests/test_geometry.py +++ b/tests/test_geometry.py @@ -150,27 +150,24 @@ def test_duckdb_spatial_join_needs_a_point_or_a_table(grid, spatial_con): 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 default OGC:CRS84 tag, though both are WGS84 - longitude and latitude. geometry_crs="EPSG:4326" matches it.""" + 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") - query = ( - "SELECT count(*) FROM t WHERE ST_Within(geometry, " - f"ST_SetCRS(ST_GeomFromText('{_POLYGON}'), 'EPSG:4326'))" - ) - 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).fetchone() - xql.register(con, "t", grid, geometry=("x", "y"), geometry_crs="EPSG:4326") + con.execute(query.format(polygon)).fetchone() + retagged = f"ST_SetCRS({polygon}, 'OGC:CRS84')" inside = (grid.y >= -28.75) & (grid.y <= -27.9) - assert con.execute(query).fetchone() == ( - int(inside.sum()) * grid.sizes["x"], - ) + want = int(inside.sum()) * grid.sizes["x"] + assert con.execute(query.format(retagged)).fetchone() == (want,) def test_geopandas_consumes_native_points(grid):