miniproj_ops/ops/
lambert_azimuthal_equal_area.rs

1//This file is licensed under EUPL v1.2 as part of the Digital Earth Viewer
2
3use crate::{ellipsoid::Ellipsoid, traits::GetterContstruct, DbContstruct, PseudoSerialize};
4
5#[derive(Copy, Clone, Debug)]
6pub struct LambertAzimuthalEqualAreaParams {
7    /// longitude of natural origin
8    lon_orig: f64,
9    /// latitude of natural origin
10    lat_orig: f64,
11    /// false easting
12    false_e: f64,
13    /// false northing
14    false_n: f64,
15}
16
17impl LambertAzimuthalEqualAreaParams {
18    pub const fn new(lon_orig: f64, lat_orig: f64, false_e: f64, false_n: f64) -> Self {
19        Self {
20            lat_orig,
21            lon_orig,
22            false_e,
23            false_n,
24        }
25    }
26
27    /// Get longitude of natural origin in radians.
28    pub fn lon_orig(&self) -> f64 {
29        self.lon_orig
30    }
31
32    /// Get latitude of natural origin in radians.
33    pub fn lat_orig(&self) -> f64 {
34        self.lat_orig
35    }
36
37    /// Get false easting.
38    pub fn false_e(&self) -> f64 {
39        self.false_e
40    }
41
42    /// Get false northing.
43    pub fn false_n(&self) -> f64 {
44        self.false_n
45    }
46}
47
48/// Lambert Azimuthal Equal Area coordinate operation (EPSG:9820)
49#[allow(non_snake_case)]
50#[derive(Copy, Clone, Debug)]
51pub struct LambertAzimuthalEqualAreaProjection {
52    pub lon_orig: f64,
53    pub false_e: f64,
54    pub false_n: f64,
55    pub ellipsoid_e: f64,
56    pub ellipsoid_e_squared: f64,
57
58    //q_O: f64,
59    pub q_P: f64,
60    pub beta_O: f64,
61    pub R_q: f64,
62    pub D: f64,
63}
64
65impl LambertAzimuthalEqualAreaProjection {
66    #[allow(non_snake_case)]
67    pub fn new(ell: &Ellipsoid, params: &LambertAzimuthalEqualAreaParams) -> Self {
68        let q_P = (1.0 - ell.e_squared())
69            * ((1.0 / (1.0 - ell.e_squared()))
70                - ((0.5 / ell.e()) * f64::ln((1.0 - ell.e()) / (1.0 + ell.e()))));
71
72        let q_O = (1.0 - ell.e_squared())
73            * ((params.lat_orig().sin()
74                / (1.0 - ell.e_squared() * params.lat_orig().sin().powi(2)))
75                - ((0.5 / ell.e())
76                    * f64::ln(
77                        (1.0 - ell.e() * params.lat_orig().sin())
78                            / (1.0 + ell.e() * params.lat_orig().sin()),
79                    )));
80
81        let beta_O = (q_O / q_P).asin();
82
83        let R_q = ell.a() * (q_P / 2.0).sqrt();
84
85        let D = ell.a()
86            * (params.lat_orig().cos()
87                / (1.0 - ell.e_squared() * params.lat_orig().sin().powi(2)).sqrt())
88            / (R_q * beta_O.cos());
89
90        Self {
91            lon_orig: params.lon_orig(),
92            false_e: params.false_e(),
93            false_n: params.false_n(),
94            ellipsoid_e: ell.e(),
95            ellipsoid_e_squared: ell.e_squared(),
96
97            q_P,
98            //q_O,
99            beta_O,
100            R_q,
101            D,
102        }
103    }
104}
105
106impl crate::traits::Projection for LambertAzimuthalEqualAreaProjection {
107    /// as per IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – March 2020
108    /// longitude & latitude in radians
109    #[allow(non_snake_case)]
110    fn rad_to_projected(&self, longitude: f64, latitude: f64) -> (f64, f64) {
111        let q = (1.0 - self.ellipsoid_e_squared)
112            * ((latitude.sin() / (1.0 - self.ellipsoid_e_squared * latitude.sin().powi(2)))
113                - ((0.5 / self.ellipsoid_e)
114                    * f64::ln(
115                        (1.0 - self.ellipsoid_e * latitude.sin())
116                            / (1.0 + self.ellipsoid_e * latitude.sin()),
117                    )));
118
119        let beta = (q / self.q_P).asin();
120
121        let B = self.R_q
122            * (2.0
123                / (1.0
124                    + self.beta_O.sin() * beta.sin()
125                    + (self.beta_O.cos() * beta.cos() * (longitude - self.lon_orig).cos())))
126            .sqrt();
127
128        (
129            self.false_e + ((B * self.D) * (beta.cos() * (longitude - self.lon_orig).sin())),
130            self.false_n
131                + (B / self.D)
132                    * ((self.beta_O.cos() * beta.sin())
133                        - (self.beta_O.sin() * beta.cos() * (longitude - self.lon_orig).cos())),
134        )
135    }
136
137    /// as per IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – March 2020
138    /// longitude & latitude in radians
139    ///
140    /// The approximation for latitude isn't very precise (6 decimal digits)
141    #[allow(non_snake_case)]
142    fn projected_to_rad(&self, easting: f64, northing: f64) -> (f64, f64) {
143        let rho = (((easting - self.false_e) / self.D).powi(2)
144            + (self.D * (northing - self.false_n)).powi(2))
145        .sqrt();
146
147        let C = 2.0 * (rho / 2.0 / self.R_q).asin();
148
149        let beta_ = ((C.cos() * self.beta_O.sin())
150            + ((self.D * (northing - self.false_n) * C.sin() * self.beta_O.cos()) / rho))
151            .asin();
152
153        (
154            self.lon_orig
155                + f64::atan2(
156                    (easting - self.false_e) * C.sin(),
157                    self.D * rho * self.beta_O.cos() * C.cos()
158                        - self.D.powi(2) * (northing - self.false_n) * self.beta_O.sin() * C.sin(),
159                ),
160            beta_
161                + ((self.ellipsoid_e_squared / 3.0
162                    + (31.0 / 180.0) * self.ellipsoid_e_squared.powi(2)
163                    + (517.0 / 5040.0) * self.ellipsoid_e_squared.powi(3))
164                    * (beta_ * 2.0).sin()
165                    + ((23.0 / 360.0) * self.ellipsoid_e_squared.powi(2)
166                        + (251.0 / 3780.0) * self.ellipsoid_e_squared.powi(3))
167                        * (beta_ * 4.0).sin()
168                    + (761.0 / 45360.0) * self.ellipsoid_e_squared.powi(3) * (beta_ + 6.0).sin()),
169        )
170    }
171}
172
173impl PseudoSerialize for LambertAzimuthalEqualAreaProjection {
174    fn to_constructed(&self) -> String {
175        format!(
176            r"LambertAzimuthalEqualAreaProjection{{
177    lon_orig: {}f64,
178    false_e: {}f64,
179    false_n: {}f64,
180    ellipsoid_e: {}f64,
181    ellipsoid_e_squared: {}f64,
182
183    q_P: {}f64,
184    beta_O: {}f64,
185    R_q: {}f64,
186    D: {}f64,
187}}",
188            self.lon_orig,
189            self.false_e,
190            self.false_n,
191            self.ellipsoid_e,
192            self.ellipsoid_e_squared,
193            self.q_P,
194            self.beta_O,
195            self.R_q,
196            self.D
197        )
198    }
199}
200
201impl DbContstruct for LambertAzimuthalEqualAreaProjection {
202    fn from_database_params(params: &[(u32, f64)], ellipsoid: &Ellipsoid) -> Self {
203        /*
204        ImplementedProjection::new(
205            9820,
206            &[8802, 8801, 8806, 8807],
207            "LambertAzimuthalEqualAreaParams",
208            "LambertAzimuthalEqualAreaProjection"
209        )
210        */
211        let params = LambertAzimuthalEqualAreaParams::new(
212            params
213                .iter()
214                .find_map(|(c, v)| if *c == 8802 { Some(*v) } else { None })
215                .unwrap(),
216            params
217                .iter()
218                .find_map(|(c, v)| if *c == 8801 { Some(*v) } else { None })
219                .unwrap(),
220            params
221                .iter()
222                .find_map(|(c, v)| if *c == 8806 { Some(*v) } else { None })
223                .unwrap(),
224            params
225                .iter()
226                .find_map(|(c, v)| if *c == 8807 { Some(*v) } else { None })
227                .unwrap(),
228        );
229        Self::new(ellipsoid, &params)
230    }
231}
232
233pub fn direct_projection(params: &[(u32, f64)], ell: Ellipsoid) -> String {
234    LambertAzimuthalEqualAreaProjection::from_database_params(params, &ell).to_constructed()
235}
236
237impl GetterContstruct for LambertAzimuthalEqualAreaProjection {
238    fn with_db_getter<G>(mut getter: G, ellipsoid: &Ellipsoid) -> Option<Self>
239    where
240        G: FnMut(u32) -> Option<f64>,
241    {
242        let params = LambertAzimuthalEqualAreaParams::new(
243            getter(8802)?,
244            getter(8801)?,
245            getter(8806)?,
246            getter(8807)?,
247        );
248        Some(Self::new(ellipsoid, &params))
249    }
250}
251
252#[cfg(test)]
253mod tests {
254
255    use crate::ellipsoid::Ellipsoid;
256    use crate::lambert_azimuthal_equal_area::*;
257    use crate::traits::*;
258
259    #[test]
260    fn lambert_azimuthal_equal_area_consistency() {
261        let ell = Ellipsoid::from_a_f_inv(6378137.0, 298.2572221);
262        let params = LambertAzimuthalEqualAreaParams::new(
263            10.0f64.to_radians(),
264            52.0f64.to_radians(),
265            4_321_000.0,
266            3_210_000.0,
267        );
268
269        let projection = LambertAzimuthalEqualAreaProjection::new(&ell, &params);
270        let easting_goal = 3962799.45;
271        let northing_goal = 2999718.85;
272        let (lon, lat) = projection.projected_to_deg(easting_goal, northing_goal);
273        let (easting, northing) = projection.deg_to_projected(lon, lat);
274        eprintln!("easting: {easting_goal} - {easting}");
275        eprintln!("northing: {northing_goal} - {northing}");
276
277        assert!((easting - easting_goal).abs() < 0.01);
278
279        assert!((northing - northing_goal).abs() < 0.05);
280    }
281}
miniproj_ops/ops/lambert_azimuthal_equal_area.rs

miniproj_ops/ops/
lambert_azimuthal_equal_area.rs
