miniproj_ops/ops/
albers_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 AlbersEqualAreaParams {
7    /// longitude of false origin
8    lon_orig: f64,
9    /// latitude of false origin
10    lat_orig: f64,
11    /// latitude of first standard parallel
12    lat_sp1: f64,
13    /// latitude of second standard parallel
14    lat_sp2: f64,
15    /// easting at false origin
16    false_e: f64,
17    /// northing at false origin
18    false_n: f64,
19}
20
21impl AlbersEqualAreaParams {
22    pub const fn new(
23        lon_orig: f64,
24        lat_orig: f64,
25        lat_sp1: f64,
26        lat_sp2: f64,
27        false_e: f64,
28        false_n: f64,
29    ) -> Self {
30        Self {
31            lon_orig,
32            lat_orig,
33            lat_sp1,
34            lat_sp2,
35            false_e,
36            false_n,
37        }
38    }
39
40    /// Get longitude of false origin in radians.
41    pub fn lon_orig(&self) -> f64 {
42        self.lon_orig
43    }
44
45    /// Get latitude of false origin in radians.
46    pub fn lat_orig(&self) -> f64 {
47        self.lat_orig
48    }
49
50    /// Get latitude of first standard parallel in radians.
51    pub fn lat_sp1(&self) -> f64 {
52        self.lat_sp1
53    }
54
55    /// Get latitude of second standard parallel in radians.
56    pub fn lat_sp2(&self) -> f64 {
57        self.lat_sp2
58    }
59
60    /// Get easting at false origin.
61    pub fn false_e(&self) -> f64 {
62        self.false_e
63    }
64
65    /// Get northing at false origin.
66    pub fn false_n(&self) -> f64 {
67        self.false_n
68    }
69}
70
71/// Lambert Azimuthal Equal Area coordinate operation (EPSG:9820)
72#[allow(non_snake_case)]
73#[derive(Copy, Clone, Debug)]
74pub struct AlbersEqualAreaProjection {
75    pub false_e: f64,
76    pub false_n: f64,
77    pub lon_orig: f64,
78    pub ellipsoid_e: f64,
79    pub ellipsoid_e_sq: f64,
80    pub ellipsoid_a: f64,
81    pub C: f64,
82    pub n: f64,
83    pub rho_O: f64,
84    pub beta_fac_sin2: f64,
85    pub beta_fac_sin4: f64,
86    pub beta_fac_sin6: f64,
87}
88
89impl AlbersEqualAreaProjection {
90    #[allow(non_snake_case)]
91    pub fn new(ell: &Ellipsoid, params: &AlbersEqualAreaParams) -> Self {
92        dbg!(ell.e());
93        dbg!(ell.e_squared());
94        let alpha_O = Self::alpha(ell.e_squared(), params.lat_orig(), ell.e());
95        dbg!(alpha_O);
96        let alpha_1 = Self::alpha(ell.e_squared(), params.lat_sp1(), ell.e());
97        dbg!(alpha_1);
98        let alpha_2 = Self::alpha(ell.e_squared(), params.lat_sp2(), ell.e());
99        dbg!(alpha_2);
100        let m1 = params.lat_sp1().cos()
101            / (1f64 - ell.e_squared() * params.lat_sp1().sin().powi(2)).sqrt();
102        dbg!(m1);
103        let m2 = params.lat_sp2().cos()
104            / (1f64 - ell.e_squared() * params.lat_sp2().sin().powi(2)).sqrt();
105        dbg!(m2);
106        let n = (m1.powi(2) - m2.powi(2)) / (alpha_2 - alpha_1);
107        dbg!(n);
108        let C = m1.powi(2) + n * alpha_1;
109        dbg!(C);
110        let rho_O = (ell.a() * (C - n * alpha_O).sqrt()) / n;
111        dbg!(rho_O);
112
113        let beta_fac_sin2 = ell.e_squared() / 3f64
114            + 31f64 * ell.e_squared().powi(2) / 180f64
115            + 517f64 * ell.e_squared().powi(3) / 5040f64;
116        let beta_fac_sin4 =
117            23f64 * ell.e_squared().powi(2) / 360f64 + 251f64 * ell.e_squared().powi(3) / 3708f64;
118        let beta_fac_sin6 = 761f64 * ell.e_squared().powi(3) / 45360f64;
119
120        Self {
121            false_e: params.false_e(),
122            false_n: params.false_n(),
123            lon_orig: params.lon_orig(),
124            ellipsoid_e: ell.e(),
125            ellipsoid_e_sq: ell.e_squared(),
126            ellipsoid_a: ell.a(),
127            n,
128            C,
129            rho_O,
130            beta_fac_sin2,
131            beta_fac_sin4,
132            beta_fac_sin6,
133        }
134    }
135
136    //#[inline]
137    fn alpha(e_sq: f64, phi: f64, e: f64) -> f64 {
138        (1f64 - e_sq)
139            * ((phi.sin() / (1f64 - e_sq * phi.sin().powi(2)))
140                - (1f64 / (2f64 * e)) * ((1f64 - e * phi.sin()) / (1f64 + e * phi.sin())).ln())
141    }
142}
143
144impl crate::traits::Projection for AlbersEqualAreaProjection {
145    /// as per IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – March 2020
146    /// longitude & latitude in radians
147    #[allow(non_snake_case)]
148    fn rad_to_projected(&self, longitude: f64, latitude: f64) -> (f64, f64) {
149        let alpha = Self::alpha(self.ellipsoid_e_sq, latitude, self.ellipsoid_e);
150        dbg!(alpha);
151        let theta = self.n * (longitude - self.lon_orig);
152        dbg!(theta);
153        let rho = (self.ellipsoid_a * (self.C - self.n * alpha).sqrt()) / self.n;
154        dbg!(rho);
155        (
156            self.false_e + (rho * theta.sin()),
157            self.false_n + self.rho_O - (rho * theta.cos()),
158        )
159    }
160
161    /// as per IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – March 2020
162    /// longitude & latitude in radians
163    ///
164    /// The approximation for latitude isn't very precise (6 decimal digits)
165    #[allow(non_snake_case)]
166    fn projected_to_rad(&self, easting: f64, northing: f64) -> (f64, f64) {
167        let theta_: f64 = ((easting - self.false_e) * self.n.signum())
168            .atan2((self.rho_O - (northing - self.false_n)) * self.n.signum());
169        dbg!(theta_);
170        let rho_ = ((easting - self.false_e).powi(2)
171            + (self.rho_O - (northing - self.false_n)).powi(2))
172        .sqrt();
173        dbg!(rho_);
174        let alpha_ = (self.C - (rho_.powi(2) * self.n.powi(2) / self.ellipsoid_a.powi(2))) / self.n;
175        dbg!(alpha_);
176        let beta_ = (alpha_
177            / (1f64
178                - ((1f64 - self.ellipsoid_e_sq) / (2f64 * self.ellipsoid_e))
179                    * ((1f64 - self.ellipsoid_e) / (1f64 + self.ellipsoid_e)).ln()))
180        .asin();
181        dbg!(beta_);
182        let lat = beta_
183            + (2f64 * beta_).sin() * self.beta_fac_sin2
184            + (4f64 * beta_).sin() * self.beta_fac_sin4
185            + (6f64 * beta_).sin() * self.beta_fac_sin6;
186        let lon = self.lon_orig + theta_ / self.n;
187        (lon, lat)
188    }
189}
190
191impl PseudoSerialize for AlbersEqualAreaProjection {
192    fn to_constructed(&self) -> String {
193        format!(
194            r"AlbersEqualAreaProjection{{
195    false_e: {}f64,
196    false_n: {}f64,
197    lon_orig: {}f64,
198    ellipsoid_e: {}f64,
199    ellipsoid_e_sq: {}f64,
200    ellipsoid_a: {}f64,
201    C: {}f64,
202    n: {}f64,
203    rho_O: {}f64,
204    beta_fac_sin2: {}f64,
205    beta_fac_sin4: {}f64,
206    beta_fac_sin6: {}f64
207}}",
208            self.false_e,
209            self.false_n,
210            self.lon_orig,
211            self.ellipsoid_e,
212            self.ellipsoid_e_sq,
213            self.ellipsoid_a,
214            self.C,
215            self.n,
216            self.rho_O,
217            self.beta_fac_sin2,
218            self.beta_fac_sin4,
219            self.beta_fac_sin6
220        )
221    }
222}
223
224impl DbContstruct for AlbersEqualAreaProjection {
225    fn from_database_params(params: &[(u32, f64)], ellipsoid: &Ellipsoid) -> Self {
226        /*
227        ImplementedProjection::new(
228            9820,
229            &[8802, 8801, 8806, 8807],
230            "AlbersEqualAreaParams",
231            "AlbersEqualAreaProjection"
232        )
233        */
234        let params = AlbersEqualAreaParams::new(
235            params
236                .iter()
237                .find_map(|(c, v)| if *c == 8822 { Some(*v) } else { None })
238                .unwrap(),
239            params
240                .iter()
241                .find_map(|(c, v)| if *c == 8821 { Some(*v) } else { None })
242                .unwrap(),
243            params
244                .iter()
245                .find_map(|(c, v)| if *c == 8823 { Some(*v) } else { None })
246                .unwrap(),
247            params
248                .iter()
249                .find_map(|(c, v)| if *c == 8824 { Some(*v) } else { None })
250                .unwrap(),
251            params
252                .iter()
253                .find_map(|(c, v)| if *c == 8826 { Some(*v) } else { None })
254                .unwrap(),
255            params
256                .iter()
257                .find_map(|(c, v)| if *c == 8827 { Some(*v) } else { None })
258                .unwrap(),
259        );
260        Self::new(ellipsoid, &params)
261    }
262}
263
264impl GetterContstruct for AlbersEqualAreaProjection {
265    fn with_db_getter<G>(mut getter: G, ellipsoid: &Ellipsoid) -> Option<Self>
266    where
267        G: FnMut(u32) -> Option<f64>,
268    {
269        let params = AlbersEqualAreaParams::new(
270            getter(8822)?,
271            getter(8821)?,
272            getter(8823)?,
273            getter(8824)?,
274            getter(8826)?,
275            getter(8827)?,
276        );
277        Some(Self::new(ellipsoid, &params))
278    }
279}
280
281pub fn direct_projection(params: &[(u32, f64)], ell: Ellipsoid) -> String {
282    AlbersEqualAreaProjection::from_database_params(params, &ell).to_constructed()
283}
284
285#[cfg(test)]
286mod tests {
287
288    use crate::albers_equal_area::*;
289    use crate::ellipsoid::Ellipsoid;
290    use crate::traits::*;
291
292    // TODO: While passing the round-trip, this test does not match what is given in the EPSG Guidance Note 7-2, May 22.
293    #[test]
294    fn albers_equal_area_consistency_north() {
295        let ell = Ellipsoid::from_a_f_inv(6378137.00, 298.2572221);
296        let params = AlbersEqualAreaParams::new(
297            -1.72787596,
298            0.48578331,
299            0.49538262,
300            0.52854388,
301            1000000.000,
302            1000000.000,
303        );
304
305        let projection = AlbersEqualAreaProjection::new(&ell, &params);
306        let easting_goal = 1466493.492;
307        let northing_goal = 702903.006;
308        let (lon, lat) = projection.projected_to_deg(easting_goal, northing_goal);
309        eprintln!("lon: {lon}");
310        eprintln!("lat: {lat}");
311        let (easting, northing) = projection.deg_to_projected(lon, lat);
312        eprintln!("easting: {easting_goal} - {easting}");
313        eprintln!("northing: {northing_goal} - {northing}");
314
315        assert!((easting - easting_goal).abs() < 0.001);
316
317        assert!((northing - northing_goal).abs() < 0.001);
318    }
319
320    #[test]
321    fn albers_equal_area_consistency_south() {
322        let ell = Ellipsoid::from_a_f_inv(6378160.0, 298.25);
323        let params = AlbersEqualAreaParams::new(
324            -60f64.to_radians(),
325            -32f64.to_radians(),
326            -5f64.to_radians(),
327            -42f64.to_radians(),
328            0.0,
329            0.0,
330        );
331
332        let projection = AlbersEqualAreaProjection::new(&ell, &params);
333        let easting_goal = 1408623.196;
334        let northing_goal = 1507641.482;
335        let (lon, lat) = projection.projected_to_deg(easting_goal, northing_goal);
336        eprintln!("lon: {lon}");
337        eprintln!("lat: {lat}");
338        let (easting, northing) = projection.deg_to_projected(lon, lat);
339        eprintln!("easting: {easting_goal} - {easting}");
340        eprintln!("northing: {northing_goal} - {northing}");
341
342        assert!((easting - easting_goal).abs() < 0.001);
343
344        assert!((northing - northing_goal).abs() < 0.001);
345    }
346}
miniproj_ops/ops/albers_equal_area.rs

miniproj_ops/ops/
albers_equal_area.rs
