Search code examples
pythonastropypyephemsatellitepyproj

Calculating zenith and elevation angles for satellite as seen from the ground


I have the latitude, longitude, and height for a satellite and for a ground observer. I am trying to calculate the satellite zenith angle and satellite azimuth angle, for the satellite, as seen from the ground.

I'm currently trying to solve this using astropy. The documentation describes how to calculate the zenith and azimuth angle for the Sun.

I tried to place the satellite straight above the observer, so that the zenith angle should be 0° (or the elevation angle 90°).

In [41]: ground = astropy.coordinates.EarthLocation(lat=3*u.deg, lon=5*u.deg)

In [42]: sat = astropy.coordinates.EarthLocation.from_geodetic(3*u.deg, 5*u.deg, 700*u.km)

However, clearly .transform_to(AltAz(...)) does not mean what I think it means, for the answer is not what I want:

In [43]: print(sat.get_itrs().transform_to(AltAz(location=ground)))
<AltAz Coordinate (obstime=None, location=(6345216.684243768, 555134.5274868822, 331574.3153428908) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0, obswl=1.0 micron): (az, alt, distance) in (deg, deg, m)
    (314.96287134, 63.2970069, 773970.24385554)>

I expected the elevation angle to be close 90°, not 63°.

There's also some routines in pyorbital but it doesn't seem to contain what I'm looking for.

I believe it should also be possible using pyephem, except that this seems to expect I get satellite information from a catalogue. I already have the lat, lon, altitude.

I've also looked at pyproj but I think this is only for objects at the geoid.

How can I calculate the altitude and azimuth for the satellite as observed from the ground?


Solution

  • You can get this from pyorbital. The relevant function is get_observer_look. For some reason, it wants the time, and it apparently requires all inputs are ndarray objects. It warns when the azimuth is ill-defined (such as when the zenith angle is zero). But it does work:

    In [50]: import pyorbital.orbital
    
    In [51]: print(pyorbital.orbital.get_observer_look(atleast_1d(1), atleast_1d(1), atleast_1d(700), datetime.datetime.now(), atleast_1d(1), atleast_1d(1), atleast_1d(1)))
    (array([261.11934085]), array([89.99999915]))
    
    In [53]: print(pyorbital.orbital.get_observer_look(atleast_1d(10), atleast_1d(0), atleast_1d(700), datetime.datetime.now(), atleast_1d(0), atleast_1d(0), atleast_1d(0.1)))
    /home/zmaw/u237009/.conda/envs/FCDR37a/lib/python3.7/site-packages/pyorbital/orbital.py:112: RuntimeWarning: divide by zero encountered in true_divide
      az_ = np.arctan(-top_e / top_s)
    (array([270.]), array([25.73173269]))
    
    In [54]: print(pyorbital.orbital.get_observer_look(atleast_1d(10), atleast_1d(0), atleast_1d(700), datetime.datetime.now(), atleast_1d(0.1), atleast_1d(0.1), atleast_1d(0.1)))
    (array([90.56944651]), array([26.03504444]))