Determining Horizontal Coordinate System data (azimuth, altitude) for geographic locations

Question

I would like to determine, specifically, azimuth and elevation data for any given location on the surface of the earth relative to another location on the surface of the earth (for my purposes, the latter will be London, UK).

Data and APIs are readily available for certain celestial objects, but I haven't had any luck finding data as required.

Can anyone enlighten me as to a method for calculating such data?

Many thanks.

Answer 1

This can be done using the Python package Skyfield fairly easily.

Here is a shell of a script; I used some dictionaries to hold data, you may want to do something else. I've stored a lot of goodies but only printed azimuth, altitude and distance in kilometers.

By default Skyfield's .altaz() method calculates atmospheric refraction for anything higher than -1 degrees (1 degree below the horizon geometrically). You can turn that off by replacing it with .altaz(pressure_mbar=0).

OUTPUT:

 place     altitude  azimuth     distance (km)
Caracas    -33.67   258.26         7075.28
Edinburgh  -2.39    338.98          533.65
Canberra  -76.25     64.93        12377.25

Here's the Python script:

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import Topos, Loader

load        = Loader('~/Documents/fishing/SkyData')  # avoids multiple copies of large files
data        = load('de421.bsp')
ts          = load.timescale()

Earth       = data['earth']

London      = Earth + Topos(latitude_degrees  = 51.51,
                            longitude_degrees = -0.1275,
                            elevation_m = 11.)

places      = {'Caracas':(10.48, -66.90, 900), 'Edinburgh':(55.95, -3.19, 47.),
               'Canberra':(-35.29, 149.13, 578.)}

now         = ts.now()

answer_dict = dict()

for name, (lat, lon, elev) in places.items():
    place     = Earth + Topos(latitude_degrees  = lat,
                              longitude_degrees = lon,
                              elevation_m       = elev)
    vector_km            = place.at(now).position.km - London.at(now).position.km
    distance_km          = np.sqrt((vector_km**2).sum())
    dic                  = dict()
    answer_dict[name]  = dic
    dic['lat']  = lat
    dic['lon']  = lon
    dic['elev'] = elev
    dic['vector_km'] = vector_km
    dic['vector_distance_km'] = distance_km
    alt, az, d    = London.at(now).observe(place).apparent().altaz()
    dic['altaz']  = alt._degrees, az._degrees
    dic['dist']   = d.km

for name, dic in answer_dict.items():
    print(name, dic['altaz'], dic['dist'])