hacked together demo

src/main.rs

@@ -0,0 +1,97 @@
+use geo::{BoundingRect, Contains, Distance, Euclidean};
+use geo::{LineString, Point, Polygon};
+use rstar::{PointDistance, RTree, RTreeObject, AABB};
+use shapefile::Reader;
+
+#[derive(Debug)]
+struct PolygonWrapper(Polygon<f64>);
+
+impl RTreeObject for PolygonWrapper {
+    type Envelope = AABB<[f64; 2]>;
+
+    fn envelope(&self) -> Self::Envelope {
+        let bbox = self.0.bounding_rect().unwrap();
+        AABB::from_corners([bbox.min().x, bbox.min().y], [bbox.max().x, bbox.max().y])
+    }
+}
+
+impl PointDistance for PolygonWrapper {
+    fn contains_point(&self, point: &[f64; 2]) -> bool {
+        self.0.contains(&Point::new(point[0], point[1]))
+    }
+
+    fn distance_2(&self, point: &[f64; 2]) -> f64 {
+        3.0
+        // Euclidean
+        //     .distance(
+        //         self.0,
+        //          &Point::new(point[0], point[1])
+        //      )
+        //     .powi(2)
+    }
+}
+
+fn load_polygons(path: &str) -> Result<Vec<Polygon<f64>>, Box<dyn std::error::Error>> {
+    let mut reader = Reader::from_path(path)?;
+    let mut geometries = Vec::new();
+
+    for shape_record in reader.iter_shapes_and_records() {
+        let (shape, _) = shape_record?;
+
+        if let shapefile::Shape::Polygon(polygon) = shape {
+            let mut polygons = Vec::new();
+            let mut current_exterior: Option<LineString<f64>> = None;
+            let mut current_holes = Vec::new();
+
+            for ring in polygon.rings() {
+                let coords: Vec<_> = ring.points().iter().map(|p| (p.x, p.y)).collect();
+                let line_string = LineString::from(coords);
+
+                if is_clockwise(&line_string) {
+                    if let Some(exterior) = current_exterior.take() {
+                        polygons.push(Polygon::new(exterior, current_holes.drain(..).collect()));
+                    }
+                    current_exterior = Some(line_string);
+                } else {
+                    current_holes.push(line_string);
+                }
+            }
+
+            if let Some(exterior) = current_exterior {
+                polygons.push(Polygon::new(exterior, current_holes));
+            }
+
+            geometries.extend(polygons);
+        }
+    }
+
+    Ok(geometries)
+}
+
+fn is_clockwise(ring: &LineString<f64>) -> bool {
+    let coords: Vec<_> = ring.coords().collect();
+    let mut sum = 0.0;
+    for window in coords.windows(2) {
+        let p1 = window[0];
+        let p2 = window[1];
+        sum += (p2.x - p1.x) * (p2.y + p1.y);
+    }
+    sum > 0.0
+}
+
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    let polygons = load_polygons("il_cd119/tl_2024_17_cd119.shp")?;
+    println!("Loaded {} polygons", polygons.len());
+
+    let wrapped: Vec<_> = polygons.into_iter().map(PolygonWrapper).collect();
+    let rtree = RTree::bulk_load(wrapped);
+
+    let query_point = Point::new(-87.630, 41.878);
+    let results = rtree.locate_all_at_point(&[query_point.x(), query_point.y()]);
+
+    for poly in results {
+        println!("Polygon contains query point: {:?}", poly.0.exterior());
+    }
+
+    Ok(())
+}