Add refraction_correction option to nasa_horizons solar position algorithm (#133)

* Add adminition for nasa_horizons

* Correct admonition formatting

* Fix indentation

* Add refraction_correction option to nasa_horizons solar position algorithm

* (auto) Paper PDF Draft

* Fix doc warning

* Fix rendering

Co-authored-by: Kevin Anderson <[email protected]>

* Apply suggestions from code review

Co-authored-by: Ioannis Sifnaios <[email protected]>

* Add refraction algorithm note and reference

* (auto) Paper PDF Draft

---------

Co-authored-by: AdamRJensen <[email protected]>
Co-authored-by: Kevin Anderson <[email protected]>
Co-authored-by: Ioannis Sifnaios <[email protected]>

Author: 3 people

paper/paper.pdf

@@ -2,12 +2,11 @@
 import pandas as pd
 import io
 
-URL = "https://ssd.jpl.nasa.gov/api/horizons.api"
+URL = 'https://ssd.jpl.nasa.gov/api/horizons.api'
 
 
-def nasa_horizons(
-    latitude, longitude, start, end, elevation=0.0, *, time_step="1h", url=URL
-):
+def nasa_horizons(latitude, longitude, start, end, elevation=0., *,
+                  time_step='1h', refraction_correction=False, url=URL):
     """
     Retrieve solar positions from NASA's Horizons web service.
 
@@ -42,30 +41,39 @@ def nasa_horizons(
         Time step size of the requested time series. '1m' for minutes,
         '1h' for hours, '1d' for days, '1mo' for months, and '1y' for years.
         The default is '1h'.
+    refraction_correction : bool, optional
+        Whether to return apparent solar position angles corrected for
+        atmospheric refraction. The defalt is False.
     url : str, optional
         API endpoint. The default is
         ``'https://ssd.jpl.nasa.gov/api/horizons.api'``.
 
     Returns
     -------
     pandas.DataFrame
-        DataFrame with the following columns (all values in degrees):
+        DataFrame with the following columns in degrees (note that all columns
+        except azimuth are prefixed with `'apparent_'` if
+        ``refraction_correction=True``.)
 
-        - uncertainty_right_ascension
-        - uncertainty_declination
         - right_ascension
         - declination
-        - apparent_right_ascesion
-        - apparent_declination
-        - apparent_azimuth
-        - apparent_elevation
+        - azimuth
+        - elevation
+
+    Notes
+    -----
+    The refraction correction algorithm is the same as
+    :py:func:`solposx.refraction.spa`. See [3]_ for more
+    information.
 
     References
     ----------
     .. [1] `NASA Horizons Systems
        <https://ssd.jpl.nasa.gov/horizons/>`_
     .. [2] `NASA Horizons API
        <https://ssd-api.jpl.nasa.gov/doc/horizons.html>`_
+    .. [3] `NASA Horizons Manual
+       <https://ssd.jpl.nasa.gov/horizons/manual.html>`_
     """
     params = {
         "MAKE_EPHEM": "YES",  # generate ephemeris
@@ -74,62 +82,75 @@ def nasa_horizons(
         "CENTER": "coord@399",  # input coordinates for Earth (399)
         "COORD_TYPE": "GEODETIC",  # latitude, longitude, elevation in degrees
         "SITE_COORD": f"{longitude},{latitude},{elevation/1000}",
-        "START_TIME": pd.Timestamp(start).strftime("%Y-%m-%d %H:%M"),
-        "STOP_TIME": pd.Timestamp(end).strftime("%Y-%m-%d %H:%M"),
+        "START_TIME": pd.Timestamp(start).strftime('%Y-%m-%d %H:%M'),
+        "STOP_TIME": pd.Timestamp(end).strftime('%Y-%m-%d %H:%M'),
         "STEP_SIZE": time_step,
         "QUANTITIES": "2,4",  # corrected angles
         "REF_SYSTEM": "ICRF",
         "CAL_FORMAT": "CAL",  # output date format
         "CAL_TYPE": "MIXED",  # Gregorian or mixed Julian/Gregorian calendar
         "TIME_DIGITS": "FRACSEC",  # output time precision
         "ANG_FORMAT": "DEG",  # output angles in degrees
-        "APPARENT": "AIRLESS",  # no refraction
-        # "RANGE_UNITS": "AU",
-        # "SUPPRESS_RANGE_RATE": "NO",
         "SKIP_DAYLT": "NO",  # include daylight periods
         "SOLAR_ELONG": "0,180",
         "EXTRA_PREC": "NO",  # toggle additional digits on some angles (RA/DEC)
         "CSV_FORMAT": "NO",
-        "OBJ_DATA": "NO",  # whether to return summary data
+        "OBJ_DATA": "NO"  # whether to return summary data
     }
 
+    if refraction_correction:
+        params["APPARENT"] = "REFRACTED"  # account for refraction
+    else:
+        params["APPARENT"] = "AIRLESS"  # do not account for refraction
+
     # manual formatting of the url as all parameters except format shall be
-    # in enclosed in single quotes
+    # enclosed in single quotes
     url_formatted = (
-        url + "?" + "format=text&" + "&".join([f"{k}='{v}'" for k, v in params.items()])
+        url
+        + '?'
+        + "format=text&"
+        + '&'.join([f"{k}='{v}'" for k, v in params.items()])
     )
 
     res = requests.get(url_formatted)
 
     fbuf = io.StringIO(res.content.decode())
 
     lines = fbuf.readlines()
-    first_line = lines.index("$$SOE\n") + 1
-    last_line = lines.index("$$EOE\n", first_line)
+    first_line = lines.index('$$SOE\n') + 1
+    last_line = lines.index('$$EOE\n', first_line)
     header_line = lines[first_line - 3]
-    data_str = [header_line] + lines[first_line:last_line]
+    data_str = [header_line] + lines[first_line: last_line]
 
-    data = pd.read_fwf(
-        io.StringIO("\n".join(data_str)), index_col=[0], na_values=["n.a."]
-    )
-    data.index = pd.to_datetime(data.index, format="%Y-%b-%d %H:%M:%S.%f")
-    data.index = data.index.tz_localize("UTC")
+    data = pd.read_fwf(io.StringIO('\n'.join(data_str)),
+                       index_col=[0], na_values=['n.a.'])
+    data.index = pd.to_datetime(data.index, format='%Y-%b-%d %H:%M:%S.%f')
+    data.index = data.index.tz_localize('UTC')
 
     # split columns as several params have a shared header name for two params
-    column_name_split_map = {
-        "R.A._(a-appar)_DEC.": ["right_ascension", "declination"],
-        "Azi____(a-app)___Elev": ["azimuth", "elevation"],
-    }
+    if refraction_correction:
+        column_name_split_map = {
+            # "r" prefix stands for "refracted" aka. refracted/apparent angles
+            'R.A._(r-appar)__DEC': ['apparent_right_ascension', 'apparent_declination'],
+            'Azi____(r-app)___Elev': ['azimuth', 'apparent_elevation'],
+        }
+
+    else:
+        column_name_split_map = {
+            # "a" prefix stands for "airless" aka. unrefracted
+            'R.A._(a-appar)_DEC.': ['right_ascension', 'declination'],
+            'Azi____(a-app)___Elev': ['azimuth', 'elevation'],
+        }
 
     for old_name, new_names in column_name_split_map.items():
-        data[new_names] = data[old_name].str.split(r"\s+", expand=True).astype(float)
+        data[new_names] = \
+            data[old_name].str.split(r'\s+', expand=True).astype(float)
 
     data = data.drop(columns=list(column_name_split_map.keys()))
 
-    data.index.name = "time"
+    data.index.name = 'time'
     try:
-        del data["Unnamed: 1"]
+        del data['Unnamed: 1']
     except KeyError:
         pass
-
     return data

src/solposx/solarposition/nasa_horizons.py

@@ -559,6 +559,7 @@ def test_delta_t_series_input():
     pd.testing.assert_frame_equal(usno_array, usno_float)
 
 
+@pytest.fixture
 def expected_nasa():
     columns = ['right_ascension', 'declination', 'azimuth', 'elevation']
     values = [
@@ -595,17 +596,64 @@ def expected_nasa():
     return data
 
 
-def test_nasa_horizons():
-    expected = expected_nasa()
+@pytest.fixture
+def expected_nasa_refracted():
+    columns = ['apparent_right_ascension', 'apparent_declination', 'azimuth', 'apparent_elevation']
+    values = [
+        [280.73492, -22.42915, 18.384905, -61.460898],
+        [280.59597, -22.49103, 43.990328, -56.428642],
+        [280.54556, -22.54743, 63.159488, -48.660005],
+        [280.55006, -22.57750, 77.862491, -39.580841],
+        [280.58537, -22.58401, 90.093838, -30.008965],
+        [280.64079, -22.57277, 101.147333, -20.430330],
+        [280.71282, -22.54831, 111.849769, -11.199780],
+        [280.80156, -22.51427, 122.773903, -2.643156],
+        [281.15991, -22.87756, 134.336441, 4.426409],
+        [281.26791, -22.95121, 146.804482, 10.455202],
+        [281.33369, -22.96770, 160.234863, 14.733079],
+        [281.39128, -22.97067, 174.399466, 16.846454],
+        [281.44734, -22.96663, 188.799258, 16.581405],
+        [281.50589, -22.95497, 202.838039, 13.965712],
+        [281.57563, -22.92612, 216.067120, 9.260021],
+        [281.71485, -22.80480, 228.325009, 2.936098],
+        [282.04557, -22.49171, 239.722349, -4.450189],
+        [282.13033, -22.51721, 250.562948, -13.182278],
+        [282.19859, -22.53231, 261.295115, -22.514892],
+        [282.24984, -22.53327, 272.536349, -32.120854],
+        [282.27937, -22.51526, 285.194755, -41.626608],
+        [282.27410, -22.47247, 300.704397, -50.493032],
+        [282.20646, -22.40488, 321.186284, -57.787601],
+        [282.04689, -22.34306, 348.166418, -61.945628],
+        [281.84142, -22.34828, 18.119982, -61.406306],
+    ]
+    data = pd.DataFrame(data=values, columns=columns)
+    data.index = pd.date_range('2020-01-01', '2020-01-02', freq='1h', tz='UTC')
+    data.index.name = 'time'
+    data.index.freq = None
+    return data
 
+
+def test_nasa_horizons(expected_nasa):
     result = nasa_horizons(
         latitude=50,
         longitude=10,
         start='2020-01-01',
         end='2020-01-02',
     )
-    for c in expected.columns:
-        pd.testing.assert_series_equal(expected[c], result[c])
+    for c in expected_nasa.columns:
+        pd.testing.assert_series_equal(expected_nasa[c], result[c])
+
+
+def test_nasa_horizons_refracted(expected_nasa_refracted):
+    result = nasa_horizons(
+        latitude=50,
+        longitude=10,
+        start='2020-01-01',
+        end='2020-01-02',
+        refraction_correction=True,
+    )
+    for c in expected_nasa_refracted.columns:
+        pd.testing.assert_series_equal(expected_nasa_refracted[c], result[c])
 
 
 def test_nasa_horizons_frequency_daily():
@@ -634,6 +682,7 @@ def test_nasa_horizons_frequency_monthly():
         start=start,
         end=end,
         time_step='1mo',  # monthly
+        refraction_correction=True,
     )
     expected_index = pd.date_range(start, end, freq='MS', tz='UTC')
     expected_index.name = 'time'