from typing import Any, Tuple
import numpy as np
import numpy.typing as npt
from pyproj import Transformer
from scipy.interpolate import splev, splprep
from .models import CoordinateListInput
EARTH_RADIUS = 6371e3 # meters
[docs]
def gps_to_cartesian(
coords: CoordinateListInput, origin_offset: bool = True
) -> Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
"""
Converts GPS coordinates (lat/lon) to Cartesian coordinates (x, y) in meters.
Uses Web Mercator projection (EPSG:3857).
Parameters
----------
coords : CoordinateListInput
GPS coordinates to convert
origin_offset : bool, optional
If True, offsets coordinates so the first point is at (0, 0) (default: True)
Returns
-------
Tuple[NDArray[float64], NDArray[float64]]
(x_m, y_m) - Arrays of x and y coordinates in meters
"""
# Convert from WGS84 (EPSG:4326) to Web Mercator (EPSG:3857)
transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
x_m, y_m = transformer.transform(coords.longitudes, coords.latitudes)
x_m = np.array(x_m, dtype=np.float64)
y_m = np.array(y_m, dtype=np.float64)
# Offset to start at origin if requested
if origin_offset:
x_m = x_m - x_m[0]
y_m = y_m - y_m[0]
return x_m, y_m
[docs]
def gps_interpolate(
coords: CoordinateListInput, smooth: float = 0
) -> Tuple[CoordinateListInput, npt.NDArray[np.float64], Any]:
"""
Cleans GPS data, then interpolates to get a new set of coordinates with at
a constant distance interval
Parameters
----------
coords : CoordinateListInput
Specifies the original coordinates, as well as the constant distance interval to interpolate at
through the distance_step field
smooth : float, optional
Smoothing parameter for spline interpolation (default: 0)
Returns
-------
Tuple[CoordinateListInput, Any, npt.NDArray[np.float64]]
(coords_new, spline, cumulative_dist_new) - A CoordinateListInput object with the same GPS data
but at constant distance intervals, the underlying scipy spline object, and a cumulative distance
array based on the new coordinates
"""
_gps_data = gps_remove_duplicates(coords)
_gps_dist = gps_calculate_distance(_gps_data)
xi = np.cumsum(_gps_dist)
# Create smoothed spline
tck, u = splprep(
[_gps_data.latitudes, _gps_data.longitudes], u=xi, s=smooth
) # 's' controls smoothness
num = int(np.ceil(xi[-1] / coords.distance_step))
cumulative_dist_new = np.linspace(xi[0], xi[-1], num)
lat_smooth, lon_smooth = splev(cumulative_dist_new, tck)
return (
CoordinateListInput(
latitudes=lat_smooth.tolist(),
longitudes=lon_smooth.tolist(),
distance_step=coords.distance_step,
sample_dist=coords.sample_dist,
),
tck,
cumulative_dist_new,
)
[docs]
def gps_calculate_distance(
coords: CoordinateListInput,
) -> npt.NDArray[np.float64]:
"""
Returns the distance between trackpoints using the haversine formula.
Parameters
----------
coords : CoordinateListInput
GPS coordinates to calculate distances for
Returns
-------
npt.NDArray[np.float64]
Index i is the distance from the i - 1 to the i-th coordinate
"""
gps_dist = np.zeros(len(coords.latitudes), dtype=np.float64)
for i in range(len(gps_dist) - 1):
lat1 = np.radians(coords.latitudes[i])
lon1 = np.radians(coords.longitudes[i])
lat2 = np.radians(coords.latitudes[i + 1])
lon2 = np.radians(coords.longitudes[i + 1])
delta_lat = lat2 - lat1
delta_lon = lon2 - lon1
c = 2.0 * np.arcsin(
np.sqrt(
np.sin(delta_lat / 2.0) ** 2
+ np.cos(lat1) * np.cos(lat2) * np.sin(delta_lon / 2.0) ** 2
)
) # haversine formula
gps_dist[i + 1] = EARTH_RADIUS * c # great-circle distance
return gps_dist
[docs]
def gps_remove_duplicates(coords: CoordinateListInput) -> CoordinateListInput:
"""
De-duplicate GPS coordinates by removing coordinates with distance of 0.0 between them
Parameters
----------
coords : CoordinateListInput
GPS coordinates to de-duplicate
Returns
-------
CoordinateListInput
De-duplicated GPS coordinates
"""
gps_dist = gps_calculate_distance(coords)
i_dist = np.insert(np.nonzero(gps_dist)[0], 0, 0) # keep gps_dist[0] = 0.0
if len(i_dist) == len(gps_dist):
return coords
return CoordinateListInput(
longitudes=[coords.longitudes[idx] for idx in i_dist],
latitudes=[coords.latitudes[idx] for idx in i_dist],
distance_step=coords.distance_step,
sample_dist=coords.sample_dist,
)