[Future Spherical Class] Add SPoint and SMultiPoint

pyresample/future/spherical/__init__.py

@@ -0,0 +1,20 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2021 Pyresample developers
+#
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser 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 Lesser General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Lesser General Public License along
+# with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""Future features that are backwards incompatible with current functionality."""
+
+from .point import SMultiPoint, SPoint  # noqa

pyresample/future/spherical/point.py

@@ -0,0 +1,72 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2013 - 2021 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 SPoint and SMultiPoint class."""
+import numpy as np
+
+from pyresample.spherical import SCoordinate
+
+
+class SPoint(SCoordinate):
+    """Spherical Point.
+
+    The ``lon`` and ``lat`` coordinates must be provided in radians.
+    """
+
+    def __init__(self, lon, lat):
+        lon = np.asarray(lon)
+        lat = np.asarray(lat)
+        if lon.size > 1 or lat.size > 1:
+            raise ValueError("Use SMultiPoint to define multiple points.")
+        super().__init__(lon, lat)
+
+    def to_shapely(self):
+        """Convert the SPoint to a shapely Point (in lon/lat degrees)."""
+        from shapely.geometry import Point
+        point = Point(*np.rad2deg(self.vertices[0]))
+        return point
+
+
+class SMultiPoint(SCoordinate):
+    """Create a SMultiPoint object."""
+
+    def __init__(self, lon, lat):
+        lon = np.asarray(lon)
+        lat = np.asarray(lat)
+        if lon.ndim == 0 or lat.ndim == 0:
+            raise ValueError("Use SPoint to define single points.")
+        super().__init__(lon, lat)
+
+    def __eq__(self, other):
+        """Check equality."""
+        return np.allclose(self.lon, other.lon) and np.allclose(self.lat, other.lat)
+
+    def __str__(self):
+        """Get simplified representation of lon/lat arrays in degrees."""
+        vertices = np.rad2deg(self.vertices)
+        return str(vertices)
+
+    def __repr__(self):
+        """Get simplified representation of lon/lat arrays in degrees."""
+        vertices = np.rad2deg(self.vertices)
+        return str(vertices)
+
+    def to_shapely(self):
+        """Convert the SMultiPoint to a shapely MultiPoint (in lon/lat degrees)."""
+        from shapely.geometry import MultiPoint
+        point = MultiPoint(np.rad2deg(self.vertices))
+        return point

pyresample/spherical.py

@@ -34,22 +34,121 @@ def _unwrap_radians(val, mod=np.pi):
     return (val + mod) % (2 * mod) - mod
 
 
+def _xyz_to_vertices(x, y, z):
+    """Create vertices array from x,y,z values or vectors.
+
+    If x, y, z are np.values, it create a 1x3 np.array.
+    If x, y, z are np.array with shape nx1, it creates a nx3 np.array.
+    """
+    if x.ndim == 0:
+        vertices = np.array([x, y, z])
+    else:
+        vertices = np.vstack([x, y, z]).T
+    return vertices
+
+
+def _ensure_is_array(arr):
+    """Ensure that a possible np.value input is converted to np.array."""
+    if arr.ndim == 0:
+        arr = np.array([arr])
+    return arr
+
+
+def _vincenty_matrix(lon, lat, lon_ref, lat_ref):
+    """Compute a distance matrix using Vincenty formula.
+
+    The result must be multiplied by Earth radius to obtain distance in m or km.
+    The returned distance matrix has shape (n x n_ref).
+    """
+    lon = _ensure_is_array(lon)
+    lat = _ensure_is_array(lat)
+    lon_ref = _ensure_is_array(lon_ref)
+    lat_ref = _ensure_is_array(lat_ref)
+    lon = lon[:, np.newaxis]
+    lat = lat[:, np.newaxis]
+    diff_lon = lon - lon_ref
+    num = ((np.cos(lat_ref) * np.sin(diff_lon)) ** 2 +
+           (np.cos(lat) * np.sin(lat_ref) -
+            np.sin(lat) * np.cos(lat_ref) * np.cos(diff_lon)) ** 2)
+    den = (np.sin(lat) * np.sin(lat_ref) +
+           np.cos(lat) * np.cos(lat_ref) * np.cos(diff_lon))
+    dist = np.arctan2(num ** .5, den)
+    return dist
+
+
+def _haversine_matrix(lon, lat, lon_ref, lat_ref):
+    """Compute a distance matrix using haversine formula.
+
+    The result must be multiplied by Earth radius to obtain distance in m or km.
+    The returned distance matrix has shape (n x n_ref).
+    """
+    lon = _ensure_is_array(lon)
+    lat = _ensure_is_array(lat)
+    lon_ref = _ensure_is_array(lon_ref)
+    lat_ref = _ensure_is_array(lat_ref)
+    lon = lon[:, np.newaxis]
+    lat = lat[:, np.newaxis]
+    diff_lon = lon - lon_ref  # n x n_ref matrix
+    diff_lat = lat - lat_ref  # n x n_ref matrix
+    a = np.sin(diff_lat / 2.0) ** 2.0 + np.cos(lat) * np.cos(lat_ref) * np.sin(diff_lon / 2.0) ** 2.0
+    dist = 2.0 * np.arcsin(a ** .5)  # equivalent of; 2.0 * np.arctan2(np.sqrt(a), np.sqrt(1-a))
+    return dist
+
+
+def check_lon_validity(lon):
+    """Check longitude validity."""
+    if np.any(np.isinf(lon)):
+        raise ValueError("Longitude values can not contain inf values.")
+
+
+def check_lat_validity(lat):
+    """Check latitude validity."""
+    if np.any(np.isinf(lat)):
+        raise ValueError("Latitude values can not contain inf values.")
+    if np.any(np.logical_or(lat > np.pi / 2, lat < -np.pi / 2)):
+        raise ValueError("Latitude values must range between [-pi/2, pi/2].")
+
+
+def check_lon_lat(lon, lat):
+    """Check and format lon/lat values/arrays."""
+    lon = np.asarray(lon, dtype=np.float64)
+    lat = np.asarray(lat, dtype=np.float64)
+    check_lon_validity(lon)
+    check_lat_validity(lat)
+    return lon, lat
+
+
 class SCoordinate(object):
     """Spherical coordinates.
 
-    The ``lon`` and ``lat`` coordinates should be provided in radians.
+    The ``lon`` and ``lat`` coordinates must be provided in radians.
 
     """
 
     def __init__(self, lon, lat):
-        if np.isfinite(lon):
-            self.lon = float(_unwrap_radians(lon))
-        else:
-            self.lon = float(lon)
+        lon, lat = check_lon_lat(lon, lat)
+        self.lon = _unwrap_radians(lon)
         self.lat = lat
 
+    @property
+    def vertices(self):
+        """Return point(s) vertices in a ndarray of shape [n,2]."""
+        # Single values
+        if self.lon.ndim == 0:
+            vertices = np.array([self.lon, self.lat])[None, :]
+        # Array values
+        else:
+            vertices = np.vstack((self.lon, self.lat)).T
+        return vertices
+
     def cross2cart(self, point):
-        """Compute the cross product, and convert to cartesian coordinates."""
+        """Compute the cross product, and convert to cartesian coordinates.
+
+        Note:
+        - the cross product of the same point gives a zero vector.
+        - the cross product between points lying at the equator gives a zero vector.
+        - the cross product between points lying at the poles.
+        """
         lat1 = self.lat
         lon1 = self.lon
         lat2 = point.lat
@@ -62,35 +161,48 @@ def cross2cart(self, point):
         g = np.cos(lat1)
         h = np.cos(lat2)
         i = np.sin(lon2 - lon1)
-        res = CCoordinate(np.array([-ad * c + be * f,
-                                    ad * f + be * c,
-                                    g * h * i]))
-
+        x = -ad * c + be * f
+        y = ad * f + be * c
+        z = g * h * i
+        vertices = _xyz_to_vertices(x, y, z)
+        res = CCoordinate(vertices)
         return res
 
     def to_cart(self):
         """Convert to cartesian."""
-        return CCoordinate(np.array([np.cos(self.lat) * np.cos(self.lon),
-                                     np.cos(self.lat) * np.sin(self.lon),
-                                     np.sin(self.lat)]))
+        x = np.cos(self.lat) * np.cos(self.lon)
+        y = np.cos(self.lat) * np.sin(self.lon)
+        z = np.sin(self.lat)
+        vertices = _xyz_to_vertices(x, y, z)
+        return CCoordinate(vertices)
 
     def distance(self, point):
-        """Get distance using Vincenty formula."""
-        dlambda = self.lon - point.lon
-        num = ((np.cos(point.lat) * np.sin(dlambda)) ** 2 +
-               (np.cos(self.lat) * np.sin(point.lat) -
-                np.sin(self.lat) * np.cos(point.lat) *
-                np.cos(dlambda)) ** 2)
-        den = (np.sin(self.lat) * np.sin(point.lat) +
-               np.cos(self.lat) * np.cos(point.lat) * np.cos(dlambda))
+        """Get distance using Vincenty formula.
 
-        return np.arctan2(num ** .5, den)
+        The result must be multiplied by Earth radius to obtain distance in m or km.
+        """
+        lat = self.lat
+        lon = self.lon
+        lon_ref = point.lon
+        lat_ref = point.lat
+        dist = _vincenty_matrix(lon, lat, lon_ref, lat_ref)
+        if dist.size == 1:  # single point case
+            dist = dist.item()
+        return dist
 
     def hdistance(self, point):
-        """Get distance using Haversine formula."""
-        return 2 * np.arcsin((np.sin((point.lat - self.lat) / 2.0) ** 2.0 +
-                              np.cos(point.lat) * np.cos(self.lat) *
-                              np.sin((point.lon - self.lon) / 2.0) ** 2.0) ** .5)
+        """Get distance using Haversine formula.
+
+        The result must be multiplied by Earth radius to obtain distance in m or km.
+        """
+        lat = self.lat
+        lon = self.lon
+        lon_ref = point.lon
+        lat_ref = point.lat
+        dist = _haversine_matrix(lon, lat, lon_ref, lat_ref)
+        if dist.size == 1:  # single point case
+            dist = dist.item()
+        return dist
 
     def __ne__(self, other):
         """Check inequality."""
@@ -112,6 +224,39 @@ def __iter__(self):
         """Get iterator over lon/lat pairs."""
         return zip([self.lon, self.lat]).__iter__()
 
+    def plot(self, ax=None, **plot_kwargs):
+        """Plot the point(s) using Cartopy.
+
+        Assume vertices to be in radians.
+        """
+        import matplotlib.pyplot as plt
+        try:
+            import cartopy.crs as ccrs
+        except ModuleNotFoundError:
+            raise ModuleNotFoundError("Install cartopy to plot spherical geometries.")
+
+        # Create figure if ax not provided
+        ax_not_provided = False
+        if ax is None:
+            ax_not_provided = True
+            proj_crs = ccrs.PlateCarree()
+            fig, ax = plt.subplots(subplot_kw=dict(projection=proj_crs))
+
+        # Plot Points
+        ax.scatter(x=np.rad2deg(self.vertices[:, 0]),
+                   y=np.rad2deg(self.vertices[:, 1]),
+                   **plot_kwargs)
+
+        # Beautify plot by default
+        if ax_not_provided:
+            ax.stock_img()
+            ax.coastlines()
+            gl = ax.gridlines(draw_labels=True, linestyle='--')
+            gl.xlabels_top = False
+            gl.ylabels_right = False
+
+        return ax
+
 
 class CCoordinate(object):
     """Cartesian coordinates."""
@@ -123,14 +268,28 @@ def norm(self):
         """Get Euclidean norm of the vector."""
         return np.sqrt(np.einsum('...i, ...i', self.cart, self.cart))
 
-    def normalize(self):
-        """Normalize the vector."""
-        self.cart /= np.sqrt(np.einsum('...i, ...i', self.cart, self.cart))
+    def normalize(self, inplace=False):
+        """Normalize the vector.
 
-        return self
+        If self.cart == [0,0,0], norm=0, and cart becomes [nan, nan, nan]:
+        Note that self.cart == [0,0,0] can occurs when computing:
+        - the cross product of the same point.
+        - the cross product between points lying at the equator.
+        - the cross product between points lying at the poles.
+        """
+        norm = self.norm()
+        norm = norm[..., np.newaxis]  # enable vectorization
+        if inplace:
+            self.cart /= norm
+            return None
+        cart = self.cart / norm
+        return CCoordinate(cart)
 
     def cross(self, point):
-        """Get cross product with another vector."""
+        """Get cross product with another vector.
+
+        The cross product of the same vector gives a zero vector.
+        """
         return CCoordinate(np.cross(self.cart, point.cart))
 
     def dot(self, point):
@@ -176,9 +335,11 @@ def __rmul__(self, other):
         return self.__mul__(other)
 
     def to_spherical(self):
-        """Convert to Spherical coordinate object."""
-        return SCoordinate(np.arctan2(self.cart[1], self.cart[0]),
-                           np.arcsin(self.cart[2]))
+        """Convert to SPoint/SMultiPoint object."""
+        # TODO: this in future should point to SPoint or SMultiPoint
+        lon = np.arctan2(self.cart[..., 1], self.cart[..., 0])
+        lat = np.arcsin(self.cart[..., 2])
+        return SCoordinate(lon, lat)
 
 
 class Arc(object):
@@ -236,6 +397,7 @@ def angle(self, other_arc):
         ub_ = a__.cross2cart(c__)
 
         val = ua_.dot(ub_) / (ua_.norm() * ub_.norm())
+
         if abs(val - 1) < EPSILON:
             angle = 0
         elif abs(val + 1) < EPSILON: