1use std::f64::consts::PI;
7
8use crate::storage::schema::Value;
9
10pub mod h3;
11
12const EARTH_RADIUS_KM: f64 = 6_371.0;
13const EARTH_RADIUS_M: f64 = 6_371_000.0;
14
15const WGS84_A: f64 = 6_378_137.0; const WGS84_F: f64 = 1.0 / 298.257_223_563; const WGS84_B: f64 = WGS84_A * (1.0 - WGS84_F); pub const RECOGNIZED_GEO_SHAPES: &str = "GEO_POINT, {lat, lon} object, or GeoJSON Point";
21
22#[inline]
25fn to_rad(deg: f64) -> f64 {
26 deg * PI / 180.0
27}
28
29#[inline]
30fn to_deg(rad: f64) -> f64 {
31 rad * 180.0 / PI
32}
33
34#[inline]
35pub fn micro_to_deg(micro: i32) -> f64 {
36 micro as f64 / 1_000_000.0
37}
38
39#[inline]
40pub fn deg_to_micro(deg: f64) -> i32 {
41 (deg * 1_000_000.0).round() as i32
42}
43
44pub fn recognize_geo_value(value: &Value) -> Option<(f64, f64)> {
46 match value {
47 Value::GeoPoint(lat_micro, lon_micro) => {
48 recognize_geo_degrees(micro_to_deg(*lat_micro), micro_to_deg(*lon_micro))
49 }
50 Value::Json(bytes) => recognize_geo_json(bytes),
51 _ => None,
52 }
53}
54
55const LAT_ALIASES: [&str; 2] = ["lat", "latitude"];
57const LON_ALIASES: [&str; 3] = ["lon", "lng", "longitude"];
59
60fn recognize_geo_aliases(lookup: impl Fn(&str) -> Option<f64>) -> Option<(f64, f64)> {
66 let lat = LAT_ALIASES.iter().find_map(|alias| lookup(alias))?;
67 let lon = LON_ALIASES.iter().find_map(|alias| lookup(alias))?;
68 recognize_geo_degrees(lat, lon)
69}
70
71pub fn recognize_geo_fields<'a>(field: impl Fn(&str) -> Option<&'a Value>) -> Option<(f64, f64)> {
73 recognize_geo_aliases(|alias| field(alias).and_then(numeric_value_to_f64))
74}
75
76fn recognize_geo_degrees(lat: f64, lon: f64) -> Option<(f64, f64)> {
77 if lat.is_finite()
78 && lon.is_finite()
79 && (-90.0..=90.0).contains(&lat)
80 && (-180.0..=180.0).contains(&lon)
81 {
82 Some((lat, lon))
83 } else {
84 None
85 }
86}
87
88fn numeric_value_to_f64(value: &Value) -> Option<f64> {
99 match value {
100 Value::Float(value) => Some(*value),
101 Value::Integer(value) => Some(*value as f64),
102 Value::UnsignedInteger(value) => Some(*value as f64),
103 Value::DecimalText(text) => text.parse::<f64>().ok(),
104 _ => None,
105 }
106}
107
108fn recognize_geo_json(bytes: &[u8]) -> Option<(f64, f64)> {
109 let json = crate::json::from_slice::<crate::json::Value>(bytes).ok()?;
110 let object = json.as_object()?;
111 if object.get("type").and_then(|value| value.as_str()) == Some("Point") {
112 let coordinates = object.get("coordinates")?.as_array()?;
113 if coordinates.len() != 2 {
114 return None;
115 }
116 let lon = &coordinates[0];
117 let lat = &coordinates[1];
118 return recognize_geo_degrees(json_number_to_f64(lat)?, json_number_to_f64(lon)?);
119 }
120 recognize_geo_aliases(|alias| object.get(alias).and_then(json_number_to_f64))
121}
122
123fn json_number_to_f64(value: &crate::json::Value) -> Option<f64> {
127 value.as_f64()
128}
129
130pub fn haversine_km(lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> f64 {
137 let dlat = to_rad(lat2 - lat1);
138 let dlon = to_rad(lon2 - lon1);
139 let lat1_r = to_rad(lat1);
140 let lat2_r = to_rad(lat2);
141
142 let a = (dlat / 2.0).sin().powi(2) + lat1_r.cos() * lat2_r.cos() * (dlon / 2.0).sin().powi(2);
143 EARTH_RADIUS_KM * 2.0 * a.sqrt().asin()
144}
145
146pub fn haversine_m(lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> f64 {
148 haversine_km(lat1, lon1, lat2, lon2) * 1000.0
149}
150
151pub fn vincenty_m(lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> f64 {
160 let u1 = ((1.0 - WGS84_F) * to_rad(lat1).tan()).atan();
161 let u2 = ((1.0 - WGS84_F) * to_rad(lat2).tan()).atan();
162 let l = to_rad(lon2 - lon1);
163
164 let sin_u1 = u1.sin();
165 let cos_u1 = u1.cos();
166 let sin_u2 = u2.sin();
167 let cos_u2 = u2.cos();
168
169 let mut lambda = l;
170 let mut prev_lambda;
171 let mut sin_sigma;
172 let mut cos_sigma;
173 let mut sigma;
174 let mut sin_alpha;
175 let mut cos2_alpha;
176 let mut cos_2sigma_m;
177
178 for _ in 0..100 {
179 let sin_lambda = lambda.sin();
180 let cos_lambda = lambda.cos();
181
182 sin_sigma = ((cos_u2 * sin_lambda).powi(2)
183 + (cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_lambda).powi(2))
184 .sqrt();
185
186 if sin_sigma == 0.0 {
187 return 0.0; }
189
190 cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos_lambda;
191 sigma = sin_sigma.atan2(cos_sigma);
192 sin_alpha = cos_u1 * cos_u2 * sin_lambda / sin_sigma;
193 cos2_alpha = 1.0 - sin_alpha.powi(2);
194
195 cos_2sigma_m = if cos2_alpha != 0.0 {
196 cos_sigma - 2.0 * sin_u1 * sin_u2 / cos2_alpha
197 } else {
198 0.0
199 };
200
201 let c = WGS84_F / 16.0 * cos2_alpha * (4.0 + WGS84_F * (4.0 - 3.0 * cos2_alpha));
202
203 prev_lambda = lambda;
204 lambda = l
205 + (1.0 - c)
206 * WGS84_F
207 * sin_alpha
208 * (sigma
209 + c * sin_sigma
210 * (cos_2sigma_m + c * cos_sigma * (-1.0 + 2.0 * cos_2sigma_m.powi(2))));
211
212 if (lambda - prev_lambda).abs() < 1e-12 {
213 let u_sq = cos2_alpha * (WGS84_A.powi(2) - WGS84_B.powi(2)) / WGS84_B.powi(2);
215 let a_coeff =
216 1.0 + u_sq / 16384.0 * (4096.0 + u_sq * (-768.0 + u_sq * (320.0 - 175.0 * u_sq)));
217 let b_coeff = u_sq / 1024.0 * (256.0 + u_sq * (-128.0 + u_sq * (74.0 - 47.0 * u_sq)));
218
219 let delta_sigma = b_coeff
220 * sin_sigma
221 * (cos_2sigma_m
222 + b_coeff / 4.0
223 * (cos_sigma * (-1.0 + 2.0 * cos_2sigma_m.powi(2))
224 - b_coeff / 6.0
225 * cos_2sigma_m
226 * (-3.0 + 4.0 * sin_sigma.powi(2))
227 * (-3.0 + 4.0 * cos_2sigma_m.powi(2))));
228
229 return WGS84_B * a_coeff * (sigma - delta_sigma);
230 }
231 }
232
233 haversine_m(lat1, lon1, lat2, lon2)
235}
236
237pub fn vincenty_km(lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> f64 {
239 vincenty_m(lat1, lon1, lat2, lon2) / 1000.0
240}
241
242pub fn bearing(lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> f64 {
249 let lat1_r = to_rad(lat1);
250 let lat2_r = to_rad(lat2);
251 let dlon = to_rad(lon2 - lon1);
252
253 let y = dlon.sin() * lat2_r.cos();
254 let x = lat1_r.cos() * lat2_r.sin() - lat1_r.sin() * lat2_r.cos() * dlon.cos();
255
256 (to_deg(y.atan2(x)) + 360.0) % 360.0
257}
258
259pub fn final_bearing(lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> f64 {
261 (bearing(lat2, lon2, lat1, lon1) + 180.0) % 360.0
262}
263
264pub fn midpoint(lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> (f64, f64) {
269 let lat1_r = to_rad(lat1);
270 let lat2_r = to_rad(lat2);
271 let dlon = to_rad(lon2 - lon1);
272
273 let bx = lat2_r.cos() * dlon.cos();
274 let by = lat2_r.cos() * dlon.sin();
275
276 let lat =
277 (lat1_r.sin() + lat2_r.sin()).atan2(((lat1_r.cos() + bx).powi(2) + by.powi(2)).sqrt());
278 let lon = to_rad(lon1) + by.atan2(lat1_r.cos() + bx);
279
280 (to_deg(lat), to_deg(lon))
281}
282
283pub fn destination(lat: f64, lon: f64, bearing_deg: f64, distance_km: f64) -> (f64, f64) {
288 let lat_r = to_rad(lat);
289 let brng_r = to_rad(bearing_deg);
290 let d = distance_km / EARTH_RADIUS_KM;
291
292 let lat2 = (lat_r.sin() * d.cos() + lat_r.cos() * d.sin() * brng_r.cos()).asin();
293 let lon2 = to_rad(lon)
294 + (brng_r.sin() * d.sin() * lat_r.cos()).atan2(d.cos() - lat_r.sin() * lat2.sin());
295
296 (to_deg(lat2), to_deg(lon2))
297}
298
299pub fn bounding_box(lat: f64, lon: f64, radius_km: f64) -> (f64, f64, f64, f64) {
307 let lat_delta = radius_km / 111.32;
308 let lon_delta = radius_km / (111.32 * to_rad(lat).cos().max(0.0001));
309
310 let min_lat = (lat - lat_delta).max(-90.0);
311 let max_lat = (lat + lat_delta).min(90.0);
312 let min_lon = lon - lon_delta;
313 let max_lon = lon + lon_delta;
314
315 (min_lat, min_lon, max_lat, max_lon)
316}
317
318pub fn polygon_area_km2(vertices: &[(f64, f64)]) -> f64 {
323 let n = vertices.len();
324 if n < 3 {
325 return 0.0;
326 }
327
328 let mut total = 0.0f64;
329 for i in 0..n {
330 let (lat1, lon1) = vertices[i];
331 let (lat2, lon2) = vertices[(i + 1) % n];
332 total += to_rad(lon2 - lon1) * (2.0 + to_rad(lat1).sin() + to_rad(lat2).sin());
333 }
334
335 (total.abs() / 2.0) * EARTH_RADIUS_KM * EARTH_RADIUS_KM
336}
337
338pub fn point_in_polygon_even_odd(lat: f64, lon: f64, vertices: &[(f64, f64)]) -> bool {
344 if vertices.len() < 3 {
345 return false;
346 }
347 let mut inside = false;
348 let mut j = vertices.len() - 1;
349 for i in 0..vertices.len() {
350 let (yi, xi) = vertices[i];
351 let (yj, xj) = vertices[j];
352 if point_on_segment(lat, lon, yi, xi, yj, xj) {
353 return true;
354 }
355 if ((yi > lat) != (yj > lat)) && (lon < (xj - xi) * (lat - yi) / (yj - yi) + xi) {
356 inside = !inside;
357 }
358 j = i;
359 }
360 inside
361}
362
363fn point_on_segment(py: f64, px: f64, ay: f64, ax: f64, by: f64, bx: f64) -> bool {
364 const EPS: f64 = 1e-12;
365 let cross = (px - ax) * (by - ay) - (py - ay) * (bx - ax);
366 if cross.abs() > EPS {
367 return false;
368 }
369 px >= ax.min(bx) - EPS
370 && px <= ax.max(bx) + EPS
371 && py >= ay.min(by) - EPS
372 && py <= ay.max(by) + EPS
373}
374
375#[cfg(test)]
378mod tests {
379 use super::*;
380 use proptest::prelude::*;
381 use std::collections::HashMap;
382
383 fn json_value(json: &str) -> Value {
384 Value::Json(json.as_bytes().to_vec())
385 }
386
387 fn fields(values: &[(&str, Value)]) -> HashMap<String, Value> {
388 values
389 .iter()
390 .map(|(key, value)| ((*key).to_string(), value.clone()))
391 .collect()
392 }
393
394 #[test]
395 fn exact_decimal_text_coordinates_are_recognized() {
396 let decimal = fields(&[
399 ("lat", Value::DecimalText("38.76000000".to_string())),
400 ("lon", Value::DecimalText("-77.15000000".to_string())),
401 ]);
402 assert_eq!(
403 recognize_geo_fields(|key| decimal.get(key)),
404 Some((38.76, -77.15))
405 );
406 }
407
408 #[test]
409 fn both_halves_of_the_seam_cover_the_same_numeric_shapes() {
410 for (text, expected) in [("38.76000000", 38.76), ("-77.15000000", -77.15)] {
415 assert_eq!(
416 numeric_value_to_f64(&Value::DecimalText(text.to_string())),
417 Some(expected),
418 "storage side must read exact decimal {text}"
419 );
420 assert_eq!(
421 json_number_to_f64(&crate::json::Value::Decimal(text.to_string())),
422 Some(expected),
423 "json side must read exact decimal {text}"
424 );
425 }
426 }
427
428 #[test]
429 fn test_haversine_paris_london() {
430 let d = haversine_km(48.8566, 2.3522, 51.5074, -0.1278);
431 assert!((d - 344.0).abs() < 5.0, "Paris-London: {d} km");
432 }
433
434 #[test]
435 fn test_haversine_zero_distance() {
436 let d = haversine_km(0.0, 0.0, 0.0, 0.0);
437 assert!(d.abs() < 0.001, "same point: {d} km");
438 }
439
440 #[test]
441 fn test_haversine_antipodal() {
442 let d = haversine_km(0.0, 0.0, 0.0, 180.0);
443 assert!((d - 20015.0).abs() < 100.0, "antipodal: {d} km");
444 }
445
446 #[test]
447 fn test_vincenty_paris_london() {
448 let d = vincenty_km(48.8566, 2.3522, 51.5074, -0.1278);
449 assert!((d - 343.5).abs() < 2.0, "Vincenty Paris-London: {d} km");
450 }
451
452 #[test]
453 fn test_vincenty_coincident() {
454 let d = vincenty_m(48.8566, 2.3522, 48.8566, 2.3522);
455 assert!(d.abs() < 0.001, "coincident: {d} m");
456 }
457
458 #[test]
459 fn test_vincenty_new_york_tokyo() {
460 let d = vincenty_km(40.7128, -74.0060, 35.6762, 139.6503);
461 assert!((d - 10838.0).abs() < 50.0, "NY-Tokyo: {d} km");
462 }
463
464 #[test]
465 fn test_bearing_north() {
466 let b = bearing(0.0, 0.0, 1.0, 0.0);
467 assert!((b - 0.0).abs() < 1.0, "north bearing: {b}°");
468 }
469
470 #[test]
471 fn test_bearing_east() {
472 let b = bearing(0.0, 0.0, 0.0, 1.0);
473 assert!((b - 90.0).abs() < 1.0, "east bearing: {b}°");
474 }
475
476 #[test]
477 fn test_midpoint_equator() {
478 let (lat, lon) = midpoint(0.0, 0.0, 0.0, 10.0);
479 assert!((lat - 0.0).abs() < 0.01, "midpoint lat: {lat}");
480 assert!((lon - 5.0).abs() < 0.01, "midpoint lon: {lon}");
481 }
482
483 #[test]
484 fn test_destination() {
485 let (lat, lon) = destination(0.0, 0.0, 0.0, 111.32);
486 assert!((lat - 1.0).abs() < 0.1, "destination lat: {lat}");
487 assert!(lon.abs() < 0.1, "destination lon: {lon}");
488 }
489
490 #[test]
491 fn test_bounding_box() {
492 let (min_lat, min_lon, max_lat, max_lon) = bounding_box(0.0, 0.0, 111.32);
493 assert!((min_lat - (-1.0)).abs() < 0.1);
494 assert!((max_lat - 1.0).abs() < 0.1);
495 assert!(min_lon < 0.0);
496 assert!(max_lon > 0.0);
497 }
498
499 #[test]
500 fn test_micro_conversion() {
501 let lat = -23.550520;
502 let micro = deg_to_micro(lat);
503 let back = micro_to_deg(micro);
504 assert!((lat - back).abs() < 0.000001);
505 }
506
507 #[test]
508 fn recognize_geo_value_accepts_geopoint_json_aliases_and_geojson_point() {
509 assert_eq!(
510 recognize_geo_value(&Value::GeoPoint(38_760_000, -77_150_000)),
511 Some((38.76, -77.15))
512 );
513 assert_eq!(
514 recognize_geo_value(&json_value(r#"{"lat":38.76,"lon":-77.15}"#)),
515 Some((38.76, -77.15))
516 );
517 assert_eq!(
518 recognize_geo_value(&json_value(r#"{"latitude":38,"lng":-77}"#)),
519 Some((38.0, -77.0))
520 );
521 assert_eq!(
522 recognize_geo_value(&json_value(r#"{"latitude":38.76,"longitude":-77.15}"#)),
523 Some((38.76, -77.15))
524 );
525 assert_eq!(
526 recognize_geo_value(&json_value(
527 r#"{"type":"Point","coordinates":[-77.15,38.76]}"#
528 )),
529 Some((38.76, -77.15))
530 );
531 }
532
533 #[test]
534 fn recognize_geojson_point_uses_lon_lat_order() {
535 assert_eq!(
536 recognize_geo_value(&json_value(
537 r#"{"type":"Point","coordinates":[103.8198,1.3521]}"#
538 )),
539 Some((1.3521, 103.8198)),
540 "GeoJSON coordinates are [longitude, latitude], not [latitude, longitude]"
541 );
542 }
543
544 #[test]
545 fn recognize_geo_fields_accepts_numeric_aliases() {
546 let row = fields(&[
547 ("latitude", Value::Integer(38)),
548 ("longitude", Value::Float(-77.15)),
549 ]);
550 assert_eq!(
551 recognize_geo_fields(|key| row.get(key)),
552 Some((38.0, -77.15))
553 );
554
555 let node = fields(&[("lat", Value::Float(38.76)), ("lng", Value::Integer(-77))]);
556 assert_eq!(
557 recognize_geo_fields(|key| node.get(key)),
558 Some((38.76, -77.0))
559 );
560 }
561
562 #[test]
563 fn recognize_geo_rejects_non_geo_shapes() {
564 for value in [
565 json_value(r#"{"lat":"38.76","lon":"-77.15"}"#),
566 json_value(r#"{"type":"Polygon","coordinates":[[-77.15,38.76]]}"#),
567 json_value(r#"{"type":"LineString","coordinates":[[-77.15,38.76]]}"#),
568 json_value(r#"{"type":"Point"}"#),
569 json_value(r#"{"type":"Point","coordinates":[]}"#),
570 json_value(r#"{"type":"Point","coordinates":[-77.15]}"#),
571 json_value(r#"{"type":"Point","coordinates":[-77.15,38.76,0]}"#),
572 json_value(r#"{"type":"Point","coordinates":["-77.15",38.76]}"#),
573 json_value(r#"{"type":"Point","coordinates":[-77.15,"38.76"]}"#),
574 json_value(r#"{"type":"Point","coordinates":[181.0,38.76]}"#),
575 json_value(r#"{"type":"Point","coordinates":[-77.15,91.0]}"#),
576 json_value(r#"{"lat":38.76}"#),
577 json_value(r#"{"lat":91.0,"lon":0.0}"#),
578 json_value(r#"{"lat":0.0,"lon":181.0}"#),
579 json_value(r#"{"lat":null,"lon":0.0}"#),
580 json_value(r#"not json"#),
581 Value::text("38.76,-77.15".to_string()),
582 ] {
583 assert_eq!(recognize_geo_value(&value), None, "{value:?}");
584 }
585
586 let string_fields = fields(&[
587 ("lat", Value::text("38.76".to_string())),
588 ("lon", Value::Float(-77.15)),
589 ]);
590 assert_eq!(recognize_geo_fields(|key| string_fields.get(key)), None);
591
592 let missing_fields = fields(&[("lat", Value::Float(38.76))]);
593 assert_eq!(recognize_geo_fields(|key| missing_fields.get(key)), None);
594
595 let non_finite_fields = fields(&[
596 ("lat", Value::Float(f64::NAN)),
597 ("lon", Value::Float(-77.15)),
598 ]);
599 assert_eq!(recognize_geo_fields(|key| non_finite_fields.get(key)), None);
600
601 let out_of_range_fields =
602 fields(&[("lat", Value::Float(38.76)), ("lon", Value::Float(-181.0))]);
603 assert_eq!(
604 recognize_geo_fields(|key| out_of_range_fields.get(key)),
605 None
606 );
607 }
608
609 proptest! {
610 #[test]
611 fn recognize_geo_json_and_field_maps_do_not_drift(
612 lat in -90.0f64..=90.0,
613 lon in -180.0f64..=180.0,
614 ) {
615 prop_assume!(lat.is_finite());
616 prop_assume!(lon.is_finite());
617 let json = json_value(&format!(r#"{{"lat":{lat},"lon":{lon}}}"#));
618 let geojson = json_value(&format!(r#"{{"type":"Point","coordinates":[{lon},{lat}]}}"#));
619 let fields = fields(&[("lat", Value::Float(lat)), ("lon", Value::Float(lon))]);
620 prop_assert_eq!(
621 recognize_geo_value(&json),
622 recognize_geo_fields(|key| fields.get(key))
623 );
624 prop_assert_eq!(
625 recognize_geo_value(&geojson),
626 recognize_geo_fields(|key| fields.get(key))
627 );
628 }
629
630 #[test]
631 fn recognize_geo_rejects_generated_out_of_range_values(
632 lat in prop_oneof![-1000.0f64..-90.000_001, 90.000_001f64..1000.0],
633 lon in -180.0f64..=180.0,
634 ) {
635 let json = json_value(&format!(r#"{{"lat":{lat},"lon":{lon}}}"#));
636 let geojson = json_value(&format!(r#"{{"type":"Point","coordinates":[{lon},{lat}]}}"#));
637 let fields = fields(&[("lat", Value::Float(lat)), ("lon", Value::Float(lon))]);
638 prop_assert_eq!(recognize_geo_value(&json), None);
639 prop_assert_eq!(recognize_geo_value(&geojson), None);
640 prop_assert_eq!(recognize_geo_fields(|key| fields.get(key)), None);
641 }
642 }
643}