Skip to content
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 12 additions & 6 deletions pyresample/future/spherical/point.py
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand All @@ -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
55 changes: 35 additions & 20 deletions pyresample/spherical.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -95,26 +97,26 @@ 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.")
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):
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


Expand All @@ -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, :]
Expand All @@ -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.

Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand Down
60 changes: 50 additions & 10 deletions pyresample/test/test_sgeom/test_point.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -39,26 +40,54 @@ 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
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))
assert shapely_point.equals_exact(spherical_point.to_shapely(), tolerance=1e-10)


class TestSMultiPoint(unittest.TestCase):
Expand All @@ -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."""
Expand All @@ -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]])
Expand All @@ -119,15 +157,15 @@ 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."""
lons = [0, np.pi]
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."""
Expand All @@ -136,29 +174,31 @@ 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."""
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)
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."""
Expand All @@ -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)