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 .
+"""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..33389cf05
--- /dev/null
+++ b/pyresample/future/spherical/point.py
@@ -0,0 +1,88 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2022 Pyresample developers
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+"""Define single and multiple points on the sphere through SPoint and SMultiPoint classes."""
+import numpy as np
+
+from pyresample.spherical import SCoordinate
+
+
+class SPoint(SCoordinate):
+ """Object representing a single point on a sphere.
+
+ 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)
+
+ @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)))
+
+ 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(*self.vertices_in_degrees[0])
+ return point
+
+
+class SMultiPoint(SCoordinate):
+ """Object representing multiple points on a sphere."""
+
+ 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)
+
+ @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)
+
+ def __str__(self):
+ """Get simplified representation of lon/lat arrays in radians."""
+ return str(self.vertices)
+
+ def __repr__(self):
+ """Get simplified representation of lon/lat arrays in radians."""
+ return str(self.vertices)
+
+ def to_shapely(self):
+ """Convert the SMultiPoint to a shapely MultiPoint (in lon/lat degrees)."""
+ from shapely.geometry import MultiPoint
+ point = MultiPoint(self.vertices_in_degrees)
+ return point
diff --git a/pyresample/spherical.py b/pyresample/spherical.py
index f6676c494..9dca27904 100644
--- a/pyresample/spherical.py
+++ b/pyresample/spherical.py
@@ -34,22 +34,128 @@ 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 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:
+ 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.asarray([arr])
+ return arr
+
+
+def _vincenty_matrix(lon, lat, lon_ref, lat_ref):
+ """Compute a distance matrix using Vincenty formula.
+
+ 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)
+ 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 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)
+ 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) 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, :]
+ # Array values
+ else:
+ 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."""
+ """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 +168,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 +231,47 @@ def __iter__(self):
"""Get iterator over lon/lat pairs."""
return zip([self.lon, self.lat]).__iter__()
+ 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.
+ """
+ import matplotlib.pyplot as plt
+ try:
+ import cartopy.crs as ccrs
+ except ModuleNotFoundError:
+ raise ModuleNotFoundError(
+ "Install cartopy to plot spherical geometries. For example, 'pip install cartopy'.")
+
+ # Create figure if ax not provided
+ if ax is None:
+ if projection_crs is None:
+ projection_crs = ccrs.PlateCarree()
+ fig, ax = plt.subplots(subplot_kw=dict(projection=projection_crs))
+
+ # Plot Points
+ vertices = self.vertices_in_degrees
+ ax.scatter(x=vertices[:, 0],
+ y=vertices[:, 1],
+ **plot_kwargs)
+
+ # 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
+
+ return ax
+
class CCoordinate(object):
"""Cartesian coordinates."""
@@ -123,14 +283,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 +350,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 +412,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 .
+"""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..4a6075e45
--- /dev/null
+++ b/pyresample/test/test_sgeom/test_point.py
@@ -0,0 +1,226 @@
+#!/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 .
+"""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):
+ """Check SPoint raises error if providing bad latitude."""
+ # 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):
+ """Check SPoint raises error if providing bad longitude."""
+ # Test inf
+ lon = np.inf
+ lat = 0
+ 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
+ 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
+ lon = 0.0
+ lat = np.pi / 2
+ spherical_point = SPoint(lon, lat)
+ shapely_point = Point(0.0, 90.0)
+ assert shapely_point.equals_exact(spherical_point.to_shapely(), tolerance=1e-10)
+
+
+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_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])
+ lats = np.array([-np.pi / 2, np.pi / 2])
+ p = SMultiPoint(lons, lats)
+ res = np.array([[0., -1.57079633],
+ [-3.14159265, 1.57079633]])
+ 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."""
+ 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.]])
+ assert 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.]])
+ assert 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)
+ assert 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)
+ assert 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)
+ 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)
+ 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)
+ 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."""
+ 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)])
+ assert 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..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)
@@ -146,12 +171,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."""