Search code examples
pythongpscoordinate-systemspyephem

Galactic coordinates for point in sky from GPS-(lon, lat, timestamp) and (azimuth, zenith)?

Given the following information: Observer is located at longitude and latitude given by GPS. Looks at an object in the sky at an azimuth and zenith on a time given by GPS.

How can I get the Galactic coordinates for the point on the sky? Can this be done with PyEphem?

I think this gets me nearly there:

from datetime import datetime
import ephem
observer = ephem.Observer()
observer.pressure = 0
observer.lon = str(4.95569830927)
observer.lat = str(52.3545603701)
observer.elevation = 56.8426188165
observer.date = ephem.Date(datetime.utcfromtimestamp(1409097608))
ra, dec = observer.radec_of(str(azimuth), str(zenith))

But how to get the galactic coordinates from this.

Solution

  • You can generate galactic coordinates using the coordinate transform objects described in the Quick Reference:

    http://rhodesmill.org/pyephem/quick.html#coordinate-conversion

    The following script takes an arbitrary azimuth and zenith angle and transforms them to galactic coordinates, and the answer agrees with the NASA calculator at http://lambda.gsfc.nasa.gov/toolbox/tb_coordconv.cfm

    from datetime import datetime
import ephem

azimuth = '45:00'
zenith = '84:00'
ninety = ephem.degrees('90:00')

observer = ephem.Observer()
observer.pressure = 0
observer.lon = str(4.95569830927)
observer.lat = str(52.3545603701)
observer.elevation = 56.8426188165
observer.date = ephem.Date(datetime.utcfromtimestamp(1409097608))
ra, dec = observer.radec_of(str(azimuth), ninety - ephem.degrees(str(zenith)))
# print ra
# print dec

e = ephem.Equatorial(ra, dec, epoch='2000')
g = ephem.Galactic(e)
print g.lon
print g.lat

    Three quick notes:

    1. Note that PyEphem always talks about “altitude and azimuth” not “zenith distance and azimuth” so the script needs to subtract the zenith distance from 90° to produce an altitude.

    2. The tiny remaining error between these results and NASA’s are, I think, inherent in the libastro library upon which PyEphem is based.

    3. If you want to try that NASA calculator yourself, note that it wants RA expressed in degrees (!), which you can force PyEphem to display with the statement print ephem.degrees(ra).