egm2008/
lib.rs

1use crate::geoid::{BAND_COUNT, BAND_SIZE, DATA};
2use thiserror::Error;
3
4mod geoid;
5
6/// Errors indicate what went wrong calculating the geoid height.
7#[derive(Error, Debug, PartialEq)]
8pub enum Error {
9    /// The supplied latitude was not in the expected range.
10    #[error("latitude out of range")]
11    LatitudeOutOfRange,
12    /// The supplied longitude was not in the expected range.
13    #[error("longitude out of range")]
14    LongitudeOutOfRange,
15}
16
17/// Get the geoid height at the specified point on the globe.
18///
19/// The latitude must be in the range [-90, 90] inclusive, and the longitude must be in the range
20/// [-180, 180] inclusive, otherwise an error will be returned.
21///
22/// The result is represented in meters above or below the WGS 84 ellipsoid.
23///
24/// # Examples
25///
26/// ```
27/// # fn main() -> Result<(), egm2008::Error> {
28///     let difference = egm2008::geoid_height(0.0, 0.0)?;
29///     let message = format!("ground level is about {difference} meters above the WGS 84 ellipsoid");
30///     assert_eq!(message, "ground level is about 17.225 meters above the WGS 84 ellipsoid");
31///     Ok(())
32/// }
33/// ```
34pub fn geoid_height(latitude: f32, longitude: f32) -> Result<f32, Error> {
35    let latitude = Latitude::new(latitude)?;
36    let longitude = Longitude::new(longitude)?;
37
38    let (band_low, band_high) = latitude.bounding_bands();
39    let (offset_low, offset_high) = longitude.bounding_offsets();
40
41    // Look up the offsets at each of the four corners of the bounding box.
42    let h00 = lookup_height(band_low, offset_low);
43    let h01 = lookup_height(band_low, offset_high);
44    let h10 = lookup_height(band_high, offset_low);
45    let h11 = lookup_height(band_high, offset_high);
46
47    // Interpolate a height for the actual lat/lon inside the bounding box.
48    let latitude_scalar = latitude.interpolation_scalar();
49    let longitude_scalar = longitude.interpolation_scalar();
50    let latitude_lo = interpolate(h00, h10, latitude_scalar);
51    let latitude_hi = interpolate(h01, h11, latitude_scalar);
52    Ok(interpolate(latitude_lo, latitude_hi, longitude_scalar))
53}
54
55/// Latitude holds an f32 representing a WGS-84 latitude in degrees.
56struct Latitude(f32);
57
58impl Latitude {
59    /// Create a new Latitude from the given f32, returning `Error::LatitudeOutOfRange` if the
60    /// provided value is not between -90 and 90.
61    fn new(latitude: f32) -> Result<Self, Error> {
62        match latitude {
63            latitude if (-90.0..=90.0).contains(&latitude) => Ok(Self(latitude)),
64            _ => Err(Error::LatitudeOutOfRange),
65        }
66    }
67
68    /// Return a translated and scaled version of the latitude for easier calculation of bands.
69    ///
70    /// The value is translated from the [-90, 90] range to the [0, 180] range, and then scaled down
71    /// to the [0, 180/SCALE] range.
72    fn translate_and_scale(&self) -> f32 {
73        (self.0 + 90.0) / geoid::SCALE
74    }
75
76    /// Get a lower and upper band that contains this latitude.
77    fn bounding_bands(&self) -> (usize, usize) {
78        let lat = self.translate_and_scale();
79        let low = lat.floor() as usize;
80        let high = (low + 1).clamp(0, BAND_COUNT - 1);
81        (low, high)
82    }
83
84    /// Figure out the proportion of the distance between the lower and upper bounds.
85    fn interpolation_scalar(&self) -> f32 {
86        let lat = self.translate_and_scale();
87        let (lo, hi) = self.bounding_bands();
88        let (lo, hi) = (lo as f32, hi as f32);
89        if lo == hi {
90            return 0.0;
91        }
92        (lat - lo) / (hi - lo)
93    }
94}
95
96/// Latitude holds an f32 representing a WGS-84 longitude in degrees.
97struct Longitude(f32);
98
99impl Longitude {
100    /// Create a new Longitude from the given f32, returning `Error::Longitude` if the
101    /// provided value is not between -180 and 180.
102    fn new(longitude: f32) -> Result<Self, Error> {
103        match longitude {
104            longitude if (-180.0..=180.0).contains(&longitude) => Ok(Self(longitude)),
105            _ => Err(Error::LongitudeOutOfRange),
106        }
107    }
108
109    /// Return a translated and scaled version of the longitude for easier calculation of offsets.
110    ///
111    /// The value is translated from the [-180, 180] range to the [0, 360] range, and then scaled
112    /// down to the [0, 360/SCALE] range.
113    fn translate_and_scale(&self) -> f32 {
114        (self.0 + 180.0) / geoid::SCALE
115    }
116
117    /// Get a lower and upper offset that contains this longitude.
118    fn bounding_offsets(&self) -> (usize, usize) {
119        let lon = self.translate_and_scale();
120        let low = lon.floor() as usize;
121        let high = (low + 1).clamp(0, BAND_SIZE - 1);
122        (low, high)
123    }
124
125    /// Figure out the proportion of the distance between the lower and upper offsets.
126    fn interpolation_scalar(&self) -> f32 {
127        let lon = self.translate_and_scale();
128        let (lo, hi) = self.bounding_offsets();
129        let (lo, hi) = (lo as f32, hi as f32);
130        if lo == hi {
131            return 0.0;
132        }
133        (lon - lo) / (hi - lo)
134    }
135}
136
137/// Look up the height for the specified band/offset.
138///
139/// # Safety
140///
141/// This could panic if an invalid band/offset are specified. However, this can never happen due to
142/// the validation we perform when constructing [`Latitude`] and [`Longitude`] and how we use the
143/// scale factor generated by the Python script.
144fn lookup_height(band: usize, offset: usize) -> f32 {
145    DATA[band * BAND_SIZE + offset]
146}
147
148/// Given known values at `a` and `b`, construct a new value for a point that is `proportion` of the
149/// way from `a` to `b`.
150///
151/// For example, if we know that F(5) = 1 and F(10) = 2, we can approximate that F(6), a value 20%
152/// of the way between 5 and 10, will be 1.2 (a value 20% of the way between 1 and 2).
153fn interpolate(a: f32, b: f32, proportion: f32) -> f32 {
154    a + ((b - a) * proportion)
155}
156
157#[cfg(test)]
158mod tests {
159    use crate::{geoid_height, interpolate, Error};
160
161    #[test]
162    fn test_interpolate() {
163        assert_eq!(interpolate(0.0, 100.0, 0.5), 50.0);
164        assert_eq!(interpolate(100.0, 0.0, 0.5), 50.0);
165        assert_eq!(interpolate(0.0, 0.0, 0.5), 0.0);
166        assert_eq!(interpolate(5.0, 10.0, 0.2), 6.0);
167    }
168
169    #[test]
170    fn fuzz_inputs() {
171        let mut latitude = -90.0;
172        while latitude <= 90.0 {
173            let mut longitude = -180.0;
174            while longitude <= 180.0 {
175                assert!(
176                    geoid_height(latitude, longitude).is_ok(),
177                    "geoid_height({}, {}) caused a crash",
178                    latitude,
179                    longitude
180                );
181                longitude += 0.09;
182            }
183            latitude += 0.11;
184        }
185    }
186
187    #[test]
188    fn test_one_invalid_input() {
189        let tests = [
190            (-91.0, 0.0, Error::LatitudeOutOfRange),
191            (91.0, 0.0, Error::LatitudeOutOfRange),
192            (0.0, -181.0, Error::LongitudeOutOfRange),
193            (0.0, 181.0, Error::LongitudeOutOfRange),
194        ];
195        for (lat, lon, expected) in tests {
196            let result = geoid_height(lat, lon);
197            assert!(
198                result.is_err(),
199                "Expected an error computing height for ({}, {})",
200                lat,
201                lon
202            );
203            if let Err(actual) = result {
204                assert_eq!(expected, actual, "Got a different error than expected");
205            }
206        }
207    }
208
209    #[test]
210    fn test_two_invalid_inputs() {
211        assert!(geoid_height(-95.0, -200.0).is_err());
212    }
213
214    #[test]
215    fn validate_known_points() {
216        let tests = [
217            (-90.0, -180.0, -30.15),
218            (-90.0, 0.0, -30.15),
219            (-90.0, 180.0, -30.15),
220            (0.0, -180.0, 21.281),
221            (0.0, 0.0, 17.225),
222            (0.0, 180.0, 21.281),
223            (90.0, -180.0, 14.899),
224            (90.0, 0.0, 14.899),
225            (90.0, 180.0, 14.899),
226            (-81.0, -135.0, -45.184),
227            (72.0, 141.0, -1.603),
228            (-87.0, 49.5, -17.6775),
229        ];
230        for (lat, lon, expected) in tests {
231            let result = geoid_height(lat, lon);
232            assert!(
233                result.is_ok(),
234                "could not calculate geoid height at ({}, {})",
235                lat,
236                lon
237            );
238            if let Ok(actual) = result {
239                assert_eq!(
240                    actual, expected,
241                    "got unexpected height at ({}, {})",
242                    lat, lon
243                );
244            }
245        }
246    }
247}
egm2008/lib.rs

egm2008/
lib.rs
