diff --git a/GeoSpatialTools/__init__.py b/GeoSpatialTools/__init__.py
index 11e1e9388c3b2a90fa76a7649632b55e77c63ef1..e3a054eb5a7cc0269477bf2e346f51f4071b7370 100644
--- a/GeoSpatialTools/__init__.py
+++ b/GeoSpatialTools/__init__.py
@@ -1,10 +1,12 @@
from .neighbours import find_nearest
from .distance_metrics import haversine
+from .great_circle import GreatCircle
from .quadtree import Ellipse, QuadTree, Record, Rectangle
from .kdtree import KDTree
__all__ = [
"Ellipse",
+ "GreatCircle",
"KDTree",
"QuadTree",
"Record",
diff --git a/GeoSpatialTools/great_circle.py b/GeoSpatialTools/great_circle.py
new file mode 100644
index 0000000000000000000000000000000000000000..b2cb6d9b8fa1430845e6f52bed196d7fd7d0aa7b
--- /dev/null
+++ b/GeoSpatialTools/great_circle.py
@@ -0,0 +1,305 @@
+"""Class for a Great Circle object."""
+
+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),
+ )
diff --git a/test/test_greatcircle.py b/test/test_greatcircle.py
new file mode 100644
index 0000000000000000000000000000000000000000..06bdab2a5120280f3a4c9f25f2cd48baa6d85cef
--- /dev/null
+++ b/test/test_greatcircle.py
@@ -0,0 +1,46 @@
+import unittest
+import numpy as np
+from GeoSpatialTools.great_circle import GreatCircle
+from GeoSpatialTools.distance_metrics import gcd_slc
+
+
+class TestGreatCircle(unittest.TestCase):
+ halifax = (-63.5728, 44.6476)
+ southampton = (-1.4049, 50.9105)
+
+ def test_constructor(self):
+ out = GreatCircle(*self.halifax, *self.southampton)
+ assert hasattr(out, "dist")
+ assert out.dist == gcd_slc(*self.southampton, *self.halifax)
+
+ def test_meridian_planes(self):
+ # TEST: That great circles from a pole stay on the same longitude
+ lon0, lat0 = 45, 23
+ GC1 = GreatCircle(0, 90, lon0, lat0)
+
+ assert GC1.dist_from_point(-lon0, lat0 + 5) > 10
+ for lat in range(lat0, 90, 2):
+ dist = GC1.dist_from_point(lon0, lat)
+ assert dist < 0.01
+
+ GC2 = GreatCircle(0, -90, lon0, -lat0)
+ assert abs(GC1.dist - GC2.dist) < 0.01
+ for lat in range(-lat0, -90, -2):
+ dist = GC2.dist_from_point(lon0, lat)
+ assert dist < 0.01
+
+ assert GC1._identical_plane(GC2)
+
+ def test_equator_vs_meridian(self):
+ # TEST: that a great circle at the equator and a line of longitude
+ # intersect with angle 90
+ GC0 = GreatCircle(-5, 0, 5, 0)
+ GC1 = GreatCircle(0, -5, 0, 5)
+ assert np.isclose([GC0.dist], [GC1.dist]).all()
+ assert GC1.dist_from_point(0, 0) < 0.01
+ # assert GC0.dist_from_point(0, 0) < 0.01
+
+ int_pts = GC0.intersection(GC1)
+ int_ang = GC0.intersection_angle(GC1)
+ assert int_pts == (0, 0)
+ assert int_ang == 90