voronoi.rs
edited
1use flatgeobuf::{FallibleStreamingIterator, FgbReader, FgbWriter, GeometryType};
2use geo::{GeometryCollection, MultiPoint, Voronoi, VoronoiClip, VoronoiParams};
3use geozero::ToGeo;
4use std::fs::File;
5use std::io::{BufReader, BufWriter};
6use std::time::Instant;
7
8fn main() -> Result<(), Box<dyn std::error::Error>> {
9 // Read input FlatGeobuf
10 // Input is EPSG:27700 (British National Grid)
11 let input_path = "data/uk_postcodes_bng.fgb";
12 let rstart = Instant::now();
13 let file = BufReader::new(File::open(input_path)?);
14 let mut reader = FgbReader::open(file)?.select_all()?;
15 let mut geoms_vec: Vec<geo::Geometry> = Vec::new();
16 while let Some(feature) = reader.next()? {
17 if let Ok(geom) = feature.to_geo() {
18 geoms_vec.push(geom);
19 }
20 }
21 let geoms = GeometryCollection::from(geoms_vec);
22 println!(
23 "Loaded {} points from {}. Took {:?}",
24 geoms.0.len(),
25 input_path,
26 rstart.elapsed()
27 );
28
29 let points: MultiPoint = geoms
30 .into_iter()
31 .map(|geom| match geom {
32 geo::Geometry::Point(point) => point,
33 _ => panic!("Expected Point geometry"),
34 })
35 .collect();
36
37 let vstart = Instant::now();
38 let voronoi_cells =
39 points.voronoi_cells_with_params(VoronoiParams::new().clip(VoronoiClip::Envelope))?;
40 println!("Voronoi calculation took {:?}", vstart.elapsed());
41
42 let geometry_collection = GeometryCollection::from_iter(voronoi_cells);
43
44 // Write output FlatGeobuf (in EPSG:27700)
45 let wstart = Instant::now();
46 let output_path = "data/voronoi_output.fgb";
47 let mut fgb = FgbWriter::create("voronoi", GeometryType::Polygon)?;
48 for polygon in geometry_collection.0 {
49 fgb.add_feature_geom(polygon, |_feat| {})?;
50 }
51 let mut fout = BufWriter::new(File::create(output_path)?);
52 fgb.write(&mut fout)?;
53 println!(
54 "Wrote Voronoi diagram to {}. Took {:?}",
55 output_path,
56 wstart.elapsed()
57 );
58
59 Ok(())
60}