Merge pull request #2625 from MillionConcepts/jplh_planetodetic_coord

improve target coordinate specification in jplhorizons

Author: bsipocz

CHANGES.rst

@@ -105,6 +105,9 @@ jplhorizons
 
 - Optional keyword arguments are now keyword only. [#1802]
 
+- Topocentric coordinates can now be specified for both center and target in observer
+  and vector queries. [#2625]
+
 jplsbdb
 ^^^^^^^
 

astroquery/jplhorizons/__init__.py

@@ -50,7 +50,10 @@ class Conf(_config.ConfigNamespace):
                    'k2': ('k2', '---'),
                    'phasecoeff': ('phasecoeff', 'mag/deg'),
                    'solar_presence': ('solar_presence', '---'),
-                   'flags': ('flags', '---'),
+                   'lunar_presence': ('lunar_presence', '---'),
+                   'interfering_body': ('interfering_body', '---'),
+                   'illumination_flag': ('illumination_flag', '---'),
+                   'nearside_flag': ('nearside_flag', '---'),
                    'R.A._(ICRF)': ('RA', 'deg'),
                    'DEC_(ICRF)': ('DEC', 'deg'),
                    'R.A.___(ICRF)': ('RA', 'deg'),

astroquery/jplhorizons/core.py

@@ -1,15 +1,15 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 
-
 # 1. standard library imports
-from numpy import nan
-from numpy import isnan
-from numpy import ndarray
 from collections import OrderedDict
+from typing import Mapping
 import warnings
 
 # 2. third party imports
 from requests.exceptions import HTTPError
+from numpy import nan
+from numpy import isnan
+from numpy import ndarray
 from astropy.table import Table, Column
 from astropy.io import ascii
 from astropy.time import Time
@@ -51,8 +51,16 @@ def __init__(self, id=None, *, location=None, epochs=None,
         Parameters
         ----------
 
-        id : str, required
-            Name, number, or designation of the object to be queried.
+        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]}.
 
         location : str or dict, optional
             Observer's location for ephemerides queries or center body name for
@@ -108,9 +116,16 @@ def __init__(self, id=None, *, location=None, epochs=None,
         """
 
         super().__init__()
-        self.id = id
-        self.location = location
-
+        # 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
         # 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)):
@@ -535,16 +550,22 @@ def ephemerides_async(self, *, airmass_lessthan=99,
 
         URL = conf.horizons_server
 
-        # check for required information
+        # check for required information and assemble commanddline stub
         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']}"
+            )
+        else:
+            commandline = str(self.id)
         if self.location is None:
             self.location = '500@399'
         if self.epochs is None:
             self.epochs = Time.now().jd
+        # expand commandline based on self.id_type
 
-        # assemble commandline based on self.id_type
-        commandline = str(self.id)
         if self.id_type in ['designation', 'name',
                             'asteroid_name', 'comet_name']:
             commandline = ({'designation': 'DES=',
@@ -580,19 +601,7 @@ def ephemerides_async(self, *, airmass_lessthan=99,
             ('EXTRA_PREC', {True: 'YES', False: 'NO'}[extra_precision])])
 
         if isinstance(self.location, dict):
-            if ('lon' not in self.location or 'lat' not in self.location or
-                    'elevation' not in self.location):
-                raise ValueError(("'location' must contain lon, lat, "
-                                  "elevation"))
-
-            if 'body' not in self.location:
-                self.location['body'] = '399'
-            request_payload['CENTER'] = 'coord@{:s}'.format(
-                str(self.location['body']))
-            request_payload['COORD_TYPE'] = 'GEODETIC'
-            request_payload['SITE_COORD'] = "'{:f},{:f},{:f}'".format(
-                self.location['lon'], self.location['lat'],
-                self.location['elevation'])
+            request_payload = dict(**request_payload, **self._location_to_params(self.location))
         else:
             request_payload['CENTER'] = "'" + str(self.location) + "'"
 
@@ -1032,17 +1041,18 @@ def vectors_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 = "g:{lon},{lat},{elevation}@{body}".format(**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=',
@@ -1060,18 +1070,12 @@ def vectors_async(self, *, get_query_payload=False,
                 commandline += ' CAP{:s};'.format(closest_apparition)
             if no_fragments:
                 commandline += ' NOFRAG;'
-
-        if isinstance(self.location, dict):
-            raise ValueError(('cannot use topographic position in state'
-                              'vectors query'))
-
-        # configure request_payload for ephemerides query
+        # configure request_payload for vectors query
         request_payload = OrderedDict([
             ('format', 'text'),
             ('EPHEM_TYPE', 'VECTORS'),
             ('OUT_UNITS', 'AU-D'),
             ('COMMAND', '"' + commandline + '"'),
-            ('CENTER', ("'" + str(self.location) + "'")),
             ('CSV_FORMAT', ('"YES"')),
             ('REF_PLANE', {'ecliptic': 'ECLIPTIC',
                            'earth': 'FRAME',
@@ -1086,7 +1090,12 @@ def vectors_async(self, *, get_query_payload=False,
             ('VEC_DELTA_T', {True: 'YES', False: 'NO'}[delta_T]),
             ('OBJ_DATA', 'YES')]
         )
-
+        if isinstance(self.location, dict):
+            request_payload = dict(
+                **request_payload, **self._location_to_params(self.location)
+            )
+        else:
+            request_payload['CENTER'] = "'" + str(self.location) + "'"
         # parse self.epochs
         if isinstance(self.epochs, (list, tuple, ndarray)):
             request_payload['TLIST'] = "\n".join([str(epoch) for epoch in
@@ -1132,6 +1141,30 @@ def vectors_async(self, *, get_query_payload=False,
         return response
 
     # ---------------------------------- parser functions
+    @staticmethod
+    def _prep_loc_dict(loc_dict, attr_name):
+        """prepare coord specification dict for 'location' or 'id'"""
+        if {'lat', 'lon', 'elevation'} - loc_dict.keys():
+            raise ValueError(
+                f"dict values for '{attr_name}' must contain 'lat', 'lon', "
+                "'elevation' (and optionally 'body')"
+            )
+        if 'body' not in loc_dict:
+            loc_dict['body'] = 399
+        return loc_dict
+
+    @staticmethod
+    def _location_to_params(loc_dict):
+        """translate a 'location' dict to a dict of request parameters"""
+        loc_dict = {
+            "CENTER": f"coord@{loc_dict['body']}",
+            "COORD_TYPE": "GEODETIC",
+            "SITE_COORD": ",".join(
+                str(float(loc_dict[k])) for k in ['lon', 'lat', 'elevation']
+            )
+        }
+        loc_dict["SITE_COORD"] = f"'{loc_dict['SITE_COORD']}'"
+        return loc_dict
 
     def _parse_result(self, response, verbose=None):
         """
@@ -1181,14 +1214,18 @@ def _parse_result(self, response, verbose=None):
         H, G = nan, nan
         M1, M2, k1, k2, phcof = nan, nan, nan, nan, nan
         headerline = []
+        centername = ''
         for idx, line in enumerate(src):
             # read in ephemerides header line; replace some field names
             if (self.query_type == 'ephemerides' and
                     "Date__(UT)__HR:MN" in line):
                 headerline = str(line).split(',')
                 headerline[2] = 'solar_presence'
-                headerline[3] = 'flags'
+                headerline[3] = "lunar_presence" if "Earth" in centername else "interfering_body"
                 headerline[-1] = '_dump'
+                if isinstance(self.id, dict) or str(self.id).startswith('g:'):
+                    headerline[4] = 'nearside_flag'
+                    headerline[5] = 'illumination_flag'
             # read in elements header line
             elif (self.query_type == 'elements' and
                   "JDTDB," in line):
@@ -1208,6 +1245,9 @@ def _parse_result(self, response, verbose=None):
             # read in targetname
             if "Target body name" in line:
                 targetname = line[18:50].strip()
+            # read in center body name
+            if "Center body name" in line:
+                centername = line[18:50].strip()
             # read in H and G (if available)
             if "rotational period in hours)" in line:
                 HGline = src[idx + 2].split('=')

astroquery/jplhorizons/tests/test_jplhorizons.py

@@ -83,7 +83,7 @@ def test_ephemerides_query(patch_request):
     assert res['targetname'] == "1 Ceres (A801 AA)"
     assert res['datetime_str'] == "2000-Jan-01 00:00:00.000"
     assert res['solar_presence'] == ""
-    assert res['flags'] == ""
+    assert res['lunar_presence'] == ""
     assert res['elongFlag'] == '/L'
     assert res['airmass'] == 999
 
@@ -256,22 +256,21 @@ def test_vectors_query_payload():
     res = jplhorizons.Horizons(id='Ceres', location='500@10',
                                epochs=2451544.5).vectors(
                                    get_query_payload=True)
-
     assert res == OrderedDict([
-        ('format', 'text'),
-        ('EPHEM_TYPE', 'VECTORS'),
-        ('OUT_UNITS', 'AU-D'),
-        ('COMMAND', '"Ceres"'),
-        ('CENTER', "'500@10'"),
-        ('CSV_FORMAT', '"YES"'),
-        ('REF_PLANE', 'ECLIPTIC'),
-        ('REF_SYSTEM', 'ICRF'),
-        ('TP_TYPE', 'ABSOLUTE'),
-        ('VEC_LABELS', 'YES'),
-        ('VEC_CORR', '"NONE"'),
-        ('VEC_DELTA_T', 'NO'),
-        ('OBJ_DATA', 'YES'),
-        ('TLIST', '2451544.5')])
+                ('format', 'text'),
+                ('EPHEM_TYPE', 'VECTORS'),
+                ('OUT_UNITS', 'AU-D'),
+                ('COMMAND', '"Ceres"'),
+                ('CSV_FORMAT', '"YES"'),
+                ('REF_PLANE', 'ECLIPTIC'),
+                ('REF_SYSTEM', 'ICRF'),
+                ('TP_TYPE', 'ABSOLUTE'),
+                ('VEC_LABELS', 'YES'),
+                ('VEC_CORR', '"NONE"'),
+                ('VEC_DELTA_T', 'NO'),
+                ('OBJ_DATA', 'YES'),
+                ('CENTER', "'500@10'"),
+                ('TLIST', '2451544.5')])
 
 
 def test_no_H(patch_request):

astroquery/jplhorizons/tests/test_jplhorizons_remote.py

@@ -1,10 +1,12 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 
-
-import pytest
 from numpy.ma import is_masked
+import numpy as np
+import pytest
 
+from astropy.coordinates import spherical_to_cartesian
 from astropy.tests.helper import assert_quantity_allclose
+import astropy.units as u
 from astropy.utils.exceptions import AstropyDeprecationWarning
 
 from ... import jplhorizons
@@ -23,7 +25,7 @@ def test_ephemerides_query(self):
         assert res['targetname'] == "1 Ceres (A801 AA)"
         assert res['datetime_str'] == "2000-Jan-01 00:00:00.000"
         assert res['solar_presence'] == ""
-        assert res['flags'] == ""
+        assert res['lunar_presence'] == ""
         assert res['elongFlag'] == '/L'
         assert res['airmass'] == 999
 
@@ -75,7 +77,7 @@ def test_ephemerides_query_two(self):
         assert res['targetname'] == "1P/Halley"
         assert res['datetime_str'] == "2080-Jan-11 09:00"
         assert res['solar_presence'] == ""
-        assert res['flags'] == "m"
+        assert res['lunar_presence'] == "m"
         assert res['elongFlag'] == '/L'
 
         for value in ['H', 'G']:
@@ -98,7 +100,7 @@ def test_ephemerides_query_three(self):
         assert res['targetname'] == "73P/Schwassmann-Wachmann 3"
         assert res['datetime_str'] == "2080-Jan-01 00:00"
         assert res['solar_presence'] == "*"
-        assert res['flags'] == "m"
+        assert res['lunar_presence'] == "m"
         assert res['elongFlag'] == '/L'
 
         for value in ['H', 'G']:
@@ -123,7 +125,7 @@ def test_ephemerides_query_four(self):
         assert res['targetname'] == "167P/CINEOS"
         assert res['datetime_str'] == "2080-Jan-01 00:00"
         assert res['solar_presence'] == "*"
-        assert res['flags'] == "m"
+        assert res['lunar_presence'] == "m"
         assert res['elongFlag'] == '/T'
 
         for value in ['H', 'G', 'M1', 'k1']:
@@ -150,7 +152,7 @@ def test_ephemerides_query_five(self):
         assert res['targetname'] == "12P/Pons-Brooks"
         assert res['datetime_str'] == "2080-Jan-01 00:00"
         assert res['solar_presence'] == "*"
-        assert res['flags'] == "m"
+        assert res['lunar_presence'] == "m"
         assert res['elongFlag'] == '/L'
 
         for value in ['H', 'G', 'phasecoeff']:
@@ -405,3 +407,35 @@ def test_ephemerides_extraprecision(self):
         vec_highprec = obj.ephemerides(extra_precision=True)
 
         assert (vec_simple['RA'][0]-vec_highprec['RA'][0]) > 1e-7
+
+    def test_geodetic_queries(self):
+        """
+        black-box test for observer and vectors queries with geodetic
+        coordinates. checks spatial sensibility.
+        """
+        phobos = {'body': 401, 'lon': -30, 'lat': -20, 'elevation': 0}
+        deimos = {'body': 402, 'lon': -10, 'lat': -40, 'elevation': 0}
+        deimos_phobos = jplhorizons.Horizons(phobos, location=deimos, epochs=2.4e6)
+        phobos_deimos = jplhorizons.Horizons(deimos, location=phobos, epochs=2.4e6)
+        pd_eph, dp_eph = phobos_deimos.ephemerides(), deimos_phobos.ephemerides()
+        dp_xyz = spherical_to_cartesian(
+            dp_eph['delta'], dp_eph['DEC'], dp_eph['RA']
+        )
+        pd_xyz = spherical_to_cartesian(
+            pd_eph['delta'], pd_eph['DEC'], pd_eph['RA']
+        )
+        elementwise = [(dp_el + pd_el) for dp_el, pd_el in zip(dp_xyz, pd_xyz)]
+        eph_offset = (sum([off ** 2 for off in elementwise]) ** 0.5).to(u.km)
+        # horizons can do better than this, but we'd have to go to a little
+        # more trouble than is necessary for a software test...
+        assert np.isclose(eph_offset.value, 2.558895)
+        # ...and vectors queries are really what you're meant to use for
+        # this sort of thing.
+        pd_vec, dp_vec = phobos_deimos.vectors(), deimos_phobos.vectors()
+        vec_offset = np.sum(
+            (
+                pd_vec.as_array(names=('x', 'y', 'z')).view('f8')
+                + dp_vec.as_array(names=('x', 'y', 'z')).view('f8')
+            ) ** 2
+        )
+        assert np.isclose(vec_offset, 0)