feat(symbology): geometric interval and standard deviation classes

ROADMAP.md

@@ -57,7 +57,7 @@ The bugs users hit, the docs and a few quick wins: good first pull requests.
 - A pass over the user guide, the comparison pages and the README against what the app does
   today (vertex editing, spatial join and kriging are listed as missing or undocumented; some
   pages still say 16 languages).
-- Styling: standard deviation, geometric interval and manual breaks _(open)_.
+- Styling: manual breaks _(open)_.
 - Selection by polygon and lasso, and a radius tool on the map (Neighborhood Explorer already
   selects by radius) _(open)_.
 
@@ -187,7 +187,8 @@ Open questions rather than plans; discuss before building.
   projections), a search box for layers, places and coordinates, bookmarks, web tile
   layers (XYZ/TMS/WMS, Mapbox Vector Tiles, MapLibre styles, PMTiles and MBTiles archives),
   a basemap gallery (OpenFreeMap, USGS Imagery), worldviews.
-- Cartography: single, unique values, graduated colours and proportional symbols; layer groups;
+- Cartography: single, unique values, graduated colours (natural breaks, equal interval, quantile,
+  geometric interval, standard deviation) and proportional symbols; layer groups;
   label classes with SQL queries, placement and annotation; custom colour ramps (and QGIS, GDAL
   and GMT palette files); masking, display filters, blend modes, aggregation, extrusion; graphics.
 - Editing: create, move, rotate, scale, merge, split, explode, vertex editing; snapping, templates, elevation,

apps/rustgis-python/python/rustgis/__init__.py

@@ -157,7 +157,8 @@ def style(self, type: str = "single", field: Optional[str] = None, **params: Any
                 (colours by class) or ``"proportional"`` (symbol sizes by value).
             field: The field the colours or sizes come from (not for ``"single"``).
             **params: More ``layer.symbology`` parameters: ``classes`` (1 to 12), ``method``
-                (``"naturalBreaks"``, ``"equalInterval"``, ``"quantile"``), ``ramp``.
+                (``"naturalBreaks"``, ``"equalInterval"``, ``"quantile"``,
+                ``"geometricInterval"``, ``"standardDeviation"``), ``ramp``.
 
         Returns:
             The layer, so calls can be chained.

crates/assistant/skills/thematic-map/SKILL.md

@@ -9,7 +9,7 @@ description: Style a layer by its data (graduated colours, categories, proportio
    field suits graduated or proportional symbols; text or a few distinct values suit unique values.
 2. Style it with `run_command` `layer.symbology`:
    - Numbers: `{"layer": L, "type": "graduated", "field": F, "classes": 5, "method": "naturalBreaks"}` (or
-     `quantile`, `equalInterval`), with an optional `"ramp"`.
+     `quantile`, `equalInterval`, `geometricInterval`, `standardDeviation`), with an optional `"ramp"`.
      For counts or sizes, `"type": "proportional"` usually reads better than colour.
    - Categories: `{"layer": L, "type": "unique", "field": F}`. More than about 12 categories is hard to read, so
      suggest grouping first.

crates/engine/src/cmd/layer.rs

@@ -109,7 +109,7 @@ pub fn commands() -> Vec<C> {
         .when(need_layer),
         C::new("layer.symbology", "Symbology", "Style › Drawing", symbology)
             .when(need_feature_layer)
-            .params(r#"{"layer"?, "type": "single"|"unique"|"graduated"|"proportional", "field"?, "classes"?: 1..12, "method"?: "naturalBreaks"|"equalInterval"|"quantile", "ramp"?: string}"#),
+            .params(r#"{"layer"?, "type": "single"|"unique"|"graduated"|"proportional", "field"?, "classes"?: 1..12, "method"?: "naturalBreaks"|"equalInterval"|"quantile"|"geometricInterval"|"standardDeviation", "ramp"?: string}"#),
         C::new("layer.ramp", "Color Scheme", "Style › Drawing", |s, p| {
             let r = s.ramp(s_req(p, "ramp")?)?;
             with(s, p, |l| match &mut l.raster {

crates/map/src/symbology.rs

@@ -311,15 +311,27 @@ pub enum ClassMethod {
     NaturalBreaks,
     EqualInterval,
     Quantile,
+    /// Class widths that grow by one ratio, for values that pile up at the low end.
+    GeometricInterval,
+    /// Classes one standard deviation wide about the mean.
+    StandardDeviation,
 }
 
 impl ClassMethod {
-    pub const ALL: [ClassMethod; 3] = [ClassMethod::NaturalBreaks, ClassMethod::EqualInterval, ClassMethod::Quantile];
+    pub const ALL: [ClassMethod; 5] = [
+        ClassMethod::NaturalBreaks,
+        ClassMethod::EqualInterval,
+        ClassMethod::Quantile,
+        ClassMethod::GeometricInterval,
+        ClassMethod::StandardDeviation,
+    ];
     pub fn label(self) -> &'static str {
         match self {
             ClassMethod::NaturalBreaks => "Natural Breaks (Jenks)",
             ClassMethod::EqualInterval => "Equal Interval",
             ClassMethod::Quantile => "Quantile",
+            ClassMethod::GeometricInterval => "Geometric Interval",
+            ClassMethod::StandardDeviation => "Standard Deviation",
         }
     }
     pub fn from_name(s: &str) -> Option<ClassMethod> {
@@ -328,6 +340,8 @@ impl ClassMethod {
             "naturalbreaks" | "natural_breaks" | "jenks" | "natural breaks" => ClassMethod::NaturalBreaks,
             "equalinterval" | "equal_interval" | "equal interval" => ClassMethod::EqualInterval,
             "quantile" => ClassMethod::Quantile,
+            "geometricinterval" | "geometric_interval" | "geometric interval" | "geometric" => ClassMethod::GeometricInterval,
+            "standarddeviation" | "standard_deviation" | "standard deviation" | "stddev" => ClassMethod::StandardDeviation,
             _ => return None,
         })
     }
@@ -1048,8 +1062,50 @@ pub fn class_breaks(sorted: &[f64], k: usize, method: ClassMethod) -> Vec<f64> {
         ClassMethod::EqualInterval => (1..=k).map(|i| lo + (hi - lo) * i as f64 / k as f64).collect::<Vec<_>>(),
         ClassMethod::Quantile => (1..=k).map(|i| sorted.get(((n * i) / k).saturating_sub(1).min(n - 1)).copied().unwrap_or(hi)).collect(),
         ClassMethod::NaturalBreaks => jenks(sorted, k),
+        ClassMethod::GeometricInterval => {
+            // Positive values: each bound is the one before times (hi / lo)^(1/k), worked out on
+            // logarithms. Values from zero or below are taken as if the minimum sat at a hundredth
+            // of the range, which makes the ratio 101^(1/k); each bound is then a share of the way
+            // from lo to hi. Neither adds large values, so the bounds stay finite and in order.
+            let range = hi - lo;
+            (1..=k)
+                .map(|i| {
+                    let f = i as f64 / k as f64;
+                    let v = if lo > 0.0 {
+                        (lo.ln() + (hi.ln() - lo.ln()) * f).exp()
+                    } else {
+                        let t = (101f64.powf(f) - 1.0) / 100.0;
+                        if range.is_finite() { lo + range * t } else { lo * (1.0 - t) + hi * t }
+                    };
+                    v.clamp(lo, hi)
+                })
+                .collect()
+        }
+        ClassMethod::StandardDeviation => {
+            // On values scaled by the largest magnitude, so the squares can't overflow.
+            let s = lo.abs().max(hi.abs());
+            let mean_s = sorted.iter().map(|v| v / s).sum::<f64>() / n as f64;
+            let mean = mean_s * s;
+            let sd = s * (sorted.iter().map(|v| (v / s - mean_s).powi(2)).sum::<f64>() / n as f64).sqrt();
+            // k − 1 inner bounds one deviation apart, centred on the mean. A bound with no value
+            // above the one before it is dropped, so no class is empty.
+            let inner = (1..k).map(|i| mean + sd * (i as f64 - k as f64 / 2.0)).filter(|&v| v > lo && v < hi);
+            let mut b = Vec::new();
+            let mut below = 0;
+            for v in inner.chain(std::iter::once(hi)) {
+                let upto = sorted.partition_point(|x| *x <= v);
+                if upto > below {
+                    b.push(v);
+                    below = upto;
+                }
+            }
+            b
+        }
     };
-    b.dedup_by(|a, b| (*a - *b).abs() < f64::EPSILON);
+    // Only bounds that are the same value: any tolerance merges distinct values somewhere (tiny
+    // ones under an absolute one, large ones under a relative one). Repeated bounds come from
+    // repeated values (quantiles, Jenks), so they are equal exactly.
+    b.dedup();
     if let Some(l) = b.last_mut() {
         *l = hi;
     }
@@ -1748,6 +1804,86 @@ mod tests {
         assert_eq!(class_breaks(&[5.0, 5.0], 4, ClassMethod::NaturalBreaks), vec![5.0]);
     }
 
+    fn close(a: &[f64], b: &[f64]) -> bool {
+        a.len() == b.len() && a.iter().zip(b).all(|(x, y)| (x - y).abs() < 1e-9 * y.abs().max(1.0))
+    }
+
+    #[test]
+    fn geometric_interval_breaks_grow_by_one_ratio() {
+        // Positive values: each bound is the one before times (hi / lo)^(1/k).
+        let v = [1.0, 5.0, 40.0, 300.0, 1000.0];
+        let b = class_breaks(&v, 3, ClassMethod::GeometricInterval);
+        assert!(close(&b, &[10.0, 100.0, 1000.0]), "{b:?}");
+        // From zero (or below), the values are shifted so the minimum sits at a hundredth of the
+        // range: the widths still grow by one ratio, and the last bound is the maximum.
+        let b = class_breaks(&[-50.0, 0.0, 50.0], 2, ClassMethod::GeometricInterval);
+        assert!(close(&b, &[-51.0 + 101f64.sqrt(), 50.0]), "{b:?}");
+        let b = class_breaks(&[0.0, 100.0], 4, ClassMethod::GeometricInterval);
+        let widths: Vec<f64> = std::iter::once(b[0]).chain(b.windows(2).map(|w| w[1] - w[0])).collect();
+        assert!(widths.windows(2).all(|w| w[1] > w[0]), "{b:?}");
+        assert_eq!(b.len(), 4);
+        assert_eq!(b.last(), Some(&100.0));
+        // Far from zero the shift is not added and taken away again, so it can't round off.
+        let b = class_breaks(&[-1e16, -9_999_999_999_999_998.0], 2, ClassMethod::GeometricInterval);
+        assert!(b.iter().all(|v| v.is_finite()), "{b:?}");
+        assert_eq!(b.last(), Some(&-9_999_999_999_999_998.0));
+        // Near the largest float nothing overflows: the bounds stay finite and in order.
+        for v in [[1e308, 1.1e308], [-1e308, 1e308], [1e-308, 1e308]] {
+            let b = class_breaks(&v, 2, ClassMethod::GeometricInterval);
+            assert!(b.len() == 2 && b.iter().all(|x| x.is_finite()) && b[0] >= v[0] && b[0] < b[1], "{v:?}: {b:?}");
+        }
+    }
+
+    #[test]
+    fn standard_deviation_breaks_are_one_deviation_wide_about_the_mean() {
+        // 0..=100: mean 50, population deviation sqrt(850).
+        let v: Vec<f64> = (0..=100).map(f64::from).collect();
+        let sd = 850f64.sqrt();
+        // An even count puts a bound on the mean; an odd one centres a class on it.
+        let b = class_breaks(&v, 4, ClassMethod::StandardDeviation);
+        assert!(close(&b, &[50.0 - sd, 50.0, 50.0 + sd, 100.0]), "{b:?}");
+        let b = class_breaks(&v, 3, ClassMethod::StandardDeviation);
+        assert!(close(&b, &[50.0 - sd / 2.0, 50.0 + sd / 2.0, 100.0]), "{b:?}");
+        // Bounds past the data are dropped rather than leaving empty classes: 8 classes of one
+        // deviation would reach ±4 deviations, but the data spans less than ±2.
+        let b = class_breaks(&v, 8, ClassMethod::StandardDeviation);
+        assert!(close(&b, &[50.0 - sd, 50.0, 50.0 + sd, 100.0]), "{b:?}");
+        // A bound with no value above the one before it goes too: mean 4, deviation sqrt(24), so
+        // (4, 4 + sqrt(24)] would hold nothing between the zeros and the tens.
+        let b = class_breaks(&[0.0, 0.0, 0.0, 10.0, 10.0], 4, ClassMethod::StandardDeviation);
+        assert!(close(&b, &[4.0, 10.0]), "{b:?}");
+        // Large values: the squared deviations would overflow without scaling.
+        let b = class_breaks(&[1e160, 2e160], 2, ClassMethod::StandardDeviation);
+        assert!(close(&b, &[1.5e160, 2e160]), "{b:?}");
+        let b = class_breaks(&[-1e308, 1e308], 2, ClassMethod::StandardDeviation);
+        assert!(close(&b, &[0.0, 1e308]), "{b:?}");
+    }
+
+    #[test]
+    fn close_values_keep_their_classes() {
+        // Bounds closer than f64::EPSILON in absolute terms are still distinct when the values
+        // themselves are that small.
+        for m in ClassMethod::ALL {
+            let b = class_breaks(&[0.0, 1e-18], 2, m);
+            assert!(b.len() == 2 && b[0] >= 0.0 && b[0] < b[1], "{m:?}: {b:?}");
+        }
+        // And close values far from zero keep theirs too.
+        let b = class_breaks(&[1e16, 1e16 + 2.0], 2, ClassMethod::Quantile);
+        assert_eq!(b, vec![1e16, 1e16 + 2.0]);
+    }
+
+    #[test]
+    fn every_method_is_listed_named_and_read_back() {
+        for m in ClassMethod::ALL {
+            let name = format!("{m:?}");
+            let camel = format!("{}{}", name[..1].to_ascii_lowercase(), &name[1..]);
+            assert_eq!(ClassMethod::from_name(&camel), Some(m), "{camel}");
+            assert_eq!(ClassMethod::from_name(&m.label().replace(" (Jenks)", "")), Some(m), "{}", m.label());
+        }
+        assert_eq!(ClassMethod::from_name("stddev"), Some(ClassMethod::StandardDeviation));
+        assert_eq!(ClassMethod::from_name("geometric"), Some(ClassMethod::GeometricInterval));
+    }
+
     #[test]
     fn graduated_symbols_by_value() {
         let mut d = Dataset::new("t", GeometryType::Polygon, vec![Field::new("v", FieldType::Double)]);

docs/user-guide/cartography.md

@@ -18,8 +18,13 @@ Select a layer in Layers and open **Style › Symbology** (or the Symbology pane
   is not fully transparent. A layout legend can leave it out with **Hide other values** (see
   [Layouts](layouts.md)).
 - **Graduated Colors**: 1 to 12 classes of a numeric field by natural breaks (Jenks), equal
-  interval or quantile. If no field is chosen, the current renderer's field is kept when it is
-  numeric; otherwise the first numeric field is used. If the layer has no numeric field, or the
+  interval, quantile, geometric interval or standard deviation. Geometric interval makes class
+  widths grow by one ratio, for values piled up at the low end; for a field that reaches zero or
+  below, the minimum is taken as if it sat at a hundredth of the range, so the widths grow by
+  101^(1/classes). Standard deviation makes classes one deviation wide about the mean, and leaves
+  out a class no value falls in, so a skewed field can get fewer classes than asked. If no field is
+  chosen, the current renderer's field is kept when it is numeric; otherwise the first numeric
+  field is used. If the layer has no numeric field, or the
   switch fails for another reason, the pane shows why and keeps the current renderer. When numeric
   fields exist, a field picker lets you try another one; the error stays visible if a retry fails.
   A problem no field can fix, such as a definition query that keeps no features, needs fixing