Skip to content

Commit

Permalink
Merge branch 'astropy:main' into issue14808
Browse files Browse the repository at this point in the history
  • Loading branch information
BillCleveland-USRA committed May 11, 2023
2 parents 9e3bfd1 + cc73b24 commit f5f428e
Show file tree
Hide file tree
Showing 10 changed files with 269 additions and 125 deletions.
96 changes: 1 addition & 95 deletions astropy/coordinates/earth.py
Expand Up @@ -8,39 +8,29 @@
import urllib.request
from warnings import warn

import erfa
import numpy as np

from astropy import constants as consts
from astropy import units as u
from astropy.units.quantity import QuantityInfoBase
from astropy.utils import data
from astropy.utils.decorators import format_doc
from astropy.utils.exceptions import AstropyUserWarning

from .angles import Angle, Latitude, Longitude
from .errors import UnknownSiteException
from .matrix_utilities import matrix_transpose
from .representation import (
BaseRepresentation,
CartesianDifferential,
CartesianRepresentation,
)
from .representation.geodetic import ELLIPSOIDS

__all__ = [
"EarthLocation",
"BaseGeodeticRepresentation",
"WGS84GeodeticRepresentation",
"WGS72GeodeticRepresentation",
"GRS80GeodeticRepresentation",
]

GeodeticLocation = collections.namedtuple("GeodeticLocation", ["lon", "lat", "height"])

ELLIPSOIDS = {}
"""Available ellipsoids (defined in erfam.h, with numbers exposed in erfa)."""
# Note: they get filled by the creation of the geodetic classes.

OMEGA_EARTH = (1.002_737_811_911_354_48 * u.cycle / u.day).to(
1 / u.s, u.dimensionless_angles()
)
Expand Down Expand Up @@ -905,87 +895,3 @@ def _to_value(self, unit, equivalencies=[]):
equivalencies = self._equivalencies
new_array = self.unit.to(unit, array_view, equivalencies=equivalencies)
return new_array.view(self.dtype).reshape(self.shape)


geodetic_base_doc = """{__doc__}
Parameters
----------
lon, lat : angle-like
The longitude and latitude of the point(s), in angular units. The
latitude should be between -90 and 90 degrees, and the longitude will
be wrapped to an angle between 0 and 360 degrees. These can also be
instances of `~astropy.coordinates.Angle` and either
`~astropy.coordinates.Longitude` not `~astropy.coordinates.Latitude`,
depending on the parameter.
height : `~astropy.units.Quantity` ['length']
The height to the point(s).
copy : bool, optional
If `True` (default), arrays will be copied. If `False`, arrays will
be references, though possibly broadcast to ensure matching shapes.
"""


@format_doc(geodetic_base_doc)
class BaseGeodeticRepresentation(BaseRepresentation):
"""Base geodetic representation."""

attr_classes = {"lon": Longitude, "lat": Latitude, "height": u.Quantity}

def __init_subclass__(cls, **kwargs):
super().__init_subclass__(**kwargs)
if "_ellipsoid" in cls.__dict__:
ELLIPSOIDS[cls._ellipsoid] = cls

def __init__(self, lon, lat=None, height=None, copy=True):
if height is None and not isinstance(lon, self.__class__):
height = 0 << u.m

super().__init__(lon, lat, height, copy=copy)
if not self.height.unit.is_equivalent(u.m):
raise u.UnitTypeError(
f"{self.__class__.__name__} requires height with units of length."
)

def to_cartesian(self):
"""
Converts WGS84 geodetic coordinates to 3D rectangular (geocentric)
cartesian coordinates.
"""
xyz = erfa.gd2gc(
getattr(erfa, self._ellipsoid), self.lon, self.lat, self.height
)
return CartesianRepresentation(xyz, xyz_axis=-1, copy=False)

@classmethod
def from_cartesian(cls, cart):
"""
Converts 3D rectangular cartesian coordinates (assumed geocentric) to
WGS84 geodetic coordinates.
"""
lon, lat, height = erfa.gc2gd(
getattr(erfa, cls._ellipsoid), cart.get_xyz(xyz_axis=-1)
)
return cls(lon, lat, height, copy=False)


@format_doc(geodetic_base_doc)
class WGS84GeodeticRepresentation(BaseGeodeticRepresentation):
"""Representation of points in WGS84 3D geodetic coordinates."""

_ellipsoid = "WGS84"


@format_doc(geodetic_base_doc)
class WGS72GeodeticRepresentation(BaseGeodeticRepresentation):
"""Representation of points in WGS72 3D geodetic coordinates."""

_ellipsoid = "WGS72"


@format_doc(geodetic_base_doc)
class GRS80GeodeticRepresentation(BaseGeodeticRepresentation):
"""Representation of points in GRS80 3D geodetic coordinates."""

_ellipsoid = "GRS80"
10 changes: 10 additions & 0 deletions astropy/coordinates/representation/__init__.py
Expand Up @@ -6,6 +6,12 @@
from .base import BaseRepresentationOrDifferential, BaseRepresentation, BaseDifferential
from .cartesian import CartesianRepresentation, CartesianDifferential
from .cylindrical import CylindricalRepresentation, CylindricalDifferential
from .geodetic import (
BaseGeodeticRepresentation,
WGS84GeodeticRepresentation,
WGS72GeodeticRepresentation,
GRS80GeodeticRepresentation,
)
from .spherical import (
SphericalRepresentation,
UnitSphericalRepresentation,
Expand Down Expand Up @@ -49,4 +55,8 @@
"RadialDifferential",
"CylindricalDifferential",
"PhysicsSphericalDifferential",
"BaseGeodeticRepresentation",
"WGS84GeodeticRepresentation",
"WGS72GeodeticRepresentation",
"GRS80GeodeticRepresentation",
]
121 changes: 121 additions & 0 deletions astropy/coordinates/representation/geodetic.py
@@ -0,0 +1,121 @@
# Licensed under a 3-clause BSD style license - see LICENSE.rst

import erfa

from astropy import units as u
from astropy.coordinates.angles import Latitude, Longitude
from astropy.utils.decorators import format_doc

from .base import BaseRepresentation
from .cartesian import CartesianRepresentation


ELLIPSOIDS = {}
"""Available ellipsoids (defined in erfam.h, with numbers exposed in erfa)."""
# Note: they get filled by the creation of the geodetic classes.


geodetic_base_doc = """{__doc__}
Parameters
----------
lon, lat : angle-like
The longitude and latitude of the point(s), in angular units. The
latitude should be between -90 and 90 degrees, and the longitude will
be wrapped to an angle between 0 and 360 degrees. These can also be
instances of `~astropy.coordinates.Angle` and either
`~astropy.coordinates.Longitude` not `~astropy.coordinates.Latitude`,
depending on the parameter.
height : `~astropy.units.Quantity` ['length']
The height to the point(s).
copy : bool, optional
If `True` (default), arrays will be copied. If `False`, arrays will
be references, though possibly broadcast to ensure matching shapes.
"""


@format_doc(geodetic_base_doc)
class BaseGeodeticRepresentation(BaseRepresentation):
"""
Base class for geodetic representations.
Subclasses need to set attributes ``_equatorial_radius`` and ``_flattening``
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`).
"""

attr_classes = {"lon": Longitude, "lat": Latitude, "height": u.Quantity}

def __init_subclass__(cls, **kwargs):
if "_ellipsoid" in cls.__dict__:
equatorial_radius, flattening = erfa.eform(getattr(erfa, cls._ellipsoid))
cls._equatorial_radius = equatorial_radius * u.m
cls._flattening = flattening * u.dimensionless_unscaled
ELLIPSOIDS[cls._ellipsoid] = cls
elif (
"_equatorial_radius" not in cls.__dict__
or "_flattening" not in cls.__dict__
):
raise AttributeError(
f"{cls.__name__} requires '_ellipsoid' or '_equatorial_radius' and '_flattening'."
)
super().__init_subclass__(**kwargs)

def __init__(self, lon, lat=None, height=None, copy=True):
if height is None and not isinstance(lon, self.__class__):
height = 0 << u.m

super().__init__(lon, lat, height, copy=copy)
if not self.height.unit.is_equivalent(u.m):
raise u.UnitTypeError(
f"{self.__class__.__name__} requires height with units of length."
)

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,
)
return CartesianRepresentation(xyz, xyz_axis=-1, copy=False)

@classmethod
def from_cartesian(cls, cart):
"""
Converts 3D rectangular cartesian coordinates (assumed geocentric) to
geodetic coordinates.
"""
lon, lat, height = erfa.gc2gde(
cls._equatorial_radius, cls._flattening, cart.get_xyz(xyz_axis=-1)
)
return cls(lon, lat, height, copy=False)


@format_doc(geodetic_base_doc)
class WGS84GeodeticRepresentation(BaseGeodeticRepresentation):
"""Representation of points in WGS84 3D geodetic coordinates."""

_ellipsoid = "WGS84"


@format_doc(geodetic_base_doc)
class WGS72GeodeticRepresentation(BaseGeodeticRepresentation):
"""Representation of points in WGS72 3D geodetic coordinates."""

_ellipsoid = "WGS72"


@format_doc(geodetic_base_doc)
class GRS80GeodeticRepresentation(BaseGeodeticRepresentation):
"""Representation of points in GRS80 3D geodetic coordinates."""

_ellipsoid = "GRS80"
3 changes: 2 additions & 1 deletion astropy/coordinates/tests/test_earth.py
Expand Up @@ -10,8 +10,9 @@
from astropy import constants
from astropy import units as u
from astropy.coordinates.angles import Latitude, Longitude
from astropy.coordinates.earth import ELLIPSOIDS, EarthLocation
from astropy.coordinates.earth import EarthLocation
from astropy.coordinates.name_resolve import NameResolveError
from astropy.coordinates.representation.geodetic import ELLIPSOIDS
from astropy.time import Time
from astropy.units import allclose as quantity_allclose
from astropy.units.tests.test_quantity_erfa_ufuncs import vvd
Expand Down

0 comments on commit f5f428e

Please sign in to comment.