census areas import (#425)

#424

Cargo.toml

@@ -31,3 +31,6 @@ members = [
 # update dependencies often).
 [profile.dev.package."*"]
 opt-level = 3
+
+[patch.crates-io]
+geojson = { git = "https://github.com/georust/geojson", branch = "mkirk/try_from" }

fifteen_min/Cargo.toml

@@ -9,8 +9,8 @@ default = ["map_gui/native", "widgetry/native-backend"]
 
 [dependencies]
 abstutil = { path = "../abstutil" }
-contour = { git = "https://github.com/dabreegster/contour-rs" }
-geojson = "0.20.1"
+contour = "0.3.0"
+geojson = "0.21.0"
 geom = { path = "../geom" }
 log = "0.4"
 map_gui = { path = "../map_gui" }

game/Cargo.toml

@@ -21,10 +21,10 @@ built = { version = "0.4.3", optional = true, features=["chrono"] }
 chrono = "0.4.15"
 collisions = { path = "../collisions" }
 colorous = "1.0.3"
-contour = { git = "https://github.com/dabreegster/contour-rs" }
+contour = "0.3.0"
 downcast-rs = "1.2.0"
 enumset = "1.0.1"
-geojson = "0.20.1"
+geojson = "0.21.0"
 geom = { path = "../geom" }
 instant = "0.1.7"
 kml = { path = "../kml" }

game/src/sandbox/gameplay/census.rs

@@ -44,7 +44,8 @@ impl State<App> for CensusGenerator {
                         config,
                         &app.primary.map,
                         &mut app.primary.current_flags.sim_flags.make_rng(),
-                    );
+                    )
+                    .expect("failed to generate scenario");
                     // TODO Do something with it -- save it, launch it in sandboxmode, display some
                     // stats about it?
                     return Transition::Pop;

geom/Cargo.toml

@@ -10,7 +10,7 @@ abstutil = { path = "../abstutil" }
 earcutr = { git = "https://github.com/donbright/earcutr" }
 geo = "0.15.0"
 geo-booleanop = "= 0.3.2"
-geojson = "0.20.1"
+geojson = "0.21.0"
 histogram = "0.6.9"
 instant = "0.1.7"
 ordered-float = { version = "2.0.0", features=["serde"] }

geom/src/polygon.rs

@@ -304,7 +304,7 @@ impl Polygon {
 
     pub fn convex_hull(list: Vec<Polygon>) -> Polygon {
         let mp: geo::MultiPolygon<f64> = list.into_iter().map(|p| to_geo(p.points())).collect();
-        from_geo(mp.convex_hull())
+        mp.convex_hull().into()
     }
 
     pub fn polylabel(&self) -> Pt2D {
@@ -476,19 +476,37 @@ fn to_geo(pts: &Vec<Pt2D>) -> geo::Polygon<f64> {
     )
 }
 
-fn from_geo(p: geo::Polygon<f64>) -> Polygon {
-    Polygon::buggy_new(
-        p.into_inner()
-            .0
-            .into_points()
-            .into_iter()
-            .map(|pt| Pt2D::new(pt.x(), pt.y()))
-            .collect(),
-    )
+impl From<geo::Polygon<f64>> for Polygon {
+    fn from(poly: geo::Polygon<f64>) -> Self {
+        let (exterior, interiors) = poly.into_inner();
+        Polygon::with_holes(
+            Ring::from(exterior),
+            interiors.into_iter().map(Ring::from).collect(),
+        )
+    }
+}
+
+impl From<Polygon> for geo::Polygon<f64> {
+    fn from(poly: Polygon) -> Self {
+        if let Some(mut rings) = poly.rings {
+            let exterior = rings.pop().expect("expected poly.rings[0] to be exterior");
+            let interiors: Vec<geo::LineString<f64>> =
+                rings.into_iter().map(geo::LineString::from).collect();
+            Self::new(exterior.into(), interiors)
+        } else {
+            let exterior_coords = poly
+                .points
+                .into_iter()
+                .map(geo::Coordinate::from)
+                .collect::<Vec<_>>();
+            let exterior = geo::LineString(exterior_coords);
+            Self::new(exterior, Vec::new())
+        }
+    }
 }
 
 fn from_multi(multi: geo::MultiPolygon<f64>) -> Vec<Polygon> {
-    multi.into_iter().map(from_geo).collect()
+    multi.into_iter().map(Polygon::from).collect()
 }
 
 fn downsize(input: Vec<usize>) -> Vec<u16> {

geom/src/pt.rs

@@ -182,3 +182,30 @@ impl HashablePt2D {
         Pt2D::new(self.x_nan.into_inner(), self.y_nan.into_inner())
     }
 }
+
+impl From<Pt2D> for geo::Coordinate<f64> {
+    fn from(pt: Pt2D) -> Self {
+        geo::Coordinate {
+            x: pt.inner_x,
+            y: pt.inner_y,
+        }
+    }
+}
+
+impl From<Pt2D> for geo::Point<f64> {
+    fn from(pt: Pt2D) -> Self {
+        geo::Point::new(pt.inner_x, pt.inner_y)
+    }
+}
+
+impl From<geo::Coordinate<f64>> for Pt2D {
+    fn from(coord: geo::Coordinate<f64>) -> Self {
+        Pt2D::new(coord.x, coord.y)
+    }
+}
+
+impl From<geo::Point<f64>> for Pt2D {
+    fn from(point: geo::Point<f64>) -> Self {
+        Pt2D::new(point.x(), point.y())
+    }
+}

geom/src/ring.rs

@@ -188,3 +188,20 @@ impl fmt::Display for Ring {
         write!(f, "])")
     }
 }
+
+impl From<Ring> for geo::LineString<f64> {
+    fn from(ring: Ring) -> Self {
+        let coords = ring
+            .pts
+            .into_iter()
+            .map(geo::Coordinate::from)
+            .collect::<Vec<_>>();
+        Self(coords)
+    }
+}
+
+impl From<geo::LineString<f64>> for Ring {
+    fn from(line_string: geo::LineString<f64>) -> Self {
+        Self::must_new(line_string.0.into_iter().map(Pt2D::from).collect())
+    }
+}

headless/Cargo.toml

@@ -6,7 +6,7 @@ edition = "2018"
 
 [dependencies]
 abstutil = { path = "../abstutil" }
-geojson = "0.20.1"
+geojson = "0.21.0"
 geom = { path = "../geom" }
 hyper = "0.13.9"
 lazy_static = "1.4.0"

importer/Cargo.toml

@@ -14,7 +14,7 @@ abstutil = { path = "../abstutil" }
 collisions = { path = "../collisions" }
 convert_osm = { path = "../convert_osm" }
 csv = "1.1.4"
-geojson = "0.20.1"
+geojson = "0.21.0"
 geom = { path = "../geom" }
 gdal = { version = "0.6.0", optional = true }
 kml = { path = "../kml" }

map_gui/Cargo.toml

@@ -16,11 +16,11 @@ release_s3 = []
 aabb-quadtree = "0.1.0"
 abstutil = { path = "../abstutil" }
 colorous = "1.0.3"
-contour = { git = "https://github.com/dabreegster/contour-rs" }
+contour = "0.3.0"
 flate2 = "1.0.19"
 futures = { version = "0.3.8", optional = true }
 futures-channel = { version = "0.3.8", optional = true }
-geojson = "0.20.1"
+geojson = "0.21.0"
 geom = { path = "../geom" }
 instant = "0.1.7"
 js-sys = { version = "0.3.45", optional = true }

popdat/Cargo.toml

@@ -6,9 +6,13 @@ edition = "2018"
 
 [dependencies]
 abstutil = { path = "../abstutil" }
+anyhow = "1.0.35"
+geo = "0.16.0"
+geojson = { version = "0.21.0", features = ["geo-types"] }
 geom = { path = "../geom" }
 log = "0.4.11"
 map_model = { path = "../map_model" }
 rand = "0.7.0"
 rand_xorshift = "0.2.0"
+serde_json = "1.0.60"
 sim = { path = "../sim" }

popdat/src/distribute_people.rs

@@ -28,7 +28,7 @@ pub fn assign_people_to_houses(
         let pct_overlap = Polygon::union_all(area.polygon.intersection(map.get_boundary_polygon()))
             .area()
             / area.polygon.area();
-        let num_residents = (pct_overlap * (area.total_population as f64)) as usize;
+        let num_residents = (pct_overlap * (area.population as f64)) as usize;
         debug!(
             "Distributing {} residents to {} buildings. {}% of this area overlapped with the map, \
              scaled residents accordingly.",

popdat/src/import_census.rs

@@ -1,11 +1,129 @@
+use std::convert::TryFrom;
+
+use geo::algorithm::intersects::Intersects;
+use geojson::GeoJson;
+
+use abstutil::Timer;
 use map_model::Map;
 
 use crate::CensusArea;
 
 impl CensusArea {
-    pub fn find_data_for_map(_map: &Map) -> Result<Vec<CensusArea>, String> {
-        // TODO importer/src/utils.rs has a download() helper that we could copy here. (And later
-        // dedupe, after deciding how this crate will integrate with the importer)
-        todo!()
+    pub fn find_data_for_map(map: &Map, timer: &mut Timer) -> anyhow::Result<Vec<CensusArea>> {
+        // TODO eventually it'd be nice to lazily download the info needed. For now we expect a
+        // prepared geojson file to exist in data/system/<city>/population_areas.geojson
+        //
+        // expected geojson formatted contents like:
+        // {
+        //     "type": "FeatureCollection",
+        //     "features": [
+        //          {
+        //              "type": "Feature",
+        //              "properties": { "population": 123 },
+        //              "geometry": {
+        //                  "type": "MultiPolygon",
+        //                  "coordinates": [ [ [ [ -73.7, 40.8 ], [ -73.7, 40.8 ], ...] ] ] ]
+        //               }
+        //          },
+        //          {
+        //              "type": "Feature",
+        //              "properties": { "population": 249 },
+        //              "geometry": {
+        //                  "type": "MultiPolygon",
+        //                  "coordinates": [ [ [ [ -73.8, 40.8 ], [ -73.8, 40.8 ], ...] ] ]
+        //               }
+        //           },
+        //          ...
+        //      ]
+        // }
+        //
+        // Note: intuitively you might expect a collection of Polygon's rather than  MultiPolygons,
+        // but anecdotally, the census data I've seen uses MultiPolygons. In practice almost
+        // all are MultiPoly's with just one element, but some truly have multiple polygons.
+        //
+        // When we implement downloading, importer/src/utils.rs has a download() helper that we
+        // could copy here. (And later dedupe, after deciding how this crate will integrate with
+        // the importer)
+        let path = abstutil::path(format!(
+            "system/{}/population_areas.geojson",
+            map.get_name().city
+        ));
+        let bytes = abstutil::slurp_file(&path).map_err(|s| anyhow!(s))?;
+        debug!("parsing geojson at path: {}", &path);
+        timer.start("parsing geojson");
+        let geojson = GeoJson::from_reader(&*bytes)?;
+        timer.stop("parsing geojson");
+        let mut results = vec![];
+        let collection = geojson::FeatureCollection::try_from(geojson)?;
+
+        let map_area = map.get_boundary_polygon();
+        let bounds = map.get_gps_bounds();
+
+        use geo::algorithm::map_coords::MapCoordsInplace;
+        let mut geo_map_area: geo::Polygon<_> = map_area.clone().into();
+        geo_map_area.map_coords_inplace(|c| {
+            let projected = geom::Pt2D::new(c.0, c.1).to_gps(bounds);
+            (projected.x(), projected.y())
+        });
+
+        debug!("collection.features: {}", &collection.features.len());
+        timer.start("converting to `CensusArea`s");
+        for feature in collection.features.into_iter() {
+            let population = match feature
+                .property("population")
+                .ok_or(anyhow!("missing 'population' property"))?
+            {
+                serde_json::Value::Number(n) => n
+                    .as_u64()
+                    .ok_or(anyhow!("unexpected population number: {:?}", n))?
+                    as usize,
+                _ => {
+                    bail!(
+                        "unexpected format for 'population': {:?}",
+                        feature.property("population")
+                    );
+                }
+            };
+
+            let geometry = feature
+                .geometry
+                .ok_or(anyhow!("geojson feature missing geometry"))?;
+            let mut multi_poly = geo::MultiPolygon::<f64>::try_from(geometry.value)?;
+            let mut geo_polygon = multi_poly
+                .0
+                .pop()
+                .ok_or(anyhow!("multipolygon was unexpectedly empty"))?;
+            if !multi_poly.0.is_empty() {
+                // I haven't looked into why this is happening - but intuitively a census area could
+                // include separate polygons - e.g. across bodies of water. In practice they are a
+                // vast minority, so we naively just take the first one for now.
+                warn!(
+                    "dropping {} polygons from census area with multiple polygons",
+                    multi_poly.0.len()
+                );
+            }
+
+            if !geo_polygon.intersects(&geo_map_area) {
+                debug!(
+                    "skipping polygon outside of map area. polygon: {:?}, map_area: {:?}",
+                    geo_polygon, geo_map_area
+                );
+                continue;
+            }
+
+            geo_polygon.map_coords_inplace(|(x, y)| {
+                let point = geom::LonLat::new(*x, *y).to_pt(bounds);
+                (point.x(), point.y())
+            });
+            let polygon = geom::Polygon::from(geo_polygon);
+            results.push(CensusArea {
+                polygon,
+                population,
+            });
+        }
+        debug!("built {} CensusAreas within map bounds", results.len());
+        timer.stop("converting to `CensusArea`s");
+
+        Ok(results)
     }
 }

popdat/src/lib.rs

@@ -17,11 +17,14 @@
 //!    different Activities throughout the day.
 //! 4) Pick specific buildings to visit to satisfy the Schedule.
 
+#[macro_use]
+extern crate anyhow;
 #[macro_use]
 extern crate log;
 
 use rand_xorshift::XorShiftRng;
 
+use abstutil::Timer;
 use geom::Polygon;
 use geom::{Distance, Time};
 use map_model::{BuildingID, Map};
@@ -38,7 +41,7 @@ mod make_person;
 /// have two overlapping areas.
 pub struct CensusArea {
     pub polygon: Polygon,
-    pub total_population: usize,
+    pub population: usize,
     // TODO Not sure what goes here, whatever census data actually has that could be useful
 }
 
@@ -105,16 +108,24 @@ pub fn generate_scenario(
 ) -> Result<Scenario, String> {
     // find_data_for_map may return an error. If so, just plumb it back to the caller using the ?
     // operator
-    let areas = CensusArea::find_data_for_map(map)?;
+    let mut timer = Timer::new("generate census scenario");
+    timer.start("building population areas for map");
+    let areas = CensusArea::find_data_for_map(map, &mut timer).map_err(|e| e.to_string())?;
+    timer.stop("building population areas for map");
+
+    timer.start("assigning people to houses");
     let people = distribute_people::assign_people_to_houses(areas, map, rng, &config);
+    timer.stop("assigning people to houses");
 
     let mut scenario = Scenario::empty(map, scenario_name);
+    timer.start("building people");
     for person in people {
         // TODO If we need to parallelize because make_person is slow, the sim crate has a fork_rng
         // method that could be useful
         scenario
             .people
             .push(make_person::make_person(person, map, rng, &config));
     }
+    timer.stop("building people");
     Ok(scenario)
 }