Retrieve additional data from Horizons API

README.md

@@ -4,3 +4,8 @@ Usage:
 time matlab -nojvm -batch "process_mhs('/work/um0878/data/mhs/metopb_mhs_2024/07', 1000);"
 ```
 
+Add additional data from Horizons to CSV file:
+
+```bash
+python add-horizons.py CSVFILE
+```

add-horizons.py

+import csv
+import datetime
+import sys
+
+from horizons import retrieve_horizon_moon
+
+
+def add_diameter_to_file(matchfile):
+    with open(matchfile, "r") as f:
+        reader = csv.reader(f)
+        matches = list(reader)
+        header = matches[0]
+        matches = matches[1:]
+        horizon_values = ["diameter", "sub-lon", "sub-lat", "phaseangle", "sundistance"]
+        headerexist = {}
+        for value in horizon_values:
+            try:
+                header.index(value)
+                headerexist[value] = True
+            except ValueError:
+                header.append(value)
+                headerexist[value] = False
+
+        for match in matches:
+            year = int(match[header.index("year")])
+            doy = int(match[header.index("day")])
+            hour = int(match[header.index("hour")])
+            minute = int(match[header.index("minute")])
+            hour = int(match[header.index("hour")])
+            lon = float(match[header.index("lon")])
+            lat = float(match[header.index("lat")])
+            height = float(match[header.index("height")])
+            date = datetime.datetime(year, 1, 1, hour, minute)
+            date = (date + datetime.timedelta(days=doy - 1)).strftime("%Y-%m-%d %H:%M")
+            horizon = retrieve_horizon_moon(date, f"{lon},{lat},{height}")
+
+            for value in horizon_values:
+                if headerexist[value]:
+                    match[header.index(value)] = horizon[value]
+                else:
+                    match.append(horizon[value])
+
+
+    with open(matchfile, "w") as f:
+        writer = csv.writer(f)
+        writer.writerow(header)
+        writer.writerows(matches)
+
+
+def main():
+    if len(sys.argv) < 2:
+        print("Usage: python add-diameter.py CSVFILE [CSVFILE ...]")
+        sys.exit(1)
+
+    for matchfile in sys.argv[1:]:
+        add_diameter_to_file(matchfile)
+
+
+if __name__ == "__main__":
+    main()

horizons.py

+import csv
+import datetime
+import sys
+from scipy.constants import speed_of_light
+
+import requests
+
+def retrieve_horizon(url):
+    print(f"Retrieving: {url}")
+    response = requests.get(url)
+    print(f"Response code: {response.status_code}")
+    s = response.content.decode().split('\n')
+
+    i = 0
+    while s[i] != '$$SOE':
+        i += 1
+
+    result = []
+    i += 1
+    while s[i] != '$$EOE':
+        result.append(s[i])
+        i += 1
+    
+    return result
+
+def retrieve_horizon_moon(datestr, site_coord):
+    date = datetime.datetime.strptime(datestr, '%Y-%m-%d %H:%M')
+    starttime = date.strftime('%Y-%m-%d %H:%M')
+    stoptime = (date + datetime.timedelta(minutes=1)).strftime('%Y-%m-%d %H:%M')
+    url = f"https://ssd.jpl.nasa.gov/api/horizons.api?format=text&COMMAND='301'&OBJ_DATA='YES'&MAKE_EPHEM='YES'&EPHEM_TYPE='OBSERVER'&CENTER='coord'&SITE_COORD='{site_coord}'&START_TIME='{starttime}'&STOP_TIME='{stoptime}'&STEP_SIZE='1%20m'&QUANTITIES='13,14,43'&RANGE_UNITS='KM'&SUPPRESS_RANGE_RATE='YES'&CSV_FORMAT='YES'"
+
+    ret = retrieve_horizon(url)
+    data = ret[0].split(",")
+    result = {}
+    result["diameter"] = data[3].strip()
+    result["sub-lon"] = data[4].strip()
+    result["sub-lat"] = data[5].strip()
+    result["phaseangle"] = data[6].strip()
+
+    url = f"https://ssd.jpl.nasa.gov/api/horizons.api?format=text&COMMAND='301'&OBJ_DATA='YES'&MAKE_EPHEM='YES'&EPHEM_TYPE='OBSERVER'&CENTER='500@10'&START_TIME='{starttime}'&STOP_TIME='{stoptime}'&STEP_SIZE='1%20m'&QUANTITIES='21'&RANGE_UNITS='KM'&SUPPRESS_RANGE_RATE='YES'&CSV_FORMAT='YES'"
+
+    ret = retrieve_horizon(url)
+    data = ret[0].split(",")
+    result["sundistance"] = float(data[3].strip()) * 60 * speed_of_light
+
+    return result