diff --git a/CHANGES.rst b/CHANGES.rst index e5351c158b..8263334903 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -100,6 +100,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 ^^^^^^^ diff --git a/astroquery/jplhorizons/__init__.py b/astroquery/jplhorizons/__init__.py index 83cb056865..5ab76be890 100644 --- a/astroquery/jplhorizons/__init__.py +++ b/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'), diff --git a/astroquery/jplhorizons/core.py b/astroquery/jplhorizons/core.py index 54ad23eab6..7861425e0d 100644 --- a/astroquery/jplhorizons/core.py +++ b/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('=') diff --git a/astroquery/jplhorizons/tests/test_jplhorizons.py b/astroquery/jplhorizons/tests/test_jplhorizons.py index b8fbb4fe66..129f13f068 100644 --- a/astroquery/jplhorizons/tests/test_jplhorizons.py +++ b/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): diff --git a/astroquery/jplhorizons/tests/test_jplhorizons_remote.py b/astroquery/jplhorizons/tests/test_jplhorizons_remote.py index 4f69be6821..07c8d3784d 100644 --- a/astroquery/jplhorizons/tests/test_jplhorizons_remote.py +++ b/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) diff --git a/docs/jplhorizons/jplhorizons.rst b/docs/jplhorizons/jplhorizons.rst index 4c8fc78bed..9a6feba1b6 100644 --- a/docs/jplhorizons/jplhorizons.rst +++ b/docs/jplhorizons/jplhorizons.rst @@ -40,26 +40,59 @@ In order to query information for a specific Solar System body, a >>> print(obj) JPLHorizons instance "Ceres"; location=568, epochs=[2458133.33546], id_type=None -``id`` refers to the target identifier and is mandatory; the exact -string will be used in the query to the Horizons system. - -``location`` means either the observer's location (e.g., Horizons ephemerides -query) or the body relative to which orbital elements are provided (e.g., -Horizons orbital elements or vectors query); the same codes as `used by Horizons -`_ are used here, which -includes `MPC Observatory codes`_. The default is ``location=None``, which uses -a geocentric location for ephemerides queries and the Sun as central body for -orbital elements and state vector queries. User-defined topocentric locations -for ephemerides queries can be provided, too, in the form of a dictionary. The -dictionary has to be formatted as follows: {``'lon'``: longitude in degrees -(East positive, West negative), ``'lat'``: latitude in degrees (North positive, -South negative), ``'elevation'``: elevation in km above the reference -ellipsoid}. In addition, ``'body'`` can be set to the Horizons body ID of the -central body if different from Earth; by default, it is assumed that this -location is on Earth if it has not been specifically set. The following example -uses the coordinates of the `Statue of Liberty +``id`` refers to the target object and is mandatory. ````str`` and ``int`` values +are valid for all query types. ``Mapping`` (e.g. ``dict``) values are valid for +observer (``ephemerides``) and vectors queries only. ``str`` or ``int`` values will be +passed directly to Horizons. See the description of the ``id_type`` argument below +for how Horizons interprets these values. See the paragraph below the description +of the ``location`` argument for valid ``dict`` formatting. + +``location`` refers to the coordinate center for the ephemeris, which has +slightly different physical interpretations depending on the query type: +observer (``ephemerides``) queries: observer location +vectors queries: coordinate origin for vectors +elements queries: relative body for orbital elements + +``str`` and ``int`` values are valid for all query types. ``Mapping`` +(e.g. ``dict``) values are valid for observer (``ephemerides``) and vectors queries only. ``str`` or ``int`` +arguments will be passed directly to Horizons. See `this section of the Horizons +manual`_ for how Horizons interprets coordinate center codes; also note that, +unlike ``id``, these include (most) `MPC Observatory codes`_. See below for valid +``dict`` formatting. The default is ``location=None``, which uses Earth body center +for observer queries, and Sun body center for orbital elements and vectors queries. + +``dict``-like arguments to ``id`` or ``location`` define a topocentric location +relative to a major body. Note that this is not possible for elements queries, +and will only work for bodies with defined rotation models (Horizons does not +have rotation models for many small or recently-discovered natural satellites). +The dictionary has to be formatted as follows: {``'lon'``: longitude in degrees, +``'lat'``: latitude in degrees (North positive, South negative), ``'elevation'``: +elevation in km above the reference ellipsoid, ``body``: Horizons ID of center body, +optional; default Earth}. + +Horizons always interprets longitude values as eastward. However, there are two +major gotchas in this: +1. For most prograde rotators, which is to say most major bodies, Horizons +interprets west-longitude as positive and east-longitude as negative. However, values +must still be entered in east-longitude, which means they must be negative; Horizons +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}``. +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. + +Here is a complete list of retrograde major bodies in Horizons: Venus (299), Arial (701), +Umbriel (702), Titania (703), Oberon (704), Miranda (705), Cordelia (706), Ophelia (707), +Bianca (708), Cressida (709), Desdemona (710), Juliet (711), Portia (712), Rosalind (713), +Belinda (714), Puck (715), Uranus (799), Triton (801). All other major bodies are prograde. + +Two examples of usage for specified topocentric coordinates follow. + +1. This observer (``ephemerides``) query uses the coordinates of the `Statue of Liberty `_ -as the observer's location: +as the observer's location, and Ceres as the target: .. doctest-remote-data:: @@ -72,6 +105,24 @@ as the observer's location: >>> print(obj) JPLHorizons instance "Ceres"; location={'lon': -74.0466891, 'lat': 40.6892534, 'elevation': 0.093}, 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 +completely-specified ephemeris data. For instance, published reduced data for the +lunar satellite Chang'e-2 include orbital height and lat/lon. Although there is no published +ephemeris for Chang'e-2, Horizons (combined with the fact that Chang'e-2 looked nadir), +can be used to compute vectors from Chang'e-2 to specific points on the lunar surface. +Here is an example of using ``jplhorizons`` to find the distance from Chang'e-2 +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} + >>> 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 + >>> print(f"{distance_km.value.data[0]:.3f}") + 181.213 ``epochs`` is either a scalar or list of Julian dates (floats or strings) in the case of discrete epochs, or, in the case of a range of epochs, a dictionary that @@ -609,3 +660,4 @@ Reference/API .. _astropy units: http://docs.astropy.org/en/stable/units/index.html .. _Definition of Observer Table Quantities: https://ssd.jpl.nasa.gov/horizons/manual.html#observer-table .. _Horizons documentation: https://ssd.jpl.nasa.gov/horizons/manual.html#observer-table +.. _this section of the Horizons manual: \ No newline at end of file