geo.py 2.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990
  1. from math import asin, cos, radians, sin, sqrt
  2. # Radius of earth in meters, [as recommended by the IUGG](ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf)
  3. MEAN_EARTH_RADIUS = 6371008.8
  4. def geo_distance(lon1: float, lat1: float, lon2: float, lat2: float) -> float:
  5. """
  6. Calculate distance between two points on Earth using Haversine formula.
  7. Args:
  8. lon1: longitude of first point
  9. lat1: latitude of first point
  10. lon2: longitude of second point
  11. lat2: latitude of second point
  12. Returns:
  13. distance in meters
  14. """
  15. # convert decimal degrees to radians
  16. lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
  17. # haversine formula
  18. dlon = lon2 - lon1
  19. dlat = lat2 - lat1
  20. a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
  21. c = 2 * asin(sqrt(a))
  22. return MEAN_EARTH_RADIUS * c
  23. def test_geo_distance() -> None:
  24. moscow = {"lon": 37.6173, "lat": 55.7558}
  25. london = {"lon": -0.1278, "lat": 51.5074}
  26. berlin = {"lon": 13.4050, "lat": 52.5200}
  27. assert geo_distance(moscow["lon"], moscow["lat"], moscow["lon"], moscow["lat"]) < 1.0
  28. assert geo_distance(moscow["lon"], moscow["lat"], london["lon"], london["lat"]) > 2400 * 1000
  29. assert geo_distance(moscow["lon"], moscow["lat"], london["lon"], london["lat"]) < 2600 * 1000
  30. assert geo_distance(moscow["lon"], moscow["lat"], berlin["lon"], berlin["lat"]) > 1600 * 1000
  31. assert geo_distance(moscow["lon"], moscow["lat"], berlin["lon"], berlin["lat"]) < 1650 * 1000
  32. def boolean_point_in_polygon(
  33. point: tuple[float, float],
  34. exterior: list[tuple[float, float]],
  35. interiors: list[list[tuple[float, float]]],
  36. ) -> bool:
  37. inside_poly = False
  38. if in_ring(point, exterior, True):
  39. in_hole = False
  40. k = 0
  41. while k < len(interiors) and not in_hole:
  42. if in_ring(point, interiors[k], False):
  43. in_hole = True
  44. k += 1
  45. if not in_hole:
  46. inside_poly = True
  47. return inside_poly
  48. def in_ring(
  49. pt: tuple[float, float], ring: list[tuple[float, float]], ignore_boundary: bool
  50. ) -> bool:
  51. is_inside = False
  52. if ring[0][0] == ring[len(ring) - 1][0] and ring[0][1] == ring[len(ring) - 1][1]:
  53. ring = ring[0 : len(ring) - 1]
  54. j = len(ring) - 1
  55. for i in range(0, len(ring)):
  56. xi = ring[i][0]
  57. yi = ring[i][1]
  58. xj = ring[j][0]
  59. yj = ring[j][1]
  60. on_boundary = (
  61. (pt[1] * (xi - xj) + yi * (xj - pt[0]) + yj * (pt[0] - xi) == 0)
  62. and ((xi - pt[0]) * (xj - pt[0]) <= 0)
  63. and ((yi - pt[1]) * (yj - pt[1]) <= 0)
  64. )
  65. if on_boundary:
  66. return not ignore_boundary
  67. intersect = ((yi > pt[1]) != (yj > pt[1])) and (
  68. pt[0] < (xj - xi) * (pt[1] - yi) / (yj - yi) + xi
  69. )
  70. if intersect:
  71. is_inside = not is_inside
  72. j = i
  73. return is_inside