Skip to main content

lora_store/
spatial.rs

1use std::fmt;
2
3/// A 2D or 3D point, either Cartesian or WGS-84 Geographic.
4///
5/// SRID distinguishes coordinate systems:
6///   - SRID 7203: Cartesian 2D
7///   - SRID 9157: Cartesian 3D
8///   - SRID 4326: WGS-84 Geographic 2D
9///   - SRID 4979: WGS-84 Geographic 3D
10///
11/// `z` is `Some` iff the point is 3D. For geographic 3D points `z` holds
12/// the `height` in metres.
13#[derive(Debug, Clone)]
14pub struct LoraPoint {
15    /// First coordinate: x (Cartesian) or longitude (Geographic)
16    pub x: f64,
17    /// Second coordinate: y (Cartesian) or latitude (Geographic)
18    pub y: f64,
19    /// Third coordinate: z (Cartesian) or height in metres (Geographic).
20    /// `None` for 2D points.
21    pub z: Option<f64>,
22    /// Spatial Reference ID.
23    pub srid: u32,
24}
25
26pub const SRID_CARTESIAN: u32 = 7203;
27pub const SRID_CARTESIAN_3D: u32 = 9157;
28pub const SRID_WGS84: u32 = 4326;
29pub const SRID_WGS84_3D: u32 = 4979;
30
31/// Canonical CRS name strings as understood by `point()`.
32pub const CRS_CARTESIAN: &str = "cartesian";
33pub const CRS_CARTESIAN_3D: &str = "cartesian-3D";
34pub const CRS_WGS84_2D: &str = "WGS-84-2D";
35pub const CRS_WGS84_3D: &str = "WGS-84-3D";
36
37impl LoraPoint {
38    pub fn cartesian(x: f64, y: f64) -> Self {
39        Self {
40            x,
41            y,
42            z: None,
43            srid: SRID_CARTESIAN,
44        }
45    }
46
47    pub fn cartesian_3d(x: f64, y: f64, z: f64) -> Self {
48        Self {
49            x,
50            y,
51            z: Some(z),
52            srid: SRID_CARTESIAN_3D,
53        }
54    }
55
56    pub fn geographic(longitude: f64, latitude: f64) -> Self {
57        Self {
58            x: longitude,
59            y: latitude,
60            z: None,
61            srid: SRID_WGS84,
62        }
63    }
64
65    pub fn geographic_3d(longitude: f64, latitude: f64, height: f64) -> Self {
66        Self {
67            x: longitude,
68            y: latitude,
69            z: Some(height),
70            srid: SRID_WGS84_3D,
71        }
72    }
73
74    /// True for WGS-84 points (2D or 3D).
75    pub fn is_geographic(&self) -> bool {
76        self.srid == SRID_WGS84 || self.srid == SRID_WGS84_3D
77    }
78
79    /// True for 3D points (Cartesian or WGS-84).
80    pub fn is_3d(&self) -> bool {
81        self.z.is_some()
82    }
83
84    /// For geographic points, y is latitude.
85    pub fn latitude(&self) -> f64 {
86        self.y
87    }
88
89    /// For geographic points, x is longitude.
90    pub fn longitude(&self) -> f64 {
91        self.x
92    }
93
94    /// For geographic 3D points, z is height in metres.
95    pub fn height(&self) -> Option<f64> {
96        if self.is_geographic() {
97            self.z
98        } else {
99            None
100        }
101    }
102
103    /// Canonical CRS name string for this point's SRID.
104    pub fn crs_name(&self) -> &'static str {
105        match self.srid {
106            SRID_CARTESIAN => CRS_CARTESIAN,
107            SRID_CARTESIAN_3D => CRS_CARTESIAN_3D,
108            SRID_WGS84 => CRS_WGS84_2D,
109            SRID_WGS84_3D => CRS_WGS84_3D,
110            _ => "unknown",
111        }
112    }
113}
114
115impl PartialEq for LoraPoint {
116    fn eq(&self, other: &Self) -> bool {
117        if self.srid != other.srid {
118            return false;
119        }
120        if (self.x - other.x).abs() >= f64::EPSILON {
121            return false;
122        }
123        if (self.y - other.y).abs() >= f64::EPSILON {
124            return false;
125        }
126        match (self.z, other.z) {
127            (None, None) => true,
128            (Some(a), Some(b)) => (a - b).abs() < f64::EPSILON,
129            _ => false,
130        }
131    }
132}
133
134impl Eq for LoraPoint {}
135
136impl fmt::Display for LoraPoint {
137    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
138        match (self.is_geographic(), self.z) {
139            (true, Some(z)) => write!(
140                f,
141                "point({{srid:{}, x:{}, y:{}, z:{}}})",
142                self.srid, self.x, self.y, z
143            ),
144            (true, None) => write!(
145                f,
146                "point({{srid:{}, x:{}, y:{}}})",
147                self.srid, self.x, self.y
148            ),
149            (false, Some(z)) => write!(f, "point({{x:{}, y:{}, z:{}}})", self.x, self.y, z),
150            (false, None) => write!(f, "point({{x:{}, y:{}}})", self.x, self.y),
151        }
152    }
153}
154
155// ---------------------------------------------------------------------------
156// CRS / SRID resolution helpers
157// ---------------------------------------------------------------------------
158
159/// Which coordinate family the caller used.
160#[derive(Debug, Clone, Copy, PartialEq, Eq)]
161pub enum PointKeyFamily {
162    /// `x`/`y` (+ optional `z`).
163    Cartesian,
164    /// `longitude`/`latitude` (+ optional `height`).
165    Geographic,
166}
167
168/// Normalise a CRS name string to its canonical SRID.
169///
170/// Case-insensitive. Accepts "WGS-84" as an alias for the 2D form.
171pub fn srid_from_crs_name(name: &str) -> Option<u32> {
172    let lower = name.to_ascii_lowercase();
173    match lower.as_str() {
174        "cartesian" => Some(SRID_CARTESIAN),
175        "cartesian-3d" => Some(SRID_CARTESIAN_3D),
176        "wgs-84" | "wgs-84-2d" => Some(SRID_WGS84),
177        "wgs-84-3d" => Some(SRID_WGS84_3D),
178        _ => None,
179    }
180}
181
182/// Recognise a bare SRID as one of the four supported values.
183pub fn srid_is_supported(srid: u32) -> bool {
184    matches!(
185        srid,
186        SRID_CARTESIAN | SRID_CARTESIAN_3D | SRID_WGS84 | SRID_WGS84_3D
187    )
188}
189
190pub fn srid_is_3d(srid: u32) -> bool {
191    matches!(srid, SRID_CARTESIAN_3D | SRID_WGS84_3D)
192}
193
194pub fn srid_is_geographic(srid: u32) -> bool {
195    matches!(srid, SRID_WGS84 | SRID_WGS84_3D)
196}
197
198/// Resolve a final SRID from the optional user-provided `crs`, `srid`, plus
199/// the detected key family and dimensionality.
200///
201/// Errors on: unknown CRS name, unsupported SRID, CRS/SRID conflict,
202/// CRS/family conflict (e.g. geographic keys with `crs: "cartesian"`),
203/// CRS/dimensionality conflict (2D coords with a 3D CRS or vice-versa).
204pub fn resolve_srid(
205    crs: Option<&str>,
206    srid: Option<i64>,
207    family: PointKeyFamily,
208    is_3d: bool,
209) -> Result<u32, String> {
210    let crs_srid = match crs {
211        Some(name) => Some(srid_from_crs_name(name).ok_or_else(|| {
212            format!(
213                "point() got unsupported crs '{name}' \
214                 (expected one of cartesian, cartesian-3D, WGS-84, WGS-84-3D)"
215            )
216        })?),
217        None => None,
218    };
219
220    let explicit_srid = match srid {
221        Some(n) => {
222            if n < 0 || n > u32::MAX as i64 {
223                return Err(format!("point() got unsupported srid {n}"));
224            }
225            let n = n as u32;
226            if !srid_is_supported(n) {
227                return Err(format!(
228                    "point() got unsupported srid {n} \
229                     (expected one of 7203, 9157, 4326, 4979)"
230                ));
231            }
232            Some(n)
233        }
234        None => None,
235    };
236
237    let resolved = match (crs_srid, explicit_srid) {
238        (Some(a), Some(b)) if a != b => {
239            return Err(format!(
240                "point() crs '{}' and srid {} do not agree",
241                crs.unwrap(),
242                b
243            ));
244        }
245        (Some(a), _) => Some(a),
246        (None, Some(b)) => Some(b),
247        (None, None) => None,
248    };
249
250    let final_srid = match resolved {
251        Some(s) => {
252            // Family agreement.
253            let srid_geo = srid_is_geographic(s);
254            let family_geo = matches!(family, PointKeyFamily::Geographic);
255            if srid_geo != family_geo {
256                return Err(format!(
257                    "point() coordinates use {} keys but crs/srid is {}",
258                    if family_geo {
259                        "geographic (longitude/latitude)"
260                    } else {
261                        "cartesian (x/y)"
262                    },
263                    if srid_geo { "geographic" } else { "cartesian" }
264                ));
265            }
266            // Dimensionality agreement.
267            if srid_is_3d(s) != is_3d {
268                return Err(format!(
269                    "point() dimensionality mismatch: {} coordinates but {} crs/srid",
270                    if is_3d { "3D" } else { "2D" },
271                    if srid_is_3d(s) { "3D" } else { "2D" }
272                ));
273            }
274            s
275        }
276        None => match (family, is_3d) {
277            (PointKeyFamily::Cartesian, false) => SRID_CARTESIAN,
278            (PointKeyFamily::Cartesian, true) => SRID_CARTESIAN_3D,
279            (PointKeyFamily::Geographic, false) => SRID_WGS84,
280            (PointKeyFamily::Geographic, true) => SRID_WGS84_3D,
281        },
282    };
283
284    Ok(final_srid)
285}
286
287// ---------------------------------------------------------------------------
288// Distance
289// ---------------------------------------------------------------------------
290
291/// Euclidean distance between two Cartesian points (2D or 3D).
292///
293/// Callers are responsible for ensuring both points share the same SRID;
294/// `point_distance` is the usual entry point.
295pub fn cartesian_distance(a: &LoraPoint, b: &LoraPoint) -> f64 {
296    let dx = a.x - b.x;
297    let dy = a.y - b.y;
298    let dz = match (a.z, b.z) {
299        (Some(za), Some(zb)) => za - zb,
300        _ => 0.0,
301    };
302    (dx * dx + dy * dy + dz * dz).sqrt()
303}
304
305/// Haversine distance in metres between two WGS-84 geographic points.
306///
307/// Height is **ignored** even for WGS-84-3D inputs — we compute the
308/// great-circle distance on the reference sphere.
309pub fn haversine_distance(a: &LoraPoint, b: &LoraPoint) -> f64 {
310    const EARTH_RADIUS_M: f64 = 6_371_000.0;
311
312    let lat1 = a.latitude().to_radians();
313    let lat2 = b.latitude().to_radians();
314    let dlat = (b.latitude() - a.latitude()).to_radians();
315    let dlon = (b.longitude() - a.longitude()).to_radians();
316
317    let half = (dlat / 2.0).sin().powi(2) + lat1.cos() * lat2.cos() * (dlon / 2.0).sin().powi(2);
318    let c = 2.0 * half.sqrt().asin();
319    EARTH_RADIUS_M * c
320}
321
322/// Distance between two points — dispatches to Euclidean or Haversine
323/// depending on SRID. Returns `None` if the SRIDs don't match, which
324/// also covers the 2D-vs-3D mismatch since the dimension is baked into
325/// the SRID (7203 vs 9157, 4326 vs 4979).
326pub fn point_distance(a: &LoraPoint, b: &LoraPoint) -> Option<f64> {
327    if a.srid != b.srid {
328        return None;
329    }
330    if a.is_geographic() {
331        Some(haversine_distance(a, b))
332    } else {
333        Some(cartesian_distance(a, b))
334    }
335}
336
337#[cfg(test)]
338mod tests {
339    use super::*;
340
341    #[test]
342    fn resolve_srid_defaults_cartesian_2d() {
343        let s = resolve_srid(None, None, PointKeyFamily::Cartesian, false).unwrap();
344        assert_eq!(s, SRID_CARTESIAN);
345    }
346
347    #[test]
348    fn resolve_srid_defaults_cartesian_3d() {
349        let s = resolve_srid(None, None, PointKeyFamily::Cartesian, true).unwrap();
350        assert_eq!(s, SRID_CARTESIAN_3D);
351    }
352
353    #[test]
354    fn resolve_srid_defaults_wgs84_2d() {
355        let s = resolve_srid(None, None, PointKeyFamily::Geographic, false).unwrap();
356        assert_eq!(s, SRID_WGS84);
357    }
358
359    #[test]
360    fn resolve_srid_defaults_wgs84_3d() {
361        let s = resolve_srid(None, None, PointKeyFamily::Geographic, true).unwrap();
362        assert_eq!(s, SRID_WGS84_3D);
363    }
364
365    #[test]
366    fn resolve_srid_accepts_crs_name_case_insensitive() {
367        let s = resolve_srid(Some("WGS-84"), None, PointKeyFamily::Geographic, false).unwrap();
368        assert_eq!(s, SRID_WGS84);
369        let s = resolve_srid(Some("wgs-84-3d"), None, PointKeyFamily::Geographic, true).unwrap();
370        assert_eq!(s, SRID_WGS84_3D);
371        let s = resolve_srid(Some("CARTESIAN"), None, PointKeyFamily::Cartesian, false).unwrap();
372        assert_eq!(s, SRID_CARTESIAN);
373    }
374
375    #[test]
376    fn resolve_srid_conflict_between_crs_and_srid() {
377        let err = resolve_srid(
378            Some("cartesian"),
379            Some(4326),
380            PointKeyFamily::Cartesian,
381            false,
382        )
383        .unwrap_err();
384        assert!(err.contains("do not agree"));
385    }
386
387    #[test]
388    fn resolve_srid_rejects_unknown_crs() {
389        let err =
390            resolve_srid(Some("mars-centric"), None, PointKeyFamily::Cartesian, false).unwrap_err();
391        assert!(err.contains("unsupported crs"));
392    }
393
394    #[test]
395    fn resolve_srid_rejects_unsupported_srid() {
396        let err = resolve_srid(None, Some(9999), PointKeyFamily::Cartesian, false).unwrap_err();
397        assert!(err.contains("unsupported srid"));
398    }
399
400    #[test]
401    fn resolve_srid_rejects_2d_crs_with_3d_coords() {
402        let err =
403            resolve_srid(Some("cartesian"), None, PointKeyFamily::Cartesian, true).unwrap_err();
404        assert!(err.contains("dimensionality"));
405    }
406
407    #[test]
408    fn resolve_srid_rejects_3d_crs_with_2d_coords() {
409        let err =
410            resolve_srid(Some("WGS-84-3D"), None, PointKeyFamily::Geographic, false).unwrap_err();
411        assert!(err.contains("dimensionality"));
412    }
413
414    #[test]
415    fn resolve_srid_rejects_family_mismatch() {
416        let err =
417            resolve_srid(Some("cartesian"), None, PointKeyFamily::Geographic, false).unwrap_err();
418        assert!(err.contains("coordinates use"));
419    }
420
421    #[test]
422    fn cartesian_3d_distance() {
423        let a = LoraPoint::cartesian_3d(0.0, 0.0, 0.0);
424        let b = LoraPoint::cartesian_3d(1.0, 2.0, 2.0);
425        // sqrt(1 + 4 + 4) = 3
426        assert!((cartesian_distance(&a, &b) - 3.0).abs() < 1e-9);
427    }
428
429    #[test]
430    fn point_distance_returns_none_on_srid_mismatch() {
431        let a = LoraPoint::cartesian(0.0, 0.0);
432        let b = LoraPoint::cartesian_3d(0.0, 0.0, 0.0);
433        assert!(point_distance(&a, &b).is_none());
434    }
435}