Blizzard
91 lines
2.9 KiB
Python
91 lines
2.9 KiB
Python
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
|