"""
GreatCircle
-----------
Constructors and methods for interacting with GreatCircle objects, including
comparisons between GreatCircle objects.
"""
import numpy as np
from .distance_metrics import bearing, gcd_slc
def cartesian_to_lonlat(
x: float,
y: float,
z: float,
to_radians: bool = False,
) -> tuple[float, float]:
"""
Get lon, and lat from cartesian coordinates.
Parameters
----------
x : float
x coordinate
y : float
y coordinate
z : float
z coordinate
to_radians : bool
Return angles in radians. Otherwise return values in degrees.
Returns
-------
(float, float)
lon, lat
"""
R = np.sqrt(x**2 + y**2 + z**2)
x /= R
y /= R
z /= R
lat = np.arcsin(z)
tmp = np.cos(lat)
sign = np.arcsin(y / tmp)
lon = np.arccos(x / tmp) * sign
if to_radians:
return lon, lat
return np.degrees(lon), np.degrees(lat)
def polar_to_cartesian(
lon: float,
lat: float,
R: float = 6371,
to_radians: bool = True,
normalised: bool = True,
) -> tuple[float, float, float]:
"""
Convert from polars coordinates to cartesian.
Get cartesian coordinates from spherical polar coordinates. Default
behaviour assumes lon and lat, so converts to radians. Set
`to_radians=False` if the coordinates are already in radians.
Parameters
----------
lon : float
Longitude.
lat : float
Latitude.
R : float
Radius of sphere.
to_radians : bool
Convert lon and lat to radians.
normalised : bool
Return normalised vector (ignore R value).
Returns
-------
(float, float, float)
x, y, z cartesian coordinates.
"""
theta = np.radians(lon) if to_radians else lon
phi = np.radians(lat) if to_radians else lat
x = np.cos(theta) * np.cos(phi)
y = np.sin(theta) * np.cos(phi)
z = np.sin(phi)
return (x, y, z) if normalised else (R * x, R * y, R * z)
class GreatCircle:
"""
A GreatCircle object for a pair of positions.
Construct a great circle path between a pair of positions.
https://www.boeing-727.com/Data/fly%20odds/distance.html
Parameters
----------
lon0 : float
Longitude of start position.
lat0 : float
Latitude of start position.
lon1 : float
Longitude of end position.
lat1 : float
Latitude of end position.
R : float
Radius of the sphere. Default is Earth radius in km (6371.0).
"""
def __init__(
self,
lon0: float,
lat0: float,
lon1: float,
lat1: float,
R: float = 6371,
) -> None:
self.lon0 = lon0
self.lat0 = lat0
self.lon1 = lon1
self.lat1 = lat1
self.R = R
self.cross_prod = _cross_lonlat(
self.lon0, self.lat0, self.lon1, self.lat1
)
self.cross_prod_dist = np.linalg.norm(self.cross_prod)
self.bearing = bearing(self.lon0, self.lat0, self.lon1, self.lat1)
self.dist = gcd_slc(self.lon0, self.lat0, self.lon1, self.lat1)
def dist_from_point(
self,
lon: float,
lat: float,
) -> float:
"""
Compute distance from the GreatCircle to a point on the sphere.
Parameters
----------
lon : float
Longitude of the position to test.
lat : float
Longitude of the position to test.
Returns
-------
float
Minimum distance between point and the GreatCircle arc.
"""
cart = polar_to_cartesian(lon, lat, normalised=True)
num = np.dot(cart, self.cross_prod)
# WARN: This can be negative - hence using abs
return np.abs(np.arcsin(num / self.cross_prod_dist) * self.R)
def _identical_plane(
self,
other: object,
epsilon: float = 0.01,
) -> bool:
"""
Identify if other GreatCircle has the same plane.
Determined by comparing the norms of the planes constructed from the
two points and the centre of the sphere.
Returns True if the planes formed by the two great circles are
parallel, i.e. the normals defining the planes have the same
direction, or the exact opposite direction. This would mean that a
GreatCircle compared with the GreatCircle with oppposite start and
end points would return True.
Parameters
----------
other : GreatCircle
Intersecting GreatCircle object
epsilon : float
Threshold for intersection
Returns
-------
bool
Indicating if the planes formed by two GreatCircle objects are the
same (or mirrored) to within a given threshold.
"""
if not isinstance(other, GreatCircle):
raise TypeError("Input 'other' is not a GreatCircle")
return bool(
np.isclose(self.cross_prod, other.cross_prod, atol=epsilon).all()
or np.isclose(
self.cross_prod, -other.cross_prod, atol=epsilon
).all()
)
def intersection(
self, other: object, epsilon: float = 0.01
) -> tuple[float, float] | None:
"""
Determine intersection position with another GreatCircle.
Determine the location at which the GreatCircle intersects another
GreatCircle arc. (To within some epsilon threshold).
Returns `None` if there is no solution - either because there is no
intersection point, or the planes generated from the arc and centre of
the sphere are identical.
Parameters
----------
other : GreatCircle
Intersecting GreatCircle object
epsilon : float
Threshold for intersection
Returns
-------
(float, float) | None
Position of intersection
"""
if not isinstance(other, GreatCircle):
raise TypeError("Input 'other' is not a GreatCircle")
if self.R != other.R:
raise ValueError("GreatCircle radius values do not match")
if self._identical_plane(other, epsilon=epsilon):
return None
plane_intersection = np.cross(self.cross_prod, other.cross_prod)
epsilon *= self.R
S1 = plane_intersection / np.linalg.norm(plane_intersection)
lon, lat = cartesian_to_lonlat(*S1)
if self.dist_from_point(lon, lat) < epsilon:
return lon, lat
S2 = -S1
lon, lat = cartesian_to_lonlat(*S2)
if self.dist_from_point(lon, lat) < epsilon:
return lon, lat
return None
def intersection_angle(
self,
other: object,
epsilon: float = 0.01,
) -> float | None:
"""
Get angle of intersection with another GreatCircle.
Get the angle of intersection with another GreatCircle arc. Returns
None if there is no intersection.
The intersection angle is computed using the normals of the planes
formed by the two intersecting great circle objects.
Parameters
----------
other : GreatCircle
Intersecting GreatCircle object
epsilon : float
Threshold for intersection
Returns
-------
float | None
Intersection angle in degrees
"""
if not isinstance(other, GreatCircle):
raise TypeError("'other' must be of type 'GreatCircle'")
# Make sure we have an intersection!
if self.intersection(other, epsilon) is None:
return None
# INFO: Want to use self.cross and other.cross which are the normals
angle = np.arccos(
np.dot(self.cross_prod, other.cross_prod)
/ (self.cross_prod_dist * other.cross_prod_dist)
)
return np.rad2deg(angle)
def _cross_lonlat(
lon0: float,
lat0: float,
lon1: float,
lat1: float,
) -> np.ndarray:
"""
Get the cross-product between two positions on a sphere.
|u_1| x |u_2|
Parameters
----------
lon0 : float
Longitude of position 0 in degrees.
lat0 : float
Latitude of position 0 in degrees.
lon1 : float
Longitude of position 1 in degrees.
lat1 : float
Latitude of position 1 in degrees.
Returns
-------
np.ndarray
Cartesian vector of the cross product of the input lon/lat positions
(assuming a sphere of radius 1).
"""
return np.cross(
polar_to_cartesian(lon0, lat0, normalised=True),
polar_to_cartesian(lon1, lat1, normalised=True),
)