feat: pyproj CRS extension — generic reproject(x, y, src_crs, dst_crs) UDF - #225
feat: pyproj CRS extension — generic reproject(x, y, src_crs, dst_crs) UDF#225Mmoncadaisla wants to merge 12 commits into
Conversation
…) UDF Adds xarray_sql/proj.py, an optional pyproj extension that registers an ST_Transform-style scalar UDF. The CRS pair is part of the query (any spelling pyproj.CRS accepts, and it may vary per row via ordinary SQL expressions) instead of being baked in at UDF registration time, which is how the geospatial benchmarks previously hard-coded it. All PROJ work runs on a dedicated pool of Python-owned worker threads: constructing a Transformer on a DataFusion runtime thread segfaults inside PROJ, while identical concurrent work on Python threads is stable (pyproj 3.7 / PROJ 9.5). Each pool thread caches one transformer per CRS pair, and pyproj releases the GIL during transforms, so the UDF is parallel across partitions — the previous single-chunk/serial-UDF workaround in the benchmarks is no longer needed. XarrayContext auto-registers reproject() when pyproj is installed (pip install xarray-sql[proj]); proj.register() is the explicit hook for plain SessionContexts or custom names. Benchmarks 07 and 09 now use the extension instead of their local hard-coded UDFs. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The y-chunking exists to force several DataFusion partitions so the reproject() UDF demonstrably runs in parallel; say so explicitly instead of hiding it behind a computed chunk size. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
|
Closing for now — opened prematurely; will resubmit once it's ready for review. |
The kernel materialized the two broadcast CRS string columns with to_pylist() — two Python object allocations per row before PROJ did any work. At benchmark scale that dominated everything: a 606M-pixel aggregate-only reprojection took 496s while the same scan without the UDF took 4.5s. Establish CRS uniqueness with pyarrow.compute.unique (vectorized C) instead. A batch with one CRS pair — the overwhelmingly common case — becomes a single vectorized PROJ call with no per-row Python. The per-row grouping path remains for genuinely varying CRS (e.g. a CASE expression picking the UTM zone). Same 606M-pixel run after: 16.0s (37.9M px/s), a 31x speedup, with bit-identical results and flat (~2GB) streaming memory. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
|
xref: wrt the tokio issue, check out #145. We may want to tackle that soon(ish) since it might be useful for proj, too. |
Maybe we should make this |
Aligns with the base + add-ons direction: geo is the umbrella extra for CRS reprojection today and polygon/vector support later, rather than being named after the first feature it happens to ship. The module stays proj.py — it does CRS reprojection specifically. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
alxmrs
left a comment
There was a problem hiding this comment.
LGTM! A few technical fixes, then I think this should be ready to land. Thanks for the amazing contribution! :)
The reproject kernel's to_numpy(zero_copy_only=False) reads as a copy but is a zero-copy view for null-free batches; the flag only permits the NaN materialization that null handling requires. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The stack-size hypothesis was wrong: the crashes come from pyproj < 3.8 leaving a dangling PJ_CONTEXT when the ephemeral Python thread state of a non-Python-created thread is torn down (pyproj#1541, unreleased). Python-owned pool threads keep contexts alive, and the per-thread cache stays worthwhile on fixed pyproj (rebuild costs 0.07-12 ms/batch vs ~10 us for the pool round-trip). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Both the constant-CRS fast path and the varying-CRS path now only build the (src, dst, mask) groups; one shared block submits to the pool and scatters results. Fast-path detection stays on the vectorized Arrow kernel. NULL-CRS handling is uniform. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-authored-by: Alex Merose <al@merose.com> Signed-off-by: Miguel Moncada Isla <48254102+Mmoncadaisla@users.noreply.github.com>
| # ] | ||
| # | ||
| # [tool.uv.sources] | ||
| # xarray-sql = { path = "../../", editable = true } |
There was a problem hiding this comment.
Should we modify this to include the geo extra? (is that a thing?)
There was a problem hiding this comment.
Yes! "xarray-sql[geo]" works there and uv resolves the extra against the [tool.uv.sources] path source (I verified it with a quick uv run --script probe: the geo extra pulls pyproj into the script's environment). Adding it
alxmrs
left a comment
There was a problem hiding this comment.
A few random nits. Thanks for this significant contribution. LGTM.
Unpack the untyped Transformer.transform result before returning (warn_return_any), and annotate the optional proj module import as ModuleType | None. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The extra was renamed proj -> geo but the prose still said 'pyproj extension' in seven places; pyproj stays mentioned as the dependency. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
What
Adds
xarray_sql/proj.py, an optional pyproj extension (pip install xarray-sql[geo]) that registers anST_Transform-style scalar UDF:The CRS pair is part of the query, not baked in at UDF registration time (which is how benchmarks 07/09 previously hard-coded it, each with its own copy of the UDF). Any CRS spelling
pyproj.CRSaccepts works (authority codes, WKT, PROJ strings), and because the arguments are ordinary SQL expressions the CRS may vary per row — e.g. aCASEexpression selecting the UTM zone from the longitude.XarrayContextauto-registersreproject()when pyproj is installed (mirroring thecftime()UDF);proj.register(ctx, name=...)is the explicit hook for plainSessionContexts or custom names. Benchmarks 07 and 09 now use the extension and their local UDF definitions are removed.The threading discovery (worth flagging)
The docstring in 07 said "PROJ's context is not thread-safe," and the benchmark serialized the UDF by forcing a single chunk. While porting this I found the real failure mode is narrower and more interesting:
Transformer.from_crs+transform()with no locking and no crashes (pyproj 3.7.2 / PROJ 9.5.1 — verified with a standalone ThreadPoolExecutor stress repro).So the extension never calls pyproj on the calling thread: all PROJ work runs on a small dedicated pool of Python-owned worker threads, each caching one transformer per CRS pair. pyproj releases the GIL during transforms, so concurrent partitions still transform in parallel — the single-chunk/serial-UDF workaround in the benchmarks is gone rather than relocated.
Why the CRS is a query argument (design context)
Passing the CRS pair as SQL arguments isn't just ergonomics — given today's DataFusion, it's the only honest place to put it. DataFusion's type system cannot yet attach metadata like a CRS to a column type: that's the extension-types gap tracked in apache/datafusion#12644, and the reason spatial support overall (apache/datafusion#7859) was blocked for a year before the GeoArrow dense-union workaround. Until columns can carry their CRS, any "implicit CRS" design would have to smuggle it through UDF registration state — which is exactly the hard-coding this PR removes.
This also means the UDF stays complementary to where the ecosystem is heading: when extension types land and GeoArrow columns know their CRS (e.g. via datafusion-contrib/datafusion-geo), a geometry-typed
ST_Transform(geom, dst)can arrive alongsidereproject()rather than replace it — raster coordinates enter xarray-sql as plain float columns either way.Semantics
{x, y}struct (always_xyorder) from one call — one PROJ transform per row.inffor out-of-domain points is normalized to NaN so results round-trip to xarray like any other missing value.pyproj.exceptions.CRSError) and fails the query loudly rather than returning wrong coordinates.Testing
tests/test_proj.py(pyproj added to thetestextra): equivalence with direct pyproj on a multi-partition grid, round-trip UTM↔lon/lat, per-row CRS viaCASE, NULL/out-of-domain handling, invalid-CRS failure, and custom-name registration on a plainSessionContext. The suite passed 20/20 consecutive runs locally (the round-trip test is the one that segfaulted reliably before the worker-pool design — it mixes two CRS pairs in one query).Caveats: benchmarks 07/09 need Earth Engine credentials I don't have wired up here — they compile and the SQL shape is covered by the tests, but a maintainer EE run would be a good final check.
Docs
docs/geospatial.md(section 6 narrative + SQL snippets), both READMEs, and the module docstring double as the extension's documentation.🤖 Generated with Claude Code