ggsql spatial: what works, what breaks, and a prepare-upstream recipe

Run against ggsql 0.5.2 (built from main, 2026-10-13), DuckDB 1.5.2 bundled, DuckDB spatial 1.5.2. All queries are in render.sh calls below; figures in out/.

How ggsql does spatial today (from the source)

The spatial layer takes a geometry aesthetic; anything the backend can ST_AsWKB works. SRID is sniffed with ST_SRID, falling back to EPSG:4326.
Reprojection is ST_Transform in the database. There is no PROJ in the Rust tree. The 21 named projections in map_projections.rs exist to produce the horizon/visible-area polygon in lon/lat, which is then ST_Transformed and used to clip layers and draw the panel. PROJECT TO crs SETTING target => 3031 is parsed back to stere so it gets the same treatment.
Non-spatial layers (path, polygon, point, text) under a map projection are densified in SQL and ST_Transformed too, but are not clipped to the horizon.
Rows come back as Arrow with a WKB column; geozero decodes; the renderer (hephaestus) draws polygons with holes natively.

So the shape of it is: Rust knows about projections, never evaluates one.

1. The WKB route (`DRAW spatial`) – the advertised path

figure	query	result
`00_baseline_robinson`	`DRAW spatial PROJECT TO robinson`	fine
`01_spatial_stere_south`	`PROJECT TO stereographic SETTING origin => (0, -90)`	fine, horizon-clipped
`02_spatial_epsg3031`	`PROJECT TO crs SETTING target => 3031`	identical to 01: 3031 is recognised
`03_spatial_epsg3031_bounds`	same + `bounds => (+/-5.5e6)`	fine; note the slit from the pole in Antarctica
`20_holes_spatial`	southern Africa, `PROJECT TO geographic`	Lesotho hole rendered correctly

The slit in 03 is in the data (Natural Earth’s Antarctica ring runs down both sides of the antimeridian to -90); every renderer draws it. See recipe step 3a.

2. The vertex-table route (`DRAW polygon PARTITION BY ... ORDER BY ...`)

This is the “group + arrange is how you construct a polygon” construction. It works today, with PARTITION BY fid, poly, ring ORDER BY ord.

figure	query	result
`10_vertices_lonlat_stere`	lon/lat vertex table, `PROJECT TO stereographic origin (0,-90)`	tiny globe: northern-hemisphere vertices project to +/-1e8 m, nothing is clipped, the panel autoscales to that
`10b_..._bounds`	same with `bounds`	the data in view is fine, but edges from unclipped vertices wrap through the panel, and the Antarctica slit edge spirals around the pole
`11_vertices_3031_cartesian`	table pre-projected upstream with PROJ, no cut, `PROJECT TO cartesian`	same explosion, now the user’s fault rather than ggsql’s
`12_vertices_3031cut_cartesian`	cut at lat -40 on the sphere (a ring around the pole, so no antimeridian seam), densified, projected upstream	clean. Slit still present (data).
`21_holes_vertices_ring`	`PARTITION BY fid, poly, ring`	hole ring is a separate polygon, so it fills over the exterior (darker green). Wrong.
`22_holes_vertices_feature`	`PARTITION BY fid ORDER BY poly, ring, ord`	the hole appears (renderer winding rule) plus a stroke along the join from exterior to interior. ggpolypath, 2016, again.

Two things fall out:

The vertex route is one fill-rule/ring setting away from being complete: a polygon layer that treats ring > 0 (or a hole aesthetic) as an interior ring, with even-odd fill and no join stroke. The renderer already does this for WKB input, so it is plumbing, not geometry.
The map coord does not horizon-clip non-spatial layers. The fix inside ggsql is the same ST_Intersection with the horizon polygon it already applies to the spatial layer. The fix outside ggsql is recipe step 3.

3. The prepare-upstream route (wk / PROJ / geos / wkpool in R)

recipe.R is the R side; prep.py / prep_geoparquet.py are the same steps in pyproj/shapely, because there was no R on the box this ran on. Same PROJ.

figure	query	result
`13_spatial_prepared_3031`	GeoParquet with CRS = EPSG:3031 (PROJJSON), `DRAW spatial target => 3031`	fails: ggsql emits `ST_GeomFromWKB(CAST(geom AS BLOB))`, and DuckDB 1.5 cannot cast a typed `GEOMETRY('<projjson>')` column to BLOB. The builtin `ggsql:world` dodges this because ggsql `ST_AsWKB`s it at registration.
`13a_..._wkb_nosrc`	`SELECT ..., ST_AsWKB(geom) AS geom FROM '...'` then `target => 3031`	SRID lost on the WKB, falls back to 4326, transforms already-projected metres as degrees: a blue disc
`13b_..._wkb_src`	same + `SETTING source => 'EPSG:3031', target => 3031`	the good one: horizon, graticule, clean Antarctica, no slit
`13c_..._cartesian`	same, `PROJECT TO cartesian`	also fine; no graticule

So the whole recipe for a user is:

SELECT name, continent, population, ST_AsWKB(geom) AS geom
FROM 'prepared_3031.parquet'
VISUALISE continent AS fill
DRAW spatial
PROJECT TO crs SETTING source => 'EPSG:3031', target => 3031

and the preparation is three lines of wk/PROJ/geos (recipe.R): intersect with a horizon ring on the sphere, densify, transform; optionally snap + buffer(0) to dissolve the antimeridian slit.

Things worth raising on posit-dev/ggsql

Bug: CAST(geom AS BLOB) fails for GeoParquet with a non-trivial CRS on DuckDB 1.5. ST_AsWKB(geom) works everywhere ggsql already uses it.
Docs: source => is the escape hatch for pre-projected WKB, and ST_AsWKB + source is the recipe. Neither is obvious from the crs page.
Docs: a polar example (target => 3031, origin => (0,-90)) and a worked “prepare upstream” example. The crs page’s region table already has Antarctica with azimuthal_equidistant; 3031 is what people actually use.
Small feature: hole support in DRAW polygon (ring aesthetic or even-odd fill, no join stroke), so the vertex-table route is complete.
Small feature: horizon-clip non-spatial layers under a map coord, as the spatial layer already is.
Observation: the Antarctica slit is in ggsql:world; a snapped buffer(0) on that one feature at build time removes it from every example.

None of 1-3 needs a new dependency on either side.

Files

render.sh NAME "QUERY" -> out/NAME.svg + .png
prep.py, prep_geoparquet.py -> the three vertex tables and the GeoParquet
recipe.R -> the same preparation in wk/PROJ/geos (unrun; check API names)
world.parquet -> copy of ggsql’s builtin dataset (Natural Earth 110m)

4. The spiral (figs 30-34): an ORDER BY bug, and a densification lesson

DRAW polygon ... PARTITION BY fid, poly, ring ORDER BY ord renders Antarctica cleanly (32). Add continent AS fill and a spiral appears at the pole (33, and 10b). Cause, in the source:

execute/layer.rs appends the layer’s ORDER BY to the layer query.
densify_edges (plot/layer/geom/mod.rs) wraps that query and numbers edges with ROW_NUMBER() OVER (PARTITION BY <partitions + discrete mappings> ORDER BY (SELECT NULL)).

So edge order within a partition is whatever the engine’s window operator returns, not the user’s ORDER BY. Changing the partition key set (by mapping a discrete aesthetic) changes that order for a few rows at the pole slit, giving the edge (-180,-90) -> (180,-84.7), which ggsql interpolates linearly in lon/lat and projects: lon sweeps 360 degrees while lat climbs 5 degrees, an Archimedean spiral. Fix: use the layer’s ORDER BY as the window order.

Independently of the bug, 30_ant_nopole shows the rule itself: an edge from lon 180 to lon -180 at constant lat becomes the whole parallel. Linear densification in source coordinates is only safe on rings already cut at the horizon and the antimeridian; do that upstream (recipe step 3) and the densify-then-project path is fine.

Add to the issue list: 7. Bug: layer ORDER BY is not propagated into densify_edges’ window; path/polygon edge order under a map coord depends on engine row order within the partition.