| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990 |
- from math import asin, cos, radians, sin, sqrt
- # Radius of earth in meters, [as recommended by the IUGG](ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf)
- MEAN_EARTH_RADIUS = 6371008.8
- def geo_distance(lon1: float, lat1: float, lon2: float, lat2: float) -> float:
- """
- Calculate distance between two points on Earth using Haversine formula.
- Args:
- lon1: longitude of first point
- lat1: latitude of first point
- lon2: longitude of second point
- lat2: latitude of second point
- Returns:
- distance in meters
- """
- # convert decimal degrees to radians
- lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
- # haversine formula
- dlon = lon2 - lon1
- dlat = lat2 - lat1
- a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
- c = 2 * asin(sqrt(a))
- return MEAN_EARTH_RADIUS * c
- def test_geo_distance() -> None:
- moscow = {"lon": 37.6173, "lat": 55.7558}
- london = {"lon": -0.1278, "lat": 51.5074}
- berlin = {"lon": 13.4050, "lat": 52.5200}
- assert geo_distance(moscow["lon"], moscow["lat"], moscow["lon"], moscow["lat"]) < 1.0
- assert geo_distance(moscow["lon"], moscow["lat"], london["lon"], london["lat"]) > 2400 * 1000
- assert geo_distance(moscow["lon"], moscow["lat"], london["lon"], london["lat"]) < 2600 * 1000
- assert geo_distance(moscow["lon"], moscow["lat"], berlin["lon"], berlin["lat"]) > 1600 * 1000
- assert geo_distance(moscow["lon"], moscow["lat"], berlin["lon"], berlin["lat"]) < 1650 * 1000
- def boolean_point_in_polygon(
- point: tuple[float, float],
- exterior: list[tuple[float, float]],
- interiors: list[list[tuple[float, float]]],
- ) -> bool:
- inside_poly = False
- if in_ring(point, exterior, True):
- in_hole = False
- k = 0
- while k < len(interiors) and not in_hole:
- if in_ring(point, interiors[k], False):
- in_hole = True
- k += 1
- if not in_hole:
- inside_poly = True
- return inside_poly
- def in_ring(
- pt: tuple[float, float], ring: list[tuple[float, float]], ignore_boundary: bool
- ) -> bool:
- is_inside = False
- if ring[0][0] == ring[len(ring) - 1][0] and ring[0][1] == ring[len(ring) - 1][1]:
- ring = ring[0 : len(ring) - 1]
- j = len(ring) - 1
- for i in range(0, len(ring)):
- xi = ring[i][0]
- yi = ring[i][1]
- xj = ring[j][0]
- yj = ring[j][1]
- on_boundary = (
- (pt[1] * (xi - xj) + yi * (xj - pt[0]) + yj * (pt[0] - xi) == 0)
- and ((xi - pt[0]) * (xj - pt[0]) <= 0)
- and ((yi - pt[1]) * (yj - pt[1]) <= 0)
- )
- if on_boundary:
- return not ignore_boundary
- intersect = ((yi > pt[1]) != (yj > pt[1])) and (
- pt[0] < (xj - xi) * (pt[1] - yi) / (yj - yi) + xi
- )
- if intersect:
- is_inside = not is_inside
- j = i
- return is_inside
|