Merge pull request #2746 from mkelley/jplhorizons-vectors-topo-queries-202305

Allow astropy Quantity in query locations

Author: bsipocz

CHANGES.rst

@@ -156,6 +156,8 @@ jplhorizons
 
 - Assign units to ``"hour_angle"``, ``"solartime"``, and ``"siderealtime"`` columns. [#2794]
 
+- Allow using units in locations specified as coordinates. [#2746]
+
 jplsbdb
 ^^^^^^^
 

astroquery/jplhorizons/core.py

@@ -13,6 +13,7 @@
 from astropy.table import Table, Column
 from astropy.io import ascii
 from astropy.time import Time
+from astropy import units as u
 from astropy.utils.exceptions import AstropyDeprecationWarning
 from astropy.utils.decorators import deprecated_renamed_argument, deprecated_attribute
 
@@ -53,14 +54,14 @@ def __init__(self, id=None, *, location=None, epochs=None,
 
         id : str or dict, required
             Name, number, or designation of target object. Uses the same codes
-            as JPL Horizons. Arbitrary topocentric coordinates can be added
-            in a dict. The dict has to be of the form
-            {``'lon'``: longitude in deg (East positive, West
-            negative), ``'lat'``: latitude in deg (North positive, South
-            negative), ``'elevation'``: elevation in km above the reference
-            ellipsoid, [``'body'``: Horizons body ID of the central body;
-            optional; if this value is not provided it is assumed that this
-            location is on Earth]}.
+            as JPL Horizons. Arbitrary topocentric coordinates can be added in a
+            dict. The dict has to be of the form {``'lon'``: longitude in deg
+            (East positive, West negative), ``'lat'``: latitude in deg (North
+            positive, South negative), ``'elevation'``: elevation in km above
+            the reference ellipsoid, [``'body'``: Horizons body ID of the
+            central body; optional; if this value is not provided it is assumed
+            that this location is on Earth]}.  Float values are assumed to have
+            units of degrees and kilometers.
 
         location : str or dict, optional
             Observer's location for ephemerides queries or center body name for
@@ -69,12 +70,13 @@ def __init__(self, id=None, *, location=None, epochs=None,
             ephemerides queries and the Sun's center for elements and vectors
             queries. Arbitrary topocentric coordinates for ephemerides queries
             can be provided in the format of a dictionary. The dictionary has to
-            be of the form {``'lon'``: longitude in deg (East positive, West
-            negative), ``'lat'``: latitude in deg (North positive, South
-            negative), ``'elevation'``: elevation in km above the reference
-            ellipsoid, [``'body'``: Horizons body ID of the central body;
-            optional; if this value is not provided it is assumed that this
-            location is on Earth]}.
+            be of the form {``'lon'``: longitude (East positive, West negative),
+            ``'lat'``: latitude (North positive, South negative),
+            ``'elevation'``: elevation above the reference ellipsoid,
+            [``'body'``: Horizons body ID of the central body; optional; if this
+            value is not provided it is assumed that this location is on
+            Earth]}.  Float values are assumed to have units of degrees and
+            kilometers.
 
         epochs : scalar, list-like, or dictionary, optional
             Either a list of epochs in JD or MJD format or a dictionary defining
@@ -117,16 +119,10 @@ def __init__(self, id=None, *, location=None, epochs=None,
         """
 
         super().__init__()
-        # check & format coordinate dictionaries for id and location; simply
-        # treat other values as given
-        if isinstance(id, Mapping):
-            self.id = self._prep_loc_dict(dict(id), "id")
-        else:
-            self.id = id
-        if isinstance(location, Mapping):
-            self.location = self._prep_loc_dict(dict(location), "location")
-        else:
-            self.location = location
+
+        self.id = id
+        self.location = location
+
         # check for epochs to be dict or list-like; else: make it a list
         if epochs is not None:
             if isinstance(epochs, (list, tuple, ndarray)):
@@ -185,6 +181,32 @@ def __str__(self):
                     str(self.epochs),
                     str(self.id_type))
 
+    @property
+    def id(self):
+        return self._id
+
+    @id.setter
+    def id(self, _id):
+        # check & format coordinate dictionaries for id; simply treat other
+        # values as given
+        if isinstance(_id, Mapping):
+            self._id = self._prep_loc_dict(dict(_id), "id")
+        else:
+            self._id = _id
+
+    @property
+    def location(self):
+        return self._location
+
+    @location.setter
+    def location(self, _location):
+        # check & format coordinate dictionaries for location; simply treat
+        # other values as given
+        if isinstance(_location, Mapping):
+            self._location = self._prep_loc_dict(dict(_location), "location")
+        else:
+            self._location = _location
+
     # ---------------------------------- query functions
 
     @deprecated_renamed_argument("get_raw_response", None, since="0.4.7",
@@ -593,10 +615,7 @@ def ephemerides_async(self, *, airmass_lessthan=99,
         if self.id is None:
             raise ValueError("'id' parameter not set. Query aborted.")
         elif isinstance(self.id, dict):
-            commandline = (
-                f"g:{self.id['lon']},{self.id['lat']},"
-                f"{self.id['elevation']}@{self.id['body']}"
-            )
+            commandline = self._format_id_coords(self.id)
         else:
             commandline = str(self.id)
         if self.location is None:
@@ -831,16 +850,20 @@ def elements_async(self, *, get_query_payload=False,
 
         URL = conf.horizons_server
 
-        # check for required information
+        # check for required information and assemble commandline stub
         if self.id is None:
             raise ValueError("'id' parameter not set. Query aborted.")
+        elif isinstance(self.id, dict):
+            commandline = self._format_id_coords(self.id)
+        else:
+            commandline = str(self.id)
+
         if self.location is None:
             self.location = '500@10'
         if self.epochs is None:
             self.epochs = Time.now().jd
 
-        # assemble commandline based on self.id_type
-        commandline = str(self.id)
+        # expand commandline based on self.id_type
         if self.id_type in ['designation', 'name',
                             'asteroid_name', 'comet_name']:
             commandline = ({'designation': 'DES=',
@@ -859,7 +882,7 @@ def elements_async(self, *, get_query_payload=False,
                 commandline += ' NOFRAG;'
 
         if isinstance(self.location, dict):
-            raise ValueError(('cannot use topographic position in orbital'
+            raise ValueError(('cannot use topographic position in orbital '
                               'elements query'))
 
         # configure request_payload for ephemerides query
@@ -1080,7 +1103,7 @@ def vectors_async(self, *, get_query_payload=False,
         if self.id is None:
             raise ValueError("'id' parameter not set. Query aborted.")
         elif isinstance(self.id, dict):
-            commandline = "g:{lon},{lat},{elevation}@{body}".format(**self.id)
+            commandline = self._format_id_coords(self.id)
         else:
             commandline = str(self.id)
         if self.location is None:
@@ -1183,20 +1206,38 @@ def _prep_loc_dict(loc_dict, attr_name):
             )
         if 'body' not in loc_dict:
             loc_dict['body'] = 399
+        # assumed units are degrees and km
+        loc_dict["lat"] = u.Quantity(loc_dict["lat"], u.deg)
+        loc_dict["lon"] = u.Quantity(loc_dict["lon"], u.deg)
+        loc_dict["elevation"] = u.Quantity(loc_dict["elevation"], u.km)
         return loc_dict
 
     @staticmethod
     def _location_to_params(loc_dict):
-        """translate a 'location' dict to a dict of request parameters"""
-        loc_dict = {
+        """translate a 'location' dict to request parameters"""
+
+        location = {
             "CENTER": f"coord@{loc_dict['body']}",
             "COORD_TYPE": "GEODETIC",
-            "SITE_COORD": ",".join(
-                str(float(loc_dict[k])) for k in ['lon', 'lat', 'elevation']
-            )
+            "SITE_COORD": "'{}'".format(str(HorizonsClass._format_site_coords(loc_dict)))
         }
-        loc_dict["SITE_COORD"] = f"'{loc_dict['SITE_COORD']}'"
-        return loc_dict
+        return location
+
+    @staticmethod
+    def _format_coords(coords):
+        """Dictionary to Horizons API formatted lon/lat/elevation coordinate triplet."""
+        return (f"{coords['lon'].to_value('deg')},{coords['lat'].to_value('deg')},"
+                f"{coords['elevation'].to_value('km')}")
+
+    @staticmethod
+    def _format_site_coords(coords):
+        """`location` dictionary to SITE_COORDS parameter."""
+        return HorizonsClass._format_coords(coords)
+
+    @staticmethod
+    def _format_id_coords(coords):
+        """`id` dictionary to COMMAND parameter's coordinate format."""
+        return f"g:{HorizonsClass._format_coords(coords)}@{coords['body']}"
 
     def _parse_result(self, response, verbose=None):
         """

astroquery/jplhorizons/tests/test_jplhorizons.py

@@ -7,6 +7,7 @@
 from numpy.ma import is_masked
 from astropy.tests.helper import assert_quantity_allclose
 from astropy.utils.exceptions import AstropyDeprecationWarning
+from astropy import units as u
 
 from astroquery.utils.mocks import MockResponse
 from ...query import AstroQuery
@@ -291,3 +292,61 @@ def test_id_type_deprecation():
 
     with pytest.warns(AstropyDeprecationWarning):
         jplhorizons.Horizons(id='Ceres', id_type='majorbody')
+
+
+def test_id_geodetic_coords():
+    """Test target based on geodetic coordinates.
+
+    From the Horizons manual:
+
+    For example, while 301 specifies the target to be the center of the Moon,
+    and Apollo-11 @ 301 specifies the Apollo-11 landing site as target, the
+    following:
+
+        g: 348.8, -43.3, 0 @ 301
+
+    specifies an ephemeris for the crater Tycho on the Moon (body 301), at
+    geodetic (planetodetic) coordinates 348.8 degrees east longitude, -43.3
+    degrees latitude (south), and zero km altitude with respect to the Moon’s
+    mean-Earth reference frame and ellipsoid surface.
+
+    """
+
+    target = {
+        "lon": 348.8 * u.deg,
+        "lat": -43.3 * u.deg,
+        "elevation": 0 * u.m,
+        "body": 301
+    }
+
+    q = jplhorizons.Horizons(id=target)
+    for payload in (q.ephemerides(get_query_payload=True),
+                    q.vectors(get_query_payload=True),
+                    q.elements(get_query_payload=True)):
+        assert payload["COMMAND"] == '"g:348.8,-43.3,0.0@301"'
+
+
+def test_location_topocentric_coords():
+    """Test location from topocentric coordinates.
+
+    Similar to `test_id_geodetic_coords`.
+
+    """
+
+    location = {
+        "lon": 348.8 * u.deg,
+        "lat": -43.3 * u.deg,
+        "elevation": 0 * u.m,
+        "body": 301
+    }
+
+    q = jplhorizons.Horizons(id=399, location=location)
+    for payload in (q.ephemerides(get_query_payload=True),
+                    q.vectors(get_query_payload=True)):
+        assert payload["CENTER"] == 'coord@301'
+        assert payload["COORD_TYPE"] == "GEODETIC"
+        assert payload["SITE_COORD"] == "'348.8,-43.3,0.0'"
+
+    # not allowed for elements
+    with pytest.raises(ValueError):
+        q.elements(get_query_payload=True)

astroquery/jplhorizons/tests/test_jplhorizons_remote.py

@@ -361,6 +361,29 @@ def test_vectors_query_raw(self):
 
         assert len(res) >= 6412
 
+    @pytest.mark.parametrize(
+        "location",
+        (
+            {"lon": 244, "lat": 33, "elevation": 1.7},
+            {"lon": (244 * u.deg).to(u.rad), "lat": (33 * u.deg).to(u.rad), "elevation": 1700 * u.m},
+        )
+    )
+    def test_vectors_query_topocentric_coordinates(self, location):
+        "Test vectors query specifying observer's longitude, latitude, and elevation"
+        q = jplhorizons.Horizons(id='Ceres',
+                                 location=location,
+                                 id_type='smallbody',
+                                 epochs=2451544.5)
+        res = q.vectors_async()
+        i = res.text.find("Center geodetic :")
+        j = res.text.find("\n", i)
+        parts = res.text[i:j].split()
+        assert parts[3:6] == ['244.0,', '33.0,', '1.7']
+
+        start = res.text.find("$$SOE")
+        end = res.text.find("$$EOE")
+        assert res.text[start:end].find("2000-Jan-01") > 0
+
     def test_unknownobject(self):
         with pytest.raises(ValueError):
             jplhorizons.Horizons(id='spamspamspameggsspam', location='500',

docs/jplhorizons/jplhorizons.rst

@@ -78,7 +78,7 @@ must still be entered in east-longitude, which means they must be negative; Hori
 will raise an error if given any positive longitude value for such bodies. Instead enter
 the west-longitude - 360. For instance, a site on Mars (id code 499) at 30 degrees
 longitude, 30 degrees latitude, 0 km elevation should be specified as
-``{'body': 499, 'elevation': 0, 'lon': -330, 'lat': 30}``.
+``{'body': 499, 'elevation': 0 * u.km, 'lon': -330 * u.deg, 'lat': 30 * u.deg}``.
 2. This does not apply to the Earth, Moon, and Sun. Although they are prograde,
 Horizons interprets east-longitude as positive and west-longitude as negative for these
 bodies.
@@ -96,14 +96,15 @@ as the observer's location, and Ceres as the target:
 
 .. doctest-remote-data::
 
-    >>> statue_of_liberty = {'lon': -74.0466891,
-    ...                      'lat': 40.6892534,
-    ...                      'elevation': 0.093}
+    >>> import astropy.units as u
+    >>> statue_of_liberty = {'lon': -74.0466891 * u.deg,
+    ...                      'lat': 40.6892534 * u.deg,
+    ...                      'elevation': 0.093 * u.km}
     >>> obj = Horizons(id='Ceres',
     ...                location=statue_of_liberty,
     ...                epochs=2458133.33546)
     >>> print(obj)
-    JPLHorizons instance "Ceres"; location={'lon': -74.0466891, 'lat': 40.6892534, 'elevation': 0.093, 'body': 399}, epochs=[2458133.33546], id_type=None
+    JPLHorizons instance "Ceres"; location={'lon': <Quantity -74.0466891 deg>, 'lat': <Quantity 40.6892534 deg>, 'elevation': <Quantity 0.093 km>, 'body': 399}, epochs=[2458133.33546], id_type=None
 
 2. Specifying topocentric coordinates for both location and observer is often
 useful when performing geometric calculations for artificial satellites without
@@ -116,8 +117,8 @@ at a particular point in time to the center of the crater Double:
 
 .. doctest-remote-data::
 
-    >>> ce_2 = {'lon': 23.522, 'lat': 0.637, 'elevation': 181.2, 'body': 301}
-    >>> double = {'lon': 23.47, 'lat': 0.67, 'elevation': 0, 'body': 301}
+    >>> ce_2 = {'lon': 23.522 * u.deg, 'lat': 0.637 * u.deg, 'elevation': 181.2 * u.km, 'body': 301}
+    >>> double = {'lon': 23.47 * u.deg, 'lat': 0.67 * u.deg, 'elevation': 0 * u.km, 'body': 301}
     >>> obj = Horizons(id=double, location=ce_2, epochs=2454483.84247)
     >>> vecs = obj.vectors()
     >>> distance_km = (vecs['x'] ** 2 + vecs['y'] ** 2 + vecs['z'] ** 2) ** 0.5 * 1.496e8