grib/projection/
lcc.rs

1#[cfg(feature = "gridpoints-proj")]
2use super::OsgeoProj;
3use super::{
4    Ellipsoid, Project,
5    helpers::{m, phi2, t},
6};
7
8const EPS10: f64 = f64::EPSILON;
9const HALF_PI: f64 = std::f64::consts::FRAC_PI_2;
10const FORTH_PI: f64 = std::f64::consts::FRAC_PI_4;
11
12/// Parameters for Lambert Conformal Conic projection.
13#[derive(Debug, PartialEq, Clone)]
14pub struct Params {
15    /// Ellipsoid definition.
16    pub ellipsoid: Ellipsoid,
17    /// Latitude of origin (in degrees).
18    pub lat_0: f64,
19    /// Central meridian (in degrees).
20    pub lon_0: f64,
21    /// First standard parallel (in degrees).
22    pub lat_1: f64,
23    /// Second standard parallel (in degrees).
24    pub lat_2: f64,
25}
26
27#[cfg(feature = "gridpoints-proj")]
28impl OsgeoProj for Params {
29    fn proj_args(&self) -> String {
30        let Self {
31            ellipsoid: Ellipsoid { a, b, .. },
32            lat_0,
33            lon_0,
34            lat_1,
35            lat_2,
36        } = self;
37        format!(
38            "+a={a} +b={b} +proj=lcc +lat_0={lat_0} +lon_0={lon_0} +lat_1={lat_1} +lat_2={lat_2}"
39        )
40    }
41}
42
43/// Lambert Conformal Conic projection.
44pub struct Projection {
45    lam0: f64,
46    a: f64,
47    e: f64,
48    e_sq: f64,
49    k0: f64,
50    n: f64,
51    c: f64,
52    rho0: f64,
53}
54
55impl Projection {
56    pub fn new(p: &Params) -> Result<Self, &'static str> {
57        let Params {
58            ellipsoid: Ellipsoid { a, e, e_sq, .. },
59            lat_0,
60            lon_0,
61            lat_1,
62            lat_2,
63        } = p;
64        let lam0 = lon_0.to_radians();
65        let phi0 = lat_0.to_radians();
66        let phi1 = lat_1.to_radians();
67        let phi2 = lat_2.to_radians();
68        let k0 = 1.;
69
70        if (phi1 + phi2).abs() < EPS10 {
71            return Err("Invalid value for lat_1 and lat_2: |lat_1 + lat_2| should be > 0");
72        }
73
74        let (sinφ1, cosφ1) = phi1.sin_cos();
75
76        if cosφ1.abs() < EPS10 || phi1.abs() >= HALF_PI {
77            return Err("Invalid value for lat_1: |lat_1| should be < 90°");
78        }
79        if phi2.cos().abs() < EPS10 || phi2.abs() >= HALF_PI {
80            return Err("Invalid value for lat_2: |lat_2| should be < 90°");
81        }
82
83        let is_secant_cone = (phi1 - phi2) >= EPS10; // otherwise, tangent cone
84        let is_ellipsoidal = *e_sq != 0.;
85        let is_phi0_almost_equal_to_half_pi = (phi0.abs() - HALF_PI).abs() < EPS10;
86
87        let context = if is_ellipsoidal {
88            let m1 = m(sinφ1, cosφ1, *e_sq);
89            let t1 = t(cosφ1, sinφ1, *e);
90            let n = if is_secant_cone {
91                n_in_secant_cone_ellipsoidal(phi2, *e, *e_sq, m1, t1)?
92            } else {
93                sinφ1
94            };
95            let c = m1 * t1.powf(-n) / n;
96            let rho0 = if is_phi0_almost_equal_to_half_pi {
97                0.
98            } else {
99                c * t(phi0.cos(), phi0.sin(), *e).powf(n)
100            };
101            Projection {
102                lam0,
103                a: *a,
104                e: *e,
105                e_sq: *e_sq,
106                k0,
107                n,
108                c,
109                rho0,
110            }
111        } else {
112            let n = if is_secant_cone {
113                n_in_secant_cone_spherical(phi1, phi2, cosφ1, phi2.cos())
114            } else {
115                sinφ1
116            };
117            if n == 0. {
118                return Err("Invalid value for lat_1 and lat_2: |lat_1 + lat_2| should be > 0");
119            }
120            let c = cosφ1 * (FORTH_PI + 0.5 * phi1).tan().powf(n) / n;
121            let rho0 = if is_phi0_almost_equal_to_half_pi {
122                0.
123            } else {
124                c * (FORTH_PI + 0.5 * phi0).tan().powf(-n)
125            };
126            Projection {
127                lam0,
128                a: *a,
129                e: *e,
130                e_sq: *e_sq,
131                k0,
132                n,
133                c,
134                rho0,
135            }
136        };
137        Ok(context)
138    }
139
140    const ERR_TRANSFORMATION_OUTSIDE_DOMAIN: &str =
141        "Coordinate transformation outside projection domain";
142
143    fn forward(&self, (lambda, phi): &(f64, f64)) -> Result<(f64, f64), &'static str> {
144        let Self {
145            e,
146            e_sq,
147            k0,
148            n,
149            c,
150            rho0,
151            ..
152        } = self;
153
154        let rho = if (phi.abs() - HALF_PI).abs() < EPS10 {
155            if phi * n <= 0. {
156                return Err(Self::ERR_TRANSFORMATION_OUTSIDE_DOMAIN);
157            }
158            0.
159        } else {
160            c * if *e_sq != 0. {
161                t(phi.cos(), phi.sin(), *e).powf(*n)
162            } else {
163                (FORTH_PI + 0.5 * phi).tan().powf(-n)
164            }
165        };
166        let lambda = lambda * n;
167        let x = k0 * (rho * lambda.sin());
168        let y = k0 * (rho0 - rho * lambda.cos());
169        Ok((x, y))
170    }
171
172    fn inverse(&self, (x, y): &(f64, f64)) -> Result<(f64, f64), &'static str> {
173        let Self {
174            e,
175            e_sq,
176            k0,
177            n,
178            c,
179            rho0,
180            ..
181        } = self;
182
183        let x = x / k0;
184        let y = rho0 - y / k0;
185
186        let rho = x.hypot(y);
187        let lp = if rho != 0. {
188            let (rho, x, y) = if *n < 0. { (-rho, -x, -y) } else { (rho, x, y) };
189
190            let phi = if *e_sq != 0. {
191                let phi = phi2((rho / c).powf(1. / n), *e).ok_or(
192                    "the inverse of the isometric latitude function could not be solved numerically"
193                )?;
194                if phi == f64::INFINITY {
195                    return Err(Self::ERR_TRANSFORMATION_OUTSIDE_DOMAIN);
196                }
197                phi
198            } else {
199                2. * ((c / rho).powf(1. / n)).atan() - HALF_PI
200            };
201            let lambda = x.atan2(y) / n;
202            (lambda, phi)
203        } else {
204            let phi = if *n > 0. { HALF_PI } else { -HALF_PI };
205            (0., phi)
206        };
207
208        Ok(lp)
209    }
210}
211
212impl Project for Projection {
213    fn forward(&self, xy: &(f64, f64)) -> Result<(f64, f64), &'static str> {
214        self.forward(xy)
215    }
216
217    fn inverse(&self, xy: &(f64, f64)) -> Result<(f64, f64), &'static str> {
218        self.inverse(xy)
219    }
220
221    fn a(&self) -> &f64 {
222        &self.a
223    }
224
225    fn lam0(&self) -> &f64 {
226        &self.lam0
227    }
228}
229
230fn n_in_secant_cone_ellipsoidal(
231    phi2: f64,
232    e: f64,
233    e_sq: f64,
234    m1: f64,
235    t1: f64,
236) -> Result<f64, &'static str> {
237    let err_message = "Invalid value for eccentricity";
238    let (sinφ2, cosφ2) = phi2.sin_cos();
239    let n = (m1 / m(sinφ2, cosφ2, e_sq)).ln();
240    if n == 0. {
241        return Err(err_message);
242    }
243    let denom = (t1 / t(cosφ2, sinφ2, e)).ln();
244    if denom == 0. {
245        return Err(err_message);
246    }
247    let n = n / denom;
248    Ok(n)
249}
250
251fn n_in_secant_cone_spherical(phi1: f64, phi2: f64, cosφ1: f64, cosφ2: f64) -> f64 {
252    (cosφ1 / cosφ2).ln() / ((FORTH_PI + 0.5 * phi2).tan() / (FORTH_PI + 0.5 * phi1).tan()).ln()
253}
254
255#[cfg(all(test, feature = "gridpoints-proj"))]
256mod tests {
257    use proj::Proj;
258
259    use super::*;
260
261    const FORWARD_TOLERANCE_METERS: f64 = 1e-7;
262    const INVERSE_TOLERANCE_RADIANS: f64 = 1e-12;
263
264    #[test]
265    fn agrees_with_proj_for_ellipsoidal_secant_cone() {
266        assert_agrees_with_proj(Params {
267            ellipsoid: Ellipsoid::from_a_and_b(6_378_137., 6_356_752.314_245),
268            lat_0: 40.,
269            lon_0: 140.,
270            lat_1: 60.,
271            lat_2: 30.,
272        });
273    }
274
275    #[test]
276    fn agrees_with_proj_for_ellipsoidal_tangent_cone() {
277        assert_agrees_with_proj(Params {
278            ellipsoid: Ellipsoid::from_a_and_b(6_378_137., 6_356_752.314_245),
279            lat_0: 35.,
280            lon_0: -100.,
281            lat_1: 45.,
282            lat_2: 45.,
283        });
284    }
285
286    #[test]
287    fn agrees_with_proj_for_spherical_secant_cone() {
288        assert_agrees_with_proj(Params {
289            ellipsoid: Ellipsoid::from_a_and_b(6_371_229., 6_371_229.),
290            lat_0: -40.,
291            lon_0: 20.,
292            lat_1: -30.,
293            lat_2: -60.,
294        });
295    }
296
297    fn assert_agrees_with_proj(params: Params) {
298        let proj = Proj::new(&params.proj_args()).unwrap();
299        let projection = Projection::new(&params).unwrap();
300        let coordinates: [(f64, f64); 4] = [(-70., -10.), (0., 0.), (25., 45.), (70., 170.)];
301
302        for (lat, lon) in coordinates {
303            let lonlat = (lon.to_radians(), lat.to_radians());
304            let expected_xy = proj.project(lonlat, false).unwrap();
305            let actual_xy = projection.project(&lonlat, false).unwrap();
306            assert_coordinates_close(actual_xy, expected_xy, FORWARD_TOLERANCE_METERS);
307
308            let expected_lonlat = proj.project(expected_xy, true).unwrap();
309            let actual_lonlat = projection.project(&expected_xy, true).unwrap();
310            assert_coordinates_close(actual_lonlat, expected_lonlat, INVERSE_TOLERANCE_RADIANS);
311        }
312    }
313
314    fn assert_coordinates_close(actual: (f64, f64), expected: (f64, f64), tolerance: f64) {
315        assert!(
316            (actual.0 - expected.0).abs() <= tolerance
317                && (actual.1 - expected.1).abs() <= tolerance,
318            "actual {actual:?} differs from expected {expected:?} by more than {tolerance}"
319        );
320    }
321}
grib/projection/lcc.rs

grib/projection/
lcc.rs
