geomorph/
datum.rs

1//! Datum conventions
2
3use crate::math;
4
5/// Holds conventional datum information
6#[derive(Debug, Clone, PartialEq)]
7pub struct Datum {
8    pub a: f64,
9    pub f: f64,
10    pub k0: f64,
11    pub e2: f64,
12    pub es: f64,
13    pub e2m: f64,
14    pub b1: f64,
15    pub a1: f64,
16    pub c: f64,
17    pub n: f64,
18    pub maxpow: usize,
19    pub alp: [f64; 7],
20    pub bet: [f64; 7],
21    pub false_easting: [f64; 4],
22    pub false_northing: [f64; 4],
23}
24
25impl Datum {
26    /// A WGS84 instance.
27    pub const WGS84: Self = Self {
28        a: 6378137.0,
29        f: 0.0033528106647474805,
30        k0: 0.9996,
31        e2: 0.0066943799901413165,
32        es: 0.08181919084262149,
33        e2m: 0.9933056200098587,
34        b1: 0.9983242984312527,
35        a1: 6367449.145823414,
36        c: 1.0033565552493153,
37        n: 0.0016792203863837047,
38        maxpow: 6,
39        alp: [
40            0.0,
41            0.0008377318206244698,
42            7.608527773572309e-7,
43            1.1976455033294525e-9,
44            2.429170607201359e-12,
45            5.711757677865804e-15,
46            1.4911177312583895e-17,
47        ],
48        bet: [
49            0.0,
50            0.0008377321640579486,
51            5.9058701522202026e-8,
52            1.6734826652839968e-10,
53            2.1647980400627056e-13,
54            3.787978046168606e-16,
55            7.2487488906941545e-19,
56        ],
57        false_easting: [2000000.0, 2000000.0, 500000.0, 500000.0],
58        false_northing: [2000000.0, 2000000.0, 10000000.0, 0.0],
59    };
60
61    /// Return a new Datum instance.
62    pub fn new(
63        a: f64,
64        f: f64,
65        k0: f64,
66        alpcoeff: &[f64],
67        betcoeff: &[f64],
68        b1coeff: &[f64],
69    ) -> Datum {
70        let e2 = f * (2.0 - f);
71        let es = if f <= 0.0 {
72            -e2.abs().sqrt()
73        } else {
74            e2.abs().sqrt()
75        };
76        let e2m: f64 = 1.0 - e2;
77        let c: f64 = e2m.sqrt() * math::eatanhe(1.0, es).exp();
78        let n: f64 = f / (2.0 - f);
79        let maxpow: usize = 6;
80
81        // TODO when supporting multiple datums, this should be `maxpow + 1`
82        let mut alp = [0.0; 7];
83        let mut bet = [0.0; 7];
84
85        let false_easting = [2000000.0, 2000000.0, 500000.0, 500000.0];
86        let false_northing = [2000000.0, 2000000.0, 10000000.0, 0.0];
87
88        let m = maxpow / 2;
89        let b1: f64 = math::polyval(m, b1coeff, n.powi(2)) / (b1coeff[m + 1] * (1.0 + n));
90        let a1: f64 = b1 * a;
91
92        let mut o: usize = 0;
93        let mut d: f64 = n;
94
95        for i in 0..maxpow {
96            let m = maxpow - i - 1;
97            alp[i + 1] = d * math::polyval(m, &alpcoeff[o..], n) / alpcoeff[o + m + 1];
98            bet[i + 1] = d * math::polyval(m, &betcoeff[o..], n) / betcoeff[o + m + 1];
99            o += m + 2;
100            d *= n;
101        }
102
103        Datum {
104            a,
105            f,
106            k0,
107            e2,
108            es,
109            e2m,
110            b1,
111            a1,
112            c,
113            n,
114            maxpow,
115            alp,
116            bet,
117            false_easting,
118            false_northing,
119        }
120    }
121}
122
123#[cfg(test)]
124mod tests {
125    use super::*;
126
127    #[test]
128    fn validate_wgs84() {
129        assert_eq!(
130            Datum::WGS84,
131            Datum::new(
132                6378137.0,
133                0.0033528106647474805,
134                0.99960000000000004,
135                &[
136                    31564.0,
137                    -66675.0,
138                    34440.0,
139                    47250.0,
140                    -100800.0,
141                    75600.0,
142                    151200.0,
143                    -1983433.0,
144                    863232.0,
145                    748608.0,
146                    -1161216.0,
147                    524160.0,
148                    1935360.0,
149                    670412.0,
150                    406647.0,
151                    -533952.0,
152                    184464.0,
153                    725760.0,
154                    6601661.0,
155                    -7732800.0,
156                    2230245.0,
157                    7257600.0,
158                    -13675556.0,
159                    3438171.0,
160                    7983360.0,
161                    212378941.0,
162                    319334400.0,
163                ],
164                &[
165                    384796.0,
166                    -382725.0,
167                    -6720.0,
168                    932400.0,
169                    -1612800.0,
170                    1209600.0,
171                    2419200.0,
172                    -1118711.0,
173                    1695744.0,
174                    -1174656.0,
175                    258048.0,
176                    80640.0,
177                    3870720.0,
178                    22276.0,
179                    -16929.0,
180                    -15984.0,
181                    12852.0,
182                    362880.0,
183                    -830251.0,
184                    -158400.0,
185                    197865.0,
186                    7257600.0,
187                    -435388.0,
188                    453717.0,
189                    15966720.0,
190                    20648693.0,
191                    638668800.0,
192                ],
193                &[1.0, 4.0, 64.0, 256.0, 256.0],
194            )
195        );
196
197        assert_eq!((Datum::WGS84.n * 100000000.0).trunc(), 167922.0);
198    }
199}
geomorph/datum.rs

geomorph/
datum.rs
