1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
"""
Functions for computing navigational information. Can be used to add
navigational information to DataFrames.
"""
from math import acos, asin, atan2, cos, sin, degrees, radians, sqrt
from typing import Tuple
def gcd_slc(
lon0: float,
lat0: float,
lon1: float,
lat1: float,
) -> float:
"""
Compute great circle distance on earth surface between two locations.
Parameters
----------
lon0 : float
Longitude of position 0
lat0 : float
Latitude of position 0
lon1 : float
Longitude of position 1
lat1 : float
Latitude of position 1
Returns
-------
dist : float
Great circle distance between position 0 and position 1.
"""
if abs(lat0 - lat1) <= 1e-6 and abs(lon0 - lon1) <= 1e-6:
return 0
r_earth = 6371
# Convert to radians
lat0, lat1, lon0, lon1 = map(radians, [lat0, lat1, lon0, lon1])
return r_earth * acos(
sin(lat0) * sin(lat1) + cos(lat0) * cos(lat1) * cos(lon1 - lon0)
)
def haversine(
lon0: float,
lat0: float,
lon1: float,
lat1: float,
) -> float:
"""
Compute Haversine distance between two points.
Parameters
----------
lon0 : float
Longitude of position 0
lat0 : float
Latitude of position 0
lon1 : float
Longitude of position 1
lat1 : float
Latitude of position 1
Returns
-------
dist : float
Haversine distance between position 0 and position 1.
"""
lat0, lat1, dlon, dlat = map(
radians, [lat0, lat1, lon1 - lon0, lat1 - lat0]
)
if abs(dlon) < 1e-6 and abs(dlat) < 1e-6:
return 0
r_earth = 6371
a = sin(dlat / 2) ** 2 + cos(lat0) * cos(lat1) * sin(dlon / 2) ** 2
c = 2 * asin(sqrt(a))
return c * r_earth
def bearing(
lon0: float,
lat0: float,
lon1: float,
lat1: float,
) -> float:
"""
Compute the bearing of a track from (lon0, lat0) to (lon1, lat1).
Duplicated from geo-py
Parameters
----------
lon0 : float,
Longitude of start point
lat0 : float,
Latitude of start point
lon1 : float,
Longitude of target point
lat1 : float,
Latitude of target point
Returns
-------
bearing : float
The bearing from point (lon0, lat0) to point (lon1, lat1) in degrees.
"""
lon0, lat0, lon1, lat1 = map(radians, [lon0, lat0, lon1, lat1])
dlon = lon1 - lon0
numerator = sin(dlon) * cos(lat1)
denominator = cos(lat0) * sin(lat1) - (sin(lat0) * cos(lat1) * cos(dlon))
theta = atan2(numerator, denominator)
theta_deg = (degrees(theta) + 360) % 360
return theta_deg
def destination(
lon: float, lat: float, bearing: float, distance: float
) -> Tuple[float, float]:
"""
Compute destination of a great circle path.
Compute the destination of a track started from 'lon', 'lat', with
'bearing'. Distance is in units of km.
Duplicated from geo-py
Parameters
----------
lon : float
Longitude of initial position
lat : float
Latitude of initial position
bearing : float
Direction of track
distance : float
Distance to travel
Returns
-------
destination : tuple[float, float]
Longitude and Latitude of final position
"""
lon, lat = radians(lon), radians(lat)
radians_bearing = radians(bearing)
r_earth = 6371
delta = distance / r_earth
lat2 = asin(
sin(lat) * cos(delta) + cos(lat) * sin(delta) * cos(radians_bearing)
)
numerator = sin(radians_bearing) * sin(delta) * cos(lat)
denominator = cos(delta) - sin(lat) * sin(lat2)
lon2 = lon + atan2(numerator, denominator)
lon2_deg = (degrees(lon2) + 540) % 360 - 180
lat2_deg = degrees(lat2)
return lon2_deg, lat2_deg
def midpoint(
lon0: float,
lat0: float,
lon1: float,
lat1: float,
) -> Tuple[float, float]:
"""
Compute the midpoint of a great circle track
Parameters
----------
lon0 : float
Longitude of position 0
lat0 : float
Latitude of position 0
lon1 : float
Longitude of position 1
lat1 : float
Latitude of position 1
Returns
-------
lon, lat
Positions of midpoint between position 0 and position 1
"""
bear = bearing(lon0, lat0, lon1, lat1)
dist = haversine(lon0, lat0, lon1, lat1)
return destination(lon0, lat0, bear, dist / 2)