Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Future Spherical Class] Add SArc #478

Open
wants to merge 14 commits into
base: main
Choose a base branch
from
100 changes: 87 additions & 13 deletions pyresample/future/spherical/arc.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@

from pyresample.spherical import Arc

from .point import SPoint
from .point import SPoint, create_spherical_point

EPSILON = 0.0000001

Expand All @@ -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:
ghiggi marked this conversation as resolved.
Show resolved Hide resolved
raise ValueError("An SArc can not be uniquely defined between antipodal points.")


class SArc(Arc):
Expand All @@ -58,15 +60,11 @@ def __hash__(self):

def is_on_equator(self):
ghiggi marked this conversation as resolved.
Show resolved Hide resolved
"""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."""
Expand Down Expand Up @@ -125,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:
Expand Down Expand Up @@ -174,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
Expand All @@ -207,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
27 changes: 27 additions & 0 deletions pyresample/future/spherical/point.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)))
Expand Down Expand Up @@ -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()
Expand All @@ -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
6 changes: 2 additions & 4 deletions pyresample/spherical.py
Original file line number Diff line number Diff line change
Expand Up @@ -366,15 +366,13 @@ 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."""
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):
Expand Down Expand Up @@ -427,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.
"""
Expand Down
Loading