[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,151 @@
+#!/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
+
+
+def _plot(vertices, ax=None, **plot_kwargs):
+    """Plot the SPoint/SMultiPoint using Cartopy.
+
+    Assume vertices to 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(vertices[:, 0]),
+               y=np.rad2deg(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
+
+
+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].")
+
+
+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.")
+        check_lon_validity(lon)
+        check_lat_validity(lat)
+
+        super().__init__(lon, lat)
+
+    @property
+    def vertices(self):
+        """Return SPoint vertices in a ndarray of shape [1,2]."""
+        return np.array([self.lon, self.lat])[None, :]
+
+    def plot(self, ax=None, **plot_kwargs):
+        """Plot the SPoint using Cartopy."""
+        ax = _plot(self.vertices, ax=ax, **plot_kwargs)
+        return ax
+
+    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 SPoint or SMultiPoint object."""
+
+    def __new__(cls, lon, lat):
+        """Create SPoint or SMultiPoint object."""
+        lon = np.asarray(lon)
+        lat = np.asarray(lat)
+
+        # If a single point, define SPoint
+        if lon.ndim == 0:
+            return SPoint(lon, lat)  # TODO: TO BE DEFINED
+        # Otherwise a SMultiPoint
+        else:
+            return super().__new__(cls)
+
+    def __init__(self, lon, lat):
+        super().__init__(lon, lat)
+
+    @property
+    def vertices(self):
+        """Return SMultiPoint vertices in a ndarray of shape [n,2]."""
+        return np.vstack((self.lon, self.lat)).T
+
+    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(np.vstack((self.lon, self.lat)).T)
+        return str(vertices)
+
+    def __repr__(self):
+        """Get simplified representation of lon/lat arrays in degrees."""
+        vertices = np.rad2deg(np.vstack((self.lon, self.lat)).T)
+        return str(vertices)
+
+    def plot(self, ax=None, **plot_kwargs):
+        """Plot the SMultiPoint using Cartopy."""
+        ax = _plot(self.vertices, ax=ax, **plot_kwargs)
+        return ax
+
+    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,6 +34,62 @@ def _unwrap_radians(val, mod=np.pi):
     return (val + mod) % (2 * mod) - mod
 
 
+def _xyz_to_vertices(x, y, z):
+    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
+
+
 class SCoordinate(object):
     """Spherical coordinates.
 
@@ -42,14 +98,16 @@ class SCoordinate(object):
     """
 
     def __init__(self, lon, lat):
-        if np.isfinite(lon):
-            self.lon = float(_unwrap_radians(lon))
-        else:
-            self.lon = float(lon)
+        lon = np.asarray(lon, dtype=np.float64)
+        lat = np.asarray(lat, dtype=np.float64)
+        self.lon = _unwrap_radians(lon)
         self.lat = lat
 
     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.
+        """
         lat1 = self.lat
         lon1 = self.lon
         lat2 = point.lat
@@ -62,35 +120,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."""
@@ -125,12 +196,18 @@ def norm(self):
 
     def normalize(self):
         """Normalize the vector."""
-        self.cart /= np.sqrt(np.einsum('...i, ...i', self.cart, self.cart))
-
-        return self
+        # Note: if self.cart == [0,0,0], norm=0 and cart becomes [nan, nan, nan]
+        norm = self.norm()
+        cart = self.cart / norm[..., np.newaxis]  # enable vectorization
+        # mask_norm_0 = norm == 0
+        # cart[mask_norm_0, ...] = 0 # This remove nan propagation
+        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 +253,16 @@ 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 SCoordinate object."""
+        if self.cart.ndim == 1:
+            lon = np.arctan2(self.cart[1], self.cart[0])
+            lat = np.arcsin(self.cart[2])
+        else:
+            lon = np.arctan2(self.cart[:, 1], self.cart[:, 0])
+            lat = np.arcsin(self.cart[:, 2])
+        # TODO: in future this should point to SPoint or SMultiPoint
+        # - This is not used within the file.
+        return SCoordinate(lon, lat)
 
 
 class Arc(object):
@@ -236,6 +320,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:

pyresample/test/test_sgeom/__init__.py

@@ -0,0 +1,18 @@
+#!/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/>.
+"""Test for spherical geometries."""