[Future Spherical Class] Add SArc

pyresample/future/spherical/__init__.py

@@ -17,4 +17,5 @@
 # with this program.  If not, see <http://www.gnu.org/licenses/>.
 """Future features that are backwards incompatible with current functionality."""
 
+from .arc import SArc  # noqa
 from .point import SMultiPoint, SPoint  # noqa

pyresample/future/spherical/arc.py

@@ -0,0 +1,283 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2022 Pyresample developers
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""Define a great-circle arc between two points on the sphere through the SArc class."""
+
+import numpy as np
+import pyproj
+from shapely.geometry import LineString
+
+from pyresample.spherical import Arc
+
+from .point import SPoint, create_spherical_point
+
+EPSILON = 0.0000001
+
+
+def _check_valid_arc(start_point, end_point):
+    """Check arc validity."""
+    if start_point == end_point:
+        raise ValueError("An SArc can not be represented by the same start and end SPoint.")
+    if start_point.is_pole() and end_point.is_pole():
+        raise ValueError("An SArc can not be uniquely defined between the two poles.")
+    if start_point.is_on_equator() and end_point.is_on_equator() and abs(start_point.lon - end_point.lon) == np.pi:
+        raise ValueError(
+            "An SArc can not be uniquely defined on the equator if start and end points are 180 degrees apart.")
+    if start_point.get_antipode() == end_point:
+        raise ValueError("An SArc can not be uniquely defined between antipodal points.")
+
+
+class SArc(Arc):
+    """Object representing a great-circle arc between two points on a sphere.
+
+    The ``start_point`` and ``end_point`` must be SPoint objects.
+    The great-circle arc is defined as the shortest track(s) between the two points.
+    Between the north and south pole there are an infinite number of great-circle arcs.
+    """
+
+    def __init__(self, start_point, end_point):
+        _check_valid_arc(start_point, end_point)
+        super().__init__(start_point, end_point)
+
+    def __hash__(self):
+        """Define SArc hash."""
+        return hash((float(self.start.lon), float(self.start.lat),
+                     float(self.end.lon), float(self.end.lat)))
+
+    def is_on_equator(self):
+        """Check if the SArc lies on the equator."""
+        return self.start.lat == 0 and self.end.lat == 0
+
+    def __eq__(self, other):
+        """Check equality."""
+        return self.start == other.start and self.end == other.end
+
+    def reverse_direction(self):
+        """Reverse SArc direction."""
+        return SArc(self.end, self.start)
+
+    @property
+    def vertices(self):
+        """Get start SPoint and end SPoint (radians) vertices array."""
+        return self.start.vertices, self.end.vertices
+
+    @property
+    def vertices_in_degrees(self):
+        """Get start SPoint and end SPoint (degrees) vertices array."""
+        return self.start.vertices_in_degrees, self.end.vertices_in_degrees
+
+    def get_next_intersection(self, other_arc):
+        """Overwrite a Arc deprecated function. Inherited from Arc class."""
+        raise ValueError("This function is deprecated.")
+
+    def intersections(self, other_arc):
+        """Overwritea Arc deprecated function. Inherited from Arc class."""
+        raise ValueError("'SArc.intersections' is deprecated. Use '_great_circle_intersections' instead.")
+
+    def _great_circle_intersections(self, other_arc):
+        """Compute the intersections points of the greats circles over which the arcs lies.
+
+        A great circle divides the sphere in two equal hemispheres.
+        """
+        end_lon = self.end.lon
+        other_end_lon = other_arc.end.lon
+
+        if self.end.lon - self.start.lon > np.pi:
+            end_lon -= 2 * np.pi
+        if other_arc.end.lon - other_arc.start.lon > np.pi:
+            other_end_lon -= 2 * np.pi
+        if self.end.lon - self.start.lon < -np.pi:
+            end_lon += 2 * np.pi
+        if other_arc.end.lon - other_arc.start.lon < -np.pi:
+            other_end_lon += 2 * np.pi
+
+        end_point = SPoint(end_lon, self.end.lat)
+        other_end_point = SPoint(other_end_lon, other_arc.end.lat)
+
+        ea_ = self.start.cross2cart(end_point).normalize()
+        eb_ = other_arc.start.cross2cart(other_end_point).normalize()
+
+        cross = ea_.cross(eb_)
+        lat = np.arctan2(cross.cart[2],
+                         np.sqrt(cross.cart[0] ** 2 + cross.cart[1] ** 2))
+        lon = np.arctan2(cross.cart[1], cross.cart[0])
+
+        return (SPoint(lon, lat), SPoint(lon + np.pi, -lat))
+
+    def intersection(self, other_arc):
+        """Overwrite a Arc deprecated function. Inherited from Arc class."""
+        raise ValueError("'SArc.intersection' is deprecated. Use 'intersection_point' instead.")
+
+    def intersection_point(self, other_arc):
+        """Compute the intersection point between two great circle arcs.
+
+        If arc and *other_arc* does not intersect, it returns None.
+        If arc and *other_arc* are the same (whatever direction), it returns None.
+        If arc and *other_arc* overlaps, within or contain, it returns None.
+        If arc and *other_arc* intersect, it returns the intersection SPoint.
+        """
+        # If same arc (same direction), return None
+        if self == other_arc:
+            return None
+
+        great_circles_intersection_spoints = self._great_circle_intersections(other_arc)
+
+        for spoint in great_circles_intersection_spoints:
+            a = self.start
+            b = self.end
+            c = other_arc.start
+            d = other_arc.end
+
+            ab_dist = a.hdistance(b)
+            cd_dist = c.hdistance(d)
+            ap_dist = a.hdistance(spoint)
+            bp_dist = b.hdistance(spoint)
+            cp_dist = c.hdistance(spoint)
+            dp_dist = d.hdistance(spoint)
+
+            if (((spoint in (a, b)) or (abs(ap_dist + bp_dist - ab_dist) < EPSILON)) and
+                    ((spoint in (c, d)) or (abs(cp_dist + dp_dist - cd_dist) < EPSILON))):
+                return spoint
+        return None
+
+    def intersects(self, other_arc):
+        """Check if the current Sarc and another SArc intersect."""
+        return bool(self.intersection_point(other_arc))
+
+    def midpoint(self, ellips='sphere'):
+        """Return the SArc midpoint SPoint."""
+        geod = pyproj.Geod(ellps=ellips)
+        lon_start = self.start.lon
+        lon_end = self.end.lon
+        lat_start = self.start.lat
+        lat_end = self.end.lat
+        lon_mid, lat_mid = geod.npts(lon_start, lat_start, lon_end, lat_end, npts=1, radians=True)[0]
+        return SPoint(lon_mid, lat_mid)
+
+    def to_shapely(self):
+        """Convert to Shapely LineString."""
+        start_coord, end_coord = self.vertices_in_degrees
+        return LineString((start_coord.tolist()[0], end_coord.tolist()[0]))
+
+    def forward_points(self, distance, ellps="sphere"):
+        """Get points at given distance(s) from SArc end point in the forward direction.
+
+        The distance value(s) must be provided in meters.
+        If the distance is positive, the point will be located outside the SArc.
+        If the distance is negative, the point will be located inside the SArc.
+        The function returns an SPoint or SMultiPoint.
+        """
+        # Define geoid
+        geod = pyproj.Geod(ellps=ellps)
+        # Retrieve forward and back azimuth
+        fwd_az, back_az, _ = geod.inv(self.start.lon, self.start.lat,
+                                      self.end.lon, self.end.lat, radians=True)
+        # Retreve forward points
+        distance = np.array(distance)
+        lon, lat, back_az = geod.fwd(np.broadcast_to(np.array(self.end.lon), distance.shape),
+                                     np.broadcast_to(np.array(self.end.lat), distance.shape),
+                                     az=np.broadcast_to(fwd_az, distance.shape),
+                                     dist=distance, radians=True)
+        p = create_spherical_point(lon, lat)
+        return p
+
+    def backward_points(self, distance, ellps="sphere"):
+        """Get points at given distance(s) from SArc start point in the backward direction.
+
+        The distance value(s) must be provided in meters.
+        If the distance is positive, the point will be located outside the SArc.
+        If the distance is negative, the point will be located inside the SArc.
+        The function returns an SPoint or SMultiPoint.
+        """
+        reverse_arc = self.reverse_direction()
+        return reverse_arc.forward_points(distance=distance, ellps=ellps)
+
+    def extend(self, distance, direction="both", ellps="sphere"):
+        """Extend the SArc of a given distance in both, forward or backward directions.
+
+        If the distance is positive, it extends the SArc.
+        If the distance is negative, it shortens the SArc.
+        """
+        valid_direction = ["both", "forward", "backward"]
+        if direction not in valid_direction:
+            raise ValueError(f"Valid direction values are: {valid_direction}")
+        if direction in ["both", "forward"]:
+            end_point = self.forward_points(distance=distance, ellps="sphere")
+        else:
+            end_point = self.end
+
+        if direction in ["both", "backward"]:
+            start_point = self.backward_points(distance=distance, ellps="sphere")
+        else:
+            start_point = self.start
+        arc = create_spherical_arcs(start_point, end_point)
+        return arc
+
+    def shorten(self, distance, direction="both", ellps="sphere"):
+        """Short the SArc of a given distance in both, forward or backward directions.
+
+        If the distance is positive, it shortens the SArc.
+        If the distance is negative, it extends the SArc.
+        """
+        return self.extend(distance=-distance, direction=direction, ellps="sphere")
+
+    # def segmentize(self, npts=0, max_distance=0, ellips='sphere'):
+    #     """Segmentize the great-circle arc.
+
+    #     It returns an SLine.
+    #     npts or max_distance are mutually exclusively. Specify one of them.
+    #     max_distance must be provided in meters.
+    #     """
+    #     if npts != 0:
+    #         npts = npts + 2  # + 2 to account for initial and terminus points
+    #     geod = pyproj.Geod(ellps=ellips)
+    #     lon_start = self.start.lon
+    #     lon_end = self.end.lon
+    #     lat_start = self.start.lat
+    #     lat_end = self.end.lat
+    #     points = geod.inv_intermediate(lon_start, lat_start, lon_end, lat_end,
+    #                                    del_s=max_distance,
+    #                                    npts=npts,
+    #                                    radians=True,
+    #                                    initial_idx=0, terminus_idx=0)
+    #     lons, lats = (points.lons, points.lats)
+    #     lons = np.asarray(lons)
+    #     lats = np.asarray(lats)
+    #     vertices = np.stack((lons, lats)).T
+    #     return SLine(vertices)
+
+    # def to_line(self):
+    #     """Convert to SLine."""
+    #     vertices = np.vstack(self.vertices)
+    #     return SLine(vertices)
+
+    # def plot(self, *args, **kwargs):
+    #     """Convert to SLine."""
+    #     self.to_line.plot(*args, **kwargs)
+
+
+def create_spherical_arcs(start_point, end_point):
+    """Create a SArc or SArcs class depending on the number of points.
+
+    If a Spoint is provided, it returns an SArc.
+    If a SMultiPoint is provided, it returns an SArcs.
+    """
+    if isinstance(start_point, SPoint) and isinstance(end_point, SPoint):
+        arc = SArc(start_point, end_point)
+    else:
+        raise NotImplementedError("SArcs class is not yet available.")
+    return arc

pyresample/future/spherical/point.py

@@ -39,13 +39,34 @@ def from_degrees(cls, lon, lat):
         """Create SPoint from lon/lat coordinates in degrees."""
         return cls(np.deg2rad(lon), np.deg2rad(lat))
 
+    def is_pole(self):
+        """Test if the point is on a pole."""
+        return self.lat in (-np.pi / 2, np.pi / 2)
+
+    def is_on_equator(self):
+        """Test if the point is on the equator."""
+        return self.lat == 0
+
+    def get_antipode(self):
+        """Return the antipode point."""
+        lat = - self.lat
+        lon = self.lon - np.pi
+        return SPoint(lon, lat)
+
     def __str__(self):
         """Get simplified representation of lon/lat arrays in radians."""
         return str((float(self.lon), float(self.lat)))
 
     def __repr__(self):
         """Get simplified representation of lon/lat arrays in radians."""
-        return str((float(self.lon), float(self.lat)))
+        coords = str((float(self.lon), float(self.lat)))
+        return f"{self.__class__.__name__}: {coords}"
+
+    def __eq__(self, other):
+        """Check equality."""
+        if self.lat == other.lat and self.is_pole():
+            return True
+        return np.allclose((self.lon, self.lat), (other.lon, other.lat))
 
     def to_shapely(self):
         """Convert the SPoint to a shapely Point (in lon/lat degrees)."""
@@ -69,20 +90,48 @@ def from_degrees(cls, lon, lat):
         """Create SMultiPoint from lon/lat coordinates in degrees."""
         return cls(np.deg2rad(lon), np.deg2rad(lat))
 
+    def get_antipodes(self):
+        """Return the antipodal points."""
+        lat = - self.lat
+        lon = self.lon - np.pi
+        return SMultiPoint(lon, lat)
+
     def __eq__(self, other):
         """Check equality."""
-        return np.allclose(self.lon, other.lon) and np.allclose(self.lat, other.lat)
+        lon1 = self.lon.copy()
+        lon2 = other.lon.copy()
+        lat1 = self.lat
+        lat2 = other.lat
+        # Set longitude 0 at the pole for array comparison
+        lon1[np.isin(lat1, [-np.pi / 2, np.pi / 2])] = 0
+        lon2[np.isin(lat2, [-np.pi / 2, np.pi / 2])] = 0
+        return np.allclose(lon1, lon2) and np.allclose(lat1, lat2)
 
     def __str__(self):
         """Get simplified representation of lon/lat arrays in radians."""
         return str(self.vertices)
 
     def __repr__(self):
         """Get simplified representation of lon/lat arrays in radians."""
-        return str(self.vertices)
+        return f"{self.__class__.__name__}:\n {self.vertices}"
 
     def to_shapely(self):
         """Convert the SMultiPoint to a shapely MultiPoint (in lon/lat degrees)."""
         from shapely.geometry import MultiPoint
         point = MultiPoint(self.vertices_in_degrees)
         return point
+
+
+def create_spherical_point(lon, lat):
+    """Create a SPoint or SMultiPoint class depending on the number of point coordinates.
+
+    If a single point coordinate is provided, it returns an SPoint.
+    If multiple point coordinates are provided, it returns an SMultiPoint.
+    """
+    lon = np.asarray(lon)
+    lat = np.asarray(lat)
+    if lon.ndim == 0:
+        p = SPoint(lon, lat)
+    else:
+        p = SMultiPoint(lon, lat)
+    return p