grib/projection/
merc.rs

1#[cfg(feature = "gridpoints-proj")]
2use super::OsgeoProj;
3use super::{
4    Ellipsoid, Project,
5    helpers::{m, sinhpsi2tanphi},
6};
7
8const HALF_PI: f64 = std::f64::consts::FRAC_PI_2;
9
10/// Parameters for Mercator projection.
11#[derive(Debug, PartialEq, Clone)]
12pub struct Params {
13    /// Ellipsoid definition.
14    pub ellipsoid: Ellipsoid,
15    /// Latitude of true scale (in degrees).
16    pub lat_ts: f64,
17    /// Central meridian (in degrees).
18    pub lon_0: f64,
19}
20
21#[cfg(feature = "gridpoints-proj")]
22impl OsgeoProj for Params {
23    fn proj_args(&self) -> String {
24        let Self {
25            ellipsoid: Ellipsoid { a, b, .. },
26            lat_ts,
27            lon_0,
28        } = self;
29        format!("+a={a} +b={b} +proj=merc +lat_ts={lat_ts} +lon_0={lon_0}")
30    }
31}
32
33/// Mercator projection.
34pub struct Projection {
35    lam0: f64,
36    e: f64,
37    e_sq: f64,
38    ak0: f64,
39}
40
41impl Projection {
42    pub fn new(p: &Params) -> Result<Self, &'static str> {
43        let Params {
44            ellipsoid: Ellipsoid { a, e, e_sq, .. },
45            lat_ts,
46            lon_0,
47        } = p;
48        let lam0 = lon_0.to_radians();
49        let phi_ts = lat_ts.to_radians().abs();
50        if phi_ts >= HALF_PI {
51            return Err("Invalid value for lat_ts: |lat_ts| should be <= 90°");
52        }
53
54        let k0 = if *e_sq == 0.0 {
55            // sphere
56            phi_ts.cos()
57        } else {
58            // ellipsoid
59            let (sinφts, cosφts) = phi_ts.sin_cos();
60            m(sinφts, cosφts, *e_sq)
61        };
62        let ak0 = a * k0;
63        let context = Projection {
64            lam0,
65            e: *e,
66            e_sq: *e_sq,
67            ak0,
68        };
69        Ok(context)
70    }
71
72    fn ellipsoidal_forward(&self, (lambda, phi): &(f64, f64)) -> Result<(f64, f64), &'static str> {
73        let &x = lambda;
74        let (sinφ, cosφ) = phi.sin_cos();
75        let y = (sinφ / cosφ).asinh() - self.e * (self.e * sinφ).atanh();
76        Ok((x, y))
77    }
78
79    fn spheroidal_forward(&self, (lambda, phi): &(f64, f64)) -> Result<(f64, f64), &'static str> {
80        let &x = lambda;
81        let y = phi.tan().asinh();
82        Ok((x, y))
83    }
84
85    fn ellipsoidal_inverse(&self, (x, y): &(f64, f64)) -> Result<(f64, f64), &'static str> {
86        let phi = sinhpsi2tanphi(y.sinh(), self.e)
87            .ok_or(
88                "the inverse of the isometric latitude function could not be solved numerically",
89            )?
90            .atan();
91        let &lambda = x;
92        Ok((lambda, phi))
93    }
94
95    fn spheroidal_inverse(&self, (x, y): &(f64, f64)) -> Result<(f64, f64), &'static str> {
96        let phi = y.sinh().atan();
97        let &lambda = x;
98        Ok((lambda, phi))
99    }
100}
101
102impl Project for Projection {
103    fn forward(&self, xy: &(f64, f64)) -> Result<(f64, f64), &'static str> {
104        if self.e_sq == 0.0 {
105            self.spheroidal_forward(xy)
106        } else {
107            self.ellipsoidal_forward(xy)
108        }
109    }
110
111    fn inverse(&self, xy: &(f64, f64)) -> Result<(f64, f64), &'static str> {
112        if self.e_sq == 0.0 {
113            self.spheroidal_inverse(xy)
114        } else {
115            self.ellipsoidal_inverse(xy)
116        }
117    }
118
119    fn a(&self) -> &f64 {
120        &self.ak0
121    }
122
123    fn lam0(&self) -> &f64 {
124        &self.lam0
125    }
126}
127
128#[cfg(all(test, feature = "gridpoints-proj"))]
129mod tests {
130    use proj::Proj;
131
132    use super::*;
133
134    const FORWARD_TOLERANCE_METERS: f64 = 1e-8;
135    const INVERSE_TOLERANCE_RADIANS: f64 = 1e-12;
136
137    #[test]
138    fn agrees_with_proj_for_ellipsoid() {
139        assert_agrees_with_proj(Params {
140            ellipsoid: Ellipsoid::from_a_and_b(6_378_137., 6_356_752.314_245),
141            lat_ts: 20.,
142            lon_0: 140.,
143        });
144    }
145
146    #[test]
147    fn agrees_with_proj_for_sphere() {
148        assert_agrees_with_proj(Params {
149            ellipsoid: Ellipsoid::from_a_and_b(6_371_229., 6_371_229.),
150            lat_ts: -15.,
151            lon_0: -30.,
152        });
153    }
154
155    fn assert_agrees_with_proj(params: Params) {
156        let proj = Proj::new(&params.proj_args()).unwrap();
157        let projection = Projection::new(&params).unwrap();
158        let coordinates: [(f64, f64); 4] = [(-10., -70.), (0., 0.), (25., 45.), (80., 170.)];
159
160        for (lat, lon) in coordinates {
161            let lonlat = (lon.to_radians(), lat.to_radians());
162            let expected_xy = proj.project(lonlat, false).unwrap();
163            let actual_xy = projection.project(&lonlat, false).unwrap();
164            assert_coordinates_close(actual_xy, expected_xy, FORWARD_TOLERANCE_METERS);
165
166            let expected_lonlat = proj.project(expected_xy, true).unwrap();
167            let actual_lonlat = projection.project(&expected_xy, true).unwrap();
168            assert_coordinates_close(actual_lonlat, expected_lonlat, INVERSE_TOLERANCE_RADIANS);
169        }
170    }
171
172    fn assert_coordinates_close(actual: (f64, f64), expected: (f64, f64), tolerance: f64) {
173        assert!(
174            (actual.0 - expected.0).abs() <= tolerance
175                && (actual.1 - expected.1).abs() <= tolerance,
176            "actual {actual:?} differs from expected {expected:?} by more than {tolerance}"
177        );
178    }
179}
grib/projection/merc.rs

grib/projection/
merc.rs
