bambam_osm/algorithm/clustering/
clustering_rtree.rs

1use crate::model::osm::graph::OsmNodeId;
2
3use super::ClusteredGeometry;
4use geo::{BoundingRect, Coord, Polygon};
5use itertools::Itertools;
6use kdam::tqdm;
7use rstar::primitives::{GeomWithData, Rectangle};
8use rstar::{RTree, RTreeObject};
9use wkt::ToWkt;
10
11pub type ClusteredIntersections = GeomWithData<Rectangle<(f32, f32)>, ClusteredGeometry>;
12
13/// build an undirected graph of node geometries that intersect spatially.
14/// clusters are represented simply, without any changes to their geometries and linear-time
15/// intersection search for clusters. this performance hit is taken to avoid any edge cases
16/// where manipulation via overlay operations might lead to representation errors.
17pub fn build(
18    geometries: &[(OsmNodeId, Polygon<f32>)],
19) -> Result<RTree<ClusteredIntersections>, String> {
20    // intersection performed via RTree.
21
22    let mut rtree: RTree<ClusteredIntersections> = RTree::new();
23    let iter = tqdm!(
24        geometries.iter(),
25        total = geometries.len(),
26        desc = "spatial intersection"
27    );
28
29    // add each geometry to the rtree.
30    for (index, polygon) in iter {
31        let rect = rect_from_geometries(&[polygon])?;
32        let query = GeomWithData::new(rect, ClusteredGeometry::new(*index, polygon.clone()));
33        let intersecting = rtree
34            .drain_in_envelope_intersecting(query.envelope())
35            .sorted_by_key(|obj| obj.data.ids())
36            .collect_vec();
37        if intersecting.is_empty() {
38            // nothing intersects with this new cluster, insert it and move on to the next row.
39            rtree.insert(query);
40        } else {
41            // prepare to merge this geometry with any intersecting geometries by union.
42            let mut new_cluster = ClusteredGeometry::new(*index, polygon.clone());
43
44            for obj in intersecting.into_iter() {
45                // it is still possible that none of the "intersecting" geometries actually
46                // truly intersect since we only compared the bounding boxes at this point.
47                if obj.data.intersects(polygon) {
48                    // merge the intersecting data. since it was drained from the rtree, we are done.
49                    new_cluster.merge_and_sort_with(&obj.data);
50                } else {
51                    // false alarm, this one doesn't actually intersect, put it back
52                    // in the tree without changes since it was drained.
53                    rtree.insert(obj);
54                }
55            }
56
57            let new_rect = rect_from_geometries(&new_cluster.polygons())?;
58            let new_obj = GeomWithData::new(new_rect, new_cluster);
59            rtree.insert(new_obj);
60        }
61    }
62    eprintln!();
63
64    Ok(rtree)
65}
66
67/// helper function to create a rectangular rtree envelope for a given geometry
68fn rect_from_geometries(ps: &[&Polygon<f32>]) -> Result<Rectangle<(f32, f32)>, String> {
69    if ps.is_empty() {
70        return Err(String::from(
71            "rect_from_geometries called with empty collection",
72        ));
73    }
74    let mut mins = vec![];
75    let mut maxs = vec![];
76    for p in ps {
77        let bbox_rect = p.bounding_rect().ok_or_else(|| {
78            format!(
79                "internal error: cannot get bounds of geometry: '{}'",
80                p.to_wkt()
81            )
82        })?;
83        mins.push(bbox_rect.min());
84        maxs.push(bbox_rect.max());
85    }
86    let min_coord = mins
87        .into_iter()
88        .min_by_key(ordering_key)
89        .ok_or_else(|| String::from("internal error: empty 'mins' collection"))?;
90    let max_coord = maxs
91        .into_iter()
92        .max_by_key(ordering_key)
93        .ok_or_else(|| String::from("internal error: empty 'maxs' collection"))?;
94    let envelope = Rectangle::from_corners(min_coord.x_y(), max_coord.x_y());
95    Ok(envelope)
96}
97
98/// called on WGS84 coordinates to create an ordering. since floating point
99/// values have no ordering in Rust, we use scaling and conversion to i64
100/// which is a feasible bijection since the max float values are +- 180.
101fn ordering_key(coord: &Coord<f32>) -> (i64, i64) {
102    let x = (coord.x * 100_000.0) as i64;
103    let y = (coord.y * 100_000.0) as i64;
104    (x, y)
105}
bambam_osm/algorithm/clustering/clustering_rtree.rs

bambam_osm/algorithm/clustering/
clustering_rtree.rs
