Skip to content

Commit

Permalink
Read optional weights per point
Browse files Browse the repository at this point in the history
  • Loading branch information
@dabreegster
dabreegster committed Sep 12, 2024
1 parent 859722a commit 343d3a6
Show file tree
Hide file tree
Showing 3 changed files with 76 additions and 43 deletions.
9 changes: 6 additions & 3 deletions od2net/src/config.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,9 +45,11 @@ pub enum ODPattern {
/// Path to a CSV file that must have 3 columns "from", "to", and "count". The first
/// two must match zone names. "count" must be an integer.
csv_path: String,
/// If a zone doesn't have any matching origin points, use the zone's centroid instead.
/// If a zone doesn't have any matching origin points, use the zone's centroid instead. The
/// centroid weight will be 1.
origin_zone_centroid_fallback: bool,
/// If a zone doesn't have any matching destination points, use the zone's centroid instead.
/// If a zone doesn't have any matching destination points, use the zone's centroid
/// instead. The centroid weight will be 1.
destination_zone_centroid_fallback: bool,
},
ZoneToPoint {
Expand All @@ -58,7 +60,8 @@ pub enum ODPattern {
csv_path: String,
/// Path to a GeoJSON file containing Points with a "name" property
destinations_path: String,
/// If a zone doesn't have any matching origin points, use the zone's centroid instead.
/// If a zone doesn't have any matching origin points, use the zone's centroid instead. The
/// centroid weight will be 1.
origin_zone_centroid_fallback: bool,
},
/// Just read GeoJSON LineStrings from this path
Expand Down
3 changes: 1 addition & 2 deletions od2net/src/network/create_from_osm.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -255,8 +255,7 @@ fn split_edges(nodes: HashMap<NodeID, Position>, ways: HashMap<WayID, Way>) -> N
}

fn calculate_length_meters(pts: &[Position]) -> f64 {
let line_string =
LineString::<f64>::from(pts.iter().map(|pt| pt.to_degrees()).collect::<Vec<_>>());
let line_string = LineString::from(pts.iter().map(|pt| pt.to_degrees()).collect::<Vec<_>>());
line_string.haversine_length()
}

Expand Down
107 changes: 69 additions & 38 deletions od2net/src/od.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,11 @@ use std::io::BufReader;

use anyhow::Result;
use fs_err::File;
use geo::{BoundingRect, Centroid, Contains, MultiPolygon};
use geo::{BoundingRect, Centroid, Contains, MultiPolygon, Point};
use geojson::{FeatureReader, Value};
use indicatif::HumanCount;
use nanorand::{Rng, WyRand};
use rstar::{RTree, AABB};
use rstar::{PointDistance, RTree, RTreeObject, AABB};
use serde::Deserialize;

use super::config::{ODPattern, Requests};
Expand Down Expand Up @@ -53,10 +53,10 @@ pub fn generate_requests(
));
for pt in origins {
requests.push(Request {
x1: pt.0,
y1: pt.1,
x2: destinations[0].0,
y2: destinations[0].1,
x1: pt.lon,
y1: pt.lat,
x2: destinations[0].lon,
y2: destinations[0].lat,
});
}
timer.stop();
Expand All @@ -70,12 +70,12 @@ pub fn generate_requests(
HumanCount(origins.len() as u64),
));
for pt in origins {
let goto = closest.nearest_neighbor(&pt).unwrap();
let goto = closest.nearest_neighbor(&(pt.lon, pt.lat)).unwrap();
requests.push(Request {
x1: pt.0,
y1: pt.1,
x2: goto.0,
y2: goto.1,
x1: pt.lon,
y1: pt.lat,
x2: goto.lon,
y2: goto.lat,
});
}
timer.stop();
Expand Down Expand Up @@ -122,10 +122,10 @@ pub fn generate_requests(
}
};
requests.push(Request {
x1: from.0,
y1: from.1,
x2: to.0,
y2: to.1,
x1: from.lon,
y1: from.lat,
x2: to.lon,
y2: to.lat,
});
}
}
Expand Down Expand Up @@ -172,8 +172,8 @@ pub fn generate_requests(
}
};
requests.push(Request {
x1: from.0,
y1: from.1,
x1: from.lon,
y1: from.lat,
x2: to.0,
y2: to.1,
});
Expand All @@ -200,16 +200,48 @@ pub fn generate_requests(
Ok(requests)
}

// TODO Use geo?
fn load_points(path: String) -> Result<Vec<(f64, f64)>> {
/// A point with an associated relative weight. Higher weights are more likely to be sampled.
#[derive(Clone, Copy)]
struct WeightedPoint {
// TODO Use geo? Could maybe just read in one big batch and have more validation
lon: f64,
lat: f64,
weight: f64,
}

// TODO Is GeomWithData simpler?
impl RTreeObject for WeightedPoint {
type Envelope = AABB<(f64, f64)>;

fn envelope(&self) -> Self::Envelope {
AABB::from_point((self.lon, self.lat))
}
}

impl PointDistance for WeightedPoint {
fn distance_2(&self, point: &(f64, f64)) -> f64 {
// Use Euclidean distance on WGS84
(self.lon - point.0).powi(2) + (self.lat - point.1).powi(2)
}
}

fn load_points(path: String) -> Result<Vec<WeightedPoint>> {
println!("Loading points from {path}");
let reader = FeatureReader::from_reader(BufReader::new(File::open(path)?));
let mut points = Vec::new();
for feature in reader.features() {
let feature = feature?;
if let Some(geometry) = feature.geometry {
if let Value::Point(pt) = geometry.value {
points.push((pt[0], pt[1]));
if let Some(ref geometry) = feature.geometry {
if let Value::Point(pt) = &geometry.value {
let weight = feature
.property("weight")
.and_then(|x| x.as_f64())
.unwrap_or(1.0);
points.push(WeightedPoint {
lon: pt[0],
lat: pt[1],
weight,
});
}
}
}
Expand Down Expand Up @@ -244,10 +276,10 @@ fn load_named_points(path: &str) -> Result<HashMap<String, (f64, f64)>> {

fn points_per_polygon(
name: &str,
points: Vec<(f64, f64)>,
polygons: &HashMap<String, MultiPolygon<f64>>,
points: Vec<WeightedPoint>,
polygons: &HashMap<String, MultiPolygon>,
use_centroids_for_empty_zones: bool,
) -> Result<HashMap<String, Vec<(f64, f64)>>> {
) -> Result<HashMap<String, Vec<WeightedPoint>>> {
let tree = RTree::bulk_load(points);

let mut empty = Vec::new();
Expand All @@ -259,7 +291,7 @@ fn points_per_polygon(
let max = bounds.max();
let envelope: AABB<(f64, f64)> = AABB::from_corners((min.x, min.y), (max.x, max.y));
for pt in tree.locate_in_envelope(&envelope) {
if polygon.contains(&geo::Point::new(pt.0, pt.1)) {
if polygon.contains(&Point::new(pt.lon, pt.lat)) {
pts_inside.push(*pt);
}
}
Expand All @@ -278,7 +310,14 @@ fn points_per_polygon(
);
for key in empty {
if let Some(centroid) = polygons[key].centroid() {
output.insert(key.clone(), vec![centroid.into()]);
output.insert(
key.clone(),
vec![WeightedPoint {
lon: centroid.x(),
lat: centroid.y(),
weight: 1.0,
}],
);
} else {
bail!("{key} had no matching {name} points, and couldn't calculate its centroid");
}
Expand All @@ -287,20 +326,12 @@ fn points_per_polygon(
Ok(output)
}

// TODO Can we use this?
/*#[derive(Deserialize)]
struct Zone {
#[serde(deserialize_with = "deserialize_geometry")]
geometry: geo_types::MultiPolygon<f64>,
name: String,
}*/

/// Extract multipolygon zones from a GeoJSON file, using the "name" property as the key in the
/// resulting map.
fn load_zones(geojson_path: &str) -> Result<HashMap<String, MultiPolygon<f64>>> {
fn load_zones(geojson_path: &str) -> Result<HashMap<String, MultiPolygon>> {
let reader = FeatureReader::from_reader(BufReader::new(File::open(geojson_path)?));

let mut zones: HashMap<String, MultiPolygon<f64>> = HashMap::new();
let mut zones: HashMap<String, MultiPolygon> = HashMap::new();
for feature in reader.features() {
let feature = feature?;
if let Some(zone_name) = feature
Expand All @@ -309,7 +340,7 @@ fn load_zones(geojson_path: &str) -> Result<HashMap<String, MultiPolygon<f64>>>
.map(|x| x.to_string())
{
let gj_geom: geojson::Geometry = feature.geometry.unwrap();
let geo_geometry: geo::Geometry<f64> = gj_geom.try_into().unwrap();
let geo_geometry: geo::Geometry = gj_geom.try_into().unwrap();
if let geo::Geometry::MultiPolygon(mp) = geo_geometry {
zones.insert(zone_name, mp);
} else if let geo::Geometry::Polygon(p) = geo_geometry {
Expand Down

0 comments on commit 343d3a6

Please sign in to comment.