geomorph 2.0.2

Simple conversion between different coordinate systems
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193

//! Latitude and longitude coordinates

use crate::datum::Datum;
use crate::math;
use crate::mgrs::Mgrs;
use crate::utm::Utm;

use std::f64::consts;
use std::fmt;

use num_complex::{Complex, Complex64};

/// Holds a pair for latitude and longitude coordinates
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Coord {
    /// Latitude: Must be contained in the interval [-90.0..90.0]
    pub lat: f64,
    /// Longitude: Must be contained in the interval [-180.0..180.0]
    pub lon: f64,
}

impl Coord {
    /// Return a new Coord instance.
    ///
    /// Expected boundaries:
    /// - -90.0 < lat < 90.0
    /// - -180.0 < lon < 180.0
    pub const fn new(lat: f64, lon: f64) -> Coord {
        Coord { lat, lon }
    }

    /// Return a new Coord instance.
    ///
    /// Latitude will be modular 90.0
    /// Longitude will be mobular 180.0
    pub fn new_checked(lat: f64, lon: f64) -> Coord {
        Coord {
            lat: lat % 90.0,
            lon: lon % 180.0,
        }
    }
}

impl fmt::Display for Coord {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        write!(f, "({}, {})", self.lat, self.lon)
    }
}

impl From<Mgrs> for Coord {
    fn from(mgrs: Mgrs) -> Self {
        let utm: Utm = mgrs.into();
        utm.into()
    }
}

impl From<Utm> for Coord {
    fn from(utm: Utm) -> Self {
        if utm.ups {
            return Coord::new(0.0, 0.0);
        }

        let easting = utm.easting;
        let northing = utm.northing;
        let north = utm.north;
        let zone = utm.zone;
        let ups = utm.ups;

        let datum = Datum::WGS84;
        let ind: usize = if ups { 0 } else { 2 } + if north { 1 } else { 0 };
        let real_east: f64 = easting - datum.false_easting[ind];
        let real_north: f64 = northing - datum.false_northing[ind];

        let lon_0: f64 = 6.0 * (zone as f64) - 183.0;
        let mut xi: f64 = real_north / (datum.a1 * datum.k0);
        let mut eta: f64 = real_east / (datum.a1 * datum.k0);

        let xisign: f64 = if xi < 0.0 { -1.0 } else { 1.0 };
        let etasign: f64 = if eta < 0.0 { -1.0 } else { 1.0 };
        xi *= xisign;
        eta *= etasign;

        let backside: bool = xi > consts::PI / 2.0;
        if backside {
            xi = consts::PI - xi;
        }

        let c0: f64 = (2.0 * xi).cos();
        let ch0: f64 = (2.0 * eta).cosh();
        let s0: f64 = (2.0 * xi).sin();
        let sh0: f64 = (2.0 * eta).sinh();

        let mut a: Complex64 = Complex::new(2.0 * c0 * ch0, -2.0 * s0 * sh0);
        let mut n = datum.maxpow;
        let mut y0: Complex64 = Complex::new(if n == 0 { -datum.bet[n] } else { 0.0 }, 0.0);
        let mut y1: Complex64 = Complex::new(0.0, 0.0);
        let mut z0: Complex64 = Complex::new(
            if n == 0 {
                -2.0 * n as f64 * datum.bet[n]
            } else {
                0.0
            },
            0.0,
        );
        let mut z1: Complex64 = Complex::new(0.0, 0.0);

        if n == 0 {
            n -= 1;
        }

        while n > 0 {
            y1 = (a * y0) - (y1) - (datum.bet[n]);
            z1 = (a * z0) - (z1) - (2.0 * (n as f64) * datum.bet[n]);
            n -= 1;
            y0 = (a * y1) - (y0) - (datum.bet[n]);
            z0 = (a * z1) - (z0) - (2.0 * (n as f64) * datum.bet[n]);
            n -= 1;
        }

        a = Complex::new(s0 * ch0, c0 * sh0);
        y1 = Complex::new(xi, eta) + a * y0;

        let xip = y1.re;
        let etap = y1.im;
        let s = etap.sinh();
        let c = xip.cos().max(0.0);
        let r = s.hypot(c);

        let mut rlat: f64;
        let mut rlon: f64;

        if r != 0.0 {
            rlon = s.atan2(c).to_degrees();
            let sxip = xip.sin();
            let tau = math::tauf(sxip / r, datum.es);
            rlat = tau.atan().to_degrees();
        } else {
            rlat = 90.0;
            rlon = 0.0;
        }

        rlat *= xisign;
        if backside {
            rlon = 180.0 - rlon;
        }
        rlon *= etasign;
        rlon = math::angle_normalize(rlon + lon_0);

        let latitude = rlat;
        let longitude = rlon;

        Coord::new(latitude, longitude)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn to_utm() {
        let lat: f64 = 55.722682;
        let lon: f64 = 37.640653;
        let coord = Coord::new(lat, lon);
        let utm: Utm = coord.into();
        assert_eq!(utm.easting.trunc(), 414617.0);
        assert_eq!(utm.northing.trunc(), 6176052.0);
        assert_eq!(utm.north, true);
        assert_eq!(utm.zone, 37);
        assert_eq!(utm.band, 'U');
        assert_eq!(coord.lat, lat);
        assert_eq!(coord.lon, lon);
    }

    #[test]
    fn from_utm() {
        let easting: f64 = 725641.61743212992;
        let northing: f64 = 4911303.2874210617;
        let north: bool = true;
        let zone: i32 = 34;
        let band: char = 'N';
        let ups: bool = false;

        let lat: f64 = 44.319940;
        let lon: f64 = 23.829616;

        let utm: Utm = Utm::new(easting, northing, north, zone, band, ups);

        let coord: Coord = utm.into();
        assert_eq!((coord.lat * 100.0).trunc(), (lat * 100.0).trunc());
        assert_eq!((coord.lon * 100.0).trunc(), (lon * 100.0).trunc());
    }
}