Skip to content

Commit

Permalink
Add planetocentric options to BasicGeodeticRepresentation. Add tests …
Browse files Browse the repository at this point in the history 
…and documentation.
  • Loading branch information
@cmarmo
cmarmo committed May 25, 2023
1 parent cf34e78 commit 65af784
Show file tree
Hide file tree
Showing 2 changed files with 225 additions and 11 deletions.
85 changes: 78 additions & 7 deletions astropy/coordinates/representation/geodetic.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
# Licensed under a 3-clause BSD style license - see LICENSE.rst

import numpy as np
import erfa

from astropy import units as u
Expand All @@ -10,6 +11,16 @@
from .cartesian import CartesianRepresentation


def _compute_longitude(longitude, positive_longitude, wrap_angle):
if positive_longitude == "east":
orient = 1
elif positive_longitude == "west":
orient = -1
return Longitude(
orient * longitude, u.deg, wrap_angle=wrap_angle * u.deg, copy=False
)


ELLIPSOIDS = {}
"""Available ellipsoids (defined in erfam.h, with numbers exposed in erfa)."""
# Note: they get filled by the creation of the geodetic classes.
Expand Down Expand Up @@ -45,9 +56,16 @@ class BaseGeodeticRepresentation(BaseRepresentation):
to quantities holding correct values (with units of length and dimensionless,
respectively), or alternatively an ``_ellipsoid`` attribute to the relevant ERFA
index (as passed in to `erfa.eform`).
Longitudes are east positive and span from -180 to 180 degrees by default.
They can be made west positive setting `_positive_longitude='west'`, or spanning
from 0 to 360 degrees setting `_wrap_angle=360`.
Planetocentric latitudes can be obtained setting `_ographic=False`.
"""

attr_classes = {"lon": Longitude, "lat": Latitude, "height": u.Quantity}
_positive_longitude = "east"
_wrap_angle = 180
_ographic = True

def __init_subclass__(cls, **kwargs):
if "_ellipsoid" in cls.__dict__:
Expand All @@ -62,6 +80,16 @@ def __init_subclass__(cls, **kwargs):
raise AttributeError(
f"{cls.__name__} requires '_ellipsoid' or '_equatorial_radius' and '_flattening'."
)
if cls._positive_longitude not in ["east", "west"]:
raise ValueError(
f"Invalid argument '{cls._positive_longitude}' for '_positive_logitude' "
"attribute: valid argument are 'east' or 'west'."
)
if cls._wrap_angle not in [180, 360]:
raise ValueError(
f"Invalid argument '{cls._wrap_angle}' for '_wrap_angle' attribute: "
"valid arguments are 180 or 360."
)
super().__init_subclass__(**kwargs)

def __init__(self, lon, lat=None, height=None, copy=True):
Expand All @@ -79,13 +107,34 @@ def to_cartesian(self):
Converts geodetic coordinates to 3D rectangular (geocentric)
cartesian coordinates.
"""
xyz = erfa.gd2gce(
self._equatorial_radius,
self._flattening,
self.lon,
self.lat,
self.height,
)
lon = _compute_longitude(self.lon, self._positive_longitude, self._wrap_angle)

if self._ographic:
xyz = erfa.gd2gce(
self._equatorial_radius,
self._flattening,
lon,
self.lat,
self.height,
)
else:
x_spheroid = self._equatorial_radius * np.cos(self.lat) * np.cos(lon)
y_spheroid = self._equatorial_radius * np.cos(self.lat) * np.sin(lon)
z_spheroid = (self._equatorial_radius * (1 - self._flattening)) * np.sin(
self.lat
)
r = (
np.sqrt(
x_spheroid * x_spheroid
+ y_spheroid * y_spheroid
+ z_spheroid * z_spheroid
)
+ self.height
)
x = r * np.cos(self.lon) * np.cos(self.lat)
y = r * np.sin(self.lon) * np.cos(self.lat)
z = r * np.sin(self.lat)
xyz = np.stack([x, y, z], axis=1) << u.m
return CartesianRepresentation(xyz, xyz_axis=-1, copy=False)

@classmethod
Expand All @@ -94,9 +143,31 @@ def from_cartesian(cls, cart):
Converts 3D rectangular cartesian coordinates (assumed geocentric) to
geodetic coordinates.
"""
# Compute geodetic/planetodetic angles
lon, lat, height = erfa.gc2gde(
cls._equatorial_radius, cls._flattening, cart.get_xyz(xyz_axis=-1)
)
lon = _compute_longitude(lon, cls._positive_longitude, cls._wrap_angle)
if not cls._ographic:
# Compute planetocentric angles
xyz = cart.get_xyz()
# Compute planetocentric latitude
p = np.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1])
d = np.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1] + xyz[2] * xyz[2])
lat = np.where(
p != 0.0,
np.arctan(xyz[2] / p),
np.sign(xyz[2]) * 0.5 * np.pi * u.radian,
)
p_spheroid = cls._equatorial_radius * np.cos(lat)
z_spheroid = (cls._equatorial_radius * (1 - cls._flattening)) * np.sin(lat)
r_spheroid = np.sqrt(p_spheroid * p_spheroid + z_spheroid * z_spheroid)
height = np.where(
p_spheroid != 0.0,
(d - r_spheroid),
(np.abs(xyz[2]) - np.abs(z_spheroid)),
)

return cls(lon, lat, height, copy=False)


Expand Down
151 changes: 147 additions & 4 deletions astropy/coordinates/tests/test_geodetic_representations.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,45 @@ class CustomGeodetic(BaseGeodeticRepresentation):
_equatorial_radius = 4000000.0 * u.m


class IAUMARS2000GeodeticRepresentationEast180(BaseGeodeticRepresentation):
_equatorial_radius = 3396190.0 * u.m
_flattening = 0.5886007555512007 * u.percent


class IAUMARS2000GeodeticRepresentationWest180(BaseGeodeticRepresentation):
_equatorial_radius = 3396190.0 * u.m
_flattening = 0.5886007555512007 * u.percent
_positive_longitude = "west"


class IAUMARS2000GeodeticRepresentationEast360(BaseGeodeticRepresentation):
_equatorial_radius = 3396190.0 * u.m
_flattening = 0.5886007555512007 * u.percent
_wrap_angle = 360


class IAUMARS2000GeodeticRepresentationWest360(BaseGeodeticRepresentation):
_equatorial_radius = 3396190.0 * u.m
_flattening = 0.5886007555512007 * u.percent
_positive_longitude = "west"
_wrap_angle = 360


class IAUMARS2000GeodeticRepresentationOcentric(BaseGeodeticRepresentation):
_equatorial_radius = 3396190.0 * u.m
_flattening = 0.5886007555512007 * u.percent
_ographic = False


@pytest.mark.parametrize(
"geodeticrepresentation", [CustomGeodetic, WGS84GeodeticRepresentation]
"geodeticrepresentation",
[
CustomGeodetic,
WGS84GeodeticRepresentation,
IAUMARS2000GeodeticRepresentationEast360,
IAUMARS2000GeodeticRepresentationWest360,
IAUMARS2000GeodeticRepresentationOcentric,
],
)
def test_cartesian_geodetic_roundtrip(geodeticrepresentation):
# Test array-valued input in the process.
Expand All @@ -48,7 +85,14 @@ def test_cartesian_geodetic_roundtrip(geodeticrepresentation):


@pytest.mark.parametrize(
"geodeticrepresentation", [CustomGeodetic, WGS84GeodeticRepresentation]
"geodeticrepresentation",
[
CustomGeodetic,
WGS84GeodeticRepresentation,
IAUMARS2000GeodeticRepresentationEast360,
IAUMARS2000GeodeticRepresentationWest360,
IAUMARS2000GeodeticRepresentationOcentric,
],
)
def test_geodetic_cartesian_roundtrip(geodeticrepresentation):
initial_geodetic = geodeticrepresentation(
Expand Down Expand Up @@ -149,7 +193,7 @@ class InvalidCustomGeodeticEllipsoid(BaseGeodeticRepresentation):
_ellipsoid = "foo"

assert "foo" not in ELLIPSOIDS
assert "customgeodeticellipsoiderror" not in REPRESENTATION_CLASSES
assert "invalidcustomgeodeticellipsoid" not in REPRESENTATION_CLASSES


def test_geodetic_subclass_missing_equatorial_radius():
Expand All @@ -159,4 +203,103 @@ def test_geodetic_subclass_missing_equatorial_radius():
class MissingCustomGeodeticAttribute(BaseGeodeticRepresentation):
_flattening = 0.075 * u.dimensionless_unscaled

assert "customgeodeticerror" not in REPRESENTATION_CLASSES
assert "missingcustomgeodeticattribute" not in REPRESENTATION_CLASSES


def test_geodetic_subclass_bad_attributes():
msg = "Invalid argument 'foo' for '_positive_logitude'"
with pytest.raises(ValueError, match=msg):

class BadPositiveLongitudeCustomGeodetic(BaseGeodeticRepresentation):
_equatorial_radius = 3000.0 * u.km
_flattening = 0.075 * u.dimensionless_unscaled
_positive_longitude = "foo"

assert "badpositivelongitudecustomgeodetic" not in REPRESENTATION_CLASSES

msg = "Invalid argument 'foo' for '_wrap_angle'"
with pytest.raises(ValueError, match=msg):

class BadWrapAngleCustomGeodetic(BaseGeodeticRepresentation):
_equatorial_radius = 3000.0 * u.km
_flattening = 0.075 * u.dimensionless_unscaled
_wrap_angle = "foo"

assert "badwrapanglecustomgeodetic" not in REPRESENTATION_CLASSES


def test_from_cartesian_positive_direction_180_360():
# Test array-valued input in the process.
initial_cartesian = CartesianRepresentation(
x=[1, 3000.0, 2600.0] * u.km,
y=[7000.0, 4.0, -5000.0] * u.km,
z=[5.0, 6000.0, 3000.0] * u.km,
)

assert_quantity_allclose(
(
IAUMARS2000GeodeticRepresentationEast180.from_representation(
initial_cartesian
).lon.value
)
* u.deg,
-(
IAUMARS2000GeodeticRepresentationWest180.from_representation(
initial_cartesian
).lon.value
)
* u.deg,
)

assert_quantity_allclose(
(
IAUMARS2000GeodeticRepresentationEast180.from_representation(
initial_cartesian
).lon.value
)
% 180
* u.deg,
(
IAUMARS2000GeodeticRepresentationEast360.from_representation(
initial_cartesian
).lon.value
)
% 180
* u.deg,
)

assert_quantity_allclose(
(
IAUMARS2000GeodeticRepresentationEast180.from_representation(
initial_cartesian
).lon.value
)
% 180
* u.deg,
(
360
- (
IAUMARS2000GeodeticRepresentationWest360.from_representation(
initial_cartesian
).lon.value
)
)
% 180
* u.deg,
)

assert_quantity_allclose(
(
IAUMARS2000GeodeticRepresentationEast360.from_representation(
initial_cartesian
).lon.value
)
* u.deg,
(
360
- IAUMARS2000GeodeticRepresentationWest360.from_representation(
initial_cartesian
).lon.value
)
* u.deg,
)

0 comments on commit 65af784

Please sign in to comment.