popdat/
import_census.rs

1use anyhow::Result;
2use geo::{BoundingRect, Intersects, MapCoordsInPlace};
3
4use geom::{GPSBounds, Polygon};
5
6use crate::CensusArea;
7
8impl CensusArea {
9    pub async fn fetch_all_for_map(
10        map_area: &Polygon,
11        bounds: &GPSBounds,
12    ) -> Result<Vec<CensusArea>> {
13        use flatgeobuf::HttpFgbReader;
14        use geozero::geo_types::GeoWriter;
15
16        let mut geo_map_area: geo::Polygon = map_area.clone().into();
17        geo_map_area.map_coords_in_place(|c| {
18            let projected = geom::Pt2D::new(c.x, c.y).to_gps(bounds);
19            (projected.x(), projected.y()).into()
20        });
21
22        let bounding_rect = geo_map_area
23            .bounding_rect()
24            .ok_or_else(|| anyhow!("missing bound rect"))?;
25
26        // See the import handbook for how to prepare this file.
27        let mut fgb = HttpFgbReader::open("https://abstreet.s3.amazonaws.com/population_areas.fgb")
28            .await?
29            .select_bbox(
30                bounding_rect.min().x,
31                bounding_rect.min().y,
32                bounding_rect.max().x,
33                bounding_rect.max().y,
34            )
35            .await?;
36
37        let mut results = vec![];
38        while let Some(feature) = fgb.next().await? {
39            use flatgeobuf::FeatureProperties;
40            // PERF TODO: how to parse into usize directly? And avoid parsing entire props dict?
41            let props = feature.properties()?;
42            if !props.contains_key("population") {
43                warn!("skipping feature with missing population");
44                continue;
45            }
46            let population: usize = props["population"].parse()?;
47            let geometry = match feature.geometry() {
48                Some(g) => g,
49                None => {
50                    warn!("skipping feature with missing geometry");
51                    continue;
52                }
53            };
54            let mut geo = GeoWriter::new();
55            geometry.process(&mut geo, flatgeobuf::GeometryType::MultiPolygon)?;
56            if let Some(geo::Geometry::MultiPolygon(multi_poly)) = geo.take_geometry() {
57                let geo_polygon = multi_poly
58                    .0
59                    .first()
60                    .ok_or_else(|| anyhow!("multipolygon was unexpectedly empty"))?;
61                if multi_poly.0.len() > 1 {
62                    warn!(
63                        "dropping {} extra polygons from census area: {:?}",
64                        multi_poly.0.len() - 1,
65                        props
66                    );
67                }
68
69                if !geo_polygon.intersects(&geo_map_area) {
70                    debug!(
71                        "skipping polygon outside of map area. polygon: {:?}, map_area: {:?}",
72                        geo_polygon, geo_map_area
73                    );
74                    continue;
75                }
76
77                let mut polygon = geo_polygon.clone();
78                polygon.map_coords_in_place(|c| geom::LonLat::new(c.x, c.y).to_pt(bounds).into());
79                results.push(CensusArea {
80                    polygon,
81                    population,
82                });
83            } else {
84                warn!("skipping unexpected geometry");
85                continue;
86            }
87        }
88
89        Ok(results)
90    }
91}
popdat/import_census.rs

popdat/
import_census.rs
