clipping roads based on original boundary polygons. buggy, so disabled.

convert_osm/src/clip.rs

@@ -0,0 +1,115 @@
+use abstutil::{retain_btreemap, Timer};
+use geom::{GPSBounds, LonLat, PolyLine, Polygon};
+use map_model::{raw_data, IntersectionType};
+use std::fs::File;
+use std::io::{BufRead, BufReader};
+
+pub fn clip_map(map: &mut raw_data::Map, path: &str, timer: &mut Timer) -> GPSBounds {
+    timer.start("clipping map to boundary");
+    let raw_boundary_pts = read_osmosis_polygon(path);
+    let mut bounds = map.get_gps_bounds();
+    for pt in &raw_boundary_pts {
+        bounds.update(*pt);
+    }
+
+    let boundary_poly = Polygon::new(&bounds.must_convert(&raw_boundary_pts));
+    let boundary_lines: Vec<PolyLine> = boundary_poly
+        .points()
+        .windows(2)
+        .map(|pair| PolyLine::new(pair.to_vec()))
+        .collect();
+
+    // This is kind of indirect and slow, but first pass -- just remove roads completely outside
+    // the boundary polygon.
+    retain_btreemap(&mut map.roads, |_, r| {
+        let center_pts = bounds.must_convert(&r.points);
+        let first_in = boundary_poly.contains_pt(center_pts[0]);
+        let last_in = boundary_poly.contains_pt(*center_pts.last().unwrap());
+        first_in || last_in
+    });
+
+    let road_ids: Vec<raw_data::StableRoadID> = map.roads.keys().cloned().collect();
+    for id in road_ids {
+        let r = &map.roads[&id];
+        let center_pts = bounds.must_convert(&r.points);
+        let first_in = boundary_poly.contains_pt(center_pts[0]);
+        let last_in = boundary_poly.contains_pt(*center_pts.last().unwrap());
+
+        if first_in && last_in {
+            continue;
+        }
+
+        let move_i = if first_in { r.i2 } else { r.i1 };
+
+        // The road crosses the boundary. But if the intersection happens to have another connected
+        // road, then just allow this exception.
+        // TODO But what about a road slightly outside the bounds that'd otherwise connect two
+        // things in bounds? Really ought to flood outwards and see if we wind up back inside.
+        if map
+            .roads
+            .values()
+            .filter(|r2| r2.i1 == move_i || r2.i1 == move_i)
+            .count()
+            > 1
+        {
+            println!(
+                "{} crosses boundary, but briefly enough to not touch it",
+                id
+            );
+            continue;
+        }
+
+        let i = map.intersections.get_mut(&move_i).unwrap();
+        i.intersection_type = IntersectionType::Border;
+
+        // Convert the road points to a PolyLine here. Loop roads were breaking!
+        let center = PolyLine::new(center_pts);
+
+        // Now trim it.
+        let mut_r = map.roads.get_mut(&id).unwrap();
+        let border_pt = boundary_lines
+            .iter()
+            .find_map(|l| center.intersection(l).map(|(pt, _)| pt))
+            .unwrap();
+        if first_in {
+            mut_r.points =
+                bounds.must_convert_back(center.get_slice_ending_at(border_pt).unwrap().points());
+            i.point = *mut_r.points.last().unwrap();
+        } else {
+            mut_r.points = bounds.must_convert_back(
+                center
+                    .reversed()
+                    .get_slice_ending_at(border_pt)
+                    .unwrap()
+                    .reversed()
+                    .points(),
+            );
+            i.point = mut_r.points[0];
+        }
+    }
+
+    timer.stop("clipping map to boundary");
+    bounds
+}
+
+fn read_osmosis_polygon(path: &str) -> Vec<LonLat> {
+    let mut pts: Vec<LonLat> = Vec::new();
+    for (idx, maybe_line) in BufReader::new(File::open(path).unwrap())
+        .lines()
+        .enumerate()
+    {
+        if idx == 0 || idx == 1 {
+            continue;
+        }
+        let line = maybe_line.unwrap();
+        if line == "END" {
+            break;
+        }
+        let parts: Vec<&str> = line.trim_start().split("    ").collect();
+        assert!(parts.len() == 2);
+        let lon = parts[0].parse::<f64>().unwrap();
+        let lat = parts[1].parse::<f64>().unwrap();
+        pts.push(LonLat::new(lon, lat));
+    }
+    pts
+}

convert_osm/src/lib.rs

@@ -1,3 +1,4 @@
+mod clip;
 mod group_parcels;
 mod neighborhoods;
 mod osm;
@@ -51,6 +52,10 @@ pub struct Flags {
     #[structopt(long = "neighborhoods")]
     pub neighborhoods: String,
 
+    /// Osmosis clippiny polgon
+    #[structopt(long = "clip")]
+    pub clip: String,
+
     /// Output .abst path
     #[structopt(long = "output")]
     pub output: String,
@@ -64,10 +69,10 @@ pub fn convert(flags: &Flags, timer: &mut abstutil::Timer) -> raw_data::Map {
     let elevation = Elevation::new(&flags.elevation).expect("loading .hgt failed");
     let mut map =
         split_ways::split_up_roads(osm::osm_to_raw_roads(&flags.osm, timer), &elevation, timer);
+    let gps_bounds = clip::clip_map(&mut map, &flags.clip, timer);
     remove_disconnected::remove_disconnected_roads(&mut map, timer);
-    timer.start("calculate GPS bounds");
-    let gps_bounds = map.get_gps_bounds();
-    timer.stop("calculate GPS bounds");
+    // TODO Shouldn't we recalculate the gps_bounds after removing stuff? Or just base it off the
+    // boundary polygon?
 
     if flags.fast_dev {
         return map;
@@ -108,13 +113,7 @@ fn use_parking_hints(
     let mut closest: FindClosest<(raw_data::StableRoadID, bool)> =
         FindClosest::new(&gps_bounds.to_bounds());
     for (id, r) in &map.roads {
-        let pts = PolyLine::new(
-            r.points
-                .iter()
-                .map(|pt| Pt2D::from_gps(*pt, gps_bounds).unwrap())
-                .collect(),
-        );
-
+        let pts = PolyLine::new(gps_bounds.must_convert(&r.points));
         closest.add((*id, true), &pts.shift_right(LANE_THICKNESS));
         closest.add((*id, false), &pts.shift_left(LANE_THICKNESS));
     }

editor/src/plugins/edit/draw_neighborhoods.rs

@@ -119,10 +119,7 @@ impl BlockingPlugin for DrawNeighborhoodState {
             DrawNeighborhoodState::MovingPoint(n, current_idx) => (&n.points, Some(*current_idx)),
         };
         let gps_bounds = ctx.map.get_gps_bounds();
-        let pts: Vec<Pt2D> = raw_pts
-            .into_iter()
-            .map(|pt| Pt2D::from_gps(*pt, gps_bounds).unwrap())
-            .collect();
+        let pts: Vec<Pt2D> = gps_bounds.must_convert(&raw_pts);
 
         if pts.len() == 2 {
             g.draw_line(

geom/src/find_closest.rs

@@ -30,10 +30,7 @@ where
     }
 
     pub fn add_gps(&mut self, key: K, raw_pts: &Vec<LonLat>, gps_bounds: &GPSBounds) {
-        let pts: Vec<Pt2D> = raw_pts
-            .iter()
-            .map(|pt| Pt2D::from_gps(*pt, gps_bounds).unwrap())
-            .collect();
+        let pts: Vec<Pt2D> = gps_bounds.must_convert(raw_pts);
         self.geometries
             .insert(key.clone(), pts_to_line_string(&pts));
 

geom/src/gps.rs

@@ -125,4 +125,14 @@ impl GPSBounds {
         b.update(self.get_max_world_pt());
         b
     }
+
+    pub fn must_convert(&self, pts: &Vec<LonLat>) -> Vec<Pt2D> {
+        pts.iter()
+            .map(|pt| Pt2D::from_gps(*pt, self).unwrap())
+            .collect()
+    }
+
+    pub fn must_convert_back(&self, pts: &Vec<Pt2D>) -> Vec<LonLat> {
+        pts.iter().map(|pt| pt.to_gps(self).unwrap()).collect()
+    }
 }

import.sh

@@ -102,7 +102,8 @@ for poly in `ls ../data/polygons/`; do
 		--parking_shapes=../data/shapes/blockface \
 		--gtfs=../data/input/google_transit_2018_18_08 \
 		--neighborhoods=../data/input/neighborhoods.geojson \
+		--clip=../data/polygons/$name.poly \
 		--output=../data/raw_maps/$name.abst
 done
 
-# To run manually: cargo run -- --osm=../data/input/montlake.osm --elevation=../data/input/N47W122.hgt --traffic_signals=../data/input/traffic_signals.kml --residential_buildings=../data/input/residential_buildings.kml --parcels=../data/shapes/parcels --parking_shapes=../data/shapes/blockface --gtfs=../data/input/google_transit_2018_18_08 --neighborhoods=../data/input/neighborhoods.geojson --output=../data/raw_maps/montlake.abst --fast_dev
+# To run manually: cargo run -- --osm=../data/input/montlake.osm --elevation=../data/input/N47W122.hgt --traffic_signals=../data/input/traffic_signals.kml --residential_buildings=../data/input/residential_buildings.kml --parcels=../data/shapes/parcels --parking_shapes=../data/shapes/blockface --gtfs=../data/input/google_transit_2018_18_08 --neighborhoods=../data/input/neighborhoods.geojson --clip=../data/polygons/montlake.poly --output=../data/raw_maps/montlake.abst --fast_dev

map_model/src/make/buildings.rs

@@ -19,13 +19,7 @@ pub fn make_all_buildings(
     timer.start_iter("get building center points", input.len());
     for b in input {
         timer.next();
-        let pts = Pt2D::approx_dedupe(
-            b.points
-                .iter()
-                .map(|coord| Pt2D::from_gps(*coord, gps_bounds).unwrap())
-                .collect(),
-            geom::EPSILON_DIST,
-        );
+        let pts = Pt2D::approx_dedupe(gps_bounds.must_convert(&b.points), geom::EPSILON_DIST);
         let center: HashablePt2D = Pt2D::center(&pts).into();
         pts_per_bldg.push(pts);
         center_per_bldg.push(center);

map_model/src/make/half_map.rs

@@ -183,11 +183,7 @@ pub fn make_half_map(
         half_map.areas.push(Area {
             id: AreaID(idx),
             area_type: a.area_type,
-            points: a
-                .points
-                .iter()
-                .map(|coord| Pt2D::from_gps(*coord, &gps_bounds).unwrap())
-                .collect(),
+            points: gps_bounds.must_convert(&a.points),
             osm_tags: a.osm_tags.clone(),
             osm_id: a.osm_id,
         });

map_model/src/make/initial/mod.rs

@@ -100,12 +100,7 @@ impl InitialMap {
                 .roads
                 .insert(*stable_id);
 
-            let original_center_pts = PolyLine::new(
-                r.points
-                    .iter()
-                    .map(|coord| Pt2D::from_gps(*coord, &gps_bounds).unwrap())
-                    .collect(),
-            );
+            let original_center_pts = PolyLine::new(gps_bounds.must_convert(&r.points));
 
             let lane_specs = lane_specs::get_lane_specs(r, *stable_id, edits);
             let mut fwd_width = Distance::ZERO;

map_model/src/make/parcels.rs

@@ -17,13 +17,7 @@ pub fn make_all_parcels(
     let mut center_per_parcel: Vec<HashablePt2D> = Vec::new();
     let mut query: HashSet<HashablePt2D> = HashSet::new();
     for p in input {
-        let pts = Pt2D::approx_dedupe(
-            p.points
-                .iter()
-                .map(|coord| Pt2D::from_gps(*coord, gps_bounds).unwrap())
-                .collect(),
-            geom::EPSILON_DIST,
-        );
+        let pts = Pt2D::approx_dedupe(gps_bounds.must_convert(&p.points), geom::EPSILON_DIST);
         let center: HashablePt2D = Pt2D::center(&pts).into();
         pts_per_parcel.push(pts);
         center_per_parcel.push(center);

synthetic/src/model.rs

@@ -422,12 +422,7 @@ impl Model {
         /*for (idx, b) in data.buildings.iter().enumerate() {
             let b = Building {
                 label: None,
-                center: Pt2D::center(
-                    &b.points
-                        .iter()
-                        .map(|pt| Pt2D::from_gps(*pt, &gps_bounds).unwrap())
-                        .collect(),
-                ),
+                center: Pt2D::center(&gps_bounds.must_convert(&b.points)),
             };
             quadtree.insert_with_box(ID::Building(idx), b.polygon().get_bounds().as_bbox());
             m.buildings.insert(idx, b);