Spatial phase II: Map projections

doc/syntax/coord/map.qmd

@@ -0,0 +1,181 @@
+---
+title: Map
+---
+
+The map coordinate system facilitates the display of geographical data.
+It can project data to various coordinate reference systems or projections.
+It is used to display choropleths and other types of maps.
+
+## Default aesthetics
+
+* **Primary**: `lon` for longitude (horizontal)
+* **Secondary**: `lat` for latitude (vertical)
+
+Users can provide their own aesthetic names if needed. 
+For example, if using `x` and `y` aesthetics:
+
+```ggsql
+PROJECT x, y TO map
+```
+
+This maps `x` to longitude and `y` to latitude. 
+This is a convenience when converting from a Cartesian coordinate system, without having to modify all the mappings.
+
+## Settings
+
+* `clip`: Should data be removed if it appears outside the bounds of the coordinate system, or has become invalid due to the projection. 
+  Defaults to `true`.
+* `crs` A string giving the target Coordinate Reference System (CRS) holding parameters of a coordinate transformation.
+  This is only used with unnamed projections and therefore mutually exclusive with the `origin` and `parallel` settings.
+  By default, matches the Spatial Reference Identifiers (SRIDs) of geometry columns.
+  Typically takes one of the following forms:
+    * A [proj string](https://proj.org/en/stable/usage/quickstart.html) like `'+proj=longlat +lon_0=90'` for 2D longitude/latitude with 90° longitude offset.
+    * An [EPSG code](https://en.wikipedia.org/wiki/EPSG_Geodetic_Parameter_Dataset) like `'EPSG:3395'` for World Mercator projection.
+* `source` A string giving the source CRS for the data to use as fallback when not explicit in geometry. 
+  Useful when using non-spatial layers that don't carry SRIDs.
+  Defaults to `'EPSG:4326'` or longitude/latitude coordinates.
+* `bounds` A numeric array of length 4 giving the bounding box of the area to view.
+  The numbers mean ('xmin', 'ymin', 'xmax', 'ymax') in that order, in the units of the `crs` argument.
+  Be mindful that your data might be in longitude/latitude degrees, but the `crs` might use meters, or vice versa.
+  Bounds can have `Inf`/`-Inf` values to indicate the visible area's extent.
+  Using `null` values indicates to fall back to the data's extent, so only partial bounds need to be indicated.
+* `origin` Projection origin location.
+  This is only used with named projections and therefore mutually exclusive with `crs`.
+  One of the following:
+    * A scalar number, setting the origin longitude or 'central meridian' in degrees between −180 and +180.
+      Corresponds to the `+lon_0` proj string parameter.
+    * An array of two numbers, where the first number is the scalar option above.
+      The second number sets the latitude of origin in degrees between −90 and +90.
+      The latter corresponds to the `+lat_0` proj string parameter.
+      If unsupported by the projection, the latitude is silently ignored.
+* `parallel` Setting for conic or cylindrical projections.
+  This is only used with named projections and therefore mutually exclusive with `crs`.
+  If unsupported by the projection, this setting is silently ignored.
+  One of the following:
+  * A scalar number, setting the standard parallel as degrees between −90 and +90.
+    Corresponds to the `+lat_1` proj string parameter.
+  * An array of two numbers, where the first number is the scalar option above.
+    The second number sets the second standard parallel as degrees between −90 and +90.
+
+## Supported projections
+
+There are two main ways to set the map projection. 
+One is via a named projection, which allows for the `origin` and `parallel` settings.
+
+```ggsql
+PROJECT TO orthographic SETTING origin => (150, 52)
+```
+
+The second option is to use the generic `PROJECT TO map` together with the `crs` setting.
+The code below is equivalent to the code above.
+
+```ggsql
+PROJECT TO map SETTING crs => '+proj=ortho +lon_0=150 +lat_0=52'
+```
+
+In this context, 'supported' means that we are reasonably confident that we can draw a serviceable map.
+Below follows a table giving an overview of the supported projections. 
+The projections are recognised by name in the 'string' column or from the 'proj string' CRS setting.
+More information about them can be found on the [PRØJ website](https://proj.org/en/stable/operations/projections/index.html).
+
+| name                                  | string                  | origin | parallel     | proj string      |
+|---------------------------------------|-------------------------|--------|--------------|------------------|
+| Geographic (unprojected)              | `geographic`            | 0      | NA           | `+proj=longlat`  |
+| Mercator                              | `mercator`              | 0      | NA           | `+proj=merc`     |
+| Orthographic                          | `orthographic`          | (0, 0) | NA           | `+proj=ortho`    |
+| Miller Cylindrical                    | `miller`                | 0      | NA           | `+proj=mill`     |
+| Equidistant Cylindrical (Plate Carée) | `equirectangular`       | 0      | NA           | `+proj=eqc`      |
+| Stereographic                         | `stereographic`         | (0, 0) | NA           | `+proj=stere`    |
+| Gnomonic                              | `gnomonic`              | (0, 0) | NA           | `+proj=gnom`     |
+| Cylindrical Equal Area                | `equal_area`            | 0      | NA           | `+proj=cea`      |
+| Mollweide                             | `mollweide`             | 0      | NA           | `+proj=moll`     |
+| Sanson-Flamsteed Sinusoidal           | `sinusoidal`            | 0      | NA           | `+proj=sinu`     |
+| Eckert IV                             | `eckert4`               | 0      | NA           | `+proj=eck4`     |
+| Natural Earth                         | `natural`               | 0      | NA           | `+proj=natearth` |
+| Winkel Tripel                         | `winkel_tripel`         | 0      | NA           | `+proj=wintri`   |
+| Albers Equal Area                     | `albers`                | (0, 0) | (29.5, 45.5) | `+proj=aea`      |
+| Lambert Conformal Conic               | `lambert_conformal`     | (0, 0) | (29.5, 45.5) | `+proj=lcc`      |
+| Lambert Azimuthal                     | `lambert`               | (0, 0) | NA           | `+proj=laea`     |
+| Azimuthal Equidistant                 | `azimuthal_equidistant` | (0, 0) | NA           | `+proj=aeqd`     |
+| Interrupted Goode Homolosine          | `igh`                   | 0      | NA           | `+proj=igh`      |
+| Robinson                              | `robinson`              | 0      | NA           | `+proj=robin`    |
+
+The way to draw an unsupported map is to use the `crs` setting.
+Unsupported maps are likely drawn without projection-appropriate clipping.
+
+## Examples
+
+Note that depending on your reader, you may need to activate modules for spatial analysis.
+
+```{ggsql}
+-- For example, for DuckDB, one could use:
+INSTALL spatial;
+LOAD spatial;
+```
+
+### Via named projection
+
+```{ggsql}
+VISUALISE continent AS fill FROM ggsql:world
+DRAW spatial
+PROJECT TO robinson
+```
+
+### Via proj string
+
+```{ggsql}
+VISUALISE continent AS fill FROM ggsql:world
+DRAW spatial
+PROJECT TO map SETTING crs => '+proj=robin'
+```
+
+### Setting the origin
+
+```{ggsql}
+VISUALISE continent AS fill FROM ggsql:world
+DRAW spatial
+PROJECT TO orthographic
+  SETTING origin => (133.77, -25.27)
+```
+
+### Zooming in
+
+The default zoom behaviour is to zoom in where the data is.
+
+```{ggsql}
+VISUALISE continent AS fill FROM ggsql:world
+DRAW spatial 
+  FILTER subregion == 'Western Africa'
+PROJECT TO orthographic
+```
+
+An alternative is to zoom using the `bounds` setting.
+Note that while `ggsql:world` is defined in degrees longitude/latitude, the `orthographic` projection uses metres as unit.
+
+```{ggsql}
+VISUALISE continent AS fill FROM ggsql:world
+DRAW spatial
+PROJECT TO orthographic
+  SETTING bounds => (-1868152, 468481, 1638882, 2917218)
+```
+
+The `bounds` arguments can take infinites to take the visible area's extreme values, or `null` to take the data's bounds.
+
+```{ggsql}
+VISUALISE continent AS fill FROM ggsql:world
+DRAW spatial 
+  FILTER subregion == 'Western Africa'
+PROJECT TO orthographic
+  SETTING bounds => (-1868152, -Inf, 1638882, null)
+```
+
+### Parallels
+
+Default conic projections aren't always kind to the southern hemisphere.
+This can be changed by tweaking the `parallel` setting.
+
+```{ggsql}
+VISUALISE continent AS fill FROM ggsql:world
+DRAW spatial 
+PROJECT TO albers SETTING parallel => (-29.5, -45.5)
+```

doc/syntax/coord/polar.qmd

@@ -16,7 +16,7 @@ Users can provide their own aesthetic names if needed. For example, if using `x`
 PROJECT y, x TO polar
 ```
 
-This maps `y` to radius and `x` to angle. This is useful when converting from a cartesian coordinate system without editing all the mappings.
+This maps `y` to radius and `x` to angle. This is useful when converting from a Cartesian coordinate system without editing all the mappings.
 
 ## Settings
 * `clip`: Should data be removed if it appears outside the bounds of the coordinate system. Defaults to `true`

src/CLAUDE.md

@@ -46,11 +46,10 @@ Grammar lives in [`/tree-sitter-ggsql/`](../tree-sitter-ggsql/) — when adding
 | `duckdb.rs` | DuckDB (in-memory or file) | `duckdb` (default) |
 | `sqlite.rs` | SQLite | `sqlite` (default) |
 | `odbc.rs` | ODBC | `odbc` (default) |
-| `snowflake.rs` | Snowflake (very small placeholder today) | — |
 | `connection.rs` | Connection-string parsing for all of the above | — |
 | `data.rs`, `spec.rs` | `Spec` type returned by `execute()`, plus DataFrame conversion | — |
 
-`SqlDialect` trait in `mod.rs` lets each driver supply its own type names and information-schema queries. PostgreSQL has a feature flag (`postgres`) but no driver file yet.
+`SqlDialect` trait in `mod.rs` lets each driver supply its own type names, information-schema queries, and spatial helper methods (`sql_st_transform`, `sql_geometry_to_wkb`, `sql_geometry_bbox`, `sql_spatial_setup`).
 
 ### `execute/`
 
@@ -97,15 +96,12 @@ Defined in `Cargo.toml`:
 | `sqlite` | ✓ | SQLite reader |
 | `odbc` | ✓ | ODBC reader |
 | `parquet` | ✓ | Parquet support in readers/data |
-| `ipc` | ✓ | Arrow IPC support |
+| `spatial` | ✓ | Spatial/geometry support (geozero for WKT↔GeoJSON) |
 | `vegalite` | ✓ | Vega-Lite writer |
 | `builtin-data` | ✓ | Bundled penguins/airquality datasets |
-| `postgres` | — | PostgreSQL reader (declared, not yet implemented) |
-| `ggplot2` | — | ggplot2 writer (declared, not yet implemented) |
-| `all-readers` | — | `duckdb` + `postgres` + `sqlite` + `odbc` |
-| `all-writers` | — | `vegalite` + `ggplot2` + `plotters` |
+| `all-readers` | — | `duckdb` + `sqlite` + `odbc` |
 
-`ggsql-wasm` builds with `default-features = false` plus `vegalite`, `sqlite`, `builtin-data`. `ggsql-jupyter` builds with `duckdb`, `sqlite`, `odbc`, `vegalite`.
+`ggsql-wasm` builds with `default-features = false` plus `vegalite`, `sqlite`, `builtin-data`. `ggsql-jupyter` builds with `duckdb`, `vegalite`.
 
 ## Testing
 

src/Cargo.toml

@@ -36,7 +36,7 @@ bytes = { workspace = true }
 adbc_core = { version = "0.23", optional = true }
 
 # Spatial
-geozero = { workspace = true, optional = true, features = ["with-wkb", "with-geojson"] }
+geozero = { workspace = true, optional = true, features = ["with-wkb", "with-wkt", "with-geojson"] }
 
 # Serialization
 serde.workspace = true

src/execute/mod.rs

@@ -1368,6 +1368,19 @@ pub fn prepare_data_with_reader(query: &str, reader: &dyn Reader) -> Result<Prep
         layer_queries.push(layer_query);
     }
 
+    // Apply projection transforms (post-stat, pre-fetch)
+    let mut project = specs[0]
+        .project
+        .take()
+        .unwrap_or_else(crate::plot::projection::Projection::cartesian);
+    project.apply_projection_transforms(
+        &specs[0].layers,
+        &mut layer_queries,
+        dialect,
+        &execute_query,
+    )?;
+    specs[0].project = Some(project);
+
     // Phase 2: Deduplicate and execute unique queries
     let mut query_to_result: HashMap<String, DataFrame> = HashMap::new();
     for (idx, q) in layer_queries.iter().enumerate() {

src/execute/scale.rs

@@ -1723,6 +1723,8 @@ mod tests {
             coord,
             aesthetics,
             properties: std::collections::HashMap::new(),
+            map_projection: None,
+            computed: std::collections::HashMap::new(),
         });
 
         // Create scale for pos2 (theta in polar) without explicit expand