umsh_regiondb/sampling.rs
1//! The sampled-dilation membership rule.
2//!
3//! Expanded coverage is not stored as geometry. A region's effective
4//! membership is *defined* as: a position belongs to a region if any point of
5//! a fixed sample pattern — the position itself, six points at half the
6//! region's expansion distance, and twelve at the full distance — lands inside
7//! the region's core. The pattern is the semantics, not an approximation of
8//! something else, which is what lets every implementation agree exactly: each
9//! computes the same nineteen spherical destinations and runs the same integer
10//! point-in-polygon test.
11//!
12//! This mirrors `tools/regiondb-build/src/regiondb_build/sampling.py`
13//! operation for operation.
14
15/// Mean Earth radius, matching the suggested-default tie-break.
16pub const EARTH_RADIUS_M: f64 = 6_371_008.8;
17
18/// Bearings of the full-distance ring, degrees clockwise from north.
19pub const FULL_RING_BEARINGS: [f64; 12] = [
20 0.0, 30.0, 60.0, 90.0, 120.0, 150.0, 180.0, 210.0, 240.0, 270.0, 300.0, 330.0,
21];
22
23/// Bearings of the half-distance ring.
24pub const HALF_RING_BEARINGS: [f64; 6] = [0.0, 60.0, 120.0, 180.0, 240.0, 300.0];
25
26/// Spherical direct problem: where `distance_m` at `bearing_deg` lands.
27pub fn destination(latitude: f64, longitude: f64, bearing_deg: f64, distance_m: f64) -> (f64, f64) {
28 let angular = distance_m / EARTH_RADIUS_M;
29 let bearing = bearing_deg.to_radians();
30 let lat1 = latitude.to_radians();
31 let sin_lat2 = lat1.sin() * angular.cos() + lat1.cos() * angular.sin() * bearing.cos();
32 let lat2 = sin_lat2.clamp(-1.0, 1.0).asin();
33 let lon2 = longitude.to_radians()
34 + (bearing.sin() * angular.sin() * lat1.cos()).atan2(angular.cos() - lat1.sin() * sin_lat2);
35 (lat2.to_degrees(), lon2.to_degrees())
36}
37
38/// Every position to test against core geometry, the position itself first.
39///
40/// The position leads so a caller can stop on a core hit, and so the
41/// core-versus-expanded distinction falls out of the same loop: a hit at
42/// index zero is core, a hit anywhere later is expanded.
43pub fn sample_positions(latitude: f64, longitude: f64, expansion_m: i64) -> Vec<(f64, f64)> {
44 let mut positions = Vec::with_capacity(19);
45 positions.push((latitude, longitude));
46 if expansion_m <= 0 {
47 return positions;
48 }
49 let full = expansion_m as f64;
50 for bearing in HALF_RING_BEARINGS {
51 positions.push(destination(latitude, longitude, bearing, full / 2.0));
52 }
53 for bearing in FULL_RING_BEARINGS {
54 positions.push(destination(latitude, longitude, bearing, full));
55 }
56 positions
57}