ggplot_rs/coord/
map.rs

1//! `coord_quickmap` / `coord_map` (#411): map coordinates for plain
2//! longitude / latitude data (no `sf` geometry).
3
4use crate::render::Rect;
5
6use super::{AxisSpan, Coord};
7
8/// Which map approximation [`CoordMap`] uses.
9#[derive(Clone, Copy, Debug, PartialEq, Eq)]
10pub enum MapProjection {
11    /// ggplot2's `coord_quickmap()`: linear axes, with the aspect ratio set so
12    /// one degree of latitude and of longitude have their true relative
13    /// length at the middle latitude (`1 / cos(lat)`).
14    Quick,
15    /// `coord_map("mercator")`: latitude stretched as Mercator
16    /// (`ln tan(π/4 + φ/2)`), equal aspect in projected units.
17    Mercator,
18}
19
20/// Map coordinates over `x` = longitude, `y` = latitude (degrees), letterboxed
21/// to the map's aspect ratio.
22pub struct CoordMap {
23    pub projection: MapProjection,
24    x_span: Option<AxisSpan>,
25    y_span: Option<AxisSpan>,
26}
27
28impl CoordMap {
29    pub fn new(projection: MapProjection) -> Self {
30        CoordMap {
31            projection,
32            x_span: None,
33            y_span: None,
34        }
35    }
36
37    fn ratio(&self) -> f64 {
38        let (Some(x), Some(y)) = (self.x_span, self.y_span) else {
39            return 1.0;
40        };
41        let dx = (x.max - x.min).abs();
42        if dx < f64::EPSILON {
43            return 1.0;
44        }
45        match self.projection {
46            MapProjection::Quick => {
47                let mid = ((y.min + y.max) / 2.0).to_radians().cos().abs().max(1e-6);
48                (y.max - y.min).abs() / dx / mid
49            }
50            MapProjection::Mercator => (mercator(y.max) - mercator(y.min)).abs() / dx,
51        }
52    }
53
54    /// Mercator-warp a normalized y (within the data span; margins linear).
55    fn warp_y(&self, n: f64) -> f64 {
56        let Some(s) = self
57            .y_span
58            .filter(|s| s.max > s.min && (s.pmax - s.pmin).abs() > 1e-12)
59        else {
60            return n;
61        };
62        if self.projection != MapProjection::Mercator {
63            return n;
64        }
65        let v = s.min + (n - s.pmin) / (s.pmax - s.pmin) * (s.max - s.min);
66        if v < s.min || v > s.max {
67            return n;
68        }
69        let (a, b) = (mercator(s.min), mercator(s.max));
70        s.pmin + (mercator(v) - a) / (b - a) * (s.pmax - s.pmin)
71    }
72}
73
74/// Mercator y in degree-like units (so it shares units with longitude).
75fn mercator(lat: f64) -> f64 {
76    let phi = lat.clamp(-85.05, 85.05).to_radians();
77    (std::f64::consts::FRAC_PI_4 + phi / 2.0)
78        .tan()
79        .ln()
80        .to_degrees()
81}
82
83impl Coord for CoordMap {
84    fn transform(&self, point: (f64, f64), plot_area: &Rect) -> (f64, f64) {
85        let r = self
86            .data_rect(plot_area)
87            .expect("coord_map always letterboxes");
88        let (nx, ny) = (point.0, self.warp_y(point.1));
89        (r.x + nx * r.width, r.y + (1.0 - ny) * r.height)
90    }
91
92    fn set_domains(&mut self, x: Option<AxisSpan>, y: Option<AxisSpan>) {
93        self.x_span = x;
94        self.y_span = y;
95    }
96
97    fn fixed_aspect(&self) -> Option<f64> {
98        Some(self.ratio())
99    }
100
101    fn data_rect(&self, plot_area: &Rect) -> Option<Rect> {
102        let data_aspect = self.ratio();
103        let pixel_aspect = plot_area.height / plot_area.width;
104        let (w, h, off_x, off_y) = if pixel_aspect > data_aspect {
105            let h = plot_area.width * data_aspect;
106            (plot_area.width, h, 0.0, (plot_area.height - h) / 2.0)
107        } else {
108            let w = plot_area.height / data_aspect;
109            (w, plot_area.height, (plot_area.width - w) / 2.0, 0.0)
110        };
111        Some(Rect {
112            x: plot_area.x + off_x,
113            y: plot_area.y + off_y,
114            width: w,
115            height: h,
116        })
117    }
118}
119
120#[cfg(test)]
121mod tests {
122    use super::*;
123
124    fn span(min: f64, max: f64) -> AxisSpan {
125        AxisSpan {
126            min,
127            max,
128            pmin: 0.0,
129            pmax: 1.0,
130        }
131    }
132
133    #[test]
134    fn quickmap_aspect_is_one_over_cos_lat() {
135        // ggplot2: coord_quickmap at 60°N squeezes longitude by cos(60°) = 0.5.
136        let mut c = CoordMap::new(MapProjection::Quick);
137        c.set_domains(Some(span(0.0, 10.0)), Some(span(55.0, 65.0)));
138        assert!((c.ratio() - 2.0).abs() < 1e-9, "{}", c.ratio());
139    }
140
141    #[test]
142    fn mercator_stretches_high_latitudes() {
143        let mut c = CoordMap::new(MapProjection::Mercator);
144        c.set_domains(Some(span(0.0, 10.0)), Some(span(0.0, 60.0)));
145        // 30°N sits below the middle of a 0–60° Mercator map.
146        let mid = c.warp_y(0.5);
147        let want = mercator(30.0) / mercator(60.0);
148        assert!((mid - want).abs() < 1e-9 && mid < 0.5, "{mid}");
149    }
150}
ggplot_rs/coord/map.rs

ggplot_rs/coord/
map.rs
