From 45e5b184a72cf4de62d355a6c4e15349d2dd550c Mon Sep 17 00:00:00 2001
From: ghiggi <gionata.ghiggi@gmail.com>
Date: Tue, 8 Nov 2022 12:07:34 +0100
Subject: [PATCH 01/20] Add SPoint and SMultiPoint

---
 pyresample/future/spherical/__init__.py  |  20 +++
 pyresample/future/spherical/point.py     | 151 +++++++++++++++++++++
 pyresample/spherical.py                  | 149 ++++++++++++++++-----
 pyresample/test/test_sgeom/__init__.py   |  18 +++
 pyresample/test/test_sgeom/test_point.py | 159 +++++++++++++++++++++++
 pyresample/test/test_spherical.py        |   6 -
 6 files changed, 465 insertions(+), 38 deletions(-)
 create mode 100644 pyresample/future/spherical/__init__.py
 create mode 100644 pyresample/future/spherical/point.py
 create mode 100644 pyresample/test/test_sgeom/__init__.py
 create mode 100644 pyresample/test/test_sgeom/test_point.py

diff --git a/pyresample/future/spherical/__init__.py b/pyresample/future/spherical/__init__.py
new file mode 100644
index 000000000..817b63fce
--- /dev/null
+++ b/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
diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py
new file mode 100644
index 000000000..b036da525
--- /dev/null
+++ b/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
diff --git a/pyresample/spherical.py b/pyresample/spherical.py
index f6676c494..08c8e5e19 100644
--- a/pyresample/spherical.py
+++ b/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:
diff --git a/pyresample/test/test_sgeom/__init__.py b/pyresample/test/test_sgeom/__init__.py
new file mode 100644
index 000000000..b496d5604
--- /dev/null
+++ b/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."""
diff --git a/pyresample/test/test_sgeom/test_point.py b/pyresample/test/test_sgeom/test_point.py
new file mode 100644
index 000000000..e9318a28e
--- /dev/null
+++ b/pyresample/test/test_sgeom/test_point.py
@@ -0,0 +1,159 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# Copyright (c) 2013-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/>.
+"""Test cases for SPoint and SMultiPoint."""
+import unittest
+
+import numpy as np
+import pytest
+
+from pyresample.future.spherical import SMultiPoint, SPoint
+
+
+class TestSPoint(unittest.TestCase):
+    """Test SPoint."""
+
+    def test_latitude_validity(self):
+        # Test latitude outside range
+        lon = 0
+        lat = np.pi
+        with pytest.raises(ValueError):
+            SPoint(lon, lat)
+        # Test inf
+        lon = 0
+        lat = np.inf
+        with pytest.raises(ValueError):
+            SPoint(lon, lat)
+
+    def test_longitude_validity(self):
+        # Test inf
+        lon = np.inf
+        lat = 0
+        with pytest.raises(ValueError):
+            SPoint(lon, lat)
+
+    def test_raise_error_if_multi_point(self):
+        lons = np.array([0, np.pi])
+        lats = np.array([-np.pi / 2, np.pi / 2])
+        with pytest.raises(ValueError):
+            SPoint(lons, lats)
+
+    def test_to_shapely(self):
+        """Test conversion to shapely."""
+        from shapely.geometry import Point
+        lon = 0.0
+        lat = np.pi / 2
+        spherical_point = SPoint(lon, lat)
+        shapely_point = Point(0.0, 90.0)
+        self.assertTrue(shapely_point.equals_exact(spherical_point.to_shapely(), tolerance=1e-10))
+
+
+class TestSMultiPoint(unittest.TestCase):
+    """Test SMultiPoint."""
+
+    def test_vertices(self):
+        """Test vertices property."""
+        lons = np.array([0, np.pi])
+        lats = np.array([-np.pi / 2, np.pi / 2])
+        p = SMultiPoint(lons, lats)
+        res = np.array([[0., -1.57079633],
+                        [-3.14159265, 1.57079633]])
+        self.assertTrue(np.allclose(p.vertices, res))
+
+    def test_distance(self):
+        """Test Vincenty formula."""
+        lons = np.array([0, np.pi])
+        lats = np.array([-np.pi / 2, np.pi / 2])
+        p1 = SMultiPoint(lons, lats)
+        lons = np.array([0, np.pi / 2, np.pi])
+        lats = np.array([-np.pi / 2, 0, np.pi / 2])
+        p2 = SMultiPoint(lons, lats)
+        d12 = p1.distance(p2)
+        d21 = p2.distance(p1)
+        self.assertEqual(d12.shape, (2, 3))
+        self.assertEqual(d21.shape, (3, 2))
+        res = np.array([[0., 1.57079633, 3.14159265],
+                        [3.14159265, 1.57079633, 0.]])
+        self.assertTrue(np.allclose(d12, res))
+        # Special case with 1 point
+        p1 = SMultiPoint(lons[0], lats[0])
+        p2 = SMultiPoint(lons[0], lats[0])
+        d12 = p1.distance(p2)
+        assert isinstance(d12, float)
+
+    def test_hdistance(self):
+        """Test Haversine formula."""
+        lons = np.array([0, np.pi])
+        lats = np.array([-np.pi / 2, np.pi / 2])
+        p1 = SMultiPoint(lons, lats)
+        lons = np.array([0, np.pi / 2, np.pi])
+        lats = np.array([-np.pi / 2, 0, np.pi / 2])
+        p2 = SMultiPoint(lons, lats)
+        d12 = p1.hdistance(p2)
+        d21 = p2.hdistance(p1)
+        self.assertEqual(d12.shape, (2, 3))
+        self.assertEqual(d21.shape, (3, 2))
+        res = np.array([[0., 1.57079633, 3.14159265],
+                        [3.14159265, 1.57079633, 0.]])
+        self.assertTrue(np.allclose(d12, res))
+
+    def test_eq(self):
+        """Check the equality."""
+        lons = [0, np.pi]
+        lats = [-np.pi / 2, np.pi / 2]
+        p = SMultiPoint(lons, lats)
+        p1 = SMultiPoint(lons, lats)
+        self.assertTrue(p == p1)
+
+    def test_eq_antimeridian(self):
+        """Check the equality with longitudes at -180/180 degrees."""
+        lons = [np.pi, np.pi]
+        lons1 = [-np.pi, -np.pi]
+        lats = [-np.pi / 2, np.pi / 2]
+        p = SMultiPoint(lons, lats)
+        p1 = SMultiPoint(lons1, lats)
+        self.assertTrue(p == p1)
+
+    def test_neq(self):
+        """Check the equality."""
+        lons = np.array([0, np.pi])
+        lats = [-np.pi / 2, np.pi / 2]
+        p = SMultiPoint(lons, lats)
+        p1 = SMultiPoint(lons + 0.1, lats)
+        self.assertTrue(p != p1)
+
+    def test_str(self):
+        """Check the string representation."""
+        lons = [0, np.pi]
+        lats = [-np.pi / 2, np.pi / 2]
+        p = SMultiPoint(lons, lats)
+        self.assertEqual(str(p), '[[   0.  -90.]\n [-180.   90.]]')
+
+    def test_repr(self):
+        """Check the representation."""
+        lons = [0, np.pi]
+        lats = [-np.pi / 2, np.pi / 2]
+        p = SMultiPoint(lons, lats)
+        self.assertEqual(repr(p), '[[   0.  -90.]\n [-180.   90.]]')
+
+    def test_to_shapely(self):
+        """Test conversion to shapely."""
+        from shapely.geometry import MultiPoint
+        lons = np.array([0.0, np.pi])
+        lats = np.array([-np.pi / 2, np.pi / 2])
+        spherical_multipoint = SMultiPoint(lons, lats)
+        shapely_multipoint = MultiPoint([(0.0, -90.0), (-180.0, 90.0)])
+        self.assertTrue(shapely_multipoint.equals_exact(spherical_multipoint.to_shapely(), tolerance=1e-10))
diff --git a/pyresample/test/test_spherical.py b/pyresample/test/test_spherical.py
index 3ef7594ba..03a96cf1c 100644
--- a/pyresample/test/test_spherical.py
+++ b/pyresample/test/test_spherical.py
@@ -146,12 +146,6 @@ def test_equality_at_antimeridian(self):
         point_end = SCoordinate(np.pi, 0)
         assert point_start == point_end
 
-    def test_equality_of_infinites(self):
-        """Test that infinite coordinates are equal."""
-        coord1 = SCoordinate(np.inf, np.inf)
-        coord2 = SCoordinate(np.inf, np.inf)
-        assert coord1 == coord2
-
 
 class TestCCoordinate(unittest.TestCase):
     """Test SCoordinates."""

From 881ed4cb24a9e2fd40462894a896f025df38b6dd Mon Sep 17 00:00:00 2001
From: Gionata Ghiggi <gionata.ghiggi@gmail.com>
Date: Tue, 15 Nov 2022 07:45:39 +0100
Subject: [PATCH 02/20] Update __str__ method

Co-authored-by: Martin Raspaud <martin.raspaud@smhi.se>
---
 pyresample/future/spherical/point.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py
index b036da525..ed0db9a6e 100644
--- a/pyresample/future/spherical/point.py
+++ b/pyresample/future/spherical/point.py
@@ -131,7 +131,7 @@ def __eq__(self, other):
 
     def __str__(self):
         """Get simplified representation of lon/lat arrays in degrees."""
-        vertices = np.rad2deg(np.vstack((self.lon, self.lat)).T)
+        vertices = np.rad2deg(self.vertices)
         return str(vertices)
 
     def __repr__(self):

From 2208589a693fd5951dbcb67d9fa9540fb5c28eb8 Mon Sep 17 00:00:00 2001
From: Gionata Ghiggi <gionata.ghiggi@gmail.com>
Date: Tue, 15 Nov 2022 07:46:10 +0100
Subject: [PATCH 03/20] Update __repr__ method

Co-authored-by: Martin Raspaud <martin.raspaud@smhi.se>
---
 pyresample/future/spherical/point.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py
index ed0db9a6e..14f7f1888 100644
--- a/pyresample/future/spherical/point.py
+++ b/pyresample/future/spherical/point.py
@@ -136,7 +136,7 @@ def __str__(self):
 
     def __repr__(self):
         """Get simplified representation of lon/lat arrays in degrees."""
-        vertices = np.rad2deg(np.vstack((self.lon, self.lat)).T)
+        vertices = np.rad2deg(self.vertices)
         return str(vertices)
 
     def plot(self, ax=None, **plot_kwargs):

From fd27e143051ea9515568e71762239e3671bb3f2e Mon Sep 17 00:00:00 2001
From: ghiggi <gionata.ghiggi@gmail.com>
Date: Tue, 15 Nov 2022 11:02:41 +0100
Subject: [PATCH 04/20] Improve docs and small refactor

---
 pyresample/future/spherical/point.py     |  87 +----------------
 pyresample/spherical.py                  | 119 +++++++++++++++++++----
 pyresample/test/test_sgeom/test_point.py |  15 ++-
 pyresample/test/test_spherical.py        |  25 +++++
 4 files changed, 141 insertions(+), 105 deletions(-)

diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py
index 14f7f1888..4b42443c3 100644
--- a/pyresample/future/spherical/point.py
+++ b/pyresample/future/spherical/point.py
@@ -21,54 +21,6 @@
 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.
 
@@ -80,21 +32,8 @@ def __init__(self, lon, lat):
         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
@@ -103,28 +42,15 @@ def to_shapely(self):
 
 
 class SMultiPoint(SCoordinate):
-    """Create SPoint or SMultiPoint object."""
+    """Create a SMultiPoint object."""
 
-    def __new__(cls, lon, lat):
-        """Create SPoint or SMultiPoint object."""
+    def __init__(self, lon, lat):
         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):
+        if lon.ndim == 0 or lat.ndim == 0:
+            raise ValueError("Use SPoint to define single points.")
         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)
@@ -139,11 +65,6 @@ def __repr__(self):
         vertices = np.rad2deg(self.vertices)
         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
diff --git a/pyresample/spherical.py b/pyresample/spherical.py
index 08c8e5e19..ff7b85e8e 100644
--- a/pyresample/spherical.py
+++ b/pyresample/spherical.py
@@ -35,6 +35,11 @@ def _unwrap_radians(val, mod=np.pi):
 
 
 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:
@@ -90,23 +95,59 @@ def _haversine_matrix(lon, lat, lon_ref, lat_ref):
     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):
-        lon = np.asarray(lon, dtype=np.float64)
-        lat = np.asarray(lat, dtype=np.float64)
+        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.
 
-        Note: the cross product of the same point gives a zero vector.
+        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
@@ -183,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."""
@@ -194,14 +268,24 @@ 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."""
-        # Note: if self.cart == [0,0,0], norm=0 and cart becomes [nan, nan, nan]
+    def normalize(self, inplace=False):
+        """Normalize the vector.
+
+        Notes:
+        - If self.cart == [0,0,0], norm=0, and cart becomes [nan, nan, nan]:
+        - 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()
-        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)
+        norm = norm[..., np.newaxis]  # enable vectorization
+        if inplace:
+            self.cart /= norm
+            return None
+        else:
+            cart = self.cart / norm
+            return CCoordinate(cart)
 
     def cross(self, point):
         """Get cross product with another vector.
@@ -253,15 +337,10 @@ def __rmul__(self, other):
         return self.__mul__(other)
 
     def to_spherical(self):
-        """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.
+        """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)
 
 
diff --git a/pyresample/test/test_sgeom/test_point.py b/pyresample/test/test_sgeom/test_point.py
index e9318a28e..8b07fcf33 100644
--- a/pyresample/test/test_sgeom/test_point.py
+++ b/pyresample/test/test_sgeom/test_point.py
@@ -64,6 +64,17 @@ def test_to_shapely(self):
 class TestSMultiPoint(unittest.TestCase):
     """Test SMultiPoint."""
 
+    def test_single_point(self):
+        """Test behaviour when providing single lon,lat values."""
+        # Single values must raise error
+        with pytest.raises(ValueError):
+            SMultiPoint(2, 1)
+        # Array values must not raise error
+        p = SMultiPoint([2], [1])
+        assert p.lon.shape == (1,)
+        assert p.lat.shape == (1,)
+        assert p.vertices.shape == (1, 2)
+
     def test_vertices(self):
         """Test vertices property."""
         lons = np.array([0, np.pi])
@@ -89,8 +100,8 @@ def test_distance(self):
                         [3.14159265, 1.57079633, 0.]])
         self.assertTrue(np.allclose(d12, res))
         # Special case with 1 point
-        p1 = SMultiPoint(lons[0], lats[0])
-        p2 = SMultiPoint(lons[0], lats[0])
+        p1 = SMultiPoint(lons[[0]], lats[[0]])
+        p2 = SMultiPoint(lons[[0]], lats[[0]])
         d12 = p1.distance(p2)
         assert isinstance(d12, float)
 
diff --git a/pyresample/test/test_spherical.py b/pyresample/test/test_spherical.py
index 03a96cf1c..70ed25bf4 100644
--- a/pyresample/test/test_spherical.py
+++ b/pyresample/test/test_spherical.py
@@ -130,6 +130,31 @@ def test_hdistance(self):
         d = SCoordinate(0, 0).hdistance(SCoordinate(1, 1))
         np.testing.assert_allclose(d, 1.2745557823062943)
 
+    def test_crosscart(self):
+        """Test cross product and conversion to cartesian."""
+        # Test crossproduct between poles
+        p = SCoordinate(np.deg2rad(50), np.deg2rad(90.0))
+        p1 = SCoordinate(np.deg2rad(50), np.deg2rad(-90.0))
+        cp = p.cross2cart(p1)
+        assert np.allclose(cp.cart, [0, 0, 0])
+
+        # Test crossproduct between points along the equator
+        p = SCoordinate(np.deg2rad(-180.0), np.deg2rad(0.0))
+        p1 = SCoordinate(np.deg2rad(0.0), np.deg2rad(0.0))
+        cp = p.cross2cart(p1)
+        assert np.allclose(cp.cart, [0, 0, 0])
+
+        # Test crossproduct between same points
+        p = SCoordinate(np.deg2rad(50), np.deg2rad(20.0))
+        p1 = SCoordinate(np.deg2rad(50), np.deg2rad(20.0))
+        cp = p.cross2cart(p1)
+        assert np.allclose(cp.cart, [0, 0, 0])
+
+        p = SCoordinate(np.deg2rad(-180), np.deg2rad(90.0))
+        p1 = SCoordinate(np.deg2rad(180), np.deg2rad(90.0))
+        p.cross2cart(p1)
+        assert np.allclose(cp.cart, [0, 0, 0])
+
     def test_str(self):
         """Check the string representation."""
         d = SCoordinate(0, 0)

From d2581ae434bc9db8902a8e82d99d56e0ff37692f Mon Sep 17 00:00:00 2001
From: ghiggi <gionata.ghiggi@gmail.com>
Date: Tue, 15 Nov 2022 11:12:44 +0100
Subject: [PATCH 05/20] Fix sphinx problem

---
 pyresample/spherical.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/pyresample/spherical.py b/pyresample/spherical.py
index ff7b85e8e..a2d1df758 100644
--- a/pyresample/spherical.py
+++ b/pyresample/spherical.py
@@ -274,9 +274,9 @@ def normalize(self, inplace=False):
         Notes:
         - If self.cart == [0,0,0], norm=0, and cart becomes [nan, nan, nan]:
         - 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.
+            - 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

From 9f9599d3e70ba828703bf0e66e52af4ae59ef4a1 Mon Sep 17 00:00:00 2001
From: ghiggi <gionata.ghiggi@gmail.com>
Date: Tue, 15 Nov 2022 11:22:22 +0100
Subject: [PATCH 06/20] Fix sphinx problem

---
 pyresample/spherical.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/pyresample/spherical.py b/pyresample/spherical.py
index a2d1df758..0d446dfa7 100644
--- a/pyresample/spherical.py
+++ b/pyresample/spherical.py
@@ -274,9 +274,9 @@ def normalize(self, inplace=False):
         Notes:
         - If self.cart == [0,0,0], norm=0, and cart becomes [nan, nan, nan]:
         - 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.
+            * 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

From 73d9a1b5718cd2f0b64e85c825b9ec9f01c024b9 Mon Sep 17 00:00:00 2001
From: ghiggi <gionata.ghiggi@gmail.com>
Date: Tue, 15 Nov 2022 11:38:34 +0100
Subject: [PATCH 07/20] Hopefully finally maybe fix sphinx building

---
 pyresample/spherical.py | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/pyresample/spherical.py b/pyresample/spherical.py
index 0d446dfa7..18a48cc2e 100644
--- a/pyresample/spherical.py
+++ b/pyresample/spherical.py
@@ -271,12 +271,11 @@ def norm(self):
     def normalize(self, inplace=False):
         """Normalize the vector.
 
-        Notes:
-        - If self.cart == [0,0,0], norm=0, and cart becomes [nan, nan, nan]:
-        - 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.
+        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

From 6d65d222853196ab2bc100b4590dc2320cb20b2e Mon Sep 17 00:00:00 2001
From: Gionata Ghiggi <gionata.ghiggi@gmail.com>
Date: Tue, 15 Nov 2022 14:07:55 +0100
Subject: [PATCH 08/20] Remove else if else in CCoordinate.normalize

Co-authored-by: Martin Raspaud <martin.raspaud@smhi.se>
---
 pyresample/spherical.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/pyresample/spherical.py b/pyresample/spherical.py
index 18a48cc2e..5b5d42ea6 100644
--- a/pyresample/spherical.py
+++ b/pyresample/spherical.py
@@ -282,9 +282,8 @@ def normalize(self, inplace=False):
         if inplace:
             self.cart /= norm
             return None
-        else:
-            cart = self.cart / norm
-            return CCoordinate(cart)
+        cart = self.cart / norm
+        return CCoordinate(cart)
 
     def cross(self, point):
         """Get cross product with another vector.

From fbe8c053308744d332255d2bd6c39693884f54c2 Mon Sep 17 00:00:00 2001
From: Gionata Ghiggi <gionata.ghiggi@gmail.com>
Date: Wed, 16 Nov 2022 07:58:45 +0100
Subject: [PATCH 09/20] Update copyright

Co-authored-by: David Hoese <david.hoese@ssec.wisc.edu>
---
 pyresample/future/spherical/point.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py
index 4b42443c3..570f617bc 100644
--- a/pyresample/future/spherical/point.py
+++ b/pyresample/future/spherical/point.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
-# Copyright (c) 2013 - 2021 Pyresample developers
+# 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

From 51b259a9b6af825211c708a696583b4ae4d5beac Mon Sep 17 00:00:00 2001
From: Gionata Ghiggi <gionata.ghiggi@gmail.com>
Date: Wed, 16 Nov 2022 08:00:04 +0100
Subject: [PATCH 10/20] Update file description

---
 pyresample/future/spherical/point.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py
index 570f617bc..60643e7ca 100644
--- a/pyresample/future/spherical/point.py
+++ b/pyresample/future/spherical/point.py
@@ -15,7 +15,7 @@
 #
 # 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."""
+"""Define single and multiple points on the sphere through SPoint and SMultiPoint classes."""
 import numpy as np
 
 from pyresample.spherical import SCoordinate

From 86ff4fb62d9ec0a2611cf6fc60b1f88a6efd37d3 Mon Sep 17 00:00:00 2001
From: Gionata Ghiggi <gionata.ghiggi@gmail.com>
Date: Wed, 16 Nov 2022 08:00:49 +0100
Subject: [PATCH 11/20] Update SPoint docstring

---
 pyresample/future/spherical/point.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py
index 60643e7ca..e1d86d7bf 100644
--- a/pyresample/future/spherical/point.py
+++ b/pyresample/future/spherical/point.py
@@ -22,7 +22,7 @@
 
 
 class SPoint(SCoordinate):
-    """Spherical Point.
+    """Define a point on the sphere.
 
     The ``lon`` and ``lat`` coordinates must be provided in radians.
     """

From 5938ceaf17ef589eea32eadaff4c3bed6186c341 Mon Sep 17 00:00:00 2001
From: Gionata Ghiggi <gionata.ghiggi@gmail.com>
Date: Wed, 16 Nov 2022 08:02:06 +0100
Subject: [PATCH 12/20] Update SMultiPoint docstring

---
 pyresample/future/spherical/point.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py
index e1d86d7bf..cbd6fd111 100644
--- a/pyresample/future/spherical/point.py
+++ b/pyresample/future/spherical/point.py
@@ -42,7 +42,7 @@ def to_shapely(self):
 
 
 class SMultiPoint(SCoordinate):
-    """Create a SMultiPoint object."""
+    """Object representing multiple points on a sphere."""
 
     def __init__(self, lon, lat):
         lon = np.asarray(lon)

From fef04b4c0f7e04e07ce84b8ed099d554b79beb35 Mon Sep 17 00:00:00 2001
From: Gionata Ghiggi <gionata.ghiggi@gmail.com>
Date: Wed, 16 Nov 2022 08:02:47 +0100
Subject: [PATCH 13/20] Update SPoint docstring

---
 pyresample/future/spherical/point.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py
index cbd6fd111..3272277be 100644
--- a/pyresample/future/spherical/point.py
+++ b/pyresample/future/spherical/point.py
@@ -22,7 +22,7 @@
 
 
 class SPoint(SCoordinate):
-    """Define a point on the sphere.
+    """Object representing a single point on a sphere.
 
     The ``lon`` and ``lat`` coordinates must be provided in radians.
     """

From ae18d6d93d7af3d43c0d8e57b64a191314854cfd Mon Sep 17 00:00:00 2001
From: Gionata Ghiggi <gionata.ghiggi@gmail.com>
Date: Wed, 16 Nov 2022 08:04:31 +0100
Subject: [PATCH 14/20] Update _xyz_to_vertices docstring

---
 pyresample/spherical.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyresample/spherical.py b/pyresample/spherical.py
index 5b5d42ea6..efe7c7683 100644
--- a/pyresample/spherical.py
+++ b/pyresample/spherical.py
@@ -37,7 +37,7 @@ def _unwrap_radians(val, mod=np.pi):
 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 scalar arrays, it creates a 1x3 np.array.
     If x, y, z are np.array with shape nx1, it creates a nx3 np.array.
     """
     if x.ndim == 0:

From f52ae78b22480f7630d8fa32213973afc20bf971 Mon Sep 17 00:00:00 2001
From: Gionata Ghiggi <gionata.ghiggi@gmail.com>
Date: Wed, 16 Nov 2022 08:06:06 +0100
Subject: [PATCH 15/20] Replace [None,:] with [np.newaxis, :]

Co-authored-by: David Hoese <david.hoese@ssec.wisc.edu>
---
 pyresample/spherical.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyresample/spherical.py b/pyresample/spherical.py
index efe7c7683..9ad08ac65 100644
--- a/pyresample/spherical.py
+++ b/pyresample/spherical.py
@@ -135,7 +135,7 @@ 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, :]
+            vertices = np.array([self.lon, self.lat])[np.newaxis, :]
         # Array values
         else:
             vertices = np.vstack((self.lon, self.lat)).T

From b2cd0ba2b53c7418d1a7efe9ec0a19b5bffaf612 Mon Sep 17 00:00:00 2001
From: Gionata Ghiggi <gionata.ghiggi@gmail.com>
Date: Wed, 16 Nov 2022 08:07:11 +0100
Subject: [PATCH 16/20] Update information for installing cartopy

Co-authored-by: David Hoese <david.hoese@ssec.wisc.edu>
---
 pyresample/spherical.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyresample/spherical.py b/pyresample/spherical.py
index 9ad08ac65..c77b7be4b 100644
--- a/pyresample/spherical.py
+++ b/pyresample/spherical.py
@@ -233,7 +233,7 @@ def plot(self, ax=None, **plot_kwargs):
         try:
             import cartopy.crs as ccrs
         except ModuleNotFoundError:
-            raise ModuleNotFoundError("Install cartopy to plot spherical geometries.")
+            raise ModuleNotFoundError("Install cartopy to plot spherical geometries. For example, 'pip install cartopy'.")
 
         # Create figure if ax not provided
         ax_not_provided = False

From 0ebbe4944dabc7758856ae293a7263a952e2b04b Mon Sep 17 00:00:00 2001
From: ghiggi <gionata.ghiggi@gmail.com>
Date: Wed, 16 Nov 2022 09:16:07 +0100
Subject: [PATCH 17/20] Address review remarks

---
 pyresample/future/spherical/point.py     | 18 ++++---
 pyresample/spherical.py                  | 55 ++++++++++++++--------
 pyresample/test/test_sgeom/test_point.py | 60 ++++++++++++++++++++----
 3 files changed, 97 insertions(+), 36 deletions(-)

diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py
index 3272277be..985c888ac 100644
--- a/pyresample/future/spherical/point.py
+++ b/pyresample/future/spherical/point.py
@@ -34,10 +34,18 @@ def __init__(self, lon, lat):
             raise ValueError("Use SMultiPoint to define multiple points.")
         super().__init__(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)))
+
     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]))
+        point = Point(*self.vertices_in_degrees[0])
         return point
 
 
@@ -57,16 +65,14 @@ def __eq__(self, other):
 
     def __str__(self):
         """Get simplified representation of lon/lat arrays in degrees."""
-        vertices = np.rad2deg(self.vertices)
-        return str(vertices)
+        return str(self.vertices)
 
     def __repr__(self):
         """Get simplified representation of lon/lat arrays in degrees."""
-        vertices = np.rad2deg(self.vertices)
-        return str(vertices)
+        return str(self.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))
+        point = MultiPoint(self.vertices_in_degrees)
         return point
diff --git a/pyresample/spherical.py b/pyresample/spherical.py
index c77b7be4b..9dca27904 100644
--- a/pyresample/spherical.py
+++ b/pyresample/spherical.py
@@ -50,14 +50,15 @@ def _xyz_to_vertices(x, y, z):
 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])
+        arr = np.asarray([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 lon/lat inputs must be provided in radians !
+    The output must be multiplied by the Earth radius to obtain the distance in m or km.
     The returned distance matrix has shape (n x n_ref).
     """
     lon = _ensure_is_array(lon)
@@ -79,7 +80,8 @@ def _vincenty_matrix(lon, lat, lon_ref, lat_ref):
 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 lon/lat inputs must be provided in radians !
+    The output must be multiplied by the Earth radius to obtain the distance in m or km.
     The returned distance matrix has shape (n x n_ref).
     """
     lon = _ensure_is_array(lon)
@@ -95,13 +97,13 @@ def _haversine_matrix(lon, lat, lon_ref, lat_ref):
     return dist
 
 
-def check_lon_validity(lon):
+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):
+def _check_lat_validity(lat):
     """Check latitude validity."""
     if np.any(np.isinf(lat)):
         raise ValueError("Latitude values can not contain inf values.")
@@ -109,12 +111,12 @@ def check_lat_validity(lat):
         raise ValueError("Latitude values must range between [-pi/2, pi/2].")
 
 
-def check_lon_lat(lon, lat):
+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)
+    _check_lon_validity(lon)
+    _check_lat_validity(lat)
     return lon, lat
 
 
@@ -126,13 +128,13 @@ class SCoordinate(object):
     """
 
     def __init__(self, lon, lat):
-        lon, lat = check_lon_lat(lon, lat)
+        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]."""
+        """Return point(s) radians vertices in a ndarray of shape [n,2]."""
         # Single values
         if self.lon.ndim == 0:
             vertices = np.array([self.lon, self.lat])[np.newaxis, :]
@@ -141,6 +143,11 @@ def vertices(self):
             vertices = np.vstack((self.lon, self.lat)).T
         return vertices
 
+    @property
+    def vertices_in_degrees(self):
+        """Return point(s) degrees vertices in a ndarray of shape [n,2]."""
+        return np.rad2deg(self.vertices)
+
     def cross2cart(self, point):
         """Compute the cross product, and convert to cartesian coordinates.
 
@@ -224,7 +231,12 @@ def __iter__(self):
         """Get iterator over lon/lat pairs."""
         return zip([self.lon, self.lat]).__iter__()
 
-    def plot(self, ax=None, **plot_kwargs):
+    def plot(self, ax=None,
+             projection_crs=None,
+             add_coastlines=True,
+             add_gridlines=True,
+             add_background=True,
+             **plot_kwargs):
         """Plot the point(s) using Cartopy.
 
         Assume vertices to be in radians.
@@ -233,24 +245,27 @@ def plot(self, ax=None, **plot_kwargs):
         try:
             import cartopy.crs as ccrs
         except ModuleNotFoundError:
-            raise ModuleNotFoundError("Install cartopy to plot spherical geometries. For example, 'pip install cartopy'.")
+            raise ModuleNotFoundError(
+                "Install cartopy to plot spherical geometries. For example, 'pip install cartopy'.")
 
         # 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))
+            if projection_crs is None:
+                projection_crs = ccrs.PlateCarree()
+            fig, ax = plt.subplots(subplot_kw=dict(projection=projection_crs))
 
         # Plot Points
-        ax.scatter(x=np.rad2deg(self.vertices[:, 0]),
-                   y=np.rad2deg(self.vertices[:, 1]),
+        vertices = self.vertices_in_degrees
+        ax.scatter(x=vertices[:, 0],
+                   y=vertices[:, 1],
                    **plot_kwargs)
 
-        # Beautify plot by default
-        if ax_not_provided:
+        # Beautify plots
+        if add_background:
             ax.stock_img()
+        if add_coastlines:
             ax.coastlines()
+        if add_gridlines():
             gl = ax.gridlines(draw_labels=True, linestyle='--')
             gl.xlabels_top = False
             gl.ylabels_right = False
diff --git a/pyresample/test/test_sgeom/test_point.py b/pyresample/test/test_sgeom/test_point.py
index 8b07fcf33..6722807e3 100644
--- a/pyresample/test/test_sgeom/test_point.py
+++ b/pyresample/test/test_sgeom/test_point.py
@@ -27,6 +27,7 @@ class TestSPoint(unittest.TestCase):
     """Test SPoint."""
 
     def test_latitude_validity(self):
+        """Check SPoint raises error if providing bad latitude."""
         # Test latitude outside range
         lon = 0
         lat = np.pi
@@ -39,18 +40,46 @@ def test_latitude_validity(self):
             SPoint(lon, lat)
 
     def test_longitude_validity(self):
+        """Check SPoint raises error if providing bad longitude."""
         # Test inf
         lon = np.inf
         lat = 0
         with pytest.raises(ValueError):
             SPoint(lon, lat)
 
+    def test_vertices(self):
+        """Test vertices property."""
+        lons = 0
+        lats = np.pi / 2
+        p = SPoint(lons, lats)
+        res = np.array([[0., 1.57079633]])
+        assert np.allclose(p.vertices, res)
+
+    def test_vertices_in_degrees(self):
+        """Test vertices_in_degrees property."""
+        lons = 0
+        lats = np.pi / 2
+        p = SPoint(lons, lats)
+        res = np.array([[0., 90.]])
+        assert np.allclose(p.vertices_in_degrees, res)
+
     def test_raise_error_if_multi_point(self):
+        """Check SPoint raises error providing multiple points."""
         lons = np.array([0, np.pi])
         lats = np.array([-np.pi / 2, np.pi / 2])
         with pytest.raises(ValueError):
             SPoint(lons, lats)
 
+    def test_str(self):
+        """Check the string representation."""
+        d = SPoint(1.0, 0.5)
+        self.assertEqual(str(d), "(1.0, 0.5)")
+
+    def test_repr(self):
+        """Check the representation."""
+        d = SPoint(1.0, 0.5)
+        self.assertEqual(repr(d), "(1.0, 0.5)")
+
     def test_to_shapely(self):
         """Test conversion to shapely."""
         from shapely.geometry import Point
@@ -58,7 +87,7 @@ def test_to_shapely(self):
         lat = np.pi / 2
         spherical_point = SPoint(lon, lat)
         shapely_point = Point(0.0, 90.0)
-        self.assertTrue(shapely_point.equals_exact(spherical_point.to_shapely(), tolerance=1e-10))
+        assert shapely_point.equals_exact(spherical_point.to_shapely(), tolerance=1e-10)
 
 
 class TestSMultiPoint(unittest.TestCase):
@@ -82,7 +111,16 @@ def test_vertices(self):
         p = SMultiPoint(lons, lats)
         res = np.array([[0., -1.57079633],
                         [-3.14159265, 1.57079633]])
-        self.assertTrue(np.allclose(p.vertices, res))
+        assert np.allclose(p.vertices, res)
+
+    def test_vertices_in_degrees(self):
+        """Test vertices_in_degrees property."""
+        lons = np.array([0, np.pi])
+        lats = np.array([-np.pi / 2, np.pi / 2])
+        p = SMultiPoint(lons, lats)
+        res = np.array([[0., -90.],
+                        [-180., 90.]])
+        assert np.allclose(p.vertices_in_degrees, res)
 
     def test_distance(self):
         """Test Vincenty formula."""
@@ -98,7 +136,7 @@ def test_distance(self):
         self.assertEqual(d21.shape, (3, 2))
         res = np.array([[0., 1.57079633, 3.14159265],
                         [3.14159265, 1.57079633, 0.]])
-        self.assertTrue(np.allclose(d12, res))
+        assert np.allclose(d12, res)
         # Special case with 1 point
         p1 = SMultiPoint(lons[[0]], lats[[0]])
         p2 = SMultiPoint(lons[[0]], lats[[0]])
@@ -119,7 +157,7 @@ def test_hdistance(self):
         self.assertEqual(d21.shape, (3, 2))
         res = np.array([[0., 1.57079633, 3.14159265],
                         [3.14159265, 1.57079633, 0.]])
-        self.assertTrue(np.allclose(d12, res))
+        assert np.allclose(d12, res)
 
     def test_eq(self):
         """Check the equality."""
@@ -127,7 +165,7 @@ def test_eq(self):
         lats = [-np.pi / 2, np.pi / 2]
         p = SMultiPoint(lons, lats)
         p1 = SMultiPoint(lons, lats)
-        self.assertTrue(p == p1)
+        assert p == p1
 
     def test_eq_antimeridian(self):
         """Check the equality with longitudes at -180/180 degrees."""
@@ -136,7 +174,7 @@ def test_eq_antimeridian(self):
         lats = [-np.pi / 2, np.pi / 2]
         p = SMultiPoint(lons, lats)
         p1 = SMultiPoint(lons1, lats)
-        self.assertTrue(p == p1)
+        assert p == p1
 
     def test_neq(self):
         """Check the equality."""
@@ -144,21 +182,23 @@ def test_neq(self):
         lats = [-np.pi / 2, np.pi / 2]
         p = SMultiPoint(lons, lats)
         p1 = SMultiPoint(lons + 0.1, lats)
-        self.assertTrue(p != p1)
+        assert p != p1
 
     def test_str(self):
         """Check the string representation."""
         lons = [0, np.pi]
         lats = [-np.pi / 2, np.pi / 2]
         p = SMultiPoint(lons, lats)
-        self.assertEqual(str(p), '[[   0.  -90.]\n [-180.   90.]]')
+        expected_str = '[[ 0.         -1.57079633]\n [-3.14159265  1.57079633]]'
+        self.assertEqual(str(p), expected_str)
 
     def test_repr(self):
         """Check the representation."""
         lons = [0, np.pi]
         lats = [-np.pi / 2, np.pi / 2]
         p = SMultiPoint(lons, lats)
-        self.assertEqual(repr(p), '[[   0.  -90.]\n [-180.   90.]]')
+        expected_repr = '[[ 0.         -1.57079633]\n [-3.14159265  1.57079633]]'
+        self.assertEqual(repr(p), expected_repr)
 
     def test_to_shapely(self):
         """Test conversion to shapely."""
@@ -167,4 +207,4 @@ def test_to_shapely(self):
         lats = np.array([-np.pi / 2, np.pi / 2])
         spherical_multipoint = SMultiPoint(lons, lats)
         shapely_multipoint = MultiPoint([(0.0, -90.0), (-180.0, 90.0)])
-        self.assertTrue(shapely_multipoint.equals_exact(spherical_multipoint.to_shapely(), tolerance=1e-10))
+        assert shapely_multipoint.equals_exact(spherical_multipoint.to_shapely(), tolerance=1e-10)

From b8f14c26c654fc270a7713b3b18c2d1411a75511 Mon Sep 17 00:00:00 2001
From: Gionata Ghiggi <gionata.ghiggi@gmail.com>
Date: Wed, 16 Nov 2022 11:46:31 +0100
Subject: [PATCH 18/20] Update repr docstring

---
 pyresample/future/spherical/point.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py
index 985c888ac..86a2866a5 100644
--- a/pyresample/future/spherical/point.py
+++ b/pyresample/future/spherical/point.py
@@ -68,7 +68,7 @@ def __str__(self):
         return str(self.vertices)
 
     def __repr__(self):
-        """Get simplified representation of lon/lat arrays in degrees."""
+        """Get simplified representation of lon/lat arrays in radians."""
         return str(self.vertices)
 
     def to_shapely(self):

From cd23b3c943ec3b11ed2c1879ffbcc5ae54566df6 Mon Sep 17 00:00:00 2001
From: Gionata Ghiggi <gionata.ghiggi@gmail.com>
Date: Wed, 16 Nov 2022 11:46:48 +0100
Subject: [PATCH 19/20] Update str docstring

---
 pyresample/future/spherical/point.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py
index 86a2866a5..f93122b16 100644
--- a/pyresample/future/spherical/point.py
+++ b/pyresample/future/spherical/point.py
@@ -64,7 +64,7 @@ def __eq__(self, other):
         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."""
+        """Get simplified representation of lon/lat arrays in radians."""
         return str(self.vertices)
 
     def __repr__(self):

From a48ac65696f7d109088e95972c96d7cd5eaff240 Mon Sep 17 00:00:00 2001
From: ghiggi <gionata.ghiggi@gmail.com>
Date: Wed, 16 Nov 2022 12:05:06 +0100
Subject: [PATCH 20/20] Add from_degrees class method

---
 pyresample/future/spherical/point.py     | 10 ++++++++++
 pyresample/test/test_sgeom/test_point.py | 16 ++++++++++++++++
 2 files changed, 26 insertions(+)

diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py
index f93122b16..33389cf05 100644
--- a/pyresample/future/spherical/point.py
+++ b/pyresample/future/spherical/point.py
@@ -34,6 +34,11 @@ def __init__(self, lon, lat):
             raise ValueError("Use SMultiPoint to define multiple points.")
         super().__init__(lon, lat)
 
+    @classmethod
+    def from_degrees(cls, lon, lat):
+        """Create SPoint from lon/lat coordinates in degrees."""
+        return cls(np.deg2rad(lon), np.deg2rad(lat))
+
     def __str__(self):
         """Get simplified representation of lon/lat arrays in radians."""
         return str((float(self.lon), float(self.lat)))
@@ -59,6 +64,11 @@ def __init__(self, lon, lat):
             raise ValueError("Use SPoint to define single points.")
         super().__init__(lon, lat)
 
+    @classmethod
+    def from_degrees(cls, lon, lat):
+        """Create SMultiPoint from lon/lat coordinates in degrees."""
+        return cls(np.deg2rad(lon), np.deg2rad(lat))
+
     def __eq__(self, other):
         """Check equality."""
         return np.allclose(self.lon, other.lon) and np.allclose(self.lat, other.lat)
diff --git a/pyresample/test/test_sgeom/test_point.py b/pyresample/test/test_sgeom/test_point.py
index 6722807e3..4a6075e45 100644
--- a/pyresample/test/test_sgeom/test_point.py
+++ b/pyresample/test/test_sgeom/test_point.py
@@ -47,6 +47,14 @@ def test_longitude_validity(self):
         with pytest.raises(ValueError):
             SPoint(lon, lat)
 
+    def test_creation_from_degrees(self):
+        """Check SPoint creation from lat/lon in degrees."""
+        lon = 0
+        lat = 20
+        p1 = SPoint.from_degrees(lon, lat)
+        p2 = SPoint(np.deg2rad(lon), np.deg2rad(lat))
+        assert p1 == p2
+
     def test_vertices(self):
         """Test vertices property."""
         lons = 0
@@ -104,6 +112,14 @@ def test_single_point(self):
         assert p.lat.shape == (1,)
         assert p.vertices.shape == (1, 2)
 
+    def test_creation_from_degrees(self):
+        """Check SMultiPoint creation from lat/lon in degrees."""
+        lon = np.array([0, 10])
+        lat = np.array([20, 30])
+        p1 = SMultiPoint.from_degrees(lon, lat)
+        p2 = SMultiPoint(np.deg2rad(lon), np.deg2rad(lat))
+        assert p1 == p2
+
     def test_vertices(self):
         """Test vertices property."""
         lons = np.array([0, np.pi])