Merge branch 'great-circle' into 'main'

Move GreatCircle class objects from PyCOADS to GeoSpatialTools

See merge request nocsurfaceprocesses/geospatialtools!15

GeoSpatialTools/__init__.py

 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",

GeoSpatialTools/great_circle.py

+"""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),
+    )

test/test_greatcircle.py

+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