nodedb_types/
geometry.rs

1// SPDX-License-Identifier: Apache-2.0
2
3//! GeoJSON-compatible geometry types.
4//!
5//! Supports Point, LineString, Polygon, MultiPoint, MultiLineString,
6//! MultiPolygon, and GeometryCollection. Stored as GeoJSON for JSON
7//! compatibility. Includes distance (Haversine), area, bearing, and
8//! centroid calculations.
9
10use serde::{Deserialize, Serialize};
11
12/// A 2D coordinate (longitude, latitude) following GeoJSON convention.
13/// Note: GeoJSON uses [lng, lat] order, NOT [lat, lng].
14#[derive(
15    Debug,
16    Clone,
17    Copy,
18    PartialEq,
19    Serialize,
20    Deserialize,
21    zerompk::ToMessagePack,
22    zerompk::FromMessagePack,
23)]
24pub struct Coord {
25    pub lng: f64,
26    pub lat: f64,
27}
28
29impl Coord {
30    pub fn new(lng: f64, lat: f64) -> Self {
31        Self { lng, lat }
32    }
33}
34
35/// GeoJSON-compatible geometry types.
36#[derive(
37    Debug,
38    Clone,
39    PartialEq,
40    Serialize,
41    Deserialize,
42    zerompk::ToMessagePack,
43    zerompk::FromMessagePack,
44)]
45#[serde(tag = "type")]
46#[non_exhaustive]
47pub enum Geometry {
48    Point {
49        coordinates: [f64; 2],
50    },
51    LineString {
52        coordinates: Vec<[f64; 2]>,
53    },
54    Polygon {
55        coordinates: Vec<Vec<[f64; 2]>>,
56    },
57    MultiPoint {
58        coordinates: Vec<[f64; 2]>,
59    },
60    MultiLineString {
61        coordinates: Vec<Vec<[f64; 2]>>,
62    },
63    MultiPolygon {
64        coordinates: Vec<Vec<Vec<[f64; 2]>>>,
65    },
66    GeometryCollection {
67        geometries: Vec<Geometry>,
68    },
69}
70
71impl Geometry {
72    /// Create a Point from (longitude, latitude).
73    pub fn point(lng: f64, lat: f64) -> Self {
74        Geometry::Point {
75            coordinates: [lng, lat],
76        }
77    }
78
79    /// Create a LineString from a series of [lng, lat] pairs.
80    pub fn line_string(coords: Vec<[f64; 2]>) -> Self {
81        Geometry::LineString {
82            coordinates: coords,
83        }
84    }
85
86    /// Create a Polygon from exterior ring (and optional holes).
87    ///
88    /// The first ring is the exterior, subsequent rings are holes.
89    /// Each ring must be a closed loop (first point == last point).
90    pub fn polygon(rings: Vec<Vec<[f64; 2]>>) -> Self {
91        Geometry::Polygon { coordinates: rings }
92    }
93
94    /// Get the type name of this geometry.
95    pub fn geometry_type(&self) -> &'static str {
96        match self {
97            Geometry::Point { .. } => "Point",
98            Geometry::LineString { .. } => "LineString",
99            Geometry::Polygon { .. } => "Polygon",
100            Geometry::MultiPoint { .. } => "MultiPoint",
101            Geometry::MultiLineString { .. } => "MultiLineString",
102            Geometry::MultiPolygon { .. } => "MultiPolygon",
103            Geometry::GeometryCollection { .. } => "GeometryCollection",
104        }
105    }
106
107    /// Compute the centroid of the geometry.
108    pub fn centroid(&self) -> Option<[f64; 2]> {
109        match self {
110            Geometry::Point { coordinates } => Some(*coordinates),
111            Geometry::LineString { coordinates } => {
112                if coordinates.is_empty() {
113                    return None;
114                }
115                let n = coordinates.len() as f64;
116                let sum_lng: f64 = coordinates.iter().map(|c| c[0]).sum();
117                let sum_lat: f64 = coordinates.iter().map(|c| c[1]).sum();
118                Some([sum_lng / n, sum_lat / n])
119            }
120            Geometry::Polygon { coordinates } => {
121                // Centroid of exterior ring.
122                coordinates.first().and_then(|ring| {
123                    if ring.is_empty() {
124                        return None;
125                    }
126                    let n = ring.len() as f64;
127                    let sum_lng: f64 = ring.iter().map(|c| c[0]).sum();
128                    let sum_lat: f64 = ring.iter().map(|c| c[1]).sum();
129                    Some([sum_lng / n, sum_lat / n])
130                })
131            }
132            Geometry::MultiPoint { coordinates } => {
133                if coordinates.is_empty() {
134                    return None;
135                }
136                let n = coordinates.len() as f64;
137                let sum_lng: f64 = coordinates.iter().map(|c| c[0]).sum();
138                let sum_lat: f64 = coordinates.iter().map(|c| c[1]).sum();
139                Some([sum_lng / n, sum_lat / n])
140            }
141            _ => None,
142        }
143    }
144}
145
146/// Parse a GeoJSON string into a [`Geometry`].
147///
148/// Shared by every storage/read path that may encounter geometry stored as a
149/// JSON string rather than a native object — SQL `ST_Point(...)` inserts
150/// serialize to a GeoJSON string, while schemaless document writes and
151/// `Value::Geometry` keep the native form. Centralizing the string-parse core
152/// here keeps the three call sites (document index build, spatial read path,
153/// columnar geometry index) from drifting on which parser they use.
154///
155/// Uses `sonic_rs` per workspace policy (never `serde_json::from_str` for
156/// runtime JSON parsing). Returns `None` on malformed/non-geometry JSON.
157pub fn from_geojson_str(s: &str) -> Option<Geometry> {
158    sonic_rs::from_str(s).ok()
159}
160
161// ── Geo math functions ──
162
163const EARTH_RADIUS_M: f64 = 6_371_000.0;
164
165/// Haversine distance between two points in meters.
166///
167/// Input: (lng1, lat1) and (lng2, lat2) in degrees.
168pub fn haversine_distance(lng1: f64, lat1: f64, lng2: f64, lat2: f64) -> f64 {
169    let lat1_r = lat1.to_radians();
170    let lat2_r = lat2.to_radians();
171    let dlat = (lat2 - lat1).to_radians();
172    let dlng = (lng2 - lng1).to_radians();
173
174    let a = (dlat / 2.0).sin().powi(2) + lat1_r.cos() * lat2_r.cos() * (dlng / 2.0).sin().powi(2);
175    let c = 2.0 * a.sqrt().atan2((1.0 - a).sqrt());
176    EARTH_RADIUS_M * c
177}
178
179/// Haversine bearing from point A to point B in degrees (0-360).
180pub fn haversine_bearing(lng1: f64, lat1: f64, lng2: f64, lat2: f64) -> f64 {
181    let lat1_r = lat1.to_radians();
182    let lat2_r = lat2.to_radians();
183    let dlng = (lng2 - lng1).to_radians();
184
185    let y = dlng.sin() * lat2_r.cos();
186    let x = lat1_r.cos() * lat2_r.sin() - lat1_r.sin() * lat2_r.cos() * dlng.cos();
187    let bearing = y.atan2(x).to_degrees();
188    (bearing + 360.0) % 360.0
189}
190
191/// Approximate area of a polygon on Earth's surface in square meters.
192///
193/// Uses the Shoelace formula on projected coordinates (simple equirectangular).
194/// Accurate for small polygons; for large polygons use spherical excess.
195pub fn polygon_area(ring: &[[f64; 2]]) -> f64 {
196    if ring.len() < 3 {
197        return 0.0;
198    }
199    let mut sum = 0.0;
200    let n = ring.len();
201    for i in 0..n {
202        let j = (i + 1) % n;
203        // Convert to approximate meters from center.
204        let lat_avg = ((ring[i][1] + ring[j][1]) / 2.0).to_radians();
205        let x1 = ring[i][0].to_radians() * EARTH_RADIUS_M * lat_avg.cos();
206        let y1 = ring[i][1].to_radians() * EARTH_RADIUS_M;
207        let x2 = ring[j][0].to_radians() * EARTH_RADIUS_M * lat_avg.cos();
208        let y2 = ring[j][1].to_radians() * EARTH_RADIUS_M;
209        sum += x1 * y2 - x2 * y1;
210    }
211    (sum / 2.0).abs()
212}
213
214/// Check if a point is inside a polygon (ray casting algorithm).
215pub fn point_in_polygon(lng: f64, lat: f64, ring: &[[f64; 2]]) -> bool {
216    let mut inside = false;
217    let n = ring.len();
218    let mut j = n.wrapping_sub(1);
219    for i in 0..n {
220        let yi = ring[i][1];
221        let yj = ring[j][1];
222        if ((yi > lat) != (yj > lat))
223            && (lng < (ring[j][0] - ring[i][0]) * (lat - yi) / (yj - yi) + ring[i][0])
224        {
225            inside = !inside;
226        }
227        j = i;
228    }
229    inside
230}
231
232#[cfg(test)]
233mod tests {
234    use super::*;
235
236    #[test]
237    fn point_creation() {
238        let p = Geometry::point(-73.9857, 40.7484);
239        assert_eq!(p.geometry_type(), "Point");
240        if let Geometry::Point { coordinates } = &p {
241            assert!((coordinates[0] - (-73.9857)).abs() < 1e-6);
242            assert!((coordinates[1] - 40.7484).abs() < 1e-6);
243        }
244    }
245
246    #[test]
247    fn haversine_nyc_to_london() {
248        // NYC: -74.006, 40.7128 → London: -0.1278, 51.5074
249        let d = haversine_distance(-74.006, 40.7128, -0.1278, 51.5074);
250        // ~5,570 km
251        assert!((d - 5_570_000.0).abs() < 50_000.0, "got {d}m");
252    }
253
254    #[test]
255    fn haversine_same_point() {
256        let d = haversine_distance(0.0, 0.0, 0.0, 0.0);
257        assert!(d.abs() < 1e-6);
258    }
259
260    #[test]
261    fn bearing_north() {
262        let b = haversine_bearing(0.0, 0.0, 0.0, 1.0);
263        assert!((b - 0.0).abs() < 1.0, "expected ~0, got {b}");
264    }
265
266    #[test]
267    fn bearing_east() {
268        let b = haversine_bearing(0.0, 0.0, 1.0, 0.0);
269        assert!((b - 90.0).abs() < 1.0, "expected ~90, got {b}");
270    }
271
272    #[test]
273    fn polygon_area_simple() {
274        // ~1 degree square near equator ≈ 111km × 111km ≈ 12,321 km²
275        let ring = vec![[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.0, 0.0]];
276        let area = polygon_area(&ring);
277        let area_km2 = area / 1e6;
278        assert!(
279            (area_km2 - 12321.0).abs() < 500.0,
280            "expected ~12321 km², got {area_km2}"
281        );
282    }
283
284    #[test]
285    fn point_in_polygon_inside() {
286        let ring = vec![
287            [0.0, 0.0],
288            [10.0, 0.0],
289            [10.0, 10.0],
290            [0.0, 10.0],
291            [0.0, 0.0],
292        ];
293        assert!(point_in_polygon(5.0, 5.0, &ring));
294        assert!(!point_in_polygon(15.0, 5.0, &ring));
295    }
296
297    #[test]
298    fn centroid_point() {
299        let p = Geometry::point(10.0, 20.0);
300        assert_eq!(p.centroid(), Some([10.0, 20.0]));
301    }
302
303    #[test]
304    fn centroid_linestring() {
305        let ls = Geometry::line_string(vec![[0.0, 0.0], [10.0, 0.0], [10.0, 10.0]]);
306        let c = ls.centroid().unwrap();
307        assert!((c[0] - 6.6667).abs() < 0.01);
308        assert!((c[1] - 3.3333).abs() < 0.01);
309    }
310
311    #[test]
312    fn geojson_serialize() {
313        let p = Geometry::point(1.0, 2.0);
314        let json = sonic_rs::to_string(&p).unwrap();
315        assert!(json.contains("\"type\":\"Point\""));
316        assert!(json.contains("\"coordinates\":[1.0,2.0]"));
317    }
318
319    #[test]
320    fn geojson_roundtrip() {
321        let original = Geometry::polygon(vec![vec![
322            [0.0, 0.0],
323            [1.0, 0.0],
324            [1.0, 1.0],
325            [0.0, 1.0],
326            [0.0, 0.0],
327        ]]);
328        let json = sonic_rs::to_string(&original).unwrap();
329        let parsed: Geometry = sonic_rs::from_str(&json).unwrap();
330        assert_eq!(original, parsed);
331    }
332}
nodedb_types/geometry.rs

nodedb_types/
geometry.rs
