Skip to content

helpers

bearing_to_direction 

bearing_to_direction(bearing: ndarray) -> ndarray

Convert bearing angle to cardinal direction name Args: bearing (np.ndarray): Bearing angle in decimal degrees

Returns:

Type Description
ndarray

np.ndarray: Cardinal direction name
Source code in src/mechaphlowers/data/geography/helpers.py 
313
314
315
316
317
318
319
320
321
322
323
324
325
326
def bearing_to_direction(
    bearing: np.ndarray,
) -> np.ndarray:
    """
    Convert bearing angle to cardinal direction name
    Args:
        bearing (np.ndarray): Bearing angle in decimal degrees

    Returns:
        np.ndarray: Cardinal direction name
    """
    directions = np.array(['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW'])
    index = np.round(bearing / 45) % 8
    return directions[index.astype(int)]

distances_to_gps 

distances_to_gps(
    lat_a: ndarray,
    lon_a: ndarray,
    x_meters: ndarray,
    y_meters: ndarray,
) -> tuple[ndarray, ndarray]

Calculate GPS coordinates of point B given point A's coordinates and x,y distances in meters.

Parameters:

Name Type Description Default

lat_a

ndarray

Latitude of point A in decimal degrees

 required

lon_a

ndarray

Longitude of point A in decimal degrees

 required

x_meters

ndarray

Distance from west to east in meters (positive = east, negative = west)

 required

y_meters

ndarray

Distance from south to north in meters (positive = north, negative = south)

 required

Returns:

Type Description
tuple[ndarray, ndarray]

tuple[np.ndarray, np.ndarray]: tuple containing the latitude and longitude of point B
Source code in src/mechaphlowers/data/geography/helpers.py 
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
def distances_to_gps(
    lat_a: np.ndarray,
    lon_a: np.ndarray,
    x_meters: np.ndarray,
    y_meters: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Calculate GPS coordinates of point B given point A's coordinates and x,y distances in meters.

    Args:
        lat_a (np.ndarray): Latitude of point A in decimal degrees
        lon_a (np.ndarray): Longitude of point A in decimal degrees
        x_meters (np.ndarray): Distance from west to east in meters (positive = east, negative = west)
        y_meters (np.ndarray): Distance from south to north in meters (positive = north, negative = south)

    Returns:
        tuple[np.ndarray, np.ndarray]: tuple containing the latitude and longitude of point B
    """
    # Convert distances to degrees
    # 1 degree of latitude is approximately 111,111 meters
    lat_change = y_meters / 111111.0

    # 1 degree of longitude varies with latitude
    # At the equator, 1 degree is about 111,111 meters
    # At other latitudes, multiply by cos(latitude)
    lon_change = x_meters / (111111.0 * np.cos(np.radians(lat_a)))

    # Calculate new coordinates
    lat_b = lat_a + lat_change
    lon_b = lon_a + lon_change

    return lat_b, lon_b

gps_to_bearing_two_points 

gps_to_bearing_two_points(
    lat1: ndarray,
    lon1: ndarray,
    lat2: ndarray,
    lon2: ndarray,
    unit: Literal['rad', 'deg'] = 'rad',
) -> ndarray

Calculate the bearing between two points Returns bearing in degrees from north Args: lat1 (np.ndarray): Latitude of point A (in radians by default) lon1 (np.ndarray): Longitude of point A (in radians by default) lat2 (np.ndarray): Latitude of point B (in radians by default) lon2 (np.ndarray): Longitude of point B (in radians by default) unit (Literal["rad", "deg"]): Select the unit to use for both inputs and output. Defaults to "rad"

Returns:

Type Description
ndarray

np.ndarray: Bearing angle in degrees from north, anti clockwise (in radians by default)
Source code in src/mechaphlowers/data/geography/helpers.py 
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
def gps_to_bearing_two_points(
    lat1: np.ndarray,
    lon1: np.ndarray,
    lat2: np.ndarray,
    lon2: np.ndarray,
    unit: Literal["rad", "deg"] = "rad",
) -> np.ndarray:
    """
    Calculate the bearing between two points
    Returns bearing in degrees from north
    Args:
        lat1 (np.ndarray): Latitude of point A (in radians by default)
        lon1 (np.ndarray): Longitude of point A (in radians by default)
        lat2 (np.ndarray): Latitude of point B (in radians by default)
        lon2 (np.ndarray): Longitude of point B (in radians by default)
        unit (Literal["rad", "deg"]): Select the unit to use for both inputs and output. Defaults to "rad"

    Returns:
        np.ndarray: Bearing angle in degrees from north, anti clockwise (in radians by default)
    """
    if unit == "deg":
        lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2])
    dlon = lon2 - lon1
    y = np.sin(dlon) * np.cos(lat2)
    x = np.cos(lat1) * np.sin(lat2) - np.sin(lat1) * np.cos(lat2) * np.cos(
        dlon
    )
    # This formula needs to be reversed to get a anti clockwise angle
    bearing = -np.arctan2(y, x)

    bearing = convert_angle_signed_to_unsigned(bearing)  # Normalize to [0, 2π)
    if unit == "deg":
        bearing = np.degrees(bearing)
    return bearing

gps_to_lambert93 

gps_to_lambert93(
    latitude: Union[float64, ndarray],
    longitude: Union[float64, ndarray],
) -> tuple[ndarray, ndarray]

Convert GPS coordinates (WGS84) to Lambert 93 coordinates.

Parameters:

Name Type Description Default

latitude

Union[float64, ndarray]

Latitude in decimal degrees (WGS84). Can be scalar or numpy array.

 required

longitude

Union[float64, ndarray]

Longitude in decimal degrees (WGS84). Can be scalar or numpy array.

 required

Returns:

Type Description
tuple[ndarray, ndarray]

tuple[np.ndarray, np.ndarray]: (X, Y) coordinates in Lambert 93 projection (in meters)
Source code in src/mechaphlowers/data/geography/helpers.py 
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
def gps_to_lambert93(
    latitude: Union[np.float64, np.ndarray],
    longitude: Union[np.float64, np.ndarray],
) -> tuple[np.ndarray, np.ndarray]:
    """
    Convert GPS coordinates (WGS84) to Lambert 93 coordinates.

    Args:
        latitude: Latitude in decimal degrees (WGS84). Can be scalar or numpy array.
        longitude: Longitude in decimal degrees (WGS84). Can be scalar or numpy array.

    Returns:
        tuple[np.ndarray, np.ndarray]: (X, Y) coordinates in Lambert 93 projection (in meters)
    """

    # Convert inputs to numpy arrays if they aren't already
    latitude = np.asarray(latitude, dtype=np.float64)
    longitude = np.asarray(longitude, dtype=np.float64)

    # Lambert 93 projection parameters
    # Semi-major axis of the GRS80 ellipsoid
    a = 6378137.0
    # Flattening of the GRS80 ellipsoid
    f = 1 / 298.257222101

    # Derived parameters
    e2 = 2 * f - f * f  # First eccentricity squared

    # Lambert 93 specific parameters
    phi0 = np.radians(46.5)  # Latitude of origin
    phi1 = np.radians(44.0)  # First standard parallel
    phi2 = np.radians(49.0)  # Second standard parallel
    lambda0 = np.radians(3.0)  # Central meridian
    X0 = 700000.0  # False easting
    Y0 = 6600000.0  # False northing

    # Convert input coordinates to radians
    phi = np.radians(latitude)
    lambda_deg = np.radians(longitude)

    # Calculate auxiliary functions
    def m(phi):
        return np.cos(phi) / np.sqrt(1 - e2 * np.sin(phi) ** 2)

    def t(phi):
        return np.tan(np.pi / 4 - phi / 2) / (
            (1 - np.sqrt(e2) * np.sin(phi)) / (1 + np.sqrt(e2) * np.sin(phi))
        ) ** (np.sqrt(e2) / 2)

    # Calculate projection constants
    m1 = m(phi1)
    m2 = m(phi2)

    t0 = t(phi0)
    t1 = t(phi1)
    t2 = t(phi2)

    # Calculate n and F
    n = (np.log(m1) - np.log(m2)) / (np.log(t1) - np.log(t2))
    F = m1 / (n * t1**n)

    # Calculate rho0 (radius at origin)
    rho0 = a * F * t0**n

    # Calculate rho and theta for the point
    t_phi = t(phi)
    rho = a * F * t_phi**n
    theta = n * (lambda_deg - lambda0)

    # Calculate Lambert 93 coordinates
    X = X0 + rho * np.sin(theta)
    Y = Y0 + rho0 - rho * np.cos(theta)

    return (X, Y)

haversine 

haversine(
    lat1: ndarray,
    lon1: ndarray,
    lat2: ndarray,
    lon2: ndarray,
    unit: Literal['rad', 'deg'] = 'rad',
) -> ndarray

Calculate the great circle distance between two points on the earth Args: lat1 (np.ndarray): Latitude of point A (in radians by default) lon1 (np.ndarray): Longitude of point A (in radians by default) lat2 (np.ndarray): Latitude of point B (in radians by default) lon2 (np.ndarray): Longitude of point B (in radians by default) unit (Literal["rad", "deg"]): Select the unit to use for angle inputs. Defaults to "rad"

Returns:

Type Description
ndarray

np.ndarray: Distance in meters
Source code in src/mechaphlowers/data/geography/helpers.py 
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
def haversine(
    lat1: np.ndarray,
    lon1: np.ndarray,
    lat2: np.ndarray,
    lon2: np.ndarray,
    unit: Literal["rad", "deg"] = "rad",
) -> np.ndarray:
    """
    Calculate the great circle distance between two points
    on the earth
    Args:
        lat1 (np.ndarray): Latitude of point A (in radians by default)
        lon1 (np.ndarray): Longitude of point A (in radians by default)
        lat2 (np.ndarray): Latitude of point B (in radians by default)
        lon2 (np.ndarray): Longitude of point B (in radians by default)
        unit (Literal["rad", "deg"]): Select the unit to use for angle inputs. Defaults to "rad"

    Returns:
        np.ndarray: Distance in meters
    """
    if unit == "deg":
        lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2])

    # Haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    hav = (
        np.sin(dlat / 2) ** 2
        + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
    )
    c = 2 * np.arcsin(np.sqrt(hav))
    return c * RADIUS_EARTH

lambert93_to_gps 

lambert93_to_gps(
    lambert_e: Union[float64, ndarray],
    lambert_n: Union[float64, ndarray],
) -> tuple[ndarray, ndarray]

Convert Lambert 93 coordinates to WGS84 (longitude, latitude)

Parameters:

Name Type Description Default

lambert_e

Union[float64, ndarray]

Lambert 93 Easting coordinate. Can be scalar or numpy array.

 required

lambert_n

Union[float64, ndarray]

Lambert 93 Northing coordinate. Can be scalar or numpy array.

 required

Returns:

Type Description
tuple[ndarray, ndarray]

tuple[np.ndarray, np.ndarray]: (latitude, longitude) in decimal degrees
Source code in src/mechaphlowers/data/geography/helpers.py 
 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
def lambert93_to_gps(
    lambert_e: Union[np.float64, np.ndarray],
    lambert_n: Union[np.float64, np.ndarray],
) -> tuple[np.ndarray, np.ndarray]:
    """
    Convert Lambert 93 coordinates to WGS84 (longitude, latitude)

    Args:
        lambert_e: Lambert 93 Easting coordinate. Can be scalar or numpy array.
        lambert_n: Lambert 93 Northing coordinate. Can be scalar or numpy array.

    Returns:
        tuple[np.ndarray, np.ndarray]: (latitude, longitude) in decimal degrees
    """

    # Convert inputs to numpy arrays if they aren't already
    lambert_e = np.asarray(lambert_e, dtype=np.float64)
    lambert_n = np.asarray(lambert_n, dtype=np.float64)

    constantes = {
        'GRS80E': 0.081819191042816,
        'LONG_0': 3,
        'XS': 700000,
        'YS': 12655612.0499,
        'n': 0.7256077650532670,
        'C': 11754255.4261,
    }

    del_x = lambert_e - constantes['XS']
    del_y = lambert_n - constantes['YS']
    gamma = np.arctan(-del_x / del_y)
    r = np.sqrt(del_x * del_x + del_y * del_y)
    latiso = np.log(constantes['C'] / r) / constantes['n']

    # Iterative calculation for sinPhiit
    sin_phi_it0 = np.tanh(
        latiso
        + constantes['GRS80E'] * np.arctanh(constantes['GRS80E'] * np.sin(1))
    )
    sin_phi_it1 = np.tanh(
        latiso
        + constantes['GRS80E'] * np.arctanh(constantes['GRS80E'] * sin_phi_it0)
    )
    sin_phi_it2 = np.tanh(
        latiso
        + constantes['GRS80E'] * np.arctanh(constantes['GRS80E'] * sin_phi_it1)
    )
    sin_phi_it3 = np.tanh(
        latiso
        + constantes['GRS80E'] * np.arctanh(constantes['GRS80E'] * sin_phi_it2)
    )
    sin_phi_it4 = np.tanh(
        latiso
        + constantes['GRS80E'] * np.arctanh(constantes['GRS80E'] * sin_phi_it3)
    )
    sin_phi_it5 = np.tanh(
        latiso
        + constantes['GRS80E'] * np.arctanh(constantes['GRS80E'] * sin_phi_it4)
    )
    sin_phi_it6 = np.tanh(
        latiso
        + constantes['GRS80E'] * np.arctanh(constantes['GRS80E'] * sin_phi_it5)
    )

    long_rad = np.arcsin(sin_phi_it6)
    lat_rad = gamma / constantes['n'] + np.deg2rad(constantes['LONG_0'])

    longitude = np.rad2deg(lat_rad)
    latitude = np.rad2deg(long_rad)

    return (latitude, longitude)

reverse_haversine 

reverse_haversine(
    lat: ndarray,
    lon: ndarray,
    bearing: ndarray,
    distance: ndarray,
) -> tuple[ndarray, ndarray]

Implementation of the reverse of Haversine formula. Takes one set of latitude/longitude as a start point, a bearing, and a distance, and returns the resultant lat/long pair.

Parameters:

Name Type Description Default

lat

ndarray

Starting latitude in decimal degrees

 required

lon

ndarray

Starting longitude in decimal degrees

 required

bearing

ndarray

Bearing in decimal degrees

 required

distance

ndarray

Distance in meters

 required

Returns:

Type Description
tuple[ndarray, ndarray]

tuple[np.ndarray, np.ndarray]: tuple containing the latitude and longitude of the result
Source code in src/mechaphlowers/data/geography/helpers.py 
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
202
203
204
205
206
207
def reverse_haversine(
    lat: np.ndarray,
    lon: np.ndarray,
    bearing: np.ndarray,
    distance: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Implementation of the reverse of Haversine formula. Takes one set of
    latitude/longitude as a start point, a bearing, and a distance, and
    returns the resultant lat/long pair.

    Args:
        lat (np.ndarray): Starting latitude in decimal degrees
        lon (np.ndarray): Starting longitude in decimal degrees
        bearing (np.ndarray): Bearing in decimal degrees
        distance (np.ndarray): Distance in meters

    Returns:
        tuple[np.ndarray, np.ndarray]: tuple containing the latitude and longitude of the result
    """

    lat1 = np.radians(lat)
    lon1 = np.radians(lon)
    angdist = distance / RADIUS_EARTH
    theta = np.radians(bearing)

    lat2 = np.degrees(
        np.arcsin(
            np.sin(lat1) * np.cos(angdist)
            + np.cos(lat1) * np.sin(angdist) * np.cos(theta)
        )
    )

    lon2 = np.degrees(
        lon1
        + np.arctan2(
            np.sin(theta) * np.sin(angdist) * np.cos(lat1),
            np.cos(angdist) - np.sin(lat1) * np.sin(np.radians(lat2)),
        )
    )

    return (lat2, lon2)

reverse_haversine_float 

reverse_haversine_float(
    lat_rad: float,
    lon_rad: float,
    bearing: float,
    distance: float,
) -> tuple[float, float]

Reverse of haversine with floats.

Returns next coordinates in lat/lon according to bearing and distance.

Parameters:

Name Type Description Default

lat_rad

float

starting latitude in radians

 required

lon_rad

float

starting longitude in radians

 required

bearing

float

bearing angle: orientation of the line in radians, anticlockwise, towards north = 0

 required

distance

float

distance with the next point

 required

Returns:

Type Description
tuple[float, float]

tuple[float, float]: (lat2, lon2): GPS coordinates of next point
Source code in src/mechaphlowers/data/geography/helpers.py 
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
def reverse_haversine_float(
    lat_rad: float,
    lon_rad: float,
    bearing: float,
    distance: float,
) -> tuple[float, float]:
    """Reverse of haversine with floats.

    Returns next coordinates in lat/lon according to bearing and distance.

    Args:
        lat_rad (float): starting latitude in radians
        lon_rad (float): starting longitude in radians
        bearing (float): bearing angle: orientation of the line in radians, anticlockwise, towards north = 0
        distance (float): distance with the next point

    Returns:
        tuple[float, float]: (lat2, lon2): GPS coordinates of next point
    """
    angdist = distance / RADIUS_EARTH

    lat_rad_2 = np.arcsin(
        np.sin(lat_rad) * np.cos(angdist)
        + np.cos(lat_rad) * np.sin(angdist) * np.cos(bearing)
    )

    lon_rad_2 = lon_rad + np.arctan2(
        -np.sin(bearing) * np.sin(angdist) * np.cos(lat_rad),
        np.cos(angdist) - np.sin(lat_rad) * np.sin(lat_rad_2),
    )
    return (lat_rad_2, lon_rad_2)

support_distances_to_gps 

support_distances_to_gps(
    support_bases_x_meters: ndarray,
    support_bases_y_meters: ndarray,
    first_support_lat: float64,
    first_support_lon: float64,
) -> tuple[ndarray, ndarray]

Parse support distances to calculate gps coordinates in the form: Args: support_bases_x_meters (np.typing.NDArray[np.float64]): Array of x distances from first support base to support bases excluding the first support base support_bases_y_meters (np.typing.NDArray[np.float64]): Array of y distances from first support base to support bases excluding the first support base first_support_lat (np.float64): Latitude of the first support base first_support_lon (np.float64): Longitude of the first support base

Returns:

Type Description
tuple[ndarray, ndarray]

tuple[np.typing.NDArray[np.float64], np.typing.NDArray[np.float64]]: tuple containing the latitude and longitude of the support bases
Source code in src/mechaphlowers/data/geography/helpers.py 
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
def support_distances_to_gps(
    support_bases_x_meters: np.ndarray,
    support_bases_y_meters: np.ndarray,
    first_support_lat: np.float64,
    first_support_lon: np.float64,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Parse support distances to calculate gps coordinates in the form:
    Args:
        support_bases_x_meters (np.typing.NDArray[np.float64]): Array of x distances from first support base to support bases excluding the first support base
        support_bases_y_meters (np.typing.NDArray[np.float64]): Array of y distances from first support base to support bases excluding the first support base
        first_support_lat (np.float64): Latitude of the first support base
        first_support_lon (np.float64): Longitude of the first support base

    Returns:
        tuple[np.typing.NDArray[np.float64], np.typing.NDArray[np.float64]]: tuple containing the latitude and longitude of the support bases
    """
    # Calculate GPS coordinates for additional support bases
    additional_lats, additional_lons = distances_to_gps(
        np.full(
            len(support_bases_x_meters), first_support_lat, dtype=np.float64
        ),
        np.full(
            len(support_bases_x_meters), first_support_lon, dtype=np.float64
        ),
        support_bases_x_meters,
        support_bases_y_meters,
    )

    # Concatenate additional coordinates with first support coordinates
    all_lats = np.concatenate(
        [np.array([first_support_lat], dtype=np.float64), additional_lats]
    )
    all_lons = np.concatenate(
        [np.array([first_support_lon], dtype=np.float64), additional_lons]
    )

    return all_lats, all_lons