51 lines
1.4 KiB
Python
51 lines
1.4 KiB
Python
from __future__ import annotations
|
|
|
|
from dataclasses import dataclass
|
|
from math import asin, cos, radians, sin, sqrt
|
|
|
|
from pyproj import Geod
|
|
|
|
WGS84_GEOD = Geod(ellps="WGS84")
|
|
EARTH_RADIUS_M = 6_371_000.0
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class GeoPoint:
|
|
lat: float
|
|
lon: float
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class PathPoint(GeoPoint):
|
|
distance_m: float
|
|
|
|
|
|
def haversine(a: GeoPoint, b: GeoPoint) -> float:
|
|
lat1 = radians(a.lat)
|
|
lat2 = radians(b.lat)
|
|
dlat = radians(b.lat - a.lat)
|
|
dlon = radians(b.lon - a.lon)
|
|
h = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
|
|
return 2 * EARTH_RADIUS_M * asin(sqrt(h))
|
|
|
|
|
|
def sample_path(start: GeoPoint, end: GeoPoint, n: int) -> list[PathPoint]:
|
|
if n < 2:
|
|
raise ValueError("n must be at least 2")
|
|
|
|
_, _, total_distance = WGS84_GEOD.inv(start.lon, start.lat, end.lon, end.lat)
|
|
inner = [] if n == 2 else WGS84_GEOD.npts(start.lon, start.lat, end.lon, end.lat, n - 2)
|
|
coords = [(start.lon, start.lat), *inner, (end.lon, end.lat)]
|
|
step = total_distance / (n - 1)
|
|
return [
|
|
PathPoint(lat=lat, lon=lon, distance_m=total_distance if i == n - 1 else i * step)
|
|
for i, (lon, lat) in enumerate(coords)
|
|
]
|
|
|
|
|
|
def linestring_geojson(points: list[PathPoint]) -> dict[str, object]:
|
|
return {
|
|
"type": "LineString",
|
|
"coordinates": [[point.lon, point.lat] for point in points],
|
|
}
|