point-on-surface, XML file outline calculations

Cargo.toml

@@ -10,6 +10,7 @@ anyhow = "1.0.98"
 clap = { version = "4.5.37", features = ["derive"] }
 crossbeam-channel = "0.5.15"
 flatgeobuf = "4.6.0"
+geo = "0.30.0"
 geo-types = "0.7.16"
 geozero = "0.14.0"
 indicatif = "0.17.11"

src/error.rs

@@ -22,6 +22,8 @@ pub enum Error {
     Projection(#[from] proj4rs::errors::Error),
     #[error("IO error: {0}")]
     FS(#[from] std::io::Error),
+    #[error("Parsing integer failed: {0}")]
+    ParseInt(#[from] std::num::ParseIntError),
 }
 
 pub type Result<T> = std::result::Result<T, Error>;

src/geo.rs

@@ -0,0 +1,28 @@
+use geo_types::{MultiPolygon, Point};
+
+pub fn point_on_surface(mp: &MultiPolygon<f64>) -> Point<f64> {
+    use geo::{
+        Triangle,
+        algorithm::{Area, Centroid, TriangulateEarcut},
+    };
+
+    // get the biggest polygon
+    let polygon = mp
+        .into_iter()
+        .max_by(|a, b| a.unsigned_area().partial_cmp(&b.unsigned_area()).unwrap())
+        .expect("MultiPolygon must have at least one Polygon");
+
+    // (1) Triangulate into a Vec<Triangle<f64>>
+    let triangles: Vec<Triangle<f64>> = polygon.earcut_triangles();
+
+    // (2) Pick the triangle with the max area
+    let largest = triangles
+        .into_iter()
+        .max_by(|a, b| a.unsigned_area().partial_cmp(&b.unsigned_area()).unwrap())
+        .expect("polygon must have at least one triangle");
+
+    // (3) Its centroid is interior
+
+
+    largest.centroid()
+}