From 5bf261f9daacefe024043551214d9898780a3b3a Mon Sep 17 00:00:00 2001 From: ghiggi Date: Thu, 17 Nov 2022 15:10:01 +0100 Subject: [PATCH 01/13] Add legacy Arc code and tests --- pyresample/future/spherical/__init__.py | 2 + pyresample/future/spherical/arc.py | 181 ++++++++++++++++++++++++ pyresample/test/test_sgeom/test_arc.py | 160 +++++++++++++++++++++ 3 files changed, 343 insertions(+) create mode 100644 pyresample/future/spherical/arc.py create mode 100644 pyresample/test/test_sgeom/test_arc.py diff --git a/pyresample/future/spherical/__init__.py b/pyresample/future/spherical/__init__.py index 817b63fce..512ee3665 100644 --- a/pyresample/future/spherical/__init__.py +++ b/pyresample/future/spherical/__init__.py @@ -18,3 +18,5 @@ """Future features that are backwards incompatible with current functionality.""" from .point import SMultiPoint, SPoint # noqa + +EPSILON = 0.0000001 diff --git a/pyresample/future/spherical/arc.py b/pyresample/future/spherical/arc.py new file mode 100644 index 000000000..90cbbaa8d --- /dev/null +++ b/pyresample/future/spherical/arc.py @@ -0,0 +1,181 @@ +#!/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 a single great-circle arc on the sphere through the SArc class.""" +import numpy as np + +from pyresample.future.spherical import EPSILON +from pyresample.spherical import SCoordinate, _unwrap_radians + + +class Arc(object): + """An arc of the great circle between two points.""" + + def __init__(self, start, end): + self.start, self.end = start, end + + def __eq__(self, other): + """Check equality.""" + if self.start == other.start and self.end == other.end: + return 1 + return 0 + + def __ne__(self, other): + """Check not equal comparison.""" + return not self.__eq__(other) + + def __str__(self): + """Get simplified representation.""" + return str(self.start) + " -> " + str(self.end) + + def __repr__(self): + """Get simplified representation.""" + return str(self.start) + " -> " + str(self.end) + + def angle(self, other_arc): + """Oriented angle between two arcs. + + Returns: + Angle in radians. A straight line will be 0. A clockwise path + will be a negative angle and counter-clockwise will be positive. + + """ + if self.start == other_arc.start: + a__ = self.start + b__ = self.end + c__ = other_arc.end + elif self.start == other_arc.end: + a__ = self.start + b__ = self.end + c__ = other_arc.start + elif self.end == other_arc.end: + a__ = self.end + b__ = self.start + c__ = other_arc.start + elif self.end == other_arc.start: + a__ = self.end + b__ = self.start + c__ = other_arc.end + else: + raise ValueError("No common point in angle computation.") + + ua_ = a__.cross2cart(b__) + ub_ = a__.cross2cart(c__) + + val = ua_.dot(ub_) / (ua_.norm() * ub_.norm()) + + if abs(val - 1) < EPSILON: + angle = 0 + elif abs(val + 1) < EPSILON: + angle = np.pi + else: + angle = np.arccos(val) + + n__ = ua_.normalize() + if n__.dot(c__.to_cart()) > 0: + return -angle + else: + return angle + + def intersections(self, other_arc): + """Give the two intersections of the greats circles defined by the current arc and *other_arc*. + + From http://williams.best.vwh.net/intersect.htm + """ + end_lon = self.end.lon + other_end_lon = other_arc.end.lon + + if self.end.lon - self.start.lon > np.pi: + end_lon -= 2 * np.pi + if other_arc.end.lon - other_arc.start.lon > np.pi: + other_end_lon -= 2 * np.pi + if self.end.lon - self.start.lon < -np.pi: + end_lon += 2 * np.pi + if other_arc.end.lon - other_arc.start.lon < -np.pi: + other_end_lon += 2 * np.pi + + end_point = SCoordinate(end_lon, self.end.lat) + other_end_point = SCoordinate(other_end_lon, other_arc.end.lat) + + ea_ = self.start.cross2cart(end_point).normalize() + eb_ = other_arc.start.cross2cart(other_end_point).normalize() + + cross = ea_.cross(eb_) + lat = np.arctan2(cross.cart[2], + np.sqrt(cross.cart[0] ** 2 + cross.cart[1] ** 2)) + lon = np.arctan2(cross.cart[1], cross.cart[0]) + + return (SCoordinate(lon, lat), + SCoordinate(_unwrap_radians(lon + np.pi), -lat)) + + def intersects(self, other_arc): + """Check if the current arc and the *other_arc* intersect. + + An arc is defined as the shortest tracks between two points. + """ + return bool(self.intersection(other_arc)) + + def intersection(self, other_arc): + """Return where, if the current arc and the *other_arc* intersect. + + None is returned if there is not intersection. An arc is defined + as the shortest tracks between two points. + """ + if self == other_arc: + return None + + for i in self.intersections(other_arc): + a__ = self.start + b__ = self.end + c__ = other_arc.start + d__ = other_arc.end + + ab_ = a__.hdistance(b__) + cd_ = c__.hdistance(d__) + + if (((i in (a__, b__)) or + (abs(a__.hdistance(i) + b__.hdistance(i) - ab_) < EPSILON)) and + ((i in (c__, d__)) or + (abs(c__.hdistance(i) + d__.hdistance(i) - cd_) < EPSILON))): + return i + return None + + def get_next_intersection(self, arcs, known_inter=None): + """Get the next intersection between the current arc and *arcs*.""" + res = [] + for arc in arcs: + inter = self.intersection(arc) + if (inter is not None and + inter != arc.end and + inter != self.end): + res.append((inter, arc)) + + def dist(args): + """Get distance key.""" + return self.start.distance(args[0]) + + take_next = False + for inter, arc in sorted(res, key=dist): + if known_inter is not None: + if known_inter == inter: + take_next = True + elif take_next: + return inter, arc + else: + return inter, arc + + return None, None diff --git a/pyresample/test/test_sgeom/test_arc.py b/pyresample/test/test_sgeom/test_arc.py new file mode 100644 index 000000000..88bd16e71 --- /dev/null +++ b/pyresample/test/test_sgeom/test_arc.py @@ -0,0 +1,160 @@ +#!/usr/bin/env python +# -*- 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 . +"""Test cases for spherical geometry.""" + +import unittest + +import numpy as np + +from pyresample.spherical import Arc, SCoordinate + + +class TestArc(unittest.TestCase): + """Test arcs.""" + + def test_eq(self): + arc1 = Arc(SCoordinate(0, 0), + SCoordinate(np.deg2rad(10), np.deg2rad(10))) + arc2 = Arc(SCoordinate(0, np.deg2rad(10)), + SCoordinate(np.deg2rad(10), 0)) + + assert not arc1.__eq__(arc2) + assert arc1 == arc1 + + def test_ne(self): + arc1 = Arc(SCoordinate(0, 0), + SCoordinate(np.deg2rad(10), np.deg2rad(10))) + arc2 = Arc(SCoordinate(0, np.deg2rad(10)), + SCoordinate(np.deg2rad(10), 0)) + + assert arc1 != arc2 + assert not arc1.__ne__(arc1) + + def test_str(self): + arc1 = Arc(SCoordinate(0, 0), + SCoordinate(np.deg2rad(10), np.deg2rad(10))) + self.assertEqual(str(arc1), str(arc1.start) + " -> " + str(arc1.end)) + self.assertEqual(repr(arc1), str(arc1.start) + " -> " + str(arc1.end)) + + def test_intersection(self): + arc1 = Arc(SCoordinate(0, 0), + SCoordinate(np.deg2rad(10), np.deg2rad(10))) + arc2 = Arc(SCoordinate(0, np.deg2rad(10)), + SCoordinate(np.deg2rad(10), 0)) + lon, lat = arc1.intersection(arc2) + + np.testing.assert_allclose(np.rad2deg(lon), 5) + self.assertEqual(np.rad2deg(lat).round(7), round(5.0575148968282093, 7)) + + arc1 = Arc(SCoordinate(0, 0), + SCoordinate(np.deg2rad(10), np.deg2rad(10))) + + assert (arc1.intersection(arc1) is None) + + arc1 = Arc(SCoordinate(np.deg2rad(24.341215776575297), + np.deg2rad(44.987819588259327)), + SCoordinate(np.deg2rad(18.842727517611817), + np.deg2rad(46.512483610284178))) + arc2 = Arc(SCoordinate(np.deg2rad(20.165961750361905), + np.deg2rad(46.177305385810541)), + SCoordinate(np.deg2rad(20.253297585831707), + np.deg2rad(50.935830837274324))) + inter = SCoordinate(np.deg2rad(20.165957021925202), + np.deg2rad(46.177022633103398)) + self.assertEqual(arc1.intersection(arc2), inter) + + arc1 = Arc(SCoordinate(np.deg2rad(-2.4982818108326734), + np.deg2rad(48.596644847869655)), + SCoordinate(np.deg2rad(-2.9571441235622835), + np.deg2rad(49.165688435261394))) + arc2 = Arc(SCoordinate(np.deg2rad(-3.4976667413531688), + np.deg2rad(48.562704872921373)), + SCoordinate(np.deg2rad(-5.893976312685715), + np.deg2rad(48.445795283217116))) + + assert (arc1.intersection(arc2) is None) + + def test_angle(self): + arc1 = Arc(SCoordinate(np.deg2rad(157.5), + np.deg2rad(89.234600944314138)), + SCoordinate(np.deg2rad(90), + np.deg2rad(89))) + arc2 = Arc(SCoordinate(np.deg2rad(157.5), + np.deg2rad(89.234600944314138)), + SCoordinate(np.deg2rad(135), + np.deg2rad(89))) + + self.assertAlmostEqual(np.rad2deg(arc1.angle(arc2)), -44.996385007218926) + + arc1 = Arc(SCoordinate(np.deg2rad(112.5), + np.deg2rad(89.234600944314138)), + SCoordinate(np.deg2rad(90), np.deg2rad(89))) + arc2 = Arc(SCoordinate(np.deg2rad(112.5), + np.deg2rad(89.234600944314138)), + SCoordinate(np.deg2rad(45), np.deg2rad(89))) + + self.assertAlmostEqual(np.rad2deg(arc1.angle(arc2)), 44.996385007218883) + + arc1 = Arc(SCoordinate(0, 0), SCoordinate(1, 0)) + self.assertAlmostEqual(arc1.angle(arc1), 0) + + arc2 = Arc(SCoordinate(1, 0), SCoordinate(0, 0)) + self.assertAlmostEqual(arc1.angle(arc2), 0) + + arc2 = Arc(SCoordinate(0, 0), SCoordinate(-1, 0)) + self.assertAlmostEqual(arc1.angle(arc2), np.pi) + + arc2 = Arc(SCoordinate(2, 0), SCoordinate(1, 0)) + self.assertAlmostEqual(arc1.angle(arc2), np.pi) + + arc2 = Arc(SCoordinate(2, 0), SCoordinate(3, 0)) + self.assertRaises(ValueError, arc1.angle, arc2) + + def test_disjoint_arcs_crossing_antimeridian(self): + import copy + arc1 = Arc(SCoordinate(*np.deg2rad((143.76, 0))), + SCoordinate(*np.deg2rad((143.95, 7.33))) + ) + arc2 = Arc(SCoordinate(*np.deg2rad((170.34, 71.36))), + SCoordinate(*np.deg2rad((-171.03, 76.75))) + ) + arc1_orig = copy.deepcopy(arc1) + arc2_orig = copy.deepcopy(arc2) + point = arc1.intersection(arc2) + # Assert original arcs are unaffected + assert arc1_orig.end.lon == arc1.end.lon + assert arc2_orig.end.lon == arc2.end.lon + # Assert disjoint arcs returns None + assert isinstance(point, type(None)) + + def test_intersecting_arcs_crossing_antimeridian(self): + import copy + arc1 = Arc(SCoordinate(*np.deg2rad((-180.0, -90.0))), + SCoordinate(*np.deg2rad((-180.0, 90.0))) + ) + arc2 = Arc(SCoordinate(*np.deg2rad((-171.03, -76.75))), + SCoordinate(*np.deg2rad((170.34, -71.36))) + ) + arc1_orig = copy.deepcopy(arc1) + arc2_orig = copy.deepcopy(arc2) + point = arc1.intersection(arc2) + # Assert original arcs are unaffected + assert arc1_orig.end.lon == arc1.end.lon + assert arc2_orig.end.lon == arc2.end.lon + # Assert intersection result + assert point == SCoordinate(*np.deg2rad((-180.0, -74.7884716))) From e1cfa666c3ae1626190a22e543f6ae5d53feb5a3 Mon Sep 17 00:00:00 2001 From: ghiggi Date: Thu, 17 Nov 2022 20:08:09 +0100 Subject: [PATCH 02/13] Add new SArc and test units --- pyresample/future/spherical/__init__.py | 3 +- pyresample/future/spherical/arc.py | 257 +++++++++------- pyresample/future/spherical/point.py | 32 +- pyresample/spherical.py | 24 +- pyresample/test/test_sgeom/test_arc.py | 370 ++++++++++++++++------- pyresample/test/test_sgeom/test_point.py | 91 ++++-- 6 files changed, 521 insertions(+), 256 deletions(-) diff --git a/pyresample/future/spherical/__init__.py b/pyresample/future/spherical/__init__.py index 512ee3665..108158df7 100644 --- a/pyresample/future/spherical/__init__.py +++ b/pyresample/future/spherical/__init__.py @@ -17,6 +17,5 @@ # with this program. If not, see . """Future features that are backwards incompatible with current functionality.""" +from .arc import SArc # noqa from .point import SMultiPoint, SPoint # noqa - -EPSILON = 0.0000001 diff --git a/pyresample/future/spherical/arc.py b/pyresample/future/spherical/arc.py index 90cbbaa8d..0ee06a88a 100644 --- a/pyresample/future/spherical/arc.py +++ b/pyresample/future/spherical/arc.py @@ -15,86 +15,84 @@ # # You should have received a copy of the GNU General Public License # along with this program. If not, see . -"""Define a single great-circle arc on the sphere through the SArc class.""" +"""Define a great-circle arc between two points on the sphere through the SArc class.""" + import numpy as np +import pyproj +from shapely.geometry import LineString + +from pyresample.future.spherical import SPoint +from pyresample.spherical import Arc + +EPSILON = 0.0000001 + + +def _check_valid_arc(start_point, end_point): + """Check arc validity.""" + if start_point == end_point: + raise ValueError("An SArc can not be represented by the same start and end SPoint.") + if start_point.is_pole() and end_point.is_pole(): + raise ValueError("An SArc can not be uniquely defined between the two poles.") + if start_point.is_on_equator and end_point.is_on_equator and abs(start_point.lon - end_point.lon) == np.pi: + raise ValueError( + "An SArc can not be uniquely defined on the equator if start and end points are 180 degrees apart.") -from pyresample.future.spherical import EPSILON -from pyresample.spherical import SCoordinate, _unwrap_radians +class SArc(Arc): + """Object representing a great-circle arc between two points on a sphere. -class Arc(object): - """An arc of the great circle between two points.""" + The ``start_point`` and ``end_point`` must be SPoint objects. + The great-circle arc is defined as the shortest track(s) between the two points. + Between the north and south pole there are an infinite number of great-circle arcs. + """ - def __init__(self, start, end): - self.start, self.end = start, end + def __init__(self, start_point, end_point): + _check_valid_arc(start_point, end_point) + super().__init__(start_point, end_point) + + def __hash__(self): + """Define SArc hash.""" + return hash((float(self.start.lon), float(self.start.lat), + float(self.end.lon), float(self.end.lat))) + + def is_on_equator(self): + """Check if the SArc lies on the equator.""" + if self.start.lat == 0 and self.end.lat == 0: + return True + return False def __eq__(self, other): """Check equality.""" if self.start == other.start and self.end == other.end: - return 1 - return 0 - - def __ne__(self, other): - """Check not equal comparison.""" - return not self.__eq__(other) - - def __str__(self): - """Get simplified representation.""" - return str(self.start) + " -> " + str(self.end) + return True + return False - def __repr__(self): - """Get simplified representation.""" - return str(self.start) + " -> " + str(self.end) + def reverse_direction(self): + """Reverse SArc direction.""" + return SArc(self.end, self.start) - def angle(self, other_arc): - """Oriented angle between two arcs. + @property + def vertices(self): + """Get start SPoint and end SPoint (radians) vertices array.""" + return self.start.vertices, self.end.vertices - Returns: - Angle in radians. A straight line will be 0. A clockwise path - will be a negative angle and counter-clockwise will be positive. + @property + def vertices_in_degrees(self): + """Get start SPoint and end SPoint (degrees) vertices array.""" + return self.start.vertices_in_degrees, self.end.vertices_in_degrees - """ - if self.start == other_arc.start: - a__ = self.start - b__ = self.end - c__ = other_arc.end - elif self.start == other_arc.end: - a__ = self.start - b__ = self.end - c__ = other_arc.start - elif self.end == other_arc.end: - a__ = self.end - b__ = self.start - c__ = other_arc.start - elif self.end == other_arc.start: - a__ = self.end - b__ = self.start - c__ = other_arc.end - else: - raise ValueError("No common point in angle computation.") - - ua_ = a__.cross2cart(b__) - ub_ = a__.cross2cart(c__) - - val = ua_.dot(ub_) / (ua_.norm() * ub_.norm()) - - if abs(val - 1) < EPSILON: - angle = 0 - elif abs(val + 1) < EPSILON: - angle = np.pi - else: - angle = np.arccos(val) - - n__ = ua_.normalize() - if n__.dot(c__.to_cart()) > 0: - return -angle - else: - return angle + def get_next_intersection(self, other_arc): + """Overwrite a Arc deprecated function. Inherited from Arc class.""" + raise ValueError("This function is deprecated.") def intersections(self, other_arc): - """Give the two intersections of the greats circles defined by the current arc and *other_arc*. + """Overwritea Arc deprecated function. Inherited from Arc class.""" + raise ValueError("'SArc.intersections' is deprecated. Use '_great_circle_intersections' instead.") + + def _great_circle_intersections(self, other_arc): + """Compute the intersections points of the greats circles over which the arcs lies. - From http://williams.best.vwh.net/intersect.htm + A great circle divides the sphere in two equal hemispheres. """ end_lon = self.end.lon other_end_lon = other_arc.end.lon @@ -108,8 +106,8 @@ def intersections(self, other_arc): if other_arc.end.lon - other_arc.start.lon < -np.pi: other_end_lon += 2 * np.pi - end_point = SCoordinate(end_lon, self.end.lat) - other_end_point = SCoordinate(other_end_lon, other_arc.end.lat) + end_point = SPoint(end_lon, self.end.lat) + other_end_point = SPoint(other_end_lon, other_arc.end.lat) ea_ = self.start.cross2cart(end_point).normalize() eb_ = other_arc.start.cross2cart(other_end_point).normalize() @@ -119,63 +117,92 @@ def intersections(self, other_arc): np.sqrt(cross.cart[0] ** 2 + cross.cart[1] ** 2)) lon = np.arctan2(cross.cart[1], cross.cart[0]) - return (SCoordinate(lon, lat), - SCoordinate(_unwrap_radians(lon + np.pi), -lat)) - - def intersects(self, other_arc): - """Check if the current arc and the *other_arc* intersect. - - An arc is defined as the shortest tracks between two points. - """ - return bool(self.intersection(other_arc)) + return (SPoint(lon, lat), SPoint(lon + np.pi, -lat)) def intersection(self, other_arc): - """Return where, if the current arc and the *other_arc* intersect. + """Overwrite a Arc deprecated function. Inherited from Arc class.""" + raise ValueError("'SArc.intersection' is deprecated. Use 'intersection_point' instead.") + + def intersection_point(self, other_arc): + """Compute the intersection point between two arcs. - None is returned if there is not intersection. An arc is defined - as the shortest tracks between two points. + If arc and *other_arc* intersect, it returns the intersection SPoint. + If arc and *other_arc* does not intersect, it returns None. + If same arc (also same direction), it returns None. """ + # If same arc (same direction), return None if self == other_arc: return None - for i in self.intersections(other_arc): - a__ = self.start - b__ = self.end - c__ = other_arc.start - d__ = other_arc.end + great_circles_intersection_spoints = self._great_circle_intersections(other_arc) - ab_ = a__.hdistance(b__) - cd_ = c__.hdistance(d__) + for spoint in great_circles_intersection_spoints: + a = self.start + b = self.end + c = other_arc.start + d = other_arc.end - if (((i in (a__, b__)) or - (abs(a__.hdistance(i) + b__.hdistance(i) - ab_) < EPSILON)) and - ((i in (c__, d__)) or - (abs(c__.hdistance(i) + d__.hdistance(i) - cd_) < EPSILON))): - return i + ab_dist = a.hdistance(b) + cd_dist = c.hdistance(d) + ap_dist = a.hdistance(spoint) + bp_dist = b.hdistance(spoint) + cp_dist = c.hdistance(spoint) + dp_dist = d.hdistance(spoint) + + if (((spoint in (a, b)) or (abs(ap_dist + bp_dist - ab_dist) < EPSILON)) and + ((spoint in (c, d)) or (abs(cp_dist + dp_dist - cd_dist) < EPSILON))): + return spoint return None - def get_next_intersection(self, arcs, known_inter=None): - """Get the next intersection between the current arc and *arcs*.""" - res = [] - for arc in arcs: - inter = self.intersection(arc) - if (inter is not None and - inter != arc.end and - inter != self.end): - res.append((inter, arc)) - - def dist(args): - """Get distance key.""" - return self.start.distance(args[0]) - - take_next = False - for inter, arc in sorted(res, key=dist): - if known_inter is not None: - if known_inter == inter: - take_next = True - elif take_next: - return inter, arc - else: - return inter, arc - - return None, None + def intersects(self, other_arc): + """Check if the current Sarc and another SArc intersect.""" + return bool(self.intersection_point(other_arc)) + + def midpoint(self, ellips='sphere'): + """Return the SArc midpoint SPoint.""" + geod = pyproj.Geod(ellps=ellips) + lon_start = self.start.lon + lon_end = self.end.lon + lat_start = self.start.lat + lat_end = self.end.lat + lon_mid, lat_mid = geod.npts(lon_start, lat_start, lon_end, lat_end, npts=1, radians=True)[0] + return SPoint(lon_mid, lat_mid) + + def to_shapely(self): + """Convert to Shapely LineString.""" + start_coord, end_coord = self.vertices_in_degrees + return LineString((start_coord.tolist()[0], end_coord.tolist()[0])) + + # def segmentize(self, npts=0, max_distance=0, ellips='sphere'): + # """Segmentize the great-circle arc. + + # It returns an SLine. + # npts or max_distance are mutually exclusively. Specify one of them. + # max_distance must be provided in kilometers. + # """ + # if npts != 0: + # npts = npts + 2 # + 2 to account for initial and terminus + # geod = pyproj.Geod(ellps=ellips) + # lon_start = self.start.lon + # lon_end = self.end.lon + # lat_start = self.start.lat + # lat_end = self.end.lat + # points = geod.inv_intermediate(lon_start, lat_start, lon_end, lat_end, + # del_s=max_distance, + # npts=npts, + # radians=True, + # initial_idx=0, terminus_idx=0) + # lons, lats = (points.lons, points.lats) + # lons = np.asarray(lons) + # lats = np.asarray(lats) + # vertices = np.stack((lons, lats)).T + # return SLine(vertices) + + # def to_line(self): + # """Convert to SLine.""" + # vertices = np.vstack(self.vertices) + # return SLine(vertices) + + # def plot(self, *args, **kwargs): + # """Convert to SLine.""" + # self.to_line.plot(*args, **kwargs) diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py index 33389cf05..dd3b84e93 100644 --- a/pyresample/future/spherical/point.py +++ b/pyresample/future/spherical/point.py @@ -39,13 +39,32 @@ def from_degrees(cls, lon, lat): """Create SPoint from lon/lat coordinates in degrees.""" return cls(np.deg2rad(lon), np.deg2rad(lat)) + def is_pole(self): + """Test if the point is on a pole.""" + if self.lat in [-np.pi / 2, np.pi / 2]: + return True + return False + + def is_on_equator(self): + """Test if the point is on the equator.""" + if self.lat == 0: + return True + return False + 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))) + coords = str((float(self.lon), float(self.lat))) + return f"{self.__class__.__name__}: {coords}" + + def __eq__(self, other): + """Check equality.""" + if self.lat == other.lat and self.is_pole(): + return True + return np.allclose((self.lon, self.lat), (other.lon, other.lat)) def to_shapely(self): """Convert the SPoint to a shapely Point (in lon/lat degrees).""" @@ -71,7 +90,14 @@ def from_degrees(cls, lon, lat): def __eq__(self, other): """Check equality.""" - return np.allclose(self.lon, other.lon) and np.allclose(self.lat, other.lat) + lon1 = self.lon.copy() + lon2 = other.lon.copy() + lat1 = self.lat + lat2 = other.lat + # Set longitude 0 at the pole for array comparison + lon1[np.isin(lat1, [-np.pi / 2, np.pi / 2])] = 0 + lon2[np.isin(lat2, [-np.pi / 2, np.pi / 2])] = 0 + return np.allclose(lon1, lon2) and np.allclose(lat1, lat2) def __str__(self): """Get simplified representation of lon/lat arrays in radians.""" @@ -79,7 +105,7 @@ def __str__(self): def __repr__(self): """Get simplified representation of lon/lat arrays in radians.""" - return str(self.vertices) + return f"{self.__class__.__name__}:\n {self.vertices}" def to_shapely(self): """Convert the SMultiPoint to a shapely MultiPoint (in lon/lat degrees).""" diff --git a/pyresample/spherical.py b/pyresample/spherical.py index 9dca27904..de76494b5 100644 --- a/pyresample/spherical.py +++ b/pyresample/spherical.py @@ -366,19 +366,19 @@ def __init__(self, start, end): def __eq__(self, other): """Check equality.""" if self.start == other.start and self.end == other.end: - return 1 - return 0 + return True + return False def __ne__(self, other): """Check not equal comparison.""" return not self.__eq__(other) def __str__(self): - """Get simplified representation.""" + """Get simplified str representation.""" return str(self.start) + " -> " + str(self.end) def __repr__(self): - """Get simplified representation.""" + """Get simplified repr representation.""" return str(self.start) + " -> " + str(self.end) def angle(self, other_arc): @@ -427,9 +427,9 @@ def angle(self, other_arc): return angle def intersections(self, other_arc): - """Give the two intersections of the greats circles defined by the current arc and *other_arc*. + """Compute the intersections points of the greats circles over which the arcs lies. - From http://williams.best.vwh.net/intersect.htm + A great circle divides the sphere in two equal hemispheres. """ end_lon = self.end.lon other_end_lon = other_arc.end.lon @@ -465,10 +465,10 @@ def intersects(self, other_arc): return bool(self.intersection(other_arc)) def intersection(self, other_arc): - """Return where, if the current arc and the *other_arc* intersect. + """Compute the intersection point between two arcs. - None is returned if there is not intersection. An arc is defined - as the shortest tracks between two points. + If arc and *other_arc* intersect, it returns the intersection SPoint. + If arc and *other_arc* does not intersect, it returns None. """ if self == other_arc: return None @@ -490,7 +490,11 @@ def intersection(self, other_arc): return None def get_next_intersection(self, arcs, known_inter=None): - """Get the next intersection between the current arc and *arcs*.""" + """Get the next intersection between the current arc and *arcs*. + + It return a tuple with the intersecting point and the arc within *arcs* + that intersect the self arc. + """ res = [] for arc in arcs: inter = self.intersection(arc) diff --git a/pyresample/test/test_sgeom/test_arc.py b/pyresample/test/test_sgeom/test_arc.py index 88bd16e71..640b86af8 100644 --- a/pyresample/test/test_sgeom/test_arc.py +++ b/pyresample/test/test_sgeom/test_arc.py @@ -1,7 +1,7 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 # -*- coding: utf-8 -*- # -# Copyright (c) 2022 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 @@ -15,146 +15,298 @@ # # You should have received a copy of the GNU General Public License # along with this program. If not, see . -"""Test cases for spherical geometry.""" +"""Define tests for the SArc class.""" +import copy import unittest import numpy as np +import pytest +from shapely.geometry import LineString -from pyresample.spherical import Arc, SCoordinate +from pyresample.future.spherical import SArc, SPoint -class TestArc(unittest.TestCase): - """Test arcs.""" +class TestSArc(unittest.TestCase): + """Test SArc class.""" - def test_eq(self): - arc1 = Arc(SCoordinate(0, 0), - SCoordinate(np.deg2rad(10), np.deg2rad(10))) - arc2 = Arc(SCoordinate(0, np.deg2rad(10)), - SCoordinate(np.deg2rad(10), 0)) + # TODO: Fixtures defined here? + equator_arc1 = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(50, 0.0)) + equator_arc2 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(0.0, 0.0)) + + # pole_to_pole_arc1 = SArc(SPoint.from_degrees(0.0, -90.0), SPoint.from_degrees(0.0, 90.0)) + # pole_to_pole_arc2 = SArc(SPoint.from_degrees(60.0, -90.0), SPoint.from_degrees(60.0, 90.0)) + + def test_unvalid_point_arc(self): + """Check raise error if start and end point are equal.""" + p = SPoint.from_degrees(0.0, 80.0) + with pytest.raises(ValueError): + SArc(p, p) + + def test_unvalid_180degree_equator_arc(self): + """Check raise error if the points lies on the equator and are 180° apart.""" + p1 = SPoint.from_degrees(0, 0) + p2 = SPoint.from_degrees(180, 0) + with pytest.raises(ValueError): + SArc(p1, p2) + + p1 = SPoint.from_degrees(-10, 0) + p2 = SPoint.from_degrees(170, 0) + with pytest.raises(ValueError): + SArc(p1, p2) + + p1 = SPoint.from_degrees(10, 0) + p2 = SPoint.from_degrees(-170, 0) + with pytest.raises(ValueError): + SArc(p1, p2) + def test_unvalid_pole_to_pole_arc(self): + """Check raise error if the points defines a pole to pole arc.""" + p1 = SPoint.from_degrees(0.0, -90.0) + p2 = SPoint.from_degrees(0.0, 90.0) + with pytest.raises(ValueError): + SArc(p1, p2) + + def test_is_on_equator_arc(self): + """Check if the arc lies on the equator.""" + equator_arc1 = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(50, 0.0)) + arc = SArc(SPoint.from_degrees(0, 10), SPoint.from_degrees(10, 0)) + assert equator_arc1.is_on_equator() + assert not arc.is_on_equator() + + def test_eq(self): + """Test SArc equality.""" + arc1 = SArc(SPoint.from_degrees(0, 0), SPoint.from_degrees(10, 10)) + arc2 = SArc(SPoint.from_degrees(0, 10), SPoint.from_degrees(10, 0)) assert not arc1.__eq__(arc2) assert arc1 == arc1 def test_ne(self): - arc1 = Arc(SCoordinate(0, 0), - SCoordinate(np.deg2rad(10), np.deg2rad(10))) - arc2 = Arc(SCoordinate(0, np.deg2rad(10)), - SCoordinate(np.deg2rad(10), 0)) - + """Test SArc disequality.""" + arc1 = SArc(SPoint.from_degrees(0, 0), SPoint.from_degrees(10, 10)) + arc2 = SArc(SPoint.from_degrees(0, 10), SPoint.from_degrees(10, 0)) assert arc1 != arc2 assert not arc1.__ne__(arc1) def test_str(self): - arc1 = Arc(SCoordinate(0, 0), - SCoordinate(np.deg2rad(10), np.deg2rad(10))) - self.assertEqual(str(arc1), str(arc1.start) + " -> " + str(arc1.end)) - self.assertEqual(repr(arc1), str(arc1.start) + " -> " + str(arc1.end)) - - def test_intersection(self): - arc1 = Arc(SCoordinate(0, 0), - SCoordinate(np.deg2rad(10), np.deg2rad(10))) - arc2 = Arc(SCoordinate(0, np.deg2rad(10)), - SCoordinate(np.deg2rad(10), 0)) - lon, lat = arc1.intersection(arc2) + """Test SArc __str__ representation.""" + arc = SArc(SPoint.from_degrees(0, 0), SPoint.from_degrees(10, 10)) + expected_str = str(arc.start) + " -> " + str(arc.end) + assert str(arc) == expected_str + + def test_repr(self): + """Test SArc __repr__ representation.""" + arc = SArc(SPoint.from_degrees(0, 0), SPoint.from_degrees(10, 10)) + expected_repr = str(arc.start) + " -> " + str(arc.end) + assert repr(arc) == expected_repr + + def test_vertices(self): + """Test vertices property.""" + start_point_vertices = np.deg2rad((-180.0, 0.0)) + end_point_vertices = np.deg2rad((-180.0, 90.0)) + arc = SArc(SPoint(*start_point_vertices), SPoint(*end_point_vertices)) + start_vertices, end_vertices = arc.vertices + assert np.allclose(start_vertices, start_point_vertices) + assert np.allclose(end_vertices, end_point_vertices) + + def test_vertices_in_degrees(self): + """Test vertices_in_degrees property.""" + start_point_vertices = np.deg2rad((-180.0, 0.0)) + end_point_vertices = np.deg2rad((-180.0, 90.0)) + arc = SArc(SPoint.from_degrees(* start_point_vertices), + SPoint.from_degrees(* end_point_vertices)) + start_vertices, end_vertices = arc.vertices_in_degrees + assert np.allclose(start_vertices, start_point_vertices) + assert np.allclose(end_vertices, end_point_vertices) + + def test_to_shapely(self): + """Test conversion to shapely.""" + start_point_vertices = np.deg2rad((-180.0, 0.0)) + end_point_vertices = np.deg2rad((-180.0, 90.0)) + arc = SArc(SPoint.from_degrees(* start_point_vertices), + SPoint.from_degrees(* end_point_vertices)) + shapely_line = arc.to_shapely() + expected_line = LineString((start_point_vertices, end_point_vertices)) + assert shapely_line.equals_exact(expected_line, tolerance=1e-10) + + def test_hash(self): + """Test arc hash.""" + arc = SArc(SPoint.from_degrees(-180.0, -90.0), SPoint.from_degrees(-180.0, 0.0)) + assert hash(arc) == -3096892178517935054 + + def test_midpoint(self): + """Test arc midpoint.""" + start_point_vertices = np.deg2rad((-180.0, 0.0)) + end_point_vertices = np.deg2rad((-180.0, 90.0)) + arc = SArc(SPoint(*start_point_vertices), SPoint(*end_point_vertices)) + midpoint = arc.midpoint() + assert isinstance(midpoint, SPoint) + assert np.allclose(midpoint.vertices_in_degrees[0], (-180, 45)) + + def test_intersection_point(self): + """Test SArc(s) intersection point.""" + arc1 = SArc(SPoint.from_degrees(0, 0), SPoint.from_degrees(10, 10)) + arc2 = SArc(SPoint.from_degrees(0, 10), SPoint.from_degrees(10, 0)) + lon, lat = arc1.intersection_point(arc2) np.testing.assert_allclose(np.rad2deg(lon), 5) - self.assertEqual(np.rad2deg(lat).round(7), round(5.0575148968282093, 7)) - - arc1 = Arc(SCoordinate(0, 0), - SCoordinate(np.deg2rad(10), np.deg2rad(10))) - - assert (arc1.intersection(arc1) is None) - - arc1 = Arc(SCoordinate(np.deg2rad(24.341215776575297), - np.deg2rad(44.987819588259327)), - SCoordinate(np.deg2rad(18.842727517611817), - np.deg2rad(46.512483610284178))) - arc2 = Arc(SCoordinate(np.deg2rad(20.165961750361905), - np.deg2rad(46.177305385810541)), - SCoordinate(np.deg2rad(20.253297585831707), - np.deg2rad(50.935830837274324))) - inter = SCoordinate(np.deg2rad(20.165957021925202), - np.deg2rad(46.177022633103398)) - self.assertEqual(arc1.intersection(arc2), inter) - - arc1 = Arc(SCoordinate(np.deg2rad(-2.4982818108326734), - np.deg2rad(48.596644847869655)), - SCoordinate(np.deg2rad(-2.9571441235622835), - np.deg2rad(49.165688435261394))) - arc2 = Arc(SCoordinate(np.deg2rad(-3.4976667413531688), - np.deg2rad(48.562704872921373)), - SCoordinate(np.deg2rad(-5.893976312685715), - np.deg2rad(48.445795283217116))) - - assert (arc1.intersection(arc2) is None) - - def test_angle(self): - arc1 = Arc(SCoordinate(np.deg2rad(157.5), - np.deg2rad(89.234600944314138)), - SCoordinate(np.deg2rad(90), - np.deg2rad(89))) - arc2 = Arc(SCoordinate(np.deg2rad(157.5), - np.deg2rad(89.234600944314138)), - SCoordinate(np.deg2rad(135), - np.deg2rad(89))) - - self.assertAlmostEqual(np.rad2deg(arc1.angle(arc2)), -44.996385007218926) - - arc1 = Arc(SCoordinate(np.deg2rad(112.5), - np.deg2rad(89.234600944314138)), - SCoordinate(np.deg2rad(90), np.deg2rad(89))) - arc2 = Arc(SCoordinate(np.deg2rad(112.5), - np.deg2rad(89.234600944314138)), - SCoordinate(np.deg2rad(45), np.deg2rad(89))) - - self.assertAlmostEqual(np.rad2deg(arc1.angle(arc2)), 44.996385007218883) - - arc1 = Arc(SCoordinate(0, 0), SCoordinate(1, 0)) - self.assertAlmostEqual(arc1.angle(arc1), 0) - - arc2 = Arc(SCoordinate(1, 0), SCoordinate(0, 0)) - self.assertAlmostEqual(arc1.angle(arc2), 0) - - arc2 = Arc(SCoordinate(0, 0), SCoordinate(-1, 0)) - self.assertAlmostEqual(arc1.angle(arc2), np.pi) - - arc2 = Arc(SCoordinate(2, 0), SCoordinate(1, 0)) - self.assertAlmostEqual(arc1.angle(arc2), np.pi) - - arc2 = Arc(SCoordinate(2, 0), SCoordinate(3, 0)) - self.assertRaises(ValueError, arc1.angle, arc2) + assert np.rad2deg(lat).round(7) == round(5.0575148968282093, 7) + + # Test when same SArc + arc1 = SArc(SPoint(0, 0), SPoint.from_degrees(10, 10)) + assert (arc1.intersection_point(arc1) is None) + + # Test intersecting SArc(s) + arc1 = SArc(SPoint.from_degrees(24.341215776575297, 44.987819588259327), + SPoint.from_degrees(18.842727517611817, 46.512483610284178)) + arc2 = SArc(SPoint.from_degrees(20.165961750361905, 46.177305385810541), + SPoint.from_degrees(20.253297585831707, 50.935830837274324)) + p = SPoint.from_degrees(20.165957021925202, 46.177022633103398) + assert arc1.intersection_point(arc2) == p + assert arc1.intersects(arc2) + + # Test non-intersecting SArc(s) + arc1 = SArc(SPoint.from_degrees(-2.4982818108326734, 48.596644847869655), + SPoint.from_degrees(-2.9571441235622835, 49.165688435261394)) + arc2 = SArc(SPoint.from_degrees(-3.4976667413531688, 48.562704872921373), + SPoint.from_degrees(-5.893976312685715, 48.445795283217116)) + assert arc1.intersection_point(arc2) is None + assert not arc1.intersects(arc2) + assert not bool(None) # this occurs within the intersects method def test_disjoint_arcs_crossing_antimeridian(self): - import copy - arc1 = Arc(SCoordinate(*np.deg2rad((143.76, 0))), - SCoordinate(*np.deg2rad((143.95, 7.33))) - ) - arc2 = Arc(SCoordinate(*np.deg2rad((170.34, 71.36))), - SCoordinate(*np.deg2rad((-171.03, 76.75))) - ) + """Test SArc(s) intersection point when disjoint arcs cross the antimeridian.""" + arc1 = SArc(SPoint.from_degrees(143.76, 0), + SPoint.from_degrees(143.95, 7.33)) + arc2 = SArc(SPoint.from_degrees(170.34, 71.36), + SPoint.from_degrees(-171.03, 76.75)) arc1_orig = copy.deepcopy(arc1) arc2_orig = copy.deepcopy(arc2) - point = arc1.intersection(arc2) + point = arc1.intersection_point(arc2) + # Assert original arcs are unaffected assert arc1_orig.end.lon == arc1.end.lon assert arc2_orig.end.lon == arc2.end.lon + # Assert disjoint arcs returns None assert isinstance(point, type(None)) def test_intersecting_arcs_crossing_antimeridian(self): - import copy - arc1 = Arc(SCoordinate(*np.deg2rad((-180.0, -90.0))), - SCoordinate(*np.deg2rad((-180.0, 90.0))) - ) - arc2 = Arc(SCoordinate(*np.deg2rad((-171.03, -76.75))), - SCoordinate(*np.deg2rad((170.34, -71.36))) - ) + """Test SArc(s) intersection point when intersecting arcs cross the antimeridian.""" + arc1 = SArc(SPoint.from_degrees(-180.0, -89.0), + SPoint.from_degrees(-180.0, 89.0)) + arc2 = SArc(SPoint.from_degrees(-171.03, -76.75), + SPoint.from_degrees(170.34, -71.36)) arc1_orig = copy.deepcopy(arc1) arc2_orig = copy.deepcopy(arc2) - point = arc1.intersection(arc2) + point = arc1.intersection_point(arc2) # Assert original arcs are unaffected assert arc1_orig.end.lon == arc1.end.lon assert arc2_orig.end.lon == arc2.end.lon # Assert intersection result - assert point == SCoordinate(*np.deg2rad((-180.0, -74.7884716))) + assert point == SPoint.from_degrees(-180.0, -74.78847163) + + def test_great_circles_equator_arcs_cases(self): + """Test behaviour when 2 arcs are lying around the equator.""" + equator_arc1 = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(50, 0.0)) + equator_arc2 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(0.0, 0.0)) + + # - If one arc is reversed, the results order is swapped + p1, p2 = equator_arc1._great_circle_intersections(equator_arc2) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + p1, p2 = equator_arc2._great_circle_intersections(equator_arc1) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + # - If both arc are reversed, results are the same + p1, p2 = equator_arc2.reverse_direction()._great_circle_intersections(equator_arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + # - Reversing the direction of one arc, lead inversion of points + p1, p2 = equator_arc2._great_circle_intersections(equator_arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-180, 0)) + assert np.allclose(p2.vertices_in_degrees, (0, 0)) + + def test_intersection_equator_arcs_case(self): + """Test intersection point with 2 arcs lying around the equator.""" + equator_arc1 = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(50, 0.0)) + equator_arc2 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(0.0, 0.0)) + equator_arc3 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(-20, 0.0)) + + # Touching arcs + p = equator_arc2.intersection_point(equator_arc1) + assert np.allclose(p.vertices_in_degrees, (0, 0)) + + p = equator_arc2.intersection_point(equator_arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (0, 0)) + + # Overlapping arcs + p = equator_arc2.intersection_point(equator_arc3) + assert isinstance(p, type(None)) + + # Disjoint arcs + p = equator_arc3.intersection_point(equator_arc1) + assert isinstance(p, type(None)) + + # ---------------------------------------------------------------------. + # TEST FOR EDGE CASE SAME ARC + + def test_intersection_same_arcs(self): + """Test behaviour when same arc.""" + nh_arc = SArc(SPoint.from_degrees(0.0, 20.0), SPoint.from_degrees(10.0, 40.0)) + # TODO BUG ! CHECK LEGAGY BEHAVIOUR + # - Define behaviour when same arc !!! + p1, p2 = nh_arc._great_circle_intersections(nh_arc) # (0,0), (-180, 0) + p = nh_arc.intersection_point(nh_arc) # None + assert isinstance(p, type(None)) + + p1, p2 = nh_arc._great_circle_intersections(nh_arc.reverse_direction()) # (0,0), (-180, 0) + p = nh_arc.intersection_point(nh_arc.reverse_direction()) # None + assert isinstance(p, type(None)) + + # ---------------------------------------------------------------------. + # TEST FOR TOUCHING CASES + + def test_touching_arc_extremes_case(self): + """Test behaviour when two arc touches at the extremes.""" + # Test arc1.end == arc2.start + arc1 = SArc(SPoint.from_degrees(0.0, 10.0), SPoint.from_degrees(0.0, 20.0)) + arc2 = SArc(SPoint.from_degrees(0.0, 20.0), SPoint.from_degrees(0.0, 30.0)) + + p = arc1.intersection_point(arc2) + assert isinstance(p, type(None)) + assert not arc1.intersects(arc2) + + p = arc2.intersection_point(arc1) + assert isinstance(p, type(None)) + + p = arc2.intersection_point(arc1.reverse_direction()) + assert isinstance(p, type(None)) + + # Test arc1.start == arc2.start + arc1 = SArc(SPoint.from_degrees(0.0, 10.0), SPoint.from_degrees(0.0, 20.0)) + arc2 = SArc(SPoint.from_degrees(0.0, 10.0), SPoint.from_degrees(0.0, -20.0)) + p = arc1.intersection_point(arc2) + assert isinstance(p, type(None)) + assert not arc1.intersects(arc2) + + # TODO: BUG Bad behaviour at the equator. + equator_arc1 = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(50, 0.0)) + arc = SArc(SPoint.from_degrees(5, 0.0), SPoint.from_degrees(0.0, 30)) + + p = equator_arc1.intersection_point(arc) + assert np.allclose(p.vertices_in_degrees, (5.0, 0)) + + p = arc.intersection_point(equator_arc1) + assert np.allclose(p.vertices_in_degrees, (5.0, 0)) + + def test_touching_arc_midpoint_case(self): + """Test behaviour when touches in the mid of the arc.""" + arc1 = SArc(SPoint.from_degrees(0.0, 10.0), SPoint.from_degrees(0.0, 20.0)) + midpoint = arc1.midpoint() # (0, 15) + arc2 = SArc(midpoint, SPoint(-np.pi / 2, midpoint.lat)) + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, (0, 15)) diff --git a/pyresample/test/test_sgeom/test_point.py b/pyresample/test/test_sgeom/test_point.py index 4a6075e45..97182bf07 100644 --- a/pyresample/test/test_sgeom/test_point.py +++ b/pyresample/test/test_sgeom/test_point.py @@ -55,6 +55,18 @@ def test_creation_from_degrees(self): p2 = SPoint(np.deg2rad(lon), np.deg2rad(lat)) assert p1 == p2 + def test_is_pole(self): + """Test is pole method.""" + assert SPoint.from_degrees(0, -90).is_pole() + assert SPoint.from_degrees(0, 90).is_pole() + assert not SPoint.from_degrees(0, 89.99).is_pole() + + def test_is_on_equator(self): + """Test is on equator method.""" + assert SPoint.from_degrees(0, 0).is_on_equator() + assert SPoint.from_degrees(10, -0).is_on_equator() + assert not SPoint.from_degrees(0, 0.001).is_on_equator() + def test_vertices(self): """Test vertices property.""" lons = 0 @@ -78,15 +90,42 @@ def test_raise_error_if_multi_point(self): with pytest.raises(ValueError): SPoint(lons, lats) + def test_eq(self): + """Check the equality.""" + p1 = SPoint.from_degrees(0, 10) + p2 = SPoint.from_degrees(0, 10) + assert p1 == p2 + + def test_pole_equality(self): + """Check pole point equality.""" + p1 = SPoint.from_degrees(0, 90) + p2 = SPoint.from_degrees(10, 90) + assert p1 == p2 + + p1 = SPoint.from_degrees(0, -90) + p2 = SPoint.from_degrees(10, -90) + assert p1 == p2 + + p1 = SPoint.from_degrees(0, -90) + p2 = SPoint.from_degrees(0, 90) + assert not p1 == p2 + assert p1 != p2 + + def test_eq_antimeridian(self): + """Check the equality of a point on the antimeridian.""" + p = SPoint.from_degrees(-180, 0) + p1 = SPoint.from_degrees(180, 0) + assert p == p1 + def test_str(self): """Check the string representation.""" d = SPoint(1.0, 0.5) - self.assertEqual(str(d), "(1.0, 0.5)") + assert 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)") + assert repr(d) == 'SPoint: (1.0, 0.5)' def test_to_shapely(self): """Test conversion to shapely.""" @@ -148,8 +187,8 @@ def test_distance(self): p2 = SMultiPoint(lons, lats) d12 = p1.distance(p2) d21 = p2.distance(p1) - self.assertEqual(d12.shape, (2, 3)) - self.assertEqual(d21.shape, (3, 2)) + assert d12.shape == (2, 3) + assert d21.shape == (3, 2) res = np.array([[0., 1.57079633, 3.14159265], [3.14159265, 1.57079633, 0.]]) assert np.allclose(d12, res) @@ -169,19 +208,22 @@ def test_hdistance(self): p2 = SMultiPoint(lons, lats) d12 = p1.hdistance(p2) d21 = p2.hdistance(p1) - self.assertEqual(d12.shape, (2, 3)) - self.assertEqual(d21.shape, (3, 2)) + assert d12.shape == (2, 3) + assert 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_pole_equality(self): + """Check pole point equality.""" + p1 = SMultiPoint.from_degrees([0, 0, 1], [0, 0, 90]) + p2 = SMultiPoint.from_degrees([0, 0, 50], [0, 0, 90]) + assert p1 == p2 + + p1 = SMultiPoint.from_degrees([0, 0, 1], [0, 0, 90]) + p2 = SMultiPoint.from_degrees([0, 0, 1], [0, 0, -90]) + assert not p1 == p2 + assert p1 != p2 def test_eq_antimeridian(self): """Check the equality with longitudes at -180/180 degrees.""" @@ -192,29 +234,44 @@ def test_eq_antimeridian(self): p1 = SMultiPoint(lons1, lats) assert p == p1 + 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_neq(self): """Check the equality.""" lons = np.array([0, np.pi]) - lats = [-np.pi / 2, np.pi / 2] + lats = [0, np.pi / 2] p = SMultiPoint(lons, lats) p1 = SMultiPoint(lons + 0.1, lats) assert p != p1 + # Equality at the pole + lons = np.array([0, np.pi]) + lats = [-np.pi / 2, np.pi / 2] + p = SMultiPoint(lons, lats) + p1 = SMultiPoint(lons + 0.1, lats) + assert not 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) + assert 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) + expected_repr = 'SMultiPoint:\n [[ 0. -1.57079633]\n [-3.14159265 1.57079633]]' + assert repr(p) == expected_repr def test_to_shapely(self): """Test conversion to shapely.""" From 9d760d60ac8edc8806a3d35758cd19a2451a0f2e Mon Sep 17 00:00:00 2001 From: ghiggi Date: Thu, 17 Nov 2022 20:15:55 +0100 Subject: [PATCH 03/13] Fix circular import --- pyresample/future/spherical/arc.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyresample/future/spherical/arc.py b/pyresample/future/spherical/arc.py index 0ee06a88a..eb5211300 100644 --- a/pyresample/future/spherical/arc.py +++ b/pyresample/future/spherical/arc.py @@ -21,9 +21,10 @@ import pyproj from shapely.geometry import LineString -from pyresample.future.spherical import SPoint from pyresample.spherical import Arc +from .point import SPoint + EPSILON = 0.0000001 From bc2382c4a96ae7d9f7d04dcf9eb5bafdcef2b55f Mon Sep 17 00:00:00 2001 From: Gionata Ghiggi Date: Fri, 18 Nov 2022 08:00:53 +0100 Subject: [PATCH 04/13] Simplify is_pole method Co-authored-by: David Hoese --- pyresample/future/spherical/point.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py index dd3b84e93..7b7a86fb5 100644 --- a/pyresample/future/spherical/point.py +++ b/pyresample/future/spherical/point.py @@ -41,9 +41,7 @@ def from_degrees(cls, lon, lat): def is_pole(self): """Test if the point is on a pole.""" - if self.lat in [-np.pi / 2, np.pi / 2]: - return True - return False + return self.lat in (-np.pi / 2, np.pi / 2) def is_on_equator(self): """Test if the point is on the equator.""" From 24af697f621e4b24ae25509839cb9a1e4366d1bc Mon Sep 17 00:00:00 2001 From: Gionata Ghiggi Date: Fri, 18 Nov 2022 08:01:27 +0100 Subject: [PATCH 05/13] Simplify __eq__ method Co-authored-by: David Hoese --- pyresample/spherical.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyresample/spherical.py b/pyresample/spherical.py index de76494b5..255f71502 100644 --- a/pyresample/spherical.py +++ b/pyresample/spherical.py @@ -365,7 +365,7 @@ def __init__(self, start, end): def __eq__(self, other): """Check equality.""" - if self.start == other.start and self.end == other.end: + return self.start == other.start and self.end == other.end return True return False From 795fd6ecdb9230f3677d0fe8c0955e3d4fe3f846 Mon Sep 17 00:00:00 2001 From: Gionata Ghiggi Date: Fri, 18 Nov 2022 08:02:34 +0100 Subject: [PATCH 06/13] Simplify is_on_equator Co-authored-by: David Hoese --- pyresample/future/spherical/point.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py index 7b7a86fb5..c2330cdaa 100644 --- a/pyresample/future/spherical/point.py +++ b/pyresample/future/spherical/point.py @@ -45,9 +45,7 @@ def is_pole(self): def is_on_equator(self): """Test if the point is on the equator.""" - if self.lat == 0: - return True - return False + return self.lat == 0 def __str__(self): """Get simplified representation of lon/lat arrays in radians.""" From b8dfbea058f4b333d21802caecf9427472c9092e Mon Sep 17 00:00:00 2001 From: Gionata Ghiggi Date: Fri, 18 Nov 2022 08:02:52 +0100 Subject: [PATCH 07/13] Drop unittest Co-authored-by: David Hoese --- pyresample/test/test_sgeom/test_arc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyresample/test/test_sgeom/test_arc.py b/pyresample/test/test_sgeom/test_arc.py index 640b86af8..3a286e815 100644 --- a/pyresample/test/test_sgeom/test_arc.py +++ b/pyresample/test/test_sgeom/test_arc.py @@ -27,7 +27,7 @@ from pyresample.future.spherical import SArc, SPoint -class TestSArc(unittest.TestCase): +class TestSArc: """Test SArc class.""" # TODO: Fixtures defined here? From 5f42d031f529d02994525d19819b3545ca1b60b6 Mon Sep 17 00:00:00 2001 From: ghiggi Date: Fri, 18 Nov 2022 08:13:53 +0100 Subject: [PATCH 08/13] Clean up code --- pyresample/future/spherical/arc.py | 8 ++------ pyresample/spherical.py | 2 -- pyresample/test/test_sgeom/test_arc.py | 1 - pyresample/test/test_sgeom/test_point.py | 5 ++--- 4 files changed, 4 insertions(+), 12 deletions(-) diff --git a/pyresample/future/spherical/arc.py b/pyresample/future/spherical/arc.py index eb5211300..58b626c40 100644 --- a/pyresample/future/spherical/arc.py +++ b/pyresample/future/spherical/arc.py @@ -58,15 +58,11 @@ def __hash__(self): def is_on_equator(self): """Check if the SArc lies on the equator.""" - if self.start.lat == 0 and self.end.lat == 0: - return True - return False + return self.start.lat == 0 and self.end.lat == 0 def __eq__(self, other): """Check equality.""" - if self.start == other.start and self.end == other.end: - return True - return False + return self.start == other.start and self.end == other.end def reverse_direction(self): """Reverse SArc direction.""" diff --git a/pyresample/spherical.py b/pyresample/spherical.py index 255f71502..c4f2d15c3 100644 --- a/pyresample/spherical.py +++ b/pyresample/spherical.py @@ -366,8 +366,6 @@ def __init__(self, start, end): def __eq__(self, other): """Check equality.""" return self.start == other.start and self.end == other.end - return True - return False def __ne__(self, other): """Check not equal comparison.""" diff --git a/pyresample/test/test_sgeom/test_arc.py b/pyresample/test/test_sgeom/test_arc.py index 3a286e815..dba7884a4 100644 --- a/pyresample/test/test_sgeom/test_arc.py +++ b/pyresample/test/test_sgeom/test_arc.py @@ -18,7 +18,6 @@ """Define tests for the SArc class.""" import copy -import unittest import numpy as np import pytest diff --git a/pyresample/test/test_sgeom/test_point.py b/pyresample/test/test_sgeom/test_point.py index 97182bf07..1bc6c8183 100644 --- a/pyresample/test/test_sgeom/test_point.py +++ b/pyresample/test/test_sgeom/test_point.py @@ -15,7 +15,6 @@ # 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 @@ -23,7 +22,7 @@ from pyresample.future.spherical import SMultiPoint, SPoint -class TestSPoint(unittest.TestCase): +class TestSPoint: """Test SPoint.""" def test_latitude_validity(self): @@ -137,7 +136,7 @@ def test_to_shapely(self): assert shapely_point.equals_exact(spherical_point.to_shapely(), tolerance=1e-10) -class TestSMultiPoint(unittest.TestCase): +class TestSMultiPoint: """Test SMultiPoint.""" def test_single_point(self): From c086c9fbf8a35036dc6cee0ecf5ed67fec1b5f8d Mon Sep 17 00:00:00 2001 From: ghiggi Date: Fri, 18 Nov 2022 17:38:46 +0100 Subject: [PATCH 09/13] Add unit test coverage --- pyresample/future/spherical/arc.py | 92 ++- pyresample/future/spherical/point.py | 27 + pyresample/spherical.py | 4 +- pyresample/test/test_sgeom/test_arc.py | 899 ++++++++++++++++++++--- pyresample/test/test_sgeom/test_point.py | 46 ++ 5 files changed, 970 insertions(+), 98 deletions(-) diff --git a/pyresample/future/spherical/arc.py b/pyresample/future/spherical/arc.py index 58b626c40..3a11848bf 100644 --- a/pyresample/future/spherical/arc.py +++ b/pyresample/future/spherical/arc.py @@ -23,7 +23,7 @@ from pyresample.spherical import Arc -from .point import SPoint +from .point import SPoint, create_spherical_point EPSILON = 0.0000001 @@ -34,9 +34,11 @@ def _check_valid_arc(start_point, end_point): raise ValueError("An SArc can not be represented by the same start and end SPoint.") if start_point.is_pole() and end_point.is_pole(): raise ValueError("An SArc can not be uniquely defined between the two poles.") - if start_point.is_on_equator and end_point.is_on_equator and abs(start_point.lon - end_point.lon) == np.pi: + if start_point.is_on_equator() and end_point.is_on_equator() and abs(start_point.lon - end_point.lon) == np.pi: raise ValueError( "An SArc can not be uniquely defined on the equator if start and end points are 180 degrees apart.") + if start_point.get_antipode() == end_point: + raise ValueError("An SArc can not be uniquely defined between antipodal points.") class SArc(Arc): @@ -121,11 +123,12 @@ def intersection(self, other_arc): raise ValueError("'SArc.intersection' is deprecated. Use 'intersection_point' instead.") def intersection_point(self, other_arc): - """Compute the intersection point between two arcs. + """Compute the intersection point between two great circle arcs. - If arc and *other_arc* intersect, it returns the intersection SPoint. If arc and *other_arc* does not intersect, it returns None. - If same arc (also same direction), it returns None. + If arc and *other_arc* are the same (whatever direction), it returns None. + If arc and *other_arc* overlaps, within or contain, it returns None. + If arc and *other_arc* intersect, it returns the intersection SPoint. """ # If same arc (same direction), return None if self == other_arc: @@ -170,15 +173,77 @@ def to_shapely(self): start_coord, end_coord = self.vertices_in_degrees return LineString((start_coord.tolist()[0], end_coord.tolist()[0])) + def forward_points(self, distance, ellps="sphere"): + """Get points at given distance(s) from SArc end point in the forward direction. + + The distance value(s) must be provided in meters. + If the distance is positive, the point will be located outside the SArc. + If the distance is negative, the point will be located inside the SArc. + The function returns an SPoint or SMultiPoint. + """ + # Define geoid + geod = pyproj.Geod(ellps=ellps) + # Retrieve forward and back azimuth + fwd_az, back_az, _ = geod.inv(self.start.lon, self.start.lat, + self.end.lon, self.end.lat, radians=True) + # Retreve forward points + distance = np.array(distance) + lon, lat, back_az = geod.fwd(np.broadcast_to(np.array(self.end.lon), distance.shape), + np.broadcast_to(np.array(self.end.lat), distance.shape), + az=np.broadcast_to(fwd_az, distance.shape), + dist=distance, radians=True) + p = create_spherical_point(lon, lat) + return p + + def backward_points(self, distance, ellps="sphere"): + """Get points at given distance(s) from SArc start point in the backward direction. + + The distance value(s) must be provided in meters. + If the distance is positive, the point will be located outside the SArc. + If the distance is negative, the point will be located inside the SArc. + The function returns an SPoint or SMultiPoint. + """ + reverse_arc = self.reverse_direction() + return reverse_arc.forward_points(distance=distance, ellps=ellps) + + def extend(self, distance, direction="both", ellps="sphere"): + """Extend the SArc of a given distance in both, forward or backward directions. + + If the distance is positive, it extends the SArc. + If the distance is negative, it shortens the SArc. + """ + valid_direction = ["both", "forward", "backward"] + if direction not in valid_direction: + raise ValueError(f"Valid direction values are: {valid_direction}") + if direction in ["both", "forward"]: + end_point = self.forward_points(distance=distance, ellps="sphere") + else: + end_point = self.end + + if direction in ["both", "backward"]: + start_point = self.backward_points(distance=distance, ellps="sphere") + else: + start_point = self.start + arc = create_spherical_arcs(start_point, end_point) + return arc + + def shorten(self, distance, direction="both", ellps="sphere"): + """Short the SArc of a given distance in both, forward or backward directions. + + If the distance is positive, it shortens the SArc. + If the distance is negative, it extends the SArc. + """ + return self.extend(distance=-distance, direction=direction, ellps="sphere") + # def segmentize(self, npts=0, max_distance=0, ellips='sphere'): # """Segmentize the great-circle arc. # It returns an SLine. # npts or max_distance are mutually exclusively. Specify one of them. - # max_distance must be provided in kilometers. + # max_distance must be provided in meters. # """ # if npts != 0: - # npts = npts + 2 # + 2 to account for initial and terminus + # npts = npts + 2 # + 2 to account for initial and terminus points # geod = pyproj.Geod(ellps=ellips) # lon_start = self.start.lon # lon_end = self.end.lon @@ -203,3 +268,16 @@ def to_shapely(self): # def plot(self, *args, **kwargs): # """Convert to SLine.""" # self.to_line.plot(*args, **kwargs) + + +def create_spherical_arcs(start_point, end_point): + """Create a SArc or SArcs class depending on the number of points. + + If a Spoint is provided, it returns an SArc. + If a SMultiPoint is provided, it returns an SArcs. + """ + if isinstance(start_point, SPoint) and isinstance(end_point, SPoint): + arc = SArc(start_point, end_point) + else: + raise NotImplementedError("SArcs class is not yet available.") + return arc diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py index c2330cdaa..c5504e662 100644 --- a/pyresample/future/spherical/point.py +++ b/pyresample/future/spherical/point.py @@ -47,6 +47,12 @@ def is_on_equator(self): """Test if the point is on the equator.""" return self.lat == 0 + def get_antipode(self): + """Return the antipode point.""" + lat = - self.lat + lon = self.lon - np.pi + return SPoint(lon, lat) + def __str__(self): """Get simplified representation of lon/lat arrays in radians.""" return str((float(self.lon), float(self.lat))) @@ -84,6 +90,12 @@ def from_degrees(cls, lon, lat): """Create SMultiPoint from lon/lat coordinates in degrees.""" return cls(np.deg2rad(lon), np.deg2rad(lat)) + def get_antipodes(self): + """Return the antipodal points.""" + lat = - self.lat + lon = self.lon - np.pi + return SMultiPoint(lon, lat) + def __eq__(self, other): """Check equality.""" lon1 = self.lon.copy() @@ -108,3 +120,18 @@ def to_shapely(self): from shapely.geometry import MultiPoint point = MultiPoint(self.vertices_in_degrees) return point + + +def create_spherical_point(lon, lat): + """Create a SPoint or SMultiPoint class depending on the number of point coordinates. + + If a single point coordinate is provided, it returns an SPoint. + If multiple point coordinates are provided, it returns an SMultiPoint. + """ + lon = np.asarray(lon) + lat = np.asarray(lat) + if lon.ndim == 0: + p = SPoint(lon, lat) + else: + p = SMultiPoint(lon, lat) + return p diff --git a/pyresample/spherical.py b/pyresample/spherical.py index c4f2d15c3..d66377bea 100644 --- a/pyresample/spherical.py +++ b/pyresample/spherical.py @@ -372,7 +372,7 @@ def __ne__(self, other): return not self.__eq__(other) def __str__(self): - """Get simplified str representation.""" + """Get simplified string representation.""" return str(self.start) + " -> " + str(self.end) def __repr__(self): @@ -425,7 +425,7 @@ def angle(self, other_arc): return angle def intersections(self, other_arc): - """Compute the intersections points of the greats circles over which the arcs lies. + """Compute the intersections points of the great circles over which the arcs lies. A great circle divides the sphere in two equal hemispheres. """ diff --git a/pyresample/test/test_sgeom/test_arc.py b/pyresample/test/test_sgeom/test_arc.py index dba7884a4..64bc48063 100644 --- a/pyresample/test/test_sgeom/test_arc.py +++ b/pyresample/test/test_sgeom/test_arc.py @@ -23,13 +23,17 @@ import pytest from shapely.geometry import LineString -from pyresample.future.spherical import SArc, SPoint +from pyresample.future.spherical.arc import SArc +from pyresample.future.spherical.point import SMultiPoint, SPoint class TestSArc: """Test SArc class.""" - # TODO: Fixtures defined here? + # TODO in future + # - Fixtures defined here or outside the class + # - function that test intersection_point in all direction + equator_arc1 = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(50, 0.0)) equator_arc2 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(0.0, 0.0)) @@ -66,6 +70,13 @@ def test_unvalid_pole_to_pole_arc(self): with pytest.raises(ValueError): SArc(p1, p2) + def test_antipode_points_arcs(self): + """Check raise error if the points are the antipodes.""" + p1 = SPoint.from_degrees(45, 45.0) + p2 = SPoint.from_degrees(-135.0, -45) + with pytest.raises(ValueError): + SArc(p1, p2) + def test_is_on_equator_arc(self): """Check if the arc lies on the equator.""" equator_arc1 = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(50, 0.0)) @@ -133,6 +144,9 @@ def test_hash(self): arc = SArc(SPoint.from_degrees(-180.0, -90.0), SPoint.from_degrees(-180.0, 0.0)) assert hash(arc) == -3096892178517935054 + # ----------------------------------------------------------------------. + # Test SArc manipulations + def test_midpoint(self): """Test arc midpoint.""" start_point_vertices = np.deg2rad((-180.0, 0.0)) @@ -142,29 +156,400 @@ def test_midpoint(self): assert isinstance(midpoint, SPoint) assert np.allclose(midpoint.vertices_in_degrees[0], (-180, 45)) - def test_intersection_point(self): - """Test SArc(s) intersection point.""" + def test_arc_forward_points(self): + """Test SArc forward points.""" + arc = SArc(SPoint.from_degrees(10, 10), SPoint.from_degrees(20, 20)) + + # forward outside arc SPoint + p = arc.forward_points(distance=100) + assert isinstance(p, SPoint) + assert np.allclose(p.vertices_in_degrees, (20.00065043, 20.00065971)) + + # forward inside arc SPoint + p = arc.forward_points(distance=-100) + assert isinstance(p, SPoint) + assert np.allclose(p.vertices_in_degrees, (19.99934958, 19.99934029)) + + # Forward outside arc SMultiPoint + p = arc.forward_points(distance=[100 * 100, 1000 * 1000]) + expected_coords = np.array([[20.06506968, 20.06595905], + [26.81520629, 26.45975315]]) + assert isinstance(p, SMultiPoint) + assert np.allclose(p.vertices_in_degrees, expected_coords) + + # Forward inside arc SPoint + p = arc.forward_points(distance=[-100 * 100, -1000 * 1000]) + expected_coords = np.array([[19.93498483, 19.93401721], + [13.73281401, 13.30073419]]) + assert isinstance(p, SMultiPoint) + assert np.allclose(p.vertices_in_degrees, expected_coords) + + def test_arc_backward_points(self): + """Test SArc backward points.""" + arc = SArc(SPoint.from_degrees(10, 10), SPoint.from_degrees(20, 20)) + + # Backward outside arc SPoint + p = arc.backward_points(distance=100) + assert isinstance(p, SPoint) + assert np.allclose(p.vertices_in_degrees, (9.99934958, 9.99936874)) + + # Backward inside arc SPoint + p = arc.backward_points(distance=-100) + assert isinstance(p, SPoint) + assert np.allclose(p.vertices_in_degrees, (10.00065043, 10.00063126)) + + def test_arc_forward_points_at_equator(self): + """Test SArc forward points at equator.""" + arc = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(0.0, 0.0)) + + # Forward outside arc SPoint + p = arc.forward_points(distance=100) + assert isinstance(p, SPoint) + assert np.allclose(p.vertices_in_degrees, (0.00089932, 0.0)) + + # Forward inside arc SPoint + p = arc.forward_points(distance=-100) + assert isinstance(p, SPoint) + assert np.allclose(p.vertices_in_degrees, (-0.00089932, 0.0)) + + def test_arc_forward_points_at_pole(self): + """Test SArc forward points at pole.""" + # North pole + arc = SArc(SPoint.from_degrees(0.0, 50.0), SPoint.from_degrees(0.0, 90)) + + # Forward outside arc SPoint + p = arc.forward_points(distance=1000 * 1000) + assert isinstance(p, SPoint) + assert np.allclose(p.vertices_in_degrees, (-180, 81.00677971)) + + # Forward inside arc SPoint + p = arc.forward_points(distance=-1000 * 1000) + assert isinstance(p, SPoint) + assert np.allclose(p.vertices_in_degrees, (0, 81.00677971)) + + # South pole + arc = SArc(SPoint.from_degrees(0.0, 50.0), SPoint.from_degrees(0.0, -90)) + + # Forward outside arc SPoint + p = arc.forward_points(distance=1000 * 1000) + assert isinstance(p, SPoint) + assert np.allclose(p.vertices_in_degrees, (-180, -81.00677971)) + + # Forward inside arc SPoint + p = arc.forward_points(distance=-1000 * 1000) + assert isinstance(p, SPoint) + assert np.allclose(p.vertices_in_degrees, (0, -81.00677971)) + + def test_arc_extend(self): + """Test SArc extend method.""" + arc = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(20.0, 0.0)) + + # Test unvalid direction + with pytest.raises(ValueError): + arc.extend(distance=0, direction="bot") + + # Test backward + extended_arc = arc.extend(distance=1000 * 1000, direction="backward") + assert extended_arc.start == SPoint.from_degrees(-8.99322029, 0) + assert extended_arc.end == arc.end + + # Test forward + extended_arc = arc.extend(distance=1000 * 1000, direction="forward") + assert extended_arc.start == arc.start + assert extended_arc.end == SPoint.from_degrees(28.99322029, 0) + + # Test both + extended_arc = arc.extend(distance=1000 * 1000, direction="both") + assert extended_arc.start == SPoint.from_degrees(-8.99322029, 0) + assert extended_arc.end == SPoint.from_degrees(28.99322029, 0) + + def test_arc_shorten(self): + """Test SArc shorten method.""" + arc = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(20.0, 0.0)) + + # Test unvalid direction + with pytest.raises(ValueError): + arc.shorten(distance=0, direction="bot") + + # Test backward + shortened_arc = arc.shorten(distance=1000 * 1000, direction="backward") + assert shortened_arc.start == SPoint.from_degrees(8.99322029, 0) + assert shortened_arc.end == arc.end + + # Test forward + shortened_arc = arc.shorten(distance=1000 * 1000, direction="forward") + assert shortened_arc.start == arc.start + assert shortened_arc.end == SPoint.from_degrees(11.00677971, 0) + + # Test both + shortened_arc = arc.shorten(distance=1000 * 1000, direction="both") + assert shortened_arc.start == SPoint.from_degrees(8.99322029, 0) + assert shortened_arc.end == SPoint.from_degrees(11.00677971, 0) + # ----------------------------------------------------------------------. + # Test great circle intersections method + + def test_great_circles_intersections_of_crossing_arcs(self): + """Test behaviour when the 2 crossing arcs are not on the same great-circle.""" + # If the two arcs cross along the arc + # - one intersection point is the crossing point + # - the other intersection point is at the antipodes + arc1 = SArc(SPoint.from_degrees(0, 0), SPoint.from_degrees(10, 10)) arc2 = SArc(SPoint.from_degrees(0, 10), SPoint.from_degrees(10, 0)) - lon, lat = arc1.intersection_point(arc2) + crossing_point = (5., 5.0575149) - np.testing.assert_allclose(np.rad2deg(lon), 5) - assert np.rad2deg(lat).round(7) == round(5.0575148968282093, 7) + # - Changing the arc used to call the method changes the results order + p1, p2 = arc1._great_circle_intersections(arc2) + assert np.allclose(p1.vertices_in_degrees, (-175., -5.0575149)) + assert np.allclose(p2.vertices_in_degrees, crossing_point) - # Test when same SArc - arc1 = SArc(SPoint(0, 0), SPoint.from_degrees(10, 10)) - assert (arc1.intersection_point(arc1) is None) + p1, p2 = arc2._great_circle_intersections(arc1) + assert np.allclose(p1.vertices_in_degrees, crossing_point) + assert np.allclose(p2.vertices_in_degrees, (-175., -5.0575149)) - # Test intersecting SArc(s) - arc1 = SArc(SPoint.from_degrees(24.341215776575297, 44.987819588259327), - SPoint.from_degrees(18.842727517611817, 46.512483610284178)) - arc2 = SArc(SPoint.from_degrees(20.165961750361905, 46.177305385810541), - SPoint.from_degrees(20.253297585831707, 50.935830837274324)) - p = SPoint.from_degrees(20.165957021925202, 46.177022633103398) - assert arc1.intersection_point(arc2) == p - assert arc1.intersects(arc2) + # - If both arc are reversed, the results are the same + p1, p2 = arc1.reverse_direction()._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-175., -5.0575149)) + assert np.allclose(p2.vertices_in_degrees, crossing_point) + + p1, p2 = arc2.reverse_direction()._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, crossing_point) + assert np.allclose(p2.vertices_in_degrees, (-175., -5.0575149)) + + # - Reversing the direction of one arc changes the order of the points + p1, p2 = arc1._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, crossing_point) + assert np.allclose(p2.vertices_in_degrees, (-175., -5.0575149)) + + p1, p2 = arc2._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-175., -5.0575149)) + assert np.allclose(p2.vertices_in_degrees, crossing_point) + + def test_great_circles_intersections_on_touching_arcs(self): + """Test behaviour when the 2 arcs not on the same great-circle but touch at the extremes.""" + # If the two arcs are touching at the arc extremes + # - one intersection point is the touching point + # - the other intersection point is at the antipodes + + touching_point = (10.0, 20.0) + arc1 = SArc(SPoint.from_degrees(10.0, 10.0), SPoint.from_degrees(*touching_point)) + arc2 = SArc(SPoint.from_degrees(50.0, 60.0), SPoint.from_degrees(*touching_point)) + + # - Changing the arc used to call the method changes the results order + p1, p2 = arc1._great_circle_intersections(arc2) + assert np.allclose(p1.vertices_in_degrees, touching_point) + assert np.allclose(p2.vertices_in_degrees, (-170, -20)) + + p1, p2 = arc2._great_circle_intersections(arc1) + assert np.allclose(p1.vertices_in_degrees, (-170, -20)) + assert np.allclose(p2.vertices_in_degrees, touching_point) + + # - If both arc are reversed, the results are the same + p1, p2 = arc1.reverse_direction()._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, touching_point) + assert np.allclose(p2.vertices_in_degrees, (-170, -20)) + + p1, p2 = arc2.reverse_direction()._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-170, -20)) + assert np.allclose(p2.vertices_in_degrees, touching_point) + + # - Reversing the direction of one arc changes the order of the points + p1, p2 = arc1._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-170, -20)) + assert np.allclose(p2.vertices_in_degrees, touching_point) + + p1, p2 = arc2._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, touching_point) + assert np.allclose(p2.vertices_in_degrees, (-170, -20)) + + def test_great_circles_intersections_of_disjoint_arcs(self): + """Test behaviour when the 2 disjoint arcs are not on the same great-circle.""" + arc1 = SArc(SPoint.from_degrees(0, 0), SPoint.from_degrees(10, 10)) # north hemisphere + arc2 = SArc(SPoint.from_degrees(-10, -10), SPoint.from_degrees(-10, -20)) # south hemisphere + + # - Changing the arc used to call the method changes the results order + p1, p2 = arc1._great_circle_intersections(arc2) + assert np.allclose(p1.vertices_in_degrees, (170., 10.)) + assert np.allclose(p2.vertices_in_degrees, (-10., -10.)) + + p1, p2 = arc2._great_circle_intersections(arc1) + assert np.allclose(p1.vertices_in_degrees, (-10., -10.)) + assert np.allclose(p2.vertices_in_degrees, (170., 10.)) + + # - If both arc are reversed, the results are the same + p1, p2 = arc1.reverse_direction()._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (170., 10.)) + assert np.allclose(p2.vertices_in_degrees, (-10., -10.)) + + p1, p2 = arc2.reverse_direction()._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-10., -10.)) + assert np.allclose(p2.vertices_in_degrees, (170., 10.)) + + # - Reversing the direction of one arc changes the order of the points + p1, p2 = arc1._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-10., -10.)) + assert np.allclose(p2.vertices_in_degrees, (170., 10.)) + + p1, p2 = arc2._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (170., 10.)) + assert np.allclose(p2.vertices_in_degrees, (-10., -10.)) + + # - with arcs on the same great circle cases + def test_great_circles_intersections_of_disjoint_arcs_on_same_great_circle(self): + """Test behaviour when the 2 disjoint arcs (same great-circle case).""" + arc1 = SArc(SPoint.from_degrees(5, 1.0), SPoint.from_degrees(5.0, 20.0)) + start_point = arc1.forward_points(distance=1000 * 1000) + end_point = arc1.forward_points(distance=3000 * 1000) + arc2 = SArc(start_point, end_point) + + p1, p2 = arc1._great_circle_intersections(arc2) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + p1, p2 = arc2._great_circle_intersections(arc1) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + # - If both arc are reversed, the results are the same + p1, p2 = arc1.reverse_direction()._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + p1, p2 = arc2.reverse_direction()._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + # - Reversing the direction of one arc changes the order of the points + p1, p2 = arc2._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-180, 0)) + assert np.allclose(p2.vertices_in_degrees, (0, 0)) + + p1, p2 = arc1._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-180, 0)) + assert np.allclose(p2.vertices_in_degrees, (0, 0)) + + def test_great_circles_intersections_of_touching_equator_arcs(self): + """Test behaviour of touching equatorial arcs (same great-circle case).""" + arc1 = SArc(SPoint.from_degrees(1.0, 0.0), SPoint.from_degrees(50, 0.0)) + arc2 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(1.0, 0.0)) + + # - Changing the arc used to call the method does not change the results order + p1, p2 = arc1._great_circle_intersections(arc2) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + p1, p2 = arc2._great_circle_intersections(arc1) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + # - If both arc are reversed, the results are the same + p1, p2 = arc1.reverse_direction()._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + p1, p2 = arc2.reverse_direction()._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + # - Reversing the direction of one arc changes the order of the points + p1, p2 = arc2._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-180, 0)) + assert np.allclose(p2.vertices_in_degrees, (0, 0)) - # Test non-intersecting SArc(s) + p1, p2 = arc1._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-180, 0)) + assert np.allclose(p2.vertices_in_degrees, (0, 0)) + + def test_great_circles_intersections_on_overlapping_arcs(self): + """Test behaviour with 2 overlapping arcs (same great-circle case).""" + arc1 = SArc(SPoint.from_degrees(5, 1.0), SPoint.from_degrees(5.0, 20.0)) + arc2 = arc1.extend(distance=1000 * 1000, direction="forward") + + # - Changing the arc used to call the method does not change the results order + p1, p2 = arc1._great_circle_intersections(arc2) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + p1, p2 = arc2._great_circle_intersections(arc1) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + # - If both arc are reversed, the results are the same + p1, p2 = arc1.reverse_direction()._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + p1, p2 = arc2.reverse_direction()._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + # - Reversing the direction of one arc changes the order of the points + p1, p2 = arc2._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-180, 0)) + assert np.allclose(p2.vertices_in_degrees, (0, 0)) + + p1, p2 = arc1._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-180, 0)) + assert np.allclose(p2.vertices_in_degrees, (0, 0)) + + def test_great_circles_intersections_on_contain_within_arcs(self): + """Test behaviour with arcs contain/within each other (same great-circle case).""" + arc1 = SArc(SPoint.from_degrees(5, 1.0), SPoint.from_degrees(5.0, 20.0)) + arc2 = arc1.extend(distance=1000 * 1000, direction="both") + + # - Changing the arc used to call the method does not change the results order + p1, p2 = arc1._great_circle_intersections(arc2) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + p1, p2 = arc2._great_circle_intersections(arc1) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + # - If both arc are reversed, the results are the same + p1, p2 = arc1.reverse_direction()._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + p1, p2 = arc2.reverse_direction()._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + # - Reversing the direction of one arc changes the order of the points + p1, p2 = arc2._great_circle_intersections(arc1.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-180, 0)) + assert np.allclose(p2.vertices_in_degrees, (0, 0)) + + p1, p2 = arc1._great_circle_intersections(arc2.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (-180, 0)) + assert np.allclose(p2.vertices_in_degrees, (0, 0)) + + def test_great_circles_intersections_with_equal_arc(self): + """Test great_circles intersection points with equal arcs.""" + arc = SArc(SPoint.from_degrees(1.0, 20.0), SPoint.from_degrees(10.0, 40.0)) + + p1, p2 = arc._great_circle_intersections(arc) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + # - Reversing both arcs does not change the results + p1, p2 = arc.reverse_direction()._great_circle_intersections(arc.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + # - Reversing the direction of one arc does not change the results + p1, p2 = arc._great_circle_intersections(arc.reverse_direction()) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + p1, p2 = arc.reverse_direction()._great_circle_intersections(arc) + assert np.allclose(p1.vertices_in_degrees, (0, 0)) + assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + + # ----------------------------------------------------------------------. + # Test intersection_point method + def test_disjoint_arcs(self): + """Test disjoint arcs behaviour.""" arc1 = SArc(SPoint.from_degrees(-2.4982818108326734, 48.596644847869655), SPoint.from_degrees(-2.9571441235622835, 49.165688435261394)) arc2 = SArc(SPoint.from_degrees(-3.4976667413531688, 48.562704872921373), @@ -173,7 +558,14 @@ def test_intersection_point(self): assert not arc1.intersects(arc2) assert not bool(None) # this occurs within the intersects method - def test_disjoint_arcs_crossing_antimeridian(self): + def test_disjoint_equator_arcs(self): + """Test behaviour with disjoint arc along the equator.""" + equator_arc1 = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(50, 0.0)) + equator_arc3 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(-20, 0.0)) + p = equator_arc3.intersection_point(equator_arc1) + assert isinstance(p, type(None)) + + def test_disjoint_arcs_crossing_the_antimeridian(self): """Test SArc(s) intersection point when disjoint arcs cross the antimeridian.""" arc1 = SArc(SPoint.from_degrees(143.76, 0), SPoint.from_degrees(143.95, 7.33)) @@ -190,8 +582,62 @@ def test_disjoint_arcs_crossing_antimeridian(self): # Assert disjoint arcs returns None assert isinstance(point, type(None)) - def test_intersecting_arcs_crossing_antimeridian(self): - """Test SArc(s) intersection point when intersecting arcs cross the antimeridian.""" + def test_crossing_arcs(self): + """Test crossing arcs behaviour.""" + arc1 = SArc(SPoint.from_degrees(0, 0), SPoint.from_degrees(10, 10)) + arc2 = SArc(SPoint.from_degrees(0, 10), SPoint.from_degrees(10, 0)) + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, (5, 5.0575148968282093)) + + arc1 = SArc(SPoint.from_degrees(24.341215776575297, 44.987819588259327), + SPoint.from_degrees(18.842727517611817, 46.512483610284178)) + arc2 = SArc(SPoint.from_degrees(20.165961750361905, 46.177305385810541), + SPoint.from_degrees(20.253297585831707, 50.935830837274324)) + + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, (20.165957021925202, 46.177022633103398)) + assert arc1.intersects(arc2) + + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, (20.165957021925202, 46.177022633103398)) + + p = arc2.intersection_point(arc1) + assert np.allclose(p.vertices_in_degrees, (20.165957021925202, 46.177022633103398)) + + p = arc2.intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (20.165957021925202, 46.177022633103398)) + + p = arc1.intersection_point(arc2.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (20.165957021925202, 46.177022633103398)) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (20.165957021925202, 46.177022633103398)) + + def test_crossing_arcs_at_equator(self): + """Test arcs crossing at the equator.""" + arc1 = SArc(SPoint.from_degrees(0, 0), SPoint.from_degrees(50, 0)) + arc2 = SArc(SPoint.from_degrees(0, -10), SPoint.from_degrees(10, 10)) + + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, (5, 0)) + + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, (5, 0)) + + p = arc2.intersection_point(arc1) + assert np.allclose(p.vertices_in_degrees, (5, 0)) + + p = arc2.intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (5, 0)) + + p = arc1.intersection_point(arc2.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (5, 0)) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (5, 0)) + + def test_crossing_arcs_intersecting_at_the_antimeridian(self): + """Test SArc(s) intersection point when intersecting at the antimeridian.""" arc1 = SArc(SPoint.from_degrees(-180.0, -89.0), SPoint.from_degrees(-180.0, 89.0)) arc2 = SArc(SPoint.from_degrees(-171.03, -76.75), @@ -205,75 +651,232 @@ def test_intersecting_arcs_crossing_antimeridian(self): # Assert intersection result assert point == SPoint.from_degrees(-180.0, -74.78847163) - def test_great_circles_equator_arcs_cases(self): - """Test behaviour when 2 arcs are lying around the equator.""" - equator_arc1 = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(50, 0.0)) - equator_arc2 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(0.0, 0.0)) + def test_touching_arcs(self): + """Test behaviour when touches not at the arc extremes.""" + # Touching point at (15, 0) + arc1 = SArc(SPoint.from_degrees(0.0, 10.0), SPoint.from_degrees(0.0, 20.0)) + midpoint = arc1.midpoint() # (0, 15) + arc2 = SArc(midpoint, SPoint(-np.pi / 2, midpoint.lat)) + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, (0, 15)) - # - If one arc is reversed, the results order is swapped - p1, p2 = equator_arc1._great_circle_intersections(equator_arc2) - assert np.allclose(p1.vertices_in_degrees, (0, 0)) - assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + def test_touching_arcs_at_the_equator(self): + """Test arcs touching at the equator in the mid of one arc.""" + # Touching point at (25, 0) + arc1 = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(50, 0.0)) + arc2 = SArc(SPoint.from_degrees(25, 0.0), SPoint.from_degrees(0.0, 30)) - p1, p2 = equator_arc2._great_circle_intersections(equator_arc1) - assert np.allclose(p1.vertices_in_degrees, (0, 0)) - assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, (25.0, 0)) - # - If both arc are reversed, results are the same - p1, p2 = equator_arc2.reverse_direction()._great_circle_intersections(equator_arc1.reverse_direction()) - assert np.allclose(p1.vertices_in_degrees, (0, 0)) - assert np.allclose(p2.vertices_in_degrees, (-180, 0)) + p = arc2.intersection_point(arc1) + assert np.allclose(p.vertices_in_degrees, (25.0, 0)) - # - Reversing the direction of one arc, lead inversion of points - p1, p2 = equator_arc2._great_circle_intersections(equator_arc1.reverse_direction()) - assert np.allclose(p1.vertices_in_degrees, (-180, 0)) - assert np.allclose(p2.vertices_in_degrees, (0, 0)) + p = arc2.intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (25.0, 0)) - def test_intersection_equator_arcs_case(self): - """Test intersection point with 2 arcs lying around the equator.""" - equator_arc1 = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(50, 0.0)) - equator_arc2 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(0.0, 0.0)) - equator_arc3 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(-20, 0.0)) + p = arc1.intersection_point(arc2.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (25.0, 0)) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (25.0, 0)) - # Touching arcs - p = equator_arc2.intersection_point(equator_arc1) - assert np.allclose(p.vertices_in_degrees, (0, 0)) + def test_touching_arcs_at_extremes(self): + """Test arcs touching at the arc extremes (not on the same great-circle).""" + # Touching point at (10.0, 20.0) + arc1 = SArc(SPoint.from_degrees(10.0, 10.0), SPoint.from_degrees(10.0, 20.0)) + arc2 = SArc(SPoint.from_degrees(50.0, 60.0), SPoint.from_degrees(10.0, 20.0)) + + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, (10.0, 20.0)) + + p = arc2.intersection_point(arc1) + assert np.allclose(p.vertices_in_degrees, (10.0, 20.0)) + + p = arc2.intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (10.0, 20.0)) - p = equator_arc2.intersection_point(equator_arc1.reverse_direction()) - assert np.allclose(p.vertices_in_degrees, (0, 0)) + p = arc1.intersection_point(arc2.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (10.0, 20.0)) - # Overlapping arcs - p = equator_arc2.intersection_point(equator_arc3) + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (10.0, 20.0)) + + def test_touching_arcs_at_extremes_which_cross_antimeridian(self): + """Test arcs crossing the antimeridian, touching at the arc extremes.""" + # Touching point at (-175., 45) + arc1 = SArc(SPoint.from_degrees(-150, 10.0), SPoint.from_degrees(-175, 45)) + arc2 = SArc(SPoint.from_degrees(150, 10.0), SPoint.from_degrees(-175, 45)) + + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, (-175, 45)) + + p = arc2.intersection_point(arc1) + assert np.allclose(p.vertices_in_degrees, (-175, 45)) + + p = arc2.intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (-175, 45)) + + p = arc1.intersection_point(arc2.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (-175, 45)) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (-175, 45)) + + def test_touching_arcs_at_extremes_at_the_pole(self): + """Test arcs touching at the arc extremes at the pole.""" + # Touching point at (xxx, 90) + arc1 = SArc(SPoint.from_degrees(0.0, 10.0), SPoint.from_degrees(0.0, 90.0)) + arc2 = SArc(SPoint.from_degrees(60.0, 10.0), SPoint.from_degrees(10.0, 90.0)) + + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, (0., 90)) + + p = arc2.intersection_point(arc1) + assert np.allclose(p.vertices_in_degrees, (0., 90)) + + p = arc2.intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (0., 90)) + + p = arc1.intersection_point(arc2.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (0., 90)) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (0., 90)) + + def test_touching_arcs_at_extremes_at_the_antimeridian(self): + """Test arcs touching at the arc extremes at the antimeridian.""" + # Touching point at (-180., 45) + arc1 = SArc(SPoint.from_degrees(-150, 10.0), SPoint.from_degrees(-180, 45)) + arc2 = SArc(SPoint.from_degrees(150, 10.0), SPoint.from_degrees(180, 45)) + + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, (-180., 45)) + + p = arc2.intersection_point(arc1) + assert np.allclose(p.vertices_in_degrees, (-180., 45)) + + p = arc2.intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (-180., 45)) + + p = arc1.intersection_point(arc2.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (-180., 45)) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (-180., 45)) + + def test_touching_arcs_at_extremes_at_the_equator(self): + """Test arcs touching at the equator at the arc extremes.""" + # Touching point at the equator (5, 0) + arc1 = SArc(SPoint.from_degrees(5.0, 0.0), SPoint.from_degrees(0.0, 20.0)) + arc2 = SArc(SPoint.from_degrees(5.0, 0.0), SPoint.from_degrees(0.0, -20.0)) + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, (5.0, 0)) + + p = arc2.intersection_point(arc1) + assert np.allclose(p.vertices_in_degrees, (5.0, 0)) + + p = arc1.intersection_point(arc2.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (5.0, 0)) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, (5.0, 0)) + + # - with arcs on the same great circle + def test_touching_arcs_at_extremes_on_same_great_circle(self): + """Test arcs touching at the arc extremes with arc on the same great circle.""" + touching_point = 5.0, 20.0 + start_point = SPoint.from_degrees(6.0, 1.0) + end_point = SPoint.from_degrees(*touching_point) + arc1 = SArc(start_point, end_point) + forward_point = arc1.forward_points(distance=3000 * 1000) + arc2 = SArc(end_point, forward_point) + + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, touching_point) + + p = arc2.intersection_point(arc1) + assert np.allclose(p.vertices_in_degrees, touching_point) + + p = arc2.intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, touching_point) + + p = arc1.intersection_point(arc2.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, touching_point) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, touching_point) + + # Define another touching arc on the same great circle + # BUG [TO SOLVE] + start_point = SPoint.from_degrees(5.0, 1.0) + end_point = SPoint.from_degrees(5.0, 20.0) + arc1 = SArc(start_point, end_point) + forward_point = arc1.forward_points(distance=3000 * 1000) + arc2 = SArc(end_point, forward_point) + + # arc1._great_circle_intersections(arc2) # (0,0) , (-180,0) + + p = arc1.intersection_point(arc2) assert isinstance(p, type(None)) + assert not arc1.intersects(arc2) - # Disjoint arcs - p = equator_arc3.intersection_point(equator_arc1) + p = arc2.intersection_point(arc1) assert isinstance(p, type(None)) - # ---------------------------------------------------------------------. - # TEST FOR EDGE CASE SAME ARC + p = arc2.intersection_point(arc1.reverse_direction()) + assert isinstance(p, type(None)) - def test_intersection_same_arcs(self): - """Test behaviour when same arc.""" - nh_arc = SArc(SPoint.from_degrees(0.0, 20.0), SPoint.from_degrees(10.0, 40.0)) - # TODO BUG ! CHECK LEGAGY BEHAVIOUR - # - Define behaviour when same arc !!! - p1, p2 = nh_arc._great_circle_intersections(nh_arc) # (0,0), (-180, 0) - p = nh_arc.intersection_point(nh_arc) # None + p = arc1.intersection_point(arc2.reverse_direction()) assert isinstance(p, type(None)) - p1, p2 = nh_arc._great_circle_intersections(nh_arc.reverse_direction()) # (0,0), (-180, 0) - p = nh_arc.intersection_point(nh_arc.reverse_direction()) # None + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) assert isinstance(p, type(None)) - # ---------------------------------------------------------------------. - # TEST FOR TOUCHING CASES + def test_touching_arcs_at_extremes_with_equator_arcs(self): + """Test touching behaviour with arc lying around the equator.""" + touching_point = (0.0, 0.0) + arc1 = SArc(SPoint.from_degrees(*touching_point), SPoint.from_degrees(50, 0.0)) + arc2 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(*touching_point)) - def test_touching_arc_extremes_case(self): - """Test behaviour when two arc touches at the extremes.""" - # Test arc1.end == arc2.start - arc1 = SArc(SPoint.from_degrees(0.0, 10.0), SPoint.from_degrees(0.0, 20.0)) - arc2 = SArc(SPoint.from_degrees(0.0, 20.0), SPoint.from_degrees(0.0, 30.0)) + p = arc2.intersection_point(arc1) + assert np.allclose(p.vertices_in_degrees, touching_point) + + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, touching_point) + + p = arc2.intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, touching_point) + + p = arc1.intersection_point(arc2.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, touching_point) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, touching_point) + + touching_point = (-180.0, 0.0) + arc1 = SArc(SPoint.from_degrees(*touching_point), SPoint.from_degrees(50, 0.0)) + arc2 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(*touching_point)) + + p = arc2.intersection_point(arc1) + assert np.allclose(p.vertices_in_degrees, touching_point) + + p = arc1.intersection_point(arc2) + assert np.allclose(p.vertices_in_degrees, touching_point) + + p = arc2.intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, touching_point) + + p = arc1.intersection_point(arc2.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, touching_point) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert np.allclose(p.vertices_in_degrees, touching_point) + + # BUG [TO SOLVE] (WHATEVER LON OUTSIDE 0/180 return None) + touching_point = (-90.0, 0.0) + arc1 = SArc(SPoint.from_degrees(*touching_point), SPoint.from_degrees(50, 0.0)) + arc2 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(*touching_point)) p = arc1.intersection_point(arc2) assert isinstance(p, type(None)) @@ -285,27 +888,145 @@ def test_touching_arc_extremes_case(self): p = arc2.intersection_point(arc1.reverse_direction()) assert isinstance(p, type(None)) - # Test arc1.start == arc2.start - arc1 = SArc(SPoint.from_degrees(0.0, 10.0), SPoint.from_degrees(0.0, 20.0)) - arc2 = SArc(SPoint.from_degrees(0.0, 10.0), SPoint.from_degrees(0.0, -20.0)) + p = arc1.intersection_point(arc2.reverse_direction()) + assert isinstance(p, type(None)) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert isinstance(p, type(None)) + + # p = arc2.intersection_point(arc1) + # assert np.allclose(p.vertices_in_degrees, touching_point) + + # p = arc1.intersection_point(arc2) + # assert np.allclose(p.vertices_in_degrees, touching_point) + + # p = arc2.intersection_point(arc1.reverse_direction()) + # assert np.allclose(p.vertices_in_degrees, touching_point) + + # p = arc1.intersection_point(arc2.reverse_direction()) + # assert np.allclose(p.vertices_in_degrees, touching_point) + + # p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + # assert np.allclose(p.vertices_in_degrees, touching_point) + + def test_touching_arcs_at_extremes_with_meridian_arcs(self): + """Test arcs touching at the arc extremes with arcs pointing vertically to the pole.""" + # BUG [TO SOLVE] + + # Touching at point at (1, 20) + arc1 = SArc(SPoint.from_degrees(1.0, 10.0), SPoint.from_degrees(1.0, 20.0)) + arc2 = SArc(SPoint.from_degrees(1.0, 20.0), SPoint.from_degrees(1.0, 30.0)) + p = arc1.intersection_point(arc2) assert isinstance(p, type(None)) assert not arc1.intersects(arc2) - # TODO: BUG Bad behaviour at the equator. - equator_arc1 = SArc(SPoint.from_degrees(0.0, 0.0), SPoint.from_degrees(50, 0.0)) - arc = SArc(SPoint.from_degrees(5, 0.0), SPoint.from_degrees(0.0, 30)) + p = arc2.intersection_point(arc1) + assert isinstance(p, type(None)) - p = equator_arc1.intersection_point(arc) - assert np.allclose(p.vertices_in_degrees, (5.0, 0)) + p = arc2.intersection_point(arc1.reverse_direction()) + assert isinstance(p, type(None)) - p = arc.intersection_point(equator_arc1) - assert np.allclose(p.vertices_in_degrees, (5.0, 0)) + p = arc1.intersection_point(arc2.reverse_direction()) + assert isinstance(p, type(None)) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert isinstance(p, type(None)) + + # Touching at point at (1, -20) + arc1 = SArc(SPoint.from_degrees(1.0, -10.0), SPoint.from_degrees(1.0, -20.0)) + arc2 = SArc(SPoint.from_degrees(1.0, -20.0), SPoint.from_degrees(1.0, -30.0)) - def test_touching_arc_midpoint_case(self): - """Test behaviour when touches in the mid of the arc.""" - arc1 = SArc(SPoint.from_degrees(0.0, 10.0), SPoint.from_degrees(0.0, 20.0)) - midpoint = arc1.midpoint() # (0, 15) - arc2 = SArc(midpoint, SPoint(-np.pi / 2, midpoint.lat)) p = arc1.intersection_point(arc2) - assert np.allclose(p.vertices_in_degrees, (0, 15)) + assert isinstance(p, type(None)) + assert not arc1.intersects(arc2) + + p = arc2.intersection_point(arc1) + assert isinstance(p, type(None)) + + p = arc2.intersection_point(arc1.reverse_direction()) + assert isinstance(p, type(None)) + + p = arc1.intersection_point(arc2.reverse_direction()) + assert isinstance(p, type(None)) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert isinstance(p, type(None)) + + def test_contain_within_arcs(self): + """Test behaviour with arc contained/within each other.""" + arc1 = SArc(SPoint.from_degrees(5, 1.0), SPoint.from_degrees(5.0, 20.0)) + arc2 = arc1.extend(distance=1000 * 1000, direction="both") + + # Arc2 contain arc1 + p = arc1.intersection_point(arc2) + assert isinstance(p, type(None)) + + # Arc1 within arc1 + p = arc2.intersection_point(arc1) + assert isinstance(p, type(None)) + + def test_overlapping_arcs(self): + """Test behaviour with arc partially overlapping.""" + arc1 = SArc(SPoint.from_degrees(5, 1.0), SPoint.from_degrees(5.0, 20.0)) + arc2 = arc1.extend(distance=1000 * 1000, direction="forward") + + p = arc1.intersection_point(arc2) + assert isinstance(p, type(None)) + + p = arc2.intersection_point(arc1) + assert isinstance(p, type(None)) + + def test_overlapping_with_meridian_arcs(self): + """Test behaviour with arc overlapping along a meridian.""" + # Overlapping between (5, 10) and (5, 20) + arc1 = SArc(SPoint.from_degrees(5, 0.0), SPoint.from_degrees(5.0, 20.0)) + arc2 = SArc(SPoint.from_degrees(5, 10.0), SPoint.from_degrees(5.0, 30.0)) + + p = arc1.intersection_point(arc2) + assert isinstance(p, type(None)) + + p = arc2.intersection_point(arc1) + assert isinstance(p, type(None)) + + p = arc2.intersection_point(arc1.reverse_direction()) + assert isinstance(p, type(None)) + + p = arc1.intersection_point(arc2.reverse_direction()) + assert isinstance(p, type(None)) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert isinstance(p, type(None)) + + def test_overlapping_equator_arcs(self): + """Test behaviour with arc overlapping along the equator.""" + arc1 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(0.0, 0.0)) + arc2 = SArc(SPoint.from_degrees(-50, 0.0), SPoint.from_degrees(-20, 0.0)) + + p = arc1.intersection_point(arc2) + assert isinstance(p, type(None)) + + p = arc2.intersection_point(arc1) + assert isinstance(p, type(None)) + + p = arc2.intersection_point(arc1.reverse_direction()) + assert isinstance(p, type(None)) + + p = arc1.intersection_point(arc2.reverse_direction()) + assert isinstance(p, type(None)) + + p = arc2.reverse_direction().intersection_point(arc1.reverse_direction()) + assert isinstance(p, type(None)) + + def test_equal_arcs(self): + """Test intersection behaviour with equal arcs.""" + arc = SArc(SPoint.from_degrees(1.0, 20.0), SPoint.from_degrees(10.0, 40.0)) + + p = arc.intersection_point(arc) + assert isinstance(p, type(None)) + + p = arc.intersection_point(arc.reverse_direction()) + assert isinstance(p, type(None)) + + p = arc.reverse_direction().intersection_point(arc) + assert isinstance(p, type(None)) diff --git a/pyresample/test/test_sgeom/test_point.py b/pyresample/test/test_sgeom/test_point.py index 1bc6c8183..e272f0423 100644 --- a/pyresample/test/test_sgeom/test_point.py +++ b/pyresample/test/test_sgeom/test_point.py @@ -20,6 +20,7 @@ import pytest from pyresample.future.spherical import SMultiPoint, SPoint +from pyresample.future.spherical.point import create_spherical_point class TestSPoint: @@ -54,6 +55,16 @@ def test_creation_from_degrees(self): p2 = SPoint(np.deg2rad(lon), np.deg2rad(lat)) assert p1 == p2 + def test_spherical_point_creation(self): + """Test SPoint or SMultiPoint factory.""" + # SPoint + p = create_spherical_point(0, np.pi / 2) + assert isinstance(p, SPoint) + + # SMultiPoint + p = create_spherical_point([0, 0], [np.pi / 2, 0]) + assert isinstance(p, SMultiPoint) + def test_is_pole(self): """Test is pole method.""" assert SPoint.from_degrees(0, -90).is_pole() @@ -66,6 +77,31 @@ def test_is_on_equator(self): assert SPoint.from_degrees(10, -0).is_on_equator() assert not SPoint.from_degrees(0, 0.001).is_on_equator() + def test_antipodal_point(self): + """Test get_antipode method.""" + antipode = SPoint.from_degrees(0, 0).get_antipode() + assert np.allclose(antipode.vertices_in_degrees, (-180, 0)) + + antipode = SPoint.from_degrees(90, 0).get_antipode() + assert np.allclose(antipode.vertices_in_degrees, (-90, 0)) + + antipode = SPoint.from_degrees(-90, 0).get_antipode() + assert np.allclose(antipode.vertices_in_degrees, (90, 0)) + + antipode = SPoint.from_degrees(180, 0).get_antipode() + assert np.allclose(antipode.vertices_in_degrees, (0, 0)) + + antipode = SPoint.from_degrees(0, 90).get_antipode() + assert antipode == SPoint.from_degrees(-5, -90) + assert np.allclose(antipode.vertices_in_degrees, (-180, -90)) + + antipode = SPoint.from_degrees(0, -90).get_antipode() + assert antipode == SPoint.from_degrees(-5, 90) + assert np.allclose(antipode.vertices_in_degrees, (-180, 90)) + + antipode = SPoint.from_degrees(45, 45).get_antipode() + assert np.allclose(antipode.vertices_in_degrees, (-135., -45.)) + def test_vertices(self): """Test vertices property.""" lons = 0 @@ -158,6 +194,16 @@ def test_creation_from_degrees(self): p2 = SMultiPoint(np.deg2rad(lon), np.deg2rad(lat)) assert p1 == p2 + def test_antipodal_points(self): + """Test get_antipodes method.""" + lon = [0, 90, -90, 180, 0, 0, 45] + lat = [0, 0, 0, 0, 90, -90, 45] + antipodal_lon = [-180, -90, 90, 0, -5, -5, -135] + antipodal_lat = [0, 0, 0, 0, -90, 90, -45] + p = SMultiPoint.from_degrees(lon, lat) + antipodal_p = SMultiPoint.from_degrees(antipodal_lon, antipodal_lat) + assert p.get_antipodes() == antipodal_p + def test_vertices(self): """Test vertices property.""" lons = np.array([0, np.pi]) From 70be9777a04fbbeac403a2009d78507b71001786 Mon Sep 17 00:00:00 2001 From: ghiggi Date: Fri, 18 Nov 2022 18:19:59 +0100 Subject: [PATCH 10/13] Optimize and clean up arc intersection point --- pyresample/future/spherical/arc.py | 52 +++++++++++++++++++++--------- pyresample/spherical.py | 14 ++++---- 2 files changed, 43 insertions(+), 23 deletions(-) diff --git a/pyresample/future/spherical/arc.py b/pyresample/future/spherical/arc.py index 3a11848bf..19b6664c5 100644 --- a/pyresample/future/spherical/arc.py +++ b/pyresample/future/spherical/arc.py @@ -129,6 +129,8 @@ def intersection_point(self, other_arc): If arc and *other_arc* are the same (whatever direction), it returns None. If arc and *other_arc* overlaps, within or contain, it returns None. If arc and *other_arc* intersect, it returns the intersection SPoint. + If arc and *other_arc* touches in 1 point, it returns the intersection SPoint. + """ # If same arc (same direction), return None if self == other_arc: @@ -136,22 +138,40 @@ def intersection_point(self, other_arc): great_circles_intersection_spoints = self._great_circle_intersections(other_arc) - for spoint in great_circles_intersection_spoints: - a = self.start - b = self.end - c = other_arc.start - d = other_arc.end - - ab_dist = a.hdistance(b) - cd_dist = c.hdistance(d) - ap_dist = a.hdistance(spoint) - bp_dist = b.hdistance(spoint) - cp_dist = c.hdistance(spoint) - dp_dist = d.hdistance(spoint) - - if (((spoint in (a, b)) or (abs(ap_dist + bp_dist - ab_dist) < EPSILON)) and - ((spoint in (c, d)) or (abs(cp_dist + dp_dist - cd_dist) < EPSILON))): - return spoint + # Define arcs start and end points + a = self.start + b = self.end + c = other_arc.start + d = other_arc.end + + # Compute arc distance + ab_dist = a.hdistance(b) + cd_dist = c.hdistance(d) + + # If a great circle intersection point lies on one of the arc, + # it is the intersection point + + for point in great_circles_intersection_spoints: + # Distance from self arc start & end points to a great circle intersection point + ap_dist = a.hdistance(point) + bp_dist = b.hdistance(point) + + # Distance from other arcs start & end points to a great circle intersection point + cp_dist = c.hdistance(point) + dp_dist = d.hdistance(point) + + # Check if point is a start/end point of the arcs + point_is_on_self_extremities = point in (a, b) + point_is_on_other_arc_extremities = point in (c, d) + + # Check if point is on the arcs segment + point_is_on_self_segment = abs(ap_dist + bp_dist - ab_dist) < EPSILON + point_is_on_other_arc_segment = abs(cp_dist + dp_dist - cd_dist) < EPSILON + + point_is_on_self = point_is_on_self_extremities or point_is_on_self_segment + point_is_on_other_arc = point_is_on_other_arc_extremities or point_is_on_other_arc_segment + if point_is_on_self and point_is_on_other_arc: + return point return None def intersects(self, other_arc): diff --git a/pyresample/spherical.py b/pyresample/spherical.py index d66377bea..1c5bcf8fe 100644 --- a/pyresample/spherical.py +++ b/pyresample/spherical.py @@ -471,15 +471,15 @@ def intersection(self, other_arc): if self == other_arc: return None - for i in self.intersections(other_arc): - a__ = self.start - b__ = self.end - c__ = other_arc.start - d__ = other_arc.end + a__ = self.start + b__ = self.end + c__ = other_arc.start + d__ = other_arc.end - ab_ = a__.hdistance(b__) - cd_ = c__.hdistance(d__) + ab_ = a__.hdistance(b__) + cd_ = c__.hdistance(d__) + for i in self.intersections(other_arc): if (((i in (a__, b__)) or (abs(a__.hdistance(i) + b__.hdistance(i) - ab_) < EPSILON)) and ((i in (c__, d__)) or From 9f77f3462cfe8fa4d7cd36b661db58478fe3a1d5 Mon Sep 17 00:00:00 2001 From: ghiggi Date: Fri, 18 Nov 2022 18:31:25 +0100 Subject: [PATCH 11/13] Add SArc.__contains__ and refactor SArc.intersection_point --- pyresample/future/spherical/arc.py | 68 +++++++++++++++--------------- 1 file changed, 34 insertions(+), 34 deletions(-) diff --git a/pyresample/future/spherical/arc.py b/pyresample/future/spherical/arc.py index 19b6664c5..eead0031a 100644 --- a/pyresample/future/spherical/arc.py +++ b/pyresample/future/spherical/arc.py @@ -41,6 +41,30 @@ def _check_valid_arc(start_point, end_point): raise ValueError("An SArc can not be uniquely defined between antipodal points.") +def _is_point_on_arc(point, arc): + """Check if the point is on the arc.""" + # Define arc start and end points + start = arc.start + end = arc.end + + # Compute arc length + arc_length = start.hdistance(end) + + # Distance from self arc start & end points to a great circle intersection point + start_to_point_dist = start.hdistance(point) + end_to_point_dist = end.hdistance(point) + + # Check if point is a start/end point of the arcs + point_is_on_arc_extremities = point in (start, end) + + # Check if point is on the arc segment + point_is_on_arc_segment = abs(start_to_point_dist + end_to_point_dist - arc_length) < EPSILON + + # Assess if point is on the arc + point_is_on_arc = point_is_on_arc_extremities or point_is_on_arc_segment + return point_is_on_arc + + class SArc(Arc): """Object representing a great-circle arc between two points on a sphere. @@ -66,6 +90,13 @@ def __eq__(self, other): """Check equality.""" return self.start == other.start and self.end == other.end + def __contains__(self, point): + """Check if a point lies on the SArc.""" + if isinstance(point, SPoint): + return _is_point_on_arc(point, arc=self) + else: + raise NotImplementedError("SArc.__contains__ currently accept only SPoint objects.") + def reverse_direction(self): """Reverse SArc direction.""" return SArc(self.end, self.start) @@ -135,42 +166,11 @@ def intersection_point(self, other_arc): # If same arc (same direction), return None if self == other_arc: return None - + # Compute the great circle intersection points great_circles_intersection_spoints = self._great_circle_intersections(other_arc) - - # Define arcs start and end points - a = self.start - b = self.end - c = other_arc.start - d = other_arc.end - - # Compute arc distance - ab_dist = a.hdistance(b) - cd_dist = c.hdistance(d) - - # If a great circle intersection point lies on one of the arc, - # it is the intersection point - + # If a great circle intersection point lies on the arcs, it is the intersection point for point in great_circles_intersection_spoints: - # Distance from self arc start & end points to a great circle intersection point - ap_dist = a.hdistance(point) - bp_dist = b.hdistance(point) - - # Distance from other arcs start & end points to a great circle intersection point - cp_dist = c.hdistance(point) - dp_dist = d.hdistance(point) - - # Check if point is a start/end point of the arcs - point_is_on_self_extremities = point in (a, b) - point_is_on_other_arc_extremities = point in (c, d) - - # Check if point is on the arcs segment - point_is_on_self_segment = abs(ap_dist + bp_dist - ab_dist) < EPSILON - point_is_on_other_arc_segment = abs(cp_dist + dp_dist - cd_dist) < EPSILON - - point_is_on_self = point_is_on_self_extremities or point_is_on_self_segment - point_is_on_other_arc = point_is_on_other_arc_extremities or point_is_on_other_arc_segment - if point_is_on_self and point_is_on_other_arc: + if point in self and point in other_arc: return point return None From 6f338fdc52649c02abce817cf787e78b18dcf06e Mon Sep 17 00:00:00 2001 From: ghiggi Date: Fri, 18 Nov 2022 18:34:12 +0100 Subject: [PATCH 12/13] Improve comments of _is_point_on_arc --- pyresample/future/spherical/arc.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyresample/future/spherical/arc.py b/pyresample/future/spherical/arc.py index eead0031a..5739e6acf 100644 --- a/pyresample/future/spherical/arc.py +++ b/pyresample/future/spherical/arc.py @@ -50,11 +50,11 @@ def _is_point_on_arc(point, arc): # Compute arc length arc_length = start.hdistance(end) - # Distance from self arc start & end points to a great circle intersection point + # Distance from arc start & end points to the input point start_to_point_dist = start.hdistance(point) end_to_point_dist = end.hdistance(point) - # Check if point is a start/end point of the arcs + # Check if point is a start or end point of the arc point_is_on_arc_extremities = point in (start, end) # Check if point is on the arc segment From 59876497f6d8ff1ca99c9254c88f868b6ef69e79 Mon Sep 17 00:00:00 2001 From: ghiggi Date: Mon, 20 Feb 2023 07:57:25 +0100 Subject: [PATCH 13/13] Add antipode property --- pyresample/future/spherical/arc.py | 15 ++++++++++++++- pyresample/future/spherical/point.py | 6 ++++-- pyresample/spherical.py | 9 +++++++-- pyresample/test/test_sgeom/test_point.py | 20 ++++++++++---------- 4 files changed, 35 insertions(+), 15 deletions(-) diff --git a/pyresample/future/spherical/arc.py b/pyresample/future/spherical/arc.py index 5739e6acf..9696552c8 100644 --- a/pyresample/future/spherical/arc.py +++ b/pyresample/future/spherical/arc.py @@ -37,7 +37,7 @@ def _check_valid_arc(start_point, end_point): if start_point.is_on_equator() and end_point.is_on_equator() and abs(start_point.lon - end_point.lon) == np.pi: raise ValueError( "An SArc can not be uniquely defined on the equator if start and end points are 180 degrees apart.") - if start_point.get_antipode() == end_point: + if start_point.antipode == end_point: raise ValueError("An SArc can not be uniquely defined between antipodal points.") @@ -178,6 +178,19 @@ def intersects(self, other_arc): """Check if the current Sarc and another SArc intersect.""" return bool(self.intersection_point(other_arc)) + # def midpoint(self): + # """Return the SArc midpoint SPoint.""" + # # Retrieve start and end point in Cartesian coordinates + # start_xyz = self.start.to_cart().cart + # end_xyz = self.end.to_cart().cart + # # Find midpoint + # midpoint_xyz = (start_xyz + end_xyz) / 2.0 + # # Normalize + # midpoint_xyz = CCoordinate(midpoint_xyz).normalize() + # # Convert back to SPoint(s) + # midpoint = midpoint_xyz.to_spherical(future=True) + # return midpoint + def midpoint(self, ellips='sphere'): """Return the SArc midpoint SPoint.""" geod = pyproj.Geod(ellps=ellips) diff --git a/pyresample/future/spherical/point.py b/pyresample/future/spherical/point.py index c5504e662..a652429bd 100644 --- a/pyresample/future/spherical/point.py +++ b/pyresample/future/spherical/point.py @@ -47,7 +47,8 @@ def is_on_equator(self): """Test if the point is on the equator.""" return self.lat == 0 - def get_antipode(self): + @property + def antipode(self): """Return the antipode point.""" lat = - self.lat lon = self.lon - np.pi @@ -90,7 +91,8 @@ def from_degrees(cls, lon, lat): """Create SMultiPoint from lon/lat coordinates in degrees.""" return cls(np.deg2rad(lon), np.deg2rad(lat)) - def get_antipodes(self): + @property + def antipodes(self): """Return the antipodal points.""" lat = - self.lat lon = self.lon - np.pi diff --git a/pyresample/spherical.py b/pyresample/spherical.py index 1c5bcf8fe..cfa8b637e 100644 --- a/pyresample/spherical.py +++ b/pyresample/spherical.py @@ -349,12 +349,17 @@ def __rmul__(self, other): """Multiply.""" return self.__mul__(other) - def to_spherical(self): + def to_spherical(self, future=False): """Convert to SPoint/SMultiPoint object.""" # TODO: this in future should point to SPoint or SMultiPoint + from pyresample.future.spherical.point import create_spherical_point + lon = np.arctan2(self.cart[..., 1], self.cart[..., 0]) lat = np.arcsin(self.cart[..., 2]) - return SCoordinate(lon, lat) + if not future: + return SCoordinate(lon, lat) + else: + return create_spherical_point(lon, lat) class Arc(object): diff --git a/pyresample/test/test_sgeom/test_point.py b/pyresample/test/test_sgeom/test_point.py index e272f0423..9f6343f41 100644 --- a/pyresample/test/test_sgeom/test_point.py +++ b/pyresample/test/test_sgeom/test_point.py @@ -78,28 +78,28 @@ def test_is_on_equator(self): assert not SPoint.from_degrees(0, 0.001).is_on_equator() def test_antipodal_point(self): - """Test get_antipode method.""" - antipode = SPoint.from_degrees(0, 0).get_antipode() + """Test antipode property.""" + antipode = SPoint.from_degrees(0, 0).antipode assert np.allclose(antipode.vertices_in_degrees, (-180, 0)) - antipode = SPoint.from_degrees(90, 0).get_antipode() + antipode = SPoint.from_degrees(90, 0).antipode assert np.allclose(antipode.vertices_in_degrees, (-90, 0)) - antipode = SPoint.from_degrees(-90, 0).get_antipode() + antipode = SPoint.from_degrees(-90, 0).antipode assert np.allclose(antipode.vertices_in_degrees, (90, 0)) - antipode = SPoint.from_degrees(180, 0).get_antipode() + antipode = SPoint.from_degrees(180, 0).antipode assert np.allclose(antipode.vertices_in_degrees, (0, 0)) - antipode = SPoint.from_degrees(0, 90).get_antipode() + antipode = SPoint.from_degrees(0, 90).antipode assert antipode == SPoint.from_degrees(-5, -90) assert np.allclose(antipode.vertices_in_degrees, (-180, -90)) - antipode = SPoint.from_degrees(0, -90).get_antipode() + antipode = SPoint.from_degrees(0, -90).antipode assert antipode == SPoint.from_degrees(-5, 90) assert np.allclose(antipode.vertices_in_degrees, (-180, 90)) - antipode = SPoint.from_degrees(45, 45).get_antipode() + antipode = SPoint.from_degrees(45, 45).antipode assert np.allclose(antipode.vertices_in_degrees, (-135., -45.)) def test_vertices(self): @@ -195,14 +195,14 @@ def test_creation_from_degrees(self): assert p1 == p2 def test_antipodal_points(self): - """Test get_antipodes method.""" + """Test antipodes property.""" lon = [0, 90, -90, 180, 0, 0, 45] lat = [0, 0, 0, 0, 90, -90, 45] antipodal_lon = [-180, -90, 90, 0, -5, -5, -135] antipodal_lat = [0, 0, 0, 0, -90, 90, -45] p = SMultiPoint.from_degrees(lon, lat) antipodal_p = SMultiPoint.from_degrees(antipodal_lon, antipodal_lat) - assert p.get_antipodes() == antipodal_p + assert p.antipodes == antipodal_p def test_vertices(self): """Test vertices property."""