From c8ba5f17788916b0571ccedfefea4806f726368c Mon Sep 17 00:00:00 2001
From: Gaute Hope <eg@gaute.vetsj.com>
Date: Sun, 7 Jul 2024 08:50:17 +0200
Subject: [PATCH] try prepared geometries

---
 Cargo.toml    |  3 +++
 src/shapes.rs | 11 +++++++----
 2 files changed, 10 insertions(+), 4 deletions(-)

diff --git a/Cargo.toml b/Cargo.toml
index 829648f..43b7c05 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -28,6 +28,7 @@ xz2 = "0.1"
 ndarray = { version = "0.15", features = [ "rayon" ] }
 wkb = "0.7.1"
 rstar = "0.12.0"
+static_cell = "2.1.0"
 
 [dev-dependencies]
 rayon = "1"
@@ -46,3 +47,5 @@ default = [ ]
 [profile.release]
 debug = true
 
+[patch.crates-io]
+geo = { git = "https://github.com/georust/geo", branch = "mkirk/prepared-geom" }
diff --git a/src/shapes.rs b/src/shapes.rs
index 013049e..3a5f83d 100644
--- a/src/shapes.rs
+++ b/src/shapes.rs
@@ -4,11 +4,13 @@ use std::io;
 use std::path::Path;
 use std::{borrow::Borrow, convert::TryInto};
 
-use geo::{point, Contains, Geometry, MultiPolygon, Point, Polygon};
+use geo::{point, Contains, Geometry, MultiPolygon, Point, Polygon, PreparedGeometry, Relate};
 use numpy::{PyArray, PyReadonlyArrayDyn};
 use rstar::{PointDistance, RTree, RTreeObject, AABB};
+use static_cell::StaticCell;
 
 pub static GSHHS_F: &str = "gshhs_f_-180.000000E-90.000000N180.000000E90.000000N.wkb.xz";
+static POLYS: StaticCell<Vec<Polygon>> = StaticCell::new();
 
 #[pyclass]
 #[derive(Clone)]
@@ -17,7 +19,7 @@ pub struct Gshhg {
 }
 
 #[derive(Clone)]
-struct PolW(Polygon);
+struct PolW(PreparedGeometry<'static>);
 
 impl RTreeObject for PolW {
     type Envelope = AABB<Point<f64>>;
@@ -33,7 +35,7 @@ impl PointDistance for PolW {
     }
 
     fn contains_point(&self, point: &Point<f64>) -> bool {
-        self.0.contains(point)
+        self.0.relate(point).is_contains()
     }
 
     fn distance_2_if_less_or_equal(&self, _point: &Point<f64>, _max_distance: f64) -> Option<f64> {
@@ -45,7 +47,8 @@ impl Gshhg {
     pub fn from_geom(geom: Geometry) -> io::Result<Gshhg> {
         let geom: MultiPolygon = geom.try_into().unwrap();
         assert!(geom.0.len() > 10);
-        let geoms = geom.0.into_iter().map(|p| PolW(p)).collect();
+        POLYS.init(geom.into_iter().collect());
+        let geoms = POLYS.into_iter().map(|p| PolW(PreparedGeometry::from(p))).collect();
 
         let geom = RTree::bulk_load(geoms);
         Ok(Gshhg { geom })