Address review: fair cold-vs-cold (.compute), to_pandas headline, doc fixes

Fairness: case 05's reference used `.load()`, which caches data in place on the
forecasts/truth objects the SQL table also reads from — so running the reference
after the SQL query (as the harness does) could serve it a warm read. Switch to
`.compute()`, which returns a fresh array and leaves the inputs lazy. Verified:
reading a window repeatedly in one process stays flat, and neither side warms the
other. The other cases either reopen their data (06) or recompute eagerly with
`chunks=None` (02-04), so they were already cold; only 05 leaked.

Case 05: render the headline RMSE-by-lead table with `to_pandas()` instead of a
hand-rolled print loop; clarify that the `chunks=` arg is the Arrow batch size,
not a filter (no data dropped).

Case 06: the docstring claimed "full ARCO-ERA5" as the dataset; the query
aggregates one day's window (the WHERE prunes the read). Make that precise.

Case 07: drop the "graduate into the package" paragraph; clarify that there is
one dataset (UTM x/y as SQL input, EE's pixelLonLat as the independent
reference), not the same image opened twice in two CRS.

Docs: correct the profiling methodology note (the caching trap is in-place
`.load()`, not `open_zarr` auto-caching); point the "earns its keep" section
forward to Results/Analysis/Conclusion; move "Running the suite" above Results;
link James Bourbeau's profile; show the to_pandas headline repr.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01AWvrZYAT2NbuETBqNAN3o9

Author: alxmrs

benchmarks/geospatial/05_forecast_skill.py

@@ -90,11 +90,15 @@ def _open(url: str) -> xr.Dataset:
 def _reference_rmse(forecasts: xr.Dataset, truth: xr.Dataset) -> xr.DataArray:
     """xarray reference: per (model, lead), align truth at valid_time, take RMSE.
 
-    The 64×32 windows are tiny, so the reference loads them and reduces in
-    memory; the SQL side above stays lazy.
+    The 64×32 windows are tiny, so the reference reads them into memory and
+    reduces there; the SQL side above stays lazy. We use ``.compute()`` rather
+    than ``.load()`` deliberately: ``.load()`` caches the data *in place* on the
+    shared ``forecasts``/``truth`` objects (which the SQL table also reads from),
+    which would let a profiled reference serve a warm read — ``.compute()``
+    returns a fresh array and leaves the inputs lazy, so each measurement is cold.
     """
-    f = forecasts[_VAR].load()
-    e = truth[_VAR].load()
+    f = forecasts[_VAR].compute()
+    e = truth[_VAR].compute()
     leads = f.prediction_timedelta.values
     per_lead = []
     for lead in leads:
@@ -141,6 +145,9 @@ def main() -> None:
     )
 
     ctx = xql.XarrayContext()
+    # chunks here is the Arrow batch (partition) size the table streams in, not a
+    # filter — no data is dropped. truth spans only the valid-time window, so
+    # time:100 makes it a single partition; forecasts stream a few inits at a time.
     ctx.from_dataset("forecasts", forecasts, chunks={"time": 6})
     ctx.from_dataset("era5", truth, chunks={"time": 100})
 
@@ -170,18 +177,17 @@ def main() -> None:
 
     show_result(got)
 
-    # Headline table: error growth with forecast horizon, both models.
-    lead_days = got["lead"].values / np.timedelta64(1, "D")
-    pangu_rmse = got.sel(model="pangu")
-    graphcast_rmse = got.sel(model="graphcast")
-    print("\n  2m-temperature RMSE (K) vs lead time — lower is better:")
-    print(f"  {'lead (days)':>12} {'Pangu':>9} {'GraphCast':>11}")
-    for i in range(0, len(lead_days), 4):
-        print(
-            f"  {lead_days[i]:>12.2f} "
-            f"{float(pangu_rmse.isel(lead=i)):>9.3f} "
-            f"{float(graphcast_rmse.isel(lead=i)):>11.3f}"
-        )
+    # Headline: error growth with forecast horizon, both models. The gridded SQL
+    # result round-trips to a pandas table directly — index is lead (in days),
+    # one column per model.
+    table = (
+        got.assign_coords(lead=got["lead"].values / np.timedelta64(1, "D"))
+        .to_pandas()
+        .T
+    )
+    table.index.name = "lead (days)"
+    print("\n  2m-temperature RMSE (K) vs lead — lower is better:\n")
+    print(table.iloc[::4].round(3).to_string())
 
 
 if __name__ == "__main__":

benchmarks/geospatial/06_zonal_vector.py

@@ -31,8 +31,10 @@
 prunes it to one day), the regions are a second SQL table, and the spatial
 relationship is an ordinary ``BETWEEN``.
 
-Dataset: full ARCO-ERA5 + a handful of continental-scale bounding boxes
-(longitudes in ERA5's 0–360°E convention).
+Dataset: the full ARCO-ERA5 archive opened *lazily* — the table spans the whole
+record, but the query aggregates only one day's window (the ``WHERE`` prunes the
+read; it is not a scan of the full archive) — plus a handful of continental-scale
+bounding boxes (longitudes in ERA5's 0–360°E convention).
 """
 
 from __future__ import annotations

benchmarks/geospatial/07_reproject_udf.py

@@ -29,16 +29,15 @@
     SELECT x, y, reproject(x, y)['lon'] AS lon, reproject(x, y)['lat'] AS lat
     FROM grid
 
-**The reference is Earth Engine itself.** We open a UTM grid through
-[Xee](https://github.com/google/Xee) carrying ``ee.Image.pixelLonLat()`` — Earth
-Engine's *own* geodesy engine computes the true lon/lat of every UTM pixel
-centre. So this case validates our PROJ-in-SQL transform against a fully
-**independent** reprojection implementation (EE), not against PROJ again. They
-agree to sub-metre precision.
-
-The UDF (built on ``pyproj``, vectorized over each Arrow batch) mirrors the
-``cftime()`` UDF already shipped in xarray-sql (see ``xarray_sql/cftime.py``);
-it could graduate into the package as ``xql.register_reproject_udf``.
+**The reference is Earth Engine itself.** There is *one* dataset: a single UTM
+grid opened through [Xee](https://github.com/google/Xee) carrying
+``ee.Image.pixelLonLat()``. Each pixel arrives with two things — its UTM ``x``/
+``y`` (the grid coordinates, our SQL input) and Earth Engine's *own* per-pixel
+``longitude``/``latitude`` (data variables, the reference). So we are not
+opening the same image twice in two CRS; we feed the UTM coordinates to the PROJ
+UDF and check the lon/lat it returns against EE's independently-computed lon/lat
+for the *same* pixels. The reference is a different geodesy engine, not PROJ
+again, and they agree to sub-metre precision.
 
 PROJ's context is not thread-safe and DataFusion evaluates projection
 expressions concurrently, so we return *both* coordinates from one

docs/geospatial.md

@@ -28,7 +28,8 @@ computes the same numbers — at ERA5's real 0.25° global resolution.
 
 The operations here aren't a set we hand-picked to suit SQL. They're taken from
 [**Large Scale Geospatial Benchmarks**](https://github.com/coiled/benchmarks/discussions/1545)
-(coiled/benchmarks #1545), a discussion James Bourbeau opened in 2024 asking the
+(coiled/benchmarks #1545), a discussion [James Bourbeau](https://github.com/jrbourbeau)
+opened in 2024 asking the
 geospatial and climate community a pointed question: what are the *end-to-end
 workflows* the Xarray/Dask ecosystem needs to handle smoothly at the
 100-terabyte scale? The replies are a representative survey of what geoscience
@@ -161,11 +162,22 @@ for both models, matches an xarray reference, and reproduces the published resul
 that GraphCast edges out Pangu at every lead — the classic "error grows with
 horizon" curve (≈0.3 K at 6 h rising to ≈2.5 K at 9 days):
 
+The result round-trips to a `pandas` table directly (`got.to_pandas()`), RMSE in
+kelvin by lead time:
+
 ```
-  lead (days)     Pangu   GraphCast
-         0.25     0.336       0.296
-         5.25     1.469       1.228
-         9.25     2.814       2.380
+model        graphcast  pangu
+lead (days)
+0.25             0.296  0.336
+1.25             0.464  0.554
+2.25             0.608  0.734
+3.25             0.780  0.936
+4.25             0.988  1.191
+5.25             1.228  1.469
+6.25             1.470  1.747
+7.25             1.763  2.096
+8.25             2.092  2.489
+9.25             2.380  2.814
 ```
 
 ## 5. Raster × vector zonal statistics is a range `JOIN`
@@ -252,18 +264,43 @@ a weight matrix. `xarray-sql` sits downstream of all that as a query front-end:
 once the data is openable as an `xarray.Dataset`, these everyday operations are
 expressible — and accessible — as SQL.
 
+That is the qualitative boundary; the rest of this page puts numbers to it. The
+**Results** below report what each operation costs in SQL versus the array
+reference, **Analysis** explains *why* the relational form is slower and where the
+time goes, and the **Conclusion** turns the whole thing into a when-to-use-which.
+
+## Running the suite
+
+```shell
+python benchmarks/geospatial/02_climatology.py   # inside the repo
+uv run benchmarks/geospatial/02_climatology.py   # standalone (PEP 723 deps)
+```
+
+Each script prints its SQL, runs the array reference, and asserts the two agree.
+See [`benchmarks/geospatial/README.md`](../benchmarks/geospatial/README.md) for
+the full list and dataset notes.
+
 ## Results
 
 Correctness is the headline, but every case is also profiled. The numbers below
 come from [`run_perf.sh`](../benchmarks/geospatial/run_perf.sh) on a Google Compute
 Engine `e2-standard-8` (8 vCPU, 32 GB) in `us-central1` — in-region with the
 ARCO-ERA5 and WeatherBench 2 buckets, so the cloud read is fast. Each case runs
 **once per fresh process**, with no warmup, repeated five times: the SQL operation
-*and* the xarray reference each pay a **cold** read on every measurement. (A warm
-in-process loop would flatter the array side — `xr.open_zarr(chunks=None)` caches
-each variable after its first read, so the reference would serve later repetitions
-from RAM while the SQL side re-reads the store. One process per repetition makes
-both sides pay the read every time.)
+*and* the xarray reference each pay a **cold** read on every measurement.
+
+Fairness here took some care, because the obvious trap is caching. A reference
+that calls `.load()` caches its data *in place* on the very object the SQL table
+also reads from, so a later read — even just running the reference after the SQL
+query in the same process — could be served warm. We close that two ways. The one
+case that loads shared objects (05, forecast skill) uses `.compute()` instead,
+which returns a fresh array and leaves the inputs lazy, caching nothing; the other
+references either reopen their data or recompute their reduction eagerly on every
+read (`chunks=None` is NumPy, not Dask, so there is no graph to keep warm). And
+`run_perf.sh` runs each case in a fresh process per repetition, ruling out any
+carryover between reps. We verified the result directly: reading a window
+repeatedly in one process stays flat, and running either side after the other
+speeds up neither — the SQL query and the reference do not warm each other.
 
 | Case | Step | median (s) | stdev (s) | min (s) | max (s) | peak (MB) |
 |---|---|--:|--:|--:|--:|--:|
@@ -376,14 +413,3 @@ scale. The point of this suite is not to crown a winner but to show that the lin
 between the two is exactly where the operation is dense versus where it is
 relational, and that for a surprising share of geoscience, the operation is
 relational.
-
-## Running the suite
-
-```shell
-python benchmarks/geospatial/02_climatology.py   # inside the repo
-uv run benchmarks/geospatial/02_climatology.py   # standalone (PEP 723 deps)
-```
-
-Each script prints its SQL, runs the array reference, and asserts the two agree.
-See [`benchmarks/geospatial/README.md`](../benchmarks/geospatial/README.md) for
-the full list and dataset notes.