1

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?

gerrit
  • 24,025
  • 17
  • 97
  • 170

2 Answers2

2

Cleaning up your inputs a bit:

from astropy import units as u
from astropy.coordinates import AltAz, EarthLocation

ground = EarthLocation.from_geodetic(lat=3*u.deg, lon=5*u.deg, height=0*u.km)
sat    = EarthLocation.from_geodetic(lat=3*u.deg, lon=5*u.deg, height=700*u.km)

result = sat.get_itrs().transform_to(AltAz(location=ground))
print(result)

Which produces the expected 90 degrees:

<AltAz Coordinate (obstime=None, location=(6345.21668424, 555.13452749, 331.57431534) km, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt, distance) in (deg, deg, km)
    (89.99999672, 89.99991123, 700.)>
JCOidl
  • 452
  • 1
  • 5
  • 12
1

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]))
gerrit
  • 24,025
  • 17
  • 97
  • 170