maidenhead/
lib.rs

1use thiserror::Error;
2
3#[derive(Error, Debug)]
4pub enum MHError {
5    #[error("Invalid grid format `{0}`")]
6    InvalidGrid(String),
7    #[error("Invalid grid length {0}, only 4/6/8/10 supported")]
8    InvalidGridLength(usize),
9    #[error("Invalid Longitude/Latitude: `{0}`/`{1}`")]
10    InvalidLongLat(f64, f64),
11    #[error("unknown error when generating grid string")]
12    Unknown,
13}
14
15// Grid squares are string representations of the latitude and longitude. A good introduction to how to calculate them is in:
16// http://www.w8bh.net/grid_squares.pdf
17//
18// The format is: FFSSssEEee
19// Field / Square / Subsquare / Extended Square / Superextended Square
20// Each covering for long/lat:
21// Field: 20 / 10 degrees for long / lat
22// Square: 2 / 1 degrees
23// Subsquare: 5 / 2.5 minutes
24// Extended Square: 30 / 15 seconds
25// Superextended Square: 1.25 / 0.625 seconds
26// Note that the enumeration begins at south pole (so 90 degrees off on latitude) and is
27// unsigned positive, so needs to be subtracted by 180 to get +/- longitude.
28
29static LONG_OFFSET: f64 = 180.0;
30static LAT_OFFSET: f64 = 90.0;
31
32static LONG_F: f64 = 20.0;
33static LAT_F: f64 = 10.0;
34static LONG_SQ: f64 = 2.0;
35static LAT_SQ: f64 = 1.0;
36static LONG_SSQ: f64 = 5.0 / 60.0;
37static LAT_SSQ: f64 = 2.5 / 60.0;
38static LONG_ESQ: f64 = 30.0 / 60.0 / 60.0;
39static LAT_ESQ: f64 = 15.0 / 60.0 / 60.0;
40static LONG_SESQ: f64 = 1.25 / 60.0 / 60.0;
41static LAT_SESQ: f64 = 0.625 / 60.0 / 60.0;
42
43static LONG_MULT: [f64; 5] = [LONG_F, LONG_SQ, LONG_SSQ, LONG_ESQ, LONG_SESQ];
44static LAT_MULT: [f64; 5] = [LAT_F, LAT_SQ, LAT_SSQ, LAT_ESQ, LAT_SESQ];
45
46pub fn grid_to_longlat(grid: &str) -> Result<(f64, f64), MHError> {
47    // Validate alpha/digit format
48    // FIXME: Actual values should be A-R 0-9 a-x 0-9 A-X
49    let d = |a: char| a.is_ascii_digit();
50    let l = |a: char| a.is_ascii_alphabetic();
51    let checks = [l, l, d, d, l, l, d, d, l, l];
52    let check = grid
53        .chars()
54        .zip(checks)
55        .map(|(c, lmb)| lmb(c))
56        .collect::<Vec<bool>>();
57
58    // If any of them are false, we've got an invalid grid string
59    if check.iter().filter(|b| !*b).count() != 0 {
60        return Err(MHError::InvalidGrid(grid.to_string()));
61    }
62
63    // Also make sure the length is even (and not 2)
64    match grid.len() {
65        4 | 6 | 8 | 10 => {}
66        l => return Err(MHError::InvalidGridLength(l)),
67    }
68
69    // Now it's just a matter of calculating the offsets from the grid
70    let vals: Vec<u32> = "AA00AA00AA"
71        .chars()
72        .zip(grid.chars())
73        .map(|(t, c)| (c.to_ascii_uppercase() as u32) - (t as u32))
74        .collect();
75
76    // And multiplying each of them with their per-unit value
77    let long: f64 = vals
78        .iter()
79        .step_by(2)
80        .zip(LONG_MULT)
81        .map(|(&v, m)| v as f64 * m)
82        .sum();
83    let lat: f64 = vals
84        .iter()
85        .skip(1)
86        .step_by(2)
87        .zip(LAT_MULT)
88        .map(|(&v, m)| v as f64 * m)
89        .sum();
90
91    // Move the returned value into the middle of the precision given.
92    // This avoids imprecision due to rounding if doing grid->longlat->grid
93    // (We do this in a unit testcase)
94    let idx = grid.len() / 2 - 1;
95    let long = long + LONG_MULT[idx] / 2.0;
96    let lat = lat + LAT_MULT[idx] / 2.0;
97
98    // Finally, adjust for origin offsets
99    Ok((long - LONG_OFFSET, lat - LAT_OFFSET))
100}
101
102pub fn longlat_to_grid(long: f64, lat: f64, precision: usize) -> Result<String, MHError> {
103    let charoff = |base: char, off: u32| std::char::from_u32(base as u32 + off);
104
105    // It only makes sense to have 4+ even number of characters in a grid square
106    match precision {
107        4 | 6 | 8 | 10 => {}
108        p => return Err(MHError::InvalidGridLength(p)),
109    }
110
111    if long > 180.0 || long < -180.0 || lat < -180.0 || lat > 180.0 {
112        return Err(MHError::InvalidLongLat(long, lat));
113    }
114
115    // Do the math to calculate each position, per the w8bh website
116    let long = long + LONG_OFFSET;
117    let lat = lat + LAT_OFFSET;
118    let mut vals = Vec::new();
119    vals.push(long / LONG_F);
120    vals.push(lat / LAT_F);
121    vals.push(long % LONG_F / LONG_SQ);
122    vals.push(lat % LAT_F / LAT_SQ);
123    vals.push(long % LONG_SQ / LONG_SSQ);
124    vals.push(lat % LAT_SQ / LAT_SSQ);
125    vals.push(long % LONG_SSQ / LONG_ESQ);
126    vals.push(lat % LAT_SSQ / LAT_ESQ);
127    vals.push(long % LONG_ESQ / LONG_SESQ);
128    vals.push(lat % LAT_ESQ / LAT_SESQ);
129
130    vals.truncate(precision);
131
132    let grid: Option<String> = "AA00aa00AA"
133        .chars()
134        .zip(vals)
135        .map(|(b, o)| charoff(b, o as u32))
136        .collect();
137    match grid {
138        Some(g) => Ok(g),
139        None => Err(MHError::Unknown),
140    }
141}
142
143// Calculate the distance between two grids, using the haversine
144// formula:
145// a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
146// c = 2 ⋅ atan2( √a, √(1−a) )
147// d = R ⋅ c
148// where:
149//  φ is latitude, λ is longitude, R is earth’s radius (mean radius = 6,371km);
150//  Bearing:
151//  θ = atan2( sin Δλ ⋅ cos φ2 , cos φ1 ⋅ sin φ2 − sin φ1 ⋅ cos φ2 ⋅ cos Δλ )
152
153pub fn grid_dist_bearing(from: &str, to: &str) -> Result<(f64, f64), MHError> {
154    static RADIUS: f64 = 6371.0;
155    let (from_long, from_lat) = grid_to_longlat(from)?;
156    let (to_long, to_lat) = grid_to_longlat(to)?;
157
158    #[allow(non_snake_case)]
159    let Δλ = (to_long - from_long).to_radians();
160    #[allow(non_snake_case)]
161    let Δφ = (to_lat - from_lat).to_radians();
162    let φ1 = from_lat.to_radians();
163    let φ2 = to_lat.to_radians();
164
165    let a: f64 = (Δφ / 2.0).sin().powi(2) + φ1.cos() * φ2.cos() * (Δλ / 2.0).sin().powi(2);
166    let c: f64 = 2.0 * (a.sqrt()).atan2((1.0 - a).sqrt());
167
168    let dist = RADIUS * c;
169    let bearing = (Δλ.sin() * φ2.cos()).atan2(φ1.cos() * φ2.sin() - φ1.sin() * φ2.cos() * Δλ.cos());
170    let bearing = (bearing.to_degrees() + 360.0) % 360.0;
171
172    Ok((dist, bearing))
173}
174
175pub fn grid_distance(from: &str, to: &str) -> Result<f64, MHError> {
176    let (dist, _) = grid_dist_bearing(from, to)?;
177    Ok(dist)
178}
179
180pub fn grid_bearing(from: &str, to: &str) -> Result<f64, MHError> {
181    let (_, bearing) = grid_dist_bearing(from, to)?;
182    Ok(bearing)
183}
184
185#[cfg(test)]
186mod tests {
187    use super::*;
188
189    // From https://stackoverflow.com/questions/30856285/assert-eq-with-floating-point-numbers-and-delta
190    macro_rules! assert_delta {
191        ($x:expr, $y:expr, $d:expr) => {
192            let x = $x as f64;
193            let y = $y as f64;
194            if !((x - y).abs() < $d || (y - x).abs() < $d) {
195                panic!();
196            }
197        };
198    }
199
200    // These values come out of the PDF referenced at the top of this file
201    static TEST_GRID: &str = &"FM18lv53SL";
202    static TEST_LONG: f64 = -77.035278;
203    static TEST_LAT: f64 = 38.889484;
204
205    fn precision_n(n: usize) {
206        let grid = longlat_to_grid(TEST_LONG, TEST_LAT, n).unwrap();
207        let mut check = String::from(TEST_GRID);
208        check.truncate(n);
209        println!("Grid ({}): {}", n, check);
210        assert_eq!(grid, check);
211    }
212
213    #[test]
214    fn precision_10() {
215        precision_n(10);
216    }
217
218    #[test]
219    fn precision_8() {
220        precision_n(8);
221    }
222
223    #[test]
224    fn precision_6() {
225        precision_n(6);
226    }
227
228    #[test]
229    fn precision_4() {
230        precision_n(4);
231    }
232
233    #[test]
234    fn precision_inval() {
235        let grid = longlat_to_grid(TEST_LONG, TEST_LAT, 5);
236        assert!(grid.is_err());
237    }
238
239    #[test]
240    fn precision_inval_lat() {
241        let grid = longlat_to_grid(TEST_LONG, 921.0, 10);
242        assert!(grid.is_err());
243    }
244
245    #[test]
246    fn precision_inval_long() {
247        let grid = longlat_to_grid(-201.0, TEST_LAT, 10);
248        assert!(grid.is_err());
249    }
250
251    fn longlat_n(n: usize) {
252        let mut grid_in = String::from(TEST_GRID);
253        grid_in.truncate(n);
254
255        let ll = grid_to_longlat(&grid_in.as_str());
256        assert!(!ll.is_err());
257
258        // Make sure it's within the margin of error of the smallest field
259        let (long, lat) = ll.unwrap();
260        assert_delta!(long, TEST_LONG, LONG_MULT[n / 2 - 1]);
261        assert_delta!(lat, TEST_LAT, LAT_MULT[n / 2 - 1]);
262
263        // Let's convert it back to grid and compare
264        let grid = longlat_to_grid(long, lat, n).unwrap();
265        assert_eq!(grid_in, grid);
266    }
267
268    #[test]
269    fn longlat10() {
270        longlat_n(10);
271    }
272
273    #[test]
274    fn longlat8() {
275        longlat_n(8);
276    }
277
278    #[test]
279    fn longlat6() {
280        longlat_n(6);
281    }
282
283    #[test]
284    fn longlat4() {
285        longlat_n(4);
286    }
287
288    #[test]
289    fn longlat_invalid() {
290        let ret = grid_to_longlat("AI021");
291        assert!(ret.is_err());
292        let ret = grid_to_longlat("AIA2");
293        assert!(ret.is_err());
294        let ret = grid_to_longlat("🤷I00");
295        assert!(ret.is_err());
296        let ret = grid_to_longlat("AA00AA00AA00");
297        assert!(ret.is_err());
298        let ret = grid_to_longlat("AA00AA00AA");
299        assert!(!ret.is_err());
300    }
301
302    #[test]
303    fn test_distance_null() {
304        let dist = grid_distance(TEST_GRID, TEST_GRID).unwrap();
305        assert_eq!(dist, 0.0);
306    }
307
308    #[test]
309    fn test_distance_home() {
310        let dist = grid_distance("CM87um", "KP04ow").unwrap();
311        let bear = grid_bearing("CM87um", "KP04ow").unwrap();
312        println!("Distance: {} Bearing: {}", dist, bear);
313        println!(
314            "from: {:?} To: {:?}",
315            grid_to_longlat("CM87um"),
316            grid_to_longlat("KP04ow")
317        );
318        assert_delta!(dist, 8189.0, 1.0);
319        assert_delta!(bear, 15.224, 0.001);
320    }
321}