census areas import #425
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -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(); | ||
| mp.convex_hull().into() | ||
| } | ||
|
|
||
| pub fn polylabel(&self) -> Pt2D { | ||
|
|
@@ -476,19 +476,37 @@ fn to_geo(pts: &Vec<Pt2D>) -> geo::Polygon<f64> { | |
| ) | ||
| } | ||
|
|
||
| impl From<geo::Polygon<f64>> for Polygon { | ||
|
There was a problem hiding this comment. Are you ok with adding the idiomatic conversation traits? It's a little easier to discover this way since I know where to look rather than custom named methods. I could add more for other geo types, but only added the minimum of what we needed for now. There was a problem hiding this comment. Yes, these're better! I probably wrote the old thing before I understood |
||
| fn from(poly: geo::Polygon<f64>) -> Self { | ||
| let (exterior, interiors) = poly.into_inner(); | ||
| Polygon::with_holes( | ||
|
There was a problem hiding this comment. I wonder why I was using |
||
| Ring::from(exterior), | ||
| interiors.into_iter().map(Ring::from).collect(), | ||
| ) | ||
| } | ||
| } | ||
|
|
||
| impl From<Polygon> for geo::Polygon<f64> { | ||
|
There was a problem hiding this comment. Woo, thanks for fixing this! When I originally wrote this, I didn't know about earcutr for triangulating holes, and I think I confused multipolygons / polygons with inner portions. There was a problem hiding this comment. I've updated these |
||
| 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 | ||
|
There was a problem hiding this comment. Note this will break for a polygon produced by I think this is fine for now. Ideally every |
||
| .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(Polygon::from).collect() | ||
| } | ||
|
|
||
| fn downsize(input: Vec<usize>) -> Vec<u16> { | ||
|
|
||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,11 +1,126 @@ | ||
| use crate::CensusArea; | ||
| use abstutil::Timer; | ||
|
There was a problem hiding this comment. nit: For consistency, we should still split the imports based on category... std, external crates, other crates in this workspace, the current crate. I guess we don't have any contributors using an IDE that forces this order, so not a huge deal either way. There was a problem hiding this comment. happy to abide |
||
| use geo::algorithm::intersects::Intersects; | ||
| use geojson::GeoJson; | ||
| use map_model::Map; | ||
| use std::convert::TryFrom; | ||
|
|
||
| impl CensusArea { | ||
| 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 | ||
|
There was a problem hiding this comment.
|
||
| // | ||
| // 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))?; | ||
|
There was a problem hiding this comment. Ah, I guess |
||
| 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); | ||
|
There was a problem hiding this comment. Is there a better way to get the map poly in wgs84? Though maybe it's irrelevant, because the other question I need to answer is: is there a blessed way to convert a poly from wgs84 to map coords? There was a problem hiding this comment.
|
||
| (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) | ||
| } | ||
| } | ||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
in general I think it's a little nicer to patch things here at the top level rather than updating across all the individual crates. It's maybe a little more convenient, but more importantly, it works for transitive deps too, which you can't conveniently edit.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
+1. I always forget about this trick