Convert everything to mercator internally, after all

the new algorithm still acting weird though

backend/Cargo.toml

@@ -15,8 +15,6 @@ geo = "0.27.0"
 geojson = { git = "https://github.com/georust/geojson", features = ["geo-types"] }
 log = "0.4.20"
 osm-reader = { git = "https://github.com/a-b-street/osm-reader" }
-rand = { version = "0.8.5", default-features = false }
-rand_xorshift = "0.3.0"
 rstar = { version = "0.11.0" }
 serde = "1.0.188"
 serde_json = "1.0.105"

backend/src/heatmap.rs

@@ -1,40 +1,11 @@
 use std::collections::HashSet;
 
-use geo::{
-    BoundingRect, DensifyHaversine, Geometry, GeometryCollection, HaversineBearing,
-    HaversineDestination, Line, LineString, Point, Rect,
-};
+use geo::{Coord, Densify, Line, LineString};
 use geojson::{Feature, FeatureCollection};
-use rand::{Rng, SeedableRng};
-use rand_xorshift::XorShiftRng;
 use rstar::{primitives::GeomWithData, RTree};
 
 use crate::{CompareRouteRequest, IntersectionID, MapModel, RoadKind};
 
-pub fn measure_randomly(map: &mut MapModel, n: usize) -> FeatureCollection {
-    // TODO Expensive
-    let bbox: Rect<f64> = map
-        .roads
-        .iter()
-        .map(|r| Geometry::LineString(r.linestring.clone()))
-        .collect::<GeometryCollection>()
-        .bounding_rect()
-        .unwrap();
-    // TODO Do this in the right coordinate space
-    let dist_away = 0.01;
-
-    let mut rng = XorShiftRng::seed_from_u64(42);
-    let mut requests = Vec::new();
-    for _ in 0..n {
-        let x1 = rng.gen_range(bbox.min().x..=bbox.max().x);
-        let y1 = rng.gen_range(bbox.min().y..=bbox.max().y);
-        let x2 = x1 + rng.gen_range(-dist_away..=dist_away);
-        let y2 = y1 + rng.gen_range(-dist_away..=dist_away);
-        requests.push(CompareRouteRequest { x1, y1, x2, y2 });
-    }
-    calculate(map, requests)
-}
-
 // Walk along severances. Every X meters, try to cross from one side to the other.
 //
 // We could focus where footways connect to severances, but that's probably a crossing. Ideally we
@@ -77,19 +48,14 @@ pub fn nearby_footway_intersections(map: &mut MapModel, dist_meters: f64) -> Fea
     for i1 in &footway_intersections {
         let i1_pt = map.intersections[i1.0].point;
         for i2 in rtree.locate_within_distance(i1_pt.into(), dist_meters) {
+            // TODO Skip trivial things connected by a road
             let i2_pt = map.intersections[i2.data.0].point;
             requests.push(CompareRouteRequest {
                 x1: i1_pt.x(),
                 y1: i1_pt.y(),
                 x2: i2_pt.x(),
                 y2: i2_pt.y(),
             });
-            if requests.len() > 100 {
-                break;
-            }
-        }
-        if requests.len() > 100 {
-            break;
         }
     }
     calculate(map, requests)
@@ -100,8 +66,14 @@ fn calculate(map: &mut MapModel, requests: Vec<CompareRouteRequest>) -> FeatureC
     let mut max_score = 0.0_f64;
     for req in requests {
         let mut f = Feature::from(geojson::Geometry::from(&LineString::new(vec![
-            (req.x1, req.y1).into(),
-            (req.x2, req.y2).into(),
+            map.mercator.to_wgs84(Coord {
+                x: req.x1,
+                y: req.y1,
+            }),
+            map.mercator.to_wgs84(Coord {
+                x: req.x2,
+                y: req.y2,
+            }),
         ])));
         if let Ok(fc) = crate::route::do_route(map, req) {
             let direct = fc
@@ -142,13 +114,22 @@ fn make_perpendicular_offsets(
     let mut output = Vec::new();
     // Using lines instead of coords so we can get the angle -- but is this hard to reason about?
     // angle_at_point instead?
-    for orig_line in linestring.densify_haversine(walk_every_m).lines() {
+    for orig_line in linestring.densify(walk_every_m).lines() {
         // TODO For the last line, use the last point too
-        let pt: Point = orig_line.start.into();
-        let angle = pt.haversine_bearing(orig_line.end.into());
-        let projected_left = pt.haversine_destination(angle - 90.0, project_away_m);
-        let projected_right = pt.haversine_destination(angle + 90.0, project_away_m);
+        let angle_degs = (orig_line.end.y - orig_line.start.y)
+            .atan2(orig_line.end.x - orig_line.start.x)
+            .to_degrees();
+        let projected_left = project_away(orig_line.start, angle_degs - 90.0, project_away_m);
+        let projected_right = project_away(orig_line.start, angle_degs + 90.0, project_away_m);
         output.push(Line::new(projected_left, projected_right));
     }
     output
 }
+
+fn project_away(pt: Coord, angle_degs: f64, dist_away_m: f64) -> Coord {
+    let (sin, cos) = angle_degs.to_radians().sin_cos();
+    Coord {
+        x: pt.x + dist_away_m * cos,
+        y: pt.y + dist_away_m * sin,
+    }
+}

backend/src/lib.rs

@@ -5,13 +5,14 @@ use std::fmt;
 use std::sync::Once;
 
 use fast_paths::{FastGraph, PathCalculator};
-use geo::{Line, LineString, Point};
+use geo::{Coord, Line, LineString, MapCoordsInPlace, Point};
 use geojson::{Feature, GeoJson, Geometry};
 use rstar::{primitives::GeomWithData, RTree};
 use serde::{Deserialize, Serialize};
 use wasm_bindgen::prelude::*;
 
 mod heatmap;
+mod mercator;
 mod node_map;
 mod route;
 mod scrape;
@@ -23,6 +24,8 @@ static START: Once = Once::new();
 pub struct MapModel {
     roads: Vec<Road>,
     intersections: Vec<Intersection>,
+    // All geometry stored in worldspace, including rtrees
+    mercator: mercator::Mercator,
     // Only snaps to walkable roads
     closest_intersection: RTree<IntersectionLocation>,
     node_map: node_map::NodeMap<IntersectionID>,
@@ -103,7 +106,7 @@ impl MapModel {
         let mut features = Vec::new();
 
         for r in &self.roads {
-            features.push(r.to_gj());
+            features.push(r.to_gj(&self.mercator));
         }
 
         let gj = GeoJson::from(features);
@@ -114,15 +117,33 @@ impl MapModel {
     #[wasm_bindgen(js_name = compareRoute)]
     pub fn compare_route(&mut self, input: JsValue) -> Result<String, JsValue> {
         let req: CompareRouteRequest = serde_wasm_bindgen::from_value(input)?;
-        let gj = route::do_route(self, req).map_err(err_to_js)?;
+        let pt1 = self.mercator.to_mercator(Coord {
+            x: req.x1,
+            y: req.y1,
+        });
+        let pt2 = self.mercator.to_mercator(Coord {
+            x: req.x2,
+            y: req.y2,
+        });
+        let gj = route::do_route(
+            self,
+            CompareRouteRequest {
+                x1: pt1.x,
+                y1: pt1.y,
+                x2: pt2.x,
+                y2: pt2.y,
+            },
+        )
+        .map_err(err_to_js)?;
         let out = serde_json::to_string(&gj).map_err(err_to_js)?;
         Ok(out)
     }
 
     #[wasm_bindgen(js_name = makeHeatmap)]
     pub fn make_heatmap(&mut self) -> Result<String, JsValue> {
-        //let samples = heatmap::along_severances(self);
-        let samples = heatmap::nearby_footway_intersections(self, 0.01);
+        let samples = heatmap::along_severances(self);
+        // TODO unit here is weird or wrong or something
+        //let samples = heatmap::nearby_footway_intersections(self, 500.0);
         let out = serde_json::to_string(&samples).map_err(err_to_js)?;
         Ok(out)
     }
@@ -140,8 +161,11 @@ impl MapModel {
 }
 
 impl Road {
-    fn to_gj(&self) -> Feature {
-        let mut f = Feature::from(Geometry::from(&self.linestring));
+    fn to_gj(&self, mercator: &mercator::Mercator) -> Feature {
+        let mut linestring = self.linestring.clone();
+        linestring.map_coords_in_place(|c| mercator.to_wgs84(c));
+
+        let mut f = Feature::from(Geometry::from(&linestring));
         f.set_property("id", self.id.0);
         f.set_property("kind", format!("{:?}", self.kind));
         f.set_property("way", self.way.to_string());
@@ -154,6 +178,8 @@ impl Road {
     }
 }
 
+// Mercator worldspace internally, but not when it comes in from the app
+// TODO only use this on the boundary
 #[derive(Deserialize)]
 pub struct CompareRouteRequest {
     x1: f64,

backend/src/mercator.rs

@@ -0,0 +1,48 @@
+use geo::{BoundingRect, Coord, HaversineLength, LineString, Rect};
+
+/// Projects WGS84 points onto a Euclidean plane, using a Mercator projection. The top-left is (0,
+/// 0) and grows to the right and down (screen-drawing order, not Cartesian), with units of meters.
+/// The accuracy of this weakens for larger areas.
+pub struct Mercator {
+    wgs84_bounds: Rect,
+    width: f64,
+    height: f64,
+}
+
+impl Mercator {
+    // TODO The API is kind of annoying, or wasteful. Do builder style.
+    /// Create a boundary covering some geometry
+    pub fn from<T: BoundingRect<f64>>(geometry: T) -> Option<Self> {
+        let wgs84_bounds = geometry.bounding_rect().into()?;
+        let width = LineString::from(vec![
+            (wgs84_bounds.min().x, wgs84_bounds.min().y),
+            (wgs84_bounds.max().x, wgs84_bounds.min().y),
+        ])
+        .haversine_length();
+        let height = LineString::from(vec![
+            (wgs84_bounds.min().x, wgs84_bounds.min().y),
+            (wgs84_bounds.min().x, wgs84_bounds.max().y),
+        ])
+        .haversine_length();
+        Some(Self {
+            wgs84_bounds,
+            width,
+            height,
+        })
+    }
+
+    pub fn to_mercator(&self, pt: Coord) -> Coord {
+        let x = self.width * (pt.x - self.wgs84_bounds.min().x) / self.wgs84_bounds.width();
+        // Invert y, so that the northernmost latitude is 0
+        let y = self.height
+            - self.height * (pt.y - self.wgs84_bounds.min().y) / self.wgs84_bounds.height();
+        Coord { x, y }
+    }
+
+    pub fn to_wgs84(&self, pt: Coord) -> Coord {
+        let x = self.wgs84_bounds.min().x + (pt.x / self.width * self.wgs84_bounds.width());
+        let y = self.wgs84_bounds.min().y
+            + (self.wgs84_bounds.height() * (self.height - pt.y) / self.height);
+        Coord { x, y }
+    }
+}

backend/src/route.rs

@@ -38,7 +38,6 @@ pub fn build_router(
     (closest_intersection, node_map, ch)
 }
 
-// TODO We may be able to override the distance function? Does it work with WGS84?
 fn build_closest_intersection(
     intersections: &Vec<Intersection>,
     node_map: &NodeMap<IntersectionID>,
@@ -74,7 +73,7 @@ pub fn do_route(map: &mut MapModel, req: CompareRouteRequest) -> Result<FeatureC
             let i1 = map.node_map.translate_id(pair[0]);
             let i2 = map.node_map.translate_id(pair[1]);
             let road = map.find_edge(i1, i2);
-            features.push(road.to_gj());
+            features.push(road.to_gj(&map.mercator));
             route_length += road.linestring.haversine_length();
         }
         let direct_length = LineString::new(vec![

backend/src/scrape.rs

@@ -1,9 +1,10 @@
 use std::collections::HashMap;
 
 use anyhow::Result;
-use geo::{Coord, LineString, Point};
+use geo::{Coord, Geometry, GeometryCollection, LineString, MapCoordsInPlace, Point};
 use osm_reader::{Element, NodeID, WayID};
 
+use crate::mercator::Mercator;
 use crate::tags::Tags;
 use crate::{Intersection, IntersectionID, MapModel, Road, RoadID, RoadKind};
 
@@ -35,13 +36,35 @@ pub fn scrape_osm(input_bytes: &[u8]) -> Result<MapModel> {
         }
     }
 
-    let (roads, intersections) = split_edges(&node_mapping, highways);
+    let (mut roads, mut intersections) = split_edges(&node_mapping, highways);
+
+    // TODO expensive
+    let collection: GeometryCollection = roads
+        .iter()
+        .map(|r| Geometry::LineString(r.linestring.clone()))
+        .chain(
+            intersections
+                .iter()
+                .map(|i| Geometry::Point(i.point.clone())),
+        )
+        .collect::<Vec<_>>()
+        .into();
+    let mercator = Mercator::from(collection).unwrap();
+    for r in &mut roads {
+        r.linestring
+            .map_coords_in_place(|pt| mercator.to_mercator(pt));
+    }
+    for i in &mut intersections {
+        i.point.map_coords_in_place(|pt| mercator.to_mercator(pt));
+    }
+
     let (closest_intersection, node_map, ch) = crate::route::build_router(&intersections, &roads);
     let path_calc = fast_paths::create_calculator(&ch);
 
     Ok(MapModel {
         roads,
         intersections,
+        mercator,
         closest_intersection,
         node_map,
         ch,