Add basic population layer from popgetter

backend/Cargo.toml

@@ -20,6 +20,7 @@ futures-util = { version ="0.3.30", default-features = false }
 geo = "0.28.0"
 geojson = { git = "https://github.com/georust/geojson", features = ["geo-types"] }
 geomedea = { git = "https://github.com/michaelkirk/geomedea", default-features = false }
+geozero = { version = "0.13.0", default-features = false, features = ["with-geo"] }
 itertools = "0.13.0"
 js-sys = "0.3.69"
 log = "0.4.20"
@@ -35,6 +36,7 @@ wasm-bindgen = "0.2.87"
 wasm-bindgen-futures = "0.4.42"
 web-sys = { version = "0.3.64", features = ["console"] }
 web-time = "1.1.0"
+flatgeobuf = "4.3.0"
 
 [target.'cfg(not(target_arch = "wasm32"))'.dependencies]
 geomedea = { git = "https://github.com/michaelkirk/geomedea" }

backend/src/graph/mod.rs

@@ -7,7 +7,7 @@ mod transit_route;
 
 use anyhow::Result;
 use enum_map::{Enum, EnumMap};
-use geo::{Coord, LineLocatePoint, LineString, Point, Polygon};
+use geo::{Coord, LineLocatePoint, LineString, MultiPolygon, Point, Polygon};
 use geojson::{Feature, GeoJson, Geometry};
 use rstar::{primitives::GeomWithData, RTree};
 use serde::{Deserialize, Serialize};
@@ -32,6 +32,8 @@ pub struct Graph {
     pub amenities: Vec<Amenity>,
 
     pub gtfs: GtfsModel,
+
+    pub zones: Vec<Zone>,
 }
 
 pub type EdgeLocation = GeomWithData<LineString, RoadID>;
@@ -184,6 +186,17 @@ impl Graph {
             intersection,
         }
     }
+
+    /// Returns a GeoJSON string
+    pub fn render_zones(&self) -> Result<String> {
+        let mut features = Vec::new();
+        for zone in &self.zones {
+            let mut f = Feature::from(Geometry::from(&self.mercator.to_wgs84(&zone.geom)));
+            f.set_property("population", zone.population);
+            features.push(f);
+        }
+        Ok(serde_json::to_string(&GeoJson::from(features))?)
+    }
 }
 
 impl Road {
@@ -238,3 +251,10 @@ pub enum PathStep {
         stop2: StopID,
     },
 }
+
+#[derive(Serialize, Deserialize)]
+pub struct Zone {
+    pub geom: MultiPolygon,
+    // TODO Later on, this could be generic or user-supplied
+    pub population: u32,
+}

backend/src/graph/scrape.rs

@@ -2,17 +2,18 @@ use std::collections::HashMap;
 
 use anyhow::Result;
 use enum_map::EnumMap;
-use geo::{Coord, EuclideanLength};
+use flatgeobuf::{FeatureProperties, FgbFeature, GeozeroGeometry, HttpFgbReader};
+use geo::{Coord, EuclideanLength, MultiPolygon};
 use muv_osm::{AccessLevel, TMode};
 use osm_reader::OsmID;
 use rstar::RTree;
-use utils::Tags;
+use utils::{Mercator, Tags};
 
 use super::amenity::Amenity;
 use super::route::Router;
 use crate::graph::{
     AmenityID, Direction, EdgeLocation, Graph, GtfsSource, Intersection, IntersectionID, Mode,
-    Road, RoadID,
+    Road, RoadID, Zone,
 };
 use crate::gtfs::{GtfsModel, StopID};
 use crate::timer::Timer;
@@ -55,6 +56,7 @@ impl Graph {
     pub async fn new(
         input_bytes: &[u8],
         gtfs_source: GtfsSource,
+        population_url: Option<String>,
         mut timer: Timer,
     ) -> Result<Graph> {
         timer.step("parse OSM and split graph");
@@ -143,6 +145,13 @@ impl Graph {
         snap_stops(&mut roads, &mut gtfs, &closest_road[Mode::Foot], &mut timer);
         timer.pop();
 
+        let zones = if let Some(url) = population_url {
+            timer.step("load population zones");
+            load_zones(url, &graph.mercator).await?
+        } else {
+            Vec::new()
+        };
+
         timer.done();
 
         Ok(Graph {
@@ -155,6 +164,7 @@ impl Graph {
 
             amenities: amenities.amenities,
             gtfs,
+            zones,
         })
     }
 }
@@ -280,3 +290,34 @@ fn snap_stops(
         }
     }
 }
+
+async fn load_zones(url: String, mercator: &Mercator) -> Result<Vec<Zone>> {
+    let bbox = mercator.wgs84_bounds;
+    let mut fgb = HttpFgbReader::open(&url)
+        .await?
+        .select_bbox(bbox.min().x, bbox.min().y, bbox.max().x, bbox.max().y)
+        .await?;
+
+    let mut zones = Vec::new();
+    while let Some(feature) = fgb.next().await? {
+        // TODO Could intersect with boundary_polygon, but some extras nearby won't hurt anything
+        let mut geom = get_multipolygon(feature)?;
+        mercator.to_mercator_in_place(&mut geom);
+        zones.push(Zone {
+            geom,
+            // TODO Re-encode as UInt
+            population: feature.property::<i64>("population")?.try_into()?,
+        });
+    }
+    Ok(zones)
+}
+
+fn get_multipolygon(f: &FgbFeature) -> Result<MultiPolygon> {
+    let mut p = geozero::geo_types::GeoWriter::new();
+    f.process_geom(&mut p)?;
+    match p.take_geometry().unwrap() {
+        geo::Geometry::Polygon(p) => Ok(MultiPolygon(vec![p])),
+        geo::Geometry::MultiPolygon(mp) => Ok(mp),
+        _ => bail!("Wrong type in fgb"),
+    }
+}

backend/src/lib.rs

@@ -12,7 +12,7 @@ use geo::Coord;
 use serde::Deserialize;
 use wasm_bindgen::prelude::*;
 
-pub use graph::{Graph, Mode};
+pub use graph::{Graph, GtfsSource, Mode};
 pub use gtfs::GtfsModel;
 pub use timer::Timer;
 
@@ -40,6 +40,7 @@ impl MapModel {
         input_bytes: &[u8],
         is_osm: bool,
         gtfs_url: Option<String>,
+        population_url: Option<String>,
         progress_cb: Option<js_sys::Function>,
     ) -> Result<MapModel, JsValue> {
         // Panics shouldn't happen, but if they do, console.log them.
@@ -54,9 +55,14 @@ impl MapModel {
             None => graph::GtfsSource::None,
         };
         let graph = if is_osm {
-            Graph::new(input_bytes, gtfs, Timer::new("build graph", progress_cb))
-                .await
-                .map_err(err_to_js)?
+            Graph::new(
+                input_bytes,
+                gtfs,
+                population_url,
+                Timer::new("build graph", progress_cb),
+            )
+            .await
+            .map_err(err_to_js)?
         } else {
             bincode::deserialize_from(input_bytes).map_err(err_to_js)?
         };
@@ -88,6 +94,11 @@ impl MapModel {
         vec![b.min().x, b.min().y, b.max().x, b.max().y]
     }
 
+    #[wasm_bindgen(js_name = renderZones)]
+    pub fn render_zones(&self) -> Result<String, JsValue> {
+        self.graph.render_zones().map_err(err_to_js)
+    }
+
     #[wasm_bindgen(js_name = isochrone)]
     pub fn isochrone(&self, input: JsValue) -> Result<String, JsValue> {
         let req: IsochroneRequest = serde_wasm_bindgen::from_value(input)?;

backend/src/main.rs

@@ -23,6 +23,9 @@ async fn main() -> Result<()> {
         std::process::exit(1);
     }
 
+    // TODO Hardcoded for now
+    let population_url = Some("https://assets.od2net.org/population.fgb".to_string());
+
     if args[1] == "graph" {
         let timer = Timer::new("build graph", None);
         let osm_bytes = std::fs::read(&args[2])?;
@@ -31,7 +34,7 @@ async fn main() -> Result<()> {
             None => GtfsSource::None,
         };
 
-        let graph = Graph::new(&osm_bytes, gtfs, timer).await?;
+        let graph = Graph::new(&osm_bytes, gtfs, population_url, timer).await?;
         let writer = BufWriter::new(File::create("graph.bin")?);
         bincode::serialize_into(writer, &graph)?;
     } else if args[1] == "gmd" {

data_prep/.gitignore

@@ -0,0 +1 @@
+*.geojson

data_prep/population.sh

@@ -0,0 +1,64 @@
+#!/bin/bash
+# This script prepares a consolidated FGB file with total population count in
+# the most granular zone for different countries, using
+# https://github.com/Urban-Analytics-Technology-Platform/popgetter. This script
+# / process might live elsewhere at some point.
+
+set -e
+set -x
+
+# You need to build https://github.com/Urban-Analytics-Technology-Platform/popgetter-cli
+POPGETTER=/home/dabreegster/Downloads/popgetter-cli/target/release/popgetter
+
+# The metrics and coordinate systems were figured out manually. In the future,
+# this'll be easier with a popgetter web UI.
+
+# England
+$POPGETTER data \
+        --force-run \
+        --output-format geojson \
+        --output-file england_raw.geojson \
+        --geometry-level oa \
+        --id 1355501cf6f3b1fa8cf6a100c98f330d51a3382ed2111fef0d2fff446608a428
+mapshaper england_raw.geojson \
+        -rename-fields population='Residence type: Total; measures: Value' \
+        -each 'delete GEO_ID' \
+        -proj init=EPSG:27700 crs=wgs84 \
+        -o england.geojson
+
+# Belgium
+$POPGETTER data \
+        --force-run \
+        --output-format geojson \
+        --output-file belgium_raw.geojson \
+        --geometry-level statistical_sector \
+        --id fcf09809889c1d9715bff5f825b0c6ed4d9286f2e2b4948839accc29c15e98c5
+mapshaper belgium_raw.geojson \
+        -rename-fields population='TOTAL' \
+        -each 'delete GEO_ID' \
+        -proj init=EPSG:3812 crs=wgs84 \
+        -o belgium.geojson
+
+# USA
+# TODO This might not be the right variable, and it's only at block_group
+$POPGETTER data \
+        --force-run \
+        --output-format geojson \
+        --output-file usa_raw.geojson \
+        --geometry-level block_group \
+        --id d23e348af6ab03265b4f258178edc6b509651095f81b965c1a62396fe463d0f6
+mapshaper usa_raw.geojson \
+        -rename-fields population='B01001_E001' \
+        -each 'delete GEO_ID' \
+        -o usa.geojson
+
+# Scotland
+# TODO
+
+# Northern Ireland
+# TODO
+
+# Merge files. You need to build https://github.com/acteng/will-it-fit/tree/main/data_prep/merge_files
+MERGER=/home/dabreegster/will-it-fit/data_prep/merge_files/target/release/merge_files
+$MERGER england.geojson belgium.geojson usa.geojson
+# Hosting: mv out.fgb ~/cloudflare_sync/population.fgb, sync it

web/src/App.svelte

@@ -29,11 +29,13 @@
     showAbout,
     routeA,
     routeB,
+    showPopulation,
   } from "./stores";
   import TitleMode from "./title/TitleMode.svelte";
   import workerWrapper from "./worker?worker";
   import { type Backend } from "./worker";
   import * as Comlink from "comlink";
+        import { PopulationLayer } from "./common";
 
   onMount(async () => {
     // If you get "import declarations may only appear at top level of a
@@ -166,6 +168,10 @@
         {:else if $mode.kind == "buffer-route"}
           <BufferRouteMode gj={$mode.gj} />
         {/if}
+
+        {#if $showPopulation}
+                <PopulationLayer />
+        {/if}
       {/if}
     </MapLibre>
   </div>

web/src/BufferRouteMode.svelte

@@ -25,7 +25,7 @@
     <GeoJSON data={gj}>
       <LineLayer
         paint={{
-          "line-width": 20,
+          "line-width": ["case", ["==", ["get", "kind"], "route"], 20, 3],
           "line-color": [
             "case",
             ["==", ["get", "kind"], "route"],

web/src/IsochroneMode.svelte

@@ -26,7 +26,7 @@
       lat: lerp(0.5, bbox[1], bbox[3]),
     };
   });
-  let contours = true;
+  let contours = false;
 
   let isochroneGj: FeatureCollection | null = null;
   let routeGj: FeatureCollection | null = null;
@@ -142,7 +142,7 @@
           id="isochrone"
           filter={isLine}
           paint={{
-            "line-width": 20,
+            "line-width": 2,
             "line-color": makeColorRamp(
               ["get", "cost_seconds"],
               limitsSeconds,

web/src/common/NavBar.svelte

@@ -1,5 +1,5 @@
 <script lang="ts">
-  import { mode } from "../stores";
+  import { mode, showPopulation } from "../stores";
 </script>
 
 <nav>
@@ -39,3 +39,6 @@
     </li>
   </ul>
 </nav>
+
+<label><input type="checkbox" bind:checked={$showPopulation} />Population</label
+>

web/src/common/PopulationLayer.svelte

@@ -0,0 +1,22 @@
+<script lang="ts">
+  import { GeoJSON, LineLayer, FillLayer, hoverStateFilter } from "svelte-maplibre";
+  import { notNull } from "svelte-utils";
+  import { Popup } from "svelte-utils/map";
+  import { backend } from "../stores";
+</script>
+
+{#await notNull($backend).renderZones() then data}
+  <GeoJSON {data} generateId>
+    <FillLayer
+      manageHoverState
+      paint={{
+        "fill-color": "red",
+        "fill-opacity": hoverStateFilter(0.5, 0.8),
+      }}
+    >
+      <Popup openOn="hover" let:props>{props.population.toLocaleString()}</Popup
+      >
+    </FillLayer>
+    <LineLayer paint={{"line-color": "black", "line-width": 1}} />
+  </GeoJSON>
+{/await}

web/src/common/index.ts

@@ -4,4 +4,5 @@ export { default as Loading } from "./Loading.svelte";
 export { default as NavBar } from "./NavBar.svelte";
 export { default as PickAmenityKinds } from "./PickAmenityKinds.svelte";
 export { default as PickTravelMode } from "./PickTravelMode.svelte";
+export { default as PopulationLayer } from "./PopulationLayer.svelte";
 export { default as StopsLayer } from "./StopsLayer.svelte";

web/src/stores.ts

@@ -28,6 +28,7 @@ export type Mode =
 export let mode: Writable<Mode> = writable({ kind: "title" });
 export let map: Writable<Map | null> = writable(null);
 export let showAbout: Writable<boolean> = writable(true);
+export let showPopulation: Writable<boolean> = writable(false);
 
 export type TravelMode = "car" | "bicycle" | "foot" | "transit";