qualia_core_db/sparql_library/
geosparql.rs

1//! GeoSPARQL geometry support: a WKT (Well-Known Text) literal parser and the
2//! real geometry predicates the SPARQL extension functions dispatch to
3//! (`geof:distance`, `geof:sfContains`, `sfWithin`, `sfIntersects`, `sfTouches`).
4//!
5//! This replaces the "Simplified" placeholders in `sparql_extensions.rs` (which
6//! returned hardcoded `true`/`false` or an arbitrary threshold) with genuine
7//! computation over parsed geometry. Distances use the haversine great-circle
8//! formula (lon/lat degrees, WGS-84 mean radius); topological predicates use
9//! planar tests (ray-casting point-in-polygon, segment intersection).
10//!
11//! Scope: `POINT`, `LINESTRING`, `POLYGON` (with holes), and their `MULTI`
12//! variants — the WKT subset GeoSPARQL data uses in practice. Z/M coordinates
13//! are parsed and ignored (2D predicates). An optional SRID/`<uri>` prefix
14//! (`geo:wktLiteral` values sometimes carry `<crs> POINT(...)`) is skipped.
15
16/// Mean Earth radius (WGS-84), metres — used by the haversine distance.
17const EARTH_RADIUS_M: f64 = 6_371_008.8;
18
19/// The GeoSPARQL extension functions this engine implements.
20#[derive(Debug, Clone, Copy, PartialEq, Eq)]
21pub enum GeoFn {
22    Distance,
23    Contains,
24    Within,
25    Intersects,
26    Touches,
27}
28
29/// Canonical function IRIs, keyed by `GeoFn`. Parser and evaluator agree on
30/// these so a `geof:<local>` call hashes to the same `Function::Custom` id both
31/// sides compute.
32pub const GEO_FUNCTION_IRIS: &[(&str, GeoFn)] = &[
33    (
34        "http://www.opengis.net/def/function/geosparql/distance",
35        GeoFn::Distance,
36    ),
37    (
38        "http://www.opengis.net/def/function/geosparql/sfContains",
39        GeoFn::Contains,
40    ),
41    (
42        "http://www.opengis.net/def/function/geosparql/sfWithin",
43        GeoFn::Within,
44    ),
45    (
46        "http://www.opengis.net/def/function/geosparql/sfIntersects",
47        GeoFn::Intersects,
48    ),
49    (
50        "http://www.opengis.net/def/function/geosparql/sfTouches",
51        GeoFn::Touches,
52    ),
53];
54
55/// Canonical IRI for a `geof:<local>` function name, if recognised. Lets the
56/// parser expand a `geof:` call to the standard IRI even when the query did not
57/// declare the prefix.
58pub fn geo_function_iri(local: &str) -> Option<&'static str> {
59    let want = match local {
60        "distance" => GeoFn::Distance,
61        "sfContains" => GeoFn::Contains,
62        "sfWithin" => GeoFn::Within,
63        "sfIntersects" => GeoFn::Intersects,
64        "sfTouches" => GeoFn::Touches,
65        _ => return None,
66    };
67    GEO_FUNCTION_IRIS
68        .iter()
69        .find(|(_, k)| *k == want)
70        .map(|(iri, _)| *iri)
71}
72
73/// Map a function-IRI hash back to a `GeoFn` (used by the evaluator to dispatch
74/// a `Function::Custom(hash)`).
75pub fn geo_fn_for_hash(hash: u64) -> Option<GeoFn> {
76    GEO_FUNCTION_IRIS
77        .iter()
78        .find(|(iri, _)| crate::lexicon::generate_60bit_token(iri.as_bytes()) == hash)
79        .map(|(_, k)| *k)
80}
81
82/// Result of a GeoSPARQL predicate: a boolean (sf* topology) or a metric.
83#[derive(Debug, Clone, Copy, PartialEq)]
84pub enum GeoValue {
85    Bool(bool),
86    /// Distance in metres.
87    Number(f64),
88}
89
90/// Evaluate a `GeoFn` over two already-parsed geometries.
91pub fn eval_geo_fn(func: GeoFn, a: &Geometry, b: &Geometry) -> GeoValue {
92    match func {
93        GeoFn::Distance => GeoValue::Number(distance_metres(a, b)),
94        GeoFn::Contains => GeoValue::Bool(contains(a, b)),
95        GeoFn::Within => GeoValue::Bool(within(a, b)),
96        GeoFn::Intersects => GeoValue::Bool(intersects(a, b)),
97        GeoFn::Touches => GeoValue::Bool(touches(a, b)),
98    }
99}
100
101/// A parsed WKT geometry. Coordinates are `(x, y)` = `(longitude, latitude)`
102/// for geographic data.
103#[derive(Debug, Clone, PartialEq)]
104pub enum Geometry {
105    Point(Coord),
106    LineString(Vec<Coord>),
107    /// Rings: `rings[0]` is the exterior ring, the rest are holes.
108    Polygon(Vec<Vec<Coord>>),
109    MultiPoint(Vec<Coord>),
110    MultiLineString(Vec<Vec<Coord>>),
111    MultiPolygon(Vec<Vec<Vec<Coord>>>),
112    GeometryCollection(Vec<Geometry>),
113}
114
115/// A 2-D coordinate `(x, y)`.
116#[derive(Debug, Clone, Copy, PartialEq)]
117pub struct Coord {
118    pub x: f64,
119    pub y: f64,
120}
121
122impl Coord {
123    fn new(x: f64, y: f64) -> Self {
124        Self { x, y }
125    }
126}
127
128/// Parse a WKT literal into a `Geometry`. Accepts an optional leading CRS URI
129/// (`<http://…/CRS84> POINT(…)`) which is skipped. Case-insensitive keywords.
130pub fn parse_wkt(input: &str) -> Result<Geometry, String> {
131    let s = input.trim();
132    // Skip an optional leading `<crs-uri>`.
133    let s = if let Some(rest) = s.strip_prefix('<') {
134        match rest.find('>') {
135            Some(i) => rest[i + 1..].trim_start(),
136            None => return Err("unterminated CRS URI in WKT".to_string()),
137        }
138    } else {
139        s
140    };
141    let mut p = WktParser {
142        bytes: s.as_bytes(),
143        pos: 0,
144        src: s,
145    };
146    let g = p.parse_geometry()?;
147    Ok(g)
148}
149
150struct WktParser<'a> {
151    bytes: &'a [u8],
152    pos: usize,
153    src: &'a str,
154}
155
156impl<'a> WktParser<'a> {
157    fn skip_ws(&mut self) {
158        while self.pos < self.bytes.len() && (self.bytes[self.pos] as char).is_ascii_whitespace() {
159            self.pos += 1;
160        }
161    }
162
163    /// Read a keyword (letters), upper-cased.
164    fn keyword(&mut self) -> String {
165        self.skip_ws();
166        let start = self.pos;
167        while self.pos < self.bytes.len() && (self.bytes[self.pos] as char).is_ascii_alphabetic() {
168            self.pos += 1;
169        }
170        self.src[start..self.pos].to_ascii_uppercase()
171    }
172
173    /// After a keyword, skip an optional `Z`/`M`/`ZM` dimensionality token.
174    fn skip_dim(&mut self) {
175        self.skip_ws();
176        let start = self.pos;
177        while self.pos < self.bytes.len() {
178            match self.bytes[self.pos] {
179                b'Z' | b'z' | b'M' | b'm' => self.pos += 1,
180                _ => break,
181            }
182        }
183        // Only consume if it was a standalone dim token (followed by ws or `(`).
184        if self.pos > start {
185            let ok = self.pos >= self.bytes.len()
186                || matches!(self.bytes[self.pos] as char, ' ' | '\t' | '(' | '\n' | '\r');
187            if !ok {
188                self.pos = start; // it was part of something else (shouldn't happen)
189            }
190        }
191    }
192
193    fn expect(&mut self, c: u8) -> Result<(), String> {
194        self.skip_ws();
195        if self.pos < self.bytes.len() && self.bytes[self.pos] == c {
196            self.pos += 1;
197            Ok(())
198        } else {
199            Err(format!(
200                "expected '{}' in WKT at byte {}",
201                c as char, self.pos
202            ))
203        }
204    }
205
206    fn peek(&mut self) -> Option<u8> {
207        self.skip_ws();
208        self.bytes.get(self.pos).copied()
209    }
210
211    fn parse_geometry(&mut self) -> Result<Geometry, String> {
212        let kw = self.keyword();
213        self.skip_dim();
214        // EMPTY geometries.
215        self.skip_ws();
216        if self.src[self.pos..]
217            .to_ascii_uppercase()
218            .starts_with("EMPTY")
219        {
220            self.pos += 5;
221            return Ok(match kw.as_str() {
222                "POINT" => Geometry::Point(Coord::new(f64::NAN, f64::NAN)),
223                "LINESTRING" => Geometry::LineString(Vec::new()),
224                "POLYGON" => Geometry::Polygon(Vec::new()),
225                _ => Geometry::GeometryCollection(Vec::new()),
226            });
227        }
228        match kw.as_str() {
229            "POINT" => {
230                self.expect(b'(')?;
231                let c = self.coord()?;
232                self.expect(b')')?;
233                Ok(Geometry::Point(c))
234            }
235            "LINESTRING" => Ok(Geometry::LineString(self.coord_list()?)),
236            "POLYGON" => Ok(Geometry::Polygon(self.ring_list()?)),
237            "MULTIPOINT" => {
238                // MULTIPOINT allows `(1 2, 3 4)` or `((1 2), (3 4))`.
239                self.expect(b'(')?;
240                let mut pts = Vec::new();
241                loop {
242                    if self.peek() == Some(b'(') {
243                        self.expect(b'(')?;
244                        pts.push(self.coord()?);
245                        self.expect(b')')?;
246                    } else {
247                        pts.push(self.coord()?);
248                    }
249                    self.skip_ws();
250                    if self.peek() == Some(b',') {
251                        self.expect(b',')?;
252                    } else {
253                        break;
254                    }
255                }
256                self.expect(b')')?;
257                Ok(Geometry::MultiPoint(pts))
258            }
259            "MULTILINESTRING" => {
260                self.expect(b'(')?;
261                let mut lines = Vec::new();
262                loop {
263                    lines.push(self.coord_list()?);
264                    if self.peek() == Some(b',') {
265                        self.expect(b',')?;
266                    } else {
267                        break;
268                    }
269                }
270                self.expect(b')')?;
271                Ok(Geometry::MultiLineString(lines))
272            }
273            "MULTIPOLYGON" => {
274                self.expect(b'(')?;
275                let mut polys = Vec::new();
276                loop {
277                    polys.push(self.ring_list()?);
278                    if self.peek() == Some(b',') {
279                        self.expect(b',')?;
280                    } else {
281                        break;
282                    }
283                }
284                self.expect(b')')?;
285                Ok(Geometry::MultiPolygon(polys))
286            }
287            "GEOMETRYCOLLECTION" => {
288                self.expect(b'(')?;
289                let mut geoms = Vec::new();
290                loop {
291                    geoms.push(self.parse_geometry()?);
292                    if self.peek() == Some(b',') {
293                        self.expect(b',')?;
294                    } else {
295                        break;
296                    }
297                }
298                self.expect(b')')?;
299                Ok(Geometry::GeometryCollection(geoms))
300            }
301            other => Err(format!("unsupported WKT geometry type '{other}'")),
302        }
303    }
304
305    /// `( x y, x y, … )`
306    fn coord_list(&mut self) -> Result<Vec<Coord>, String> {
307        self.expect(b'(')?;
308        let mut out = Vec::new();
309        loop {
310            out.push(self.coord()?);
311            self.skip_ws();
312            if self.peek() == Some(b',') {
313                self.expect(b',')?;
314            } else {
315                break;
316            }
317        }
318        self.expect(b')')?;
319        Ok(out)
320    }
321
322    /// `( (ring), (hole), … )`
323    fn ring_list(&mut self) -> Result<Vec<Vec<Coord>>, String> {
324        self.expect(b'(')?;
325        let mut rings = Vec::new();
326        loop {
327            rings.push(self.coord_list()?);
328            if self.peek() == Some(b',') {
329                self.expect(b',')?;
330            } else {
331                break;
332            }
333        }
334        self.expect(b')')?;
335        Ok(rings)
336    }
337
338    /// `x y [z] [m]` — a single coordinate (extra ordinates ignored).
339    fn coord(&mut self) -> Result<Coord, String> {
340        let x = self.number()?;
341        let y = self.number()?;
342        // Consume any extra Z/M ordinates.
343        while let Some(b) = self.peek() {
344            if (b as char).is_ascii_digit() || b == b'-' || b == b'+' || b == b'.' {
345                let _ = self.number()?;
346            } else {
347                break;
348            }
349        }
350        Ok(Coord::new(x, y))
351    }
352
353    fn number(&mut self) -> Result<f64, String> {
354        self.skip_ws();
355        let start = self.pos;
356        while self.pos < self.bytes.len() {
357            let c = self.bytes[self.pos] as char;
358            if c.is_ascii_digit() || matches!(c, '-' | '+' | '.' | 'e' | 'E') {
359                self.pos += 1;
360            } else {
361                break;
362            }
363        }
364        if self.pos == start {
365            return Err(format!("expected a number in WKT at byte {}", self.pos));
366        }
367        self.src[start..self.pos]
368            .parse::<f64>()
369            .map_err(|_| format!("invalid number '{}' in WKT", &self.src[start..self.pos]))
370    }
371}
372
373// ─── Predicates ──────────────────────────────────────────────────────────────
374
375/// Great-circle distance in metres between two geometries' representative
376/// points (haversine). For non-point geometries the centroid of the coordinate
377/// set is used — a documented approximation adequate for `geof:distance`.
378pub fn distance_metres(a: &Geometry, b: &Geometry) -> f64 {
379    let pa = representative_point(a);
380    let pb = representative_point(b);
381    haversine(pa, pb)
382}
383
384fn haversine(a: Coord, b: Coord) -> f64 {
385    let lat1 = a.y.to_radians();
386    let lat2 = b.y.to_radians();
387    let dlat = (b.y - a.y).to_radians();
388    let dlon = (b.x - a.x).to_radians();
389    let h = (dlat / 2.0).sin().powi(2) + lat1.cos() * lat2.cos() * (dlon / 2.0).sin().powi(2);
390    2.0 * EARTH_RADIUS_M * h.sqrt().asin()
391}
392
393fn representative_point(g: &Geometry) -> Coord {
394    match g {
395        Geometry::Point(c) => *c,
396        Geometry::LineString(cs) | Geometry::MultiPoint(cs) => centroid(cs),
397        Geometry::Polygon(rings) => rings
398            .first()
399            .map(|r| centroid(r))
400            .unwrap_or(Coord::new(0.0, 0.0)),
401        Geometry::MultiLineString(ls) => {
402            centroid(&ls.iter().flatten().copied().collect::<Vec<_>>())
403        }
404        Geometry::MultiPolygon(ps) => centroid(
405            &ps.iter()
406                .filter_map(|p| p.first())
407                .flatten()
408                .copied()
409                .collect::<Vec<_>>(),
410        ),
411        Geometry::GeometryCollection(gs) => gs
412            .first()
413            .map(representative_point)
414            .unwrap_or(Coord::new(0.0, 0.0)),
415    }
416}
417
418fn centroid(cs: &[Coord]) -> Coord {
419    if cs.is_empty() {
420        return Coord::new(0.0, 0.0);
421    }
422    let (sx, sy) = cs
423        .iter()
424        .fold((0.0, 0.0), |(ax, ay), c| (ax + c.x, ay + c.y));
425    Coord::new(sx / cs.len() as f64, sy / cs.len() as f64)
426}
427
428/// `geof:sfContains` — does `a` contain `b`? Implemented for the common case of
429/// a polygon containing a point / all points of another geometry.
430pub fn contains(a: &Geometry, b: &Geometry) -> bool {
431    match a {
432        Geometry::Polygon(rings) => all_points(b).iter().all(|p| point_in_polygon(*p, rings)),
433        Geometry::MultiPolygon(polys) => all_points(b)
434            .iter()
435            .all(|p| polys.iter().any(|r| point_in_polygon(*p, r))),
436        _ => false,
437    }
438}
439
440/// `geof:sfWithin` — `a` within `b` ≡ `b` contains `a`.
441pub fn within(a: &Geometry, b: &Geometry) -> bool {
442    contains(b, a)
443}
444
445/// `geof:sfIntersects` — do `a` and `b` share any point? Covers point-in-polygon,
446/// shared vertices, and segment crossings between line/polygon boundaries.
447pub fn intersects(a: &Geometry, b: &Geometry) -> bool {
448    // Any point of one inside a polygon of the other.
449    if contains(a, b) || contains(b, a) {
450        return true;
451    }
452    // Any shared point.
453    let pa = all_points(a);
454    let pb = all_points(b);
455    for x in &pa {
456        for y in &pb {
457            if coords_eq(*x, *y) {
458                return true;
459            }
460        }
461    }
462    // Boundary segment crossings.
463    let sa = segments(a);
464    let sb = segments(b);
465    for (p1, p2) in &sa {
466        for (q1, q2) in &sb {
467            if segments_intersect(*p1, *p2, *q1, *q2) {
468                return true;
469            }
470        }
471    }
472    false
473}
474
475/// `geof:sfTouches` — geometries share a boundary point but no interior. A
476/// pragmatic test: they intersect, but neither contains an interior point of the
477/// other (approximated as: they intersect and no vertex of one is strictly
478/// inside a polygon of the other).
479pub fn touches(a: &Geometry, b: &Geometry) -> bool {
480    if !intersects(a, b) {
481        return false;
482    }
483    // Interiors must be disjoint. Vertex-only tests miss overlaps whose shared
484    // region has corners on both boundaries (e.g. two squares offset by half a
485    // side); edge midpoints catch those — an overlapping square's edge midpoint
486    // lands strictly inside the other, while a merely-touching one's does not.
487    let interior = probe_points(b).iter().any(|p| strictly_inside(*p, a))
488        || probe_points(a).iter().any(|p| strictly_inside(*p, b));
489    !interior
490}
491
492/// Vertices plus edge midpoints — used to detect interior overlap for `touches`.
493fn probe_points(g: &Geometry) -> Vec<Coord> {
494    let mut pts = all_points(g);
495    for (p, q) in segments(g) {
496        pts.push(Coord::new((p.x + q.x) / 2.0, (p.y + q.y) / 2.0));
497    }
498    pts
499}
500
501fn strictly_inside(p: Coord, g: &Geometry) -> bool {
502    match g {
503        Geometry::Polygon(rings) => point_in_polygon(p, rings) && !on_boundary(p, rings),
504        Geometry::MultiPolygon(polys) => polys
505            .iter()
506            .any(|r| point_in_polygon(p, r) && !on_boundary(p, r)),
507        _ => false,
508    }
509}
510
511fn on_boundary(p: Coord, rings: &[Vec<Coord>]) -> bool {
512    for ring in rings {
513        for w in ring.windows(2) {
514            if point_on_segment(p, w[0], w[1]) {
515                return true;
516            }
517        }
518    }
519    false
520}
521
522fn all_points(g: &Geometry) -> Vec<Coord> {
523    match g {
524        Geometry::Point(c) => vec![*c],
525        Geometry::LineString(cs) | Geometry::MultiPoint(cs) => cs.clone(),
526        Geometry::Polygon(rings) => rings.iter().flatten().copied().collect(),
527        Geometry::MultiLineString(ls) => ls.iter().flatten().copied().collect(),
528        Geometry::MultiPolygon(ps) => ps.iter().flatten().flatten().copied().collect(),
529        Geometry::GeometryCollection(gs) => gs.iter().flat_map(all_points).collect(),
530    }
531}
532
533fn segments(g: &Geometry) -> Vec<(Coord, Coord)> {
534    let mut out = Vec::new();
535    let ring_segs = |cs: &[Coord], out: &mut Vec<(Coord, Coord)>| {
536        for w in cs.windows(2) {
537            out.push((w[0], w[1]));
538        }
539    };
540    match g {
541        Geometry::LineString(cs) => ring_segs(cs, &mut out),
542        Geometry::Polygon(rings) => {
543            for r in rings {
544                ring_segs(r, &mut out);
545            }
546        }
547        Geometry::MultiLineString(ls) => {
548            for l in ls {
549                ring_segs(l, &mut out);
550            }
551        }
552        Geometry::MultiPolygon(ps) => {
553            for p in ps {
554                for r in p {
555                    ring_segs(r, &mut out);
556                }
557            }
558        }
559        _ => {}
560    }
561    out
562}
563
564/// Ray-casting point-in-polygon (with holes): inside the exterior ring and not
565/// inside any hole. Points on the boundary count as inside.
566fn point_in_polygon(p: Coord, rings: &[Vec<Coord>]) -> bool {
567    let Some(exterior) = rings.first() else {
568        return false;
569    };
570    if on_boundary(p, rings) {
571        return true;
572    }
573    if !point_in_ring(p, exterior) {
574        return false;
575    }
576    // Inside a hole → not contained.
577    for hole in &rings[1..] {
578        if point_in_ring(p, hole) {
579            return false;
580        }
581    }
582    true
583}
584
585fn point_in_ring(p: Coord, ring: &[Coord]) -> bool {
586    let n = ring.len();
587    if n < 3 {
588        return false;
589    }
590    let mut inside = false;
591    let mut j = n - 1;
592    for i in 0..n {
593        let vi = ring[i];
594        let vj = ring[j];
595        if (vi.y > p.y) != (vj.y > p.y) {
596            let x_int = (vj.x - vi.x) * (p.y - vi.y) / (vj.y - vi.y) + vi.x;
597            if p.x < x_int {
598                inside = !inside;
599            }
600        }
601        j = i;
602    }
603    inside
604}
605
606fn coords_eq(a: Coord, b: Coord) -> bool {
607    (a.x - b.x).abs() < 1e-9 && (a.y - b.y).abs() < 1e-9
608}
609
610fn point_on_segment(p: Coord, a: Coord, b: Coord) -> bool {
611    let cross = (b.x - a.x) * (p.y - a.y) - (b.y - a.y) * (p.x - a.x);
612    if cross.abs() > 1e-9 {
613        return false;
614    }
615    let within_x = p.x >= a.x.min(b.x) - 1e-9 && p.x <= a.x.max(b.x) + 1e-9;
616    let within_y = p.y >= a.y.min(b.y) - 1e-9 && p.y <= a.y.max(b.y) + 1e-9;
617    within_x && within_y
618}
619
620fn orient(a: Coord, b: Coord, c: Coord) -> f64 {
621    (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x)
622}
623
624fn segments_intersect(p1: Coord, p2: Coord, q1: Coord, q2: Coord) -> bool {
625    let d1 = orient(q1, q2, p1);
626    let d2 = orient(q1, q2, p2);
627    let d3 = orient(p1, p2, q1);
628    let d4 = orient(p1, p2, q2);
629    if ((d1 > 0.0) != (d2 > 0.0)) && ((d3 > 0.0) != (d4 > 0.0)) {
630        return true;
631    }
632    // Collinear boundary touches.
633    point_on_segment(p1, q1, q2)
634        || point_on_segment(p2, q1, q2)
635        || point_on_segment(q1, p1, p2)
636        || point_on_segment(q2, p1, p2)
637}
638
639#[cfg(test)]
640mod tests {
641    use super::*;
642
643    #[test]
644    fn parse_point() {
645        assert_eq!(
646            parse_wkt("POINT(1.5 -2.25)").unwrap(),
647            Geometry::Point(Coord::new(1.5, -2.25))
648        );
649    }
650
651    #[test]
652    fn parse_point_z_ignored() {
653        assert_eq!(
654            parse_wkt("POINT Z (1 2 3)").unwrap(),
655            Geometry::Point(Coord::new(1.0, 2.0))
656        );
657    }
658
659    #[test]
660    fn parse_polygon_with_hole() {
661        let g = parse_wkt("POLYGON((0 0, 4 0, 4 4, 0 4, 0 0),(1 1, 2 1, 2 2, 1 2, 1 1))").unwrap();
662        match g {
663            Geometry::Polygon(rings) => {
664                assert_eq!(rings.len(), 2);
665                assert_eq!(rings[0].len(), 5);
666            }
667            other => panic!("expected polygon, got {other:?}"),
668        }
669    }
670
671    #[test]
672    fn parse_skips_crs_prefix() {
673        let g = parse_wkt("<http://www.opengis.net/def/crs/OGC/1.3/CRS84> POINT(10 20)").unwrap();
674        assert_eq!(g, Geometry::Point(Coord::new(10.0, 20.0)));
675    }
676
677    #[test]
678    fn haversine_known_distance() {
679        // London (-0.1276, 51.5074) → Paris (2.3522, 48.8566) ≈ 343 km.
680        let d = distance_metres(
681            &Geometry::Point(Coord::new(-0.1276, 51.5074)),
682            &Geometry::Point(Coord::new(2.3522, 48.8566)),
683        );
684        assert!((d - 343_556.0).abs() < 2_000.0, "distance was {d}");
685    }
686
687    #[test]
688    fn contains_point_in_square() {
689        let square = parse_wkt("POLYGON((0 0, 4 0, 4 4, 0 4, 0 0))").unwrap();
690        assert!(contains(&square, &Geometry::Point(Coord::new(2.0, 2.0))));
691        assert!(!contains(&square, &Geometry::Point(Coord::new(5.0, 5.0))));
692    }
693
694    #[test]
695    fn hole_excludes_point() {
696        let g = parse_wkt("POLYGON((0 0, 4 0, 4 4, 0 4, 0 0),(1 1, 3 1, 3 3, 1 3, 1 1))").unwrap();
697        assert!(
698            !contains(&g, &Geometry::Point(Coord::new(2.0, 2.0))),
699            "in hole"
700        );
701        assert!(
702            contains(&g, &Geometry::Point(Coord::new(0.5, 0.5))),
703            "outside hole"
704        );
705    }
706
707    #[test]
708    fn within_is_contains_flipped() {
709        let square = parse_wkt("POLYGON((0 0, 4 0, 4 4, 0 4, 0 0))").unwrap();
710        let pt = Geometry::Point(Coord::new(2.0, 2.0));
711        assert!(within(&pt, &square));
712        assert!(!within(&square, &pt));
713    }
714
715    #[test]
716    fn intersecting_lines() {
717        let a = parse_wkt("LINESTRING(0 0, 4 4)").unwrap();
718        let b = parse_wkt("LINESTRING(0 4, 4 0)").unwrap();
719        assert!(intersects(&a, &b));
720        let c = parse_wkt("LINESTRING(0 1, 4 5)").unwrap();
721        assert!(!intersects(&a, &c), "parallel, should not intersect");
722    }
723
724    #[test]
725    fn touching_squares() {
726        // Two unit squares sharing the edge x=1.
727        let a = parse_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))").unwrap();
728        let b = parse_wkt("POLYGON((1 0, 2 0, 2 1, 1 1, 1 0))").unwrap();
729        assert!(touches(&a, &b), "edge-sharing squares should touch");
730        // Overlapping squares intersect but do NOT merely touch.
731        let c = parse_wkt("POLYGON((0.5 0, 1.5 0, 1.5 1, 0.5 1, 0.5 0))").unwrap();
732        assert!(intersects(&a, &c));
733        assert!(!touches(&a, &c), "overlap is not a touch");
734    }
735}
qualia_core_db/sparql_library/geosparql.rs

qualia_core_db/sparql_library/
geosparql.rs
