Zenith verification calculation: lat, long, time to zenith and vice versa

Astropy
,
1

I’ve posted this to stackoverflow, but if anyone here has an answer, I’d be truly grateful:

python - Astropy zenith verification calculation: lat, long, time to zenith and vice versa - Stack Overflow

I am trying to calculate the zenith equatorial coordinates (right ascension & declination) for a given latitude, longitude, and observing time. Using astropy in python, I have the following:

from datetime import datetime
from astropy.coordinates import EarthLocation, SkyCoord, AltAz
from astropy.time import Time
from astropy.io import fits
from astropy import units as u


# coordinates on Earth
lat, long = 40, 120
lat, long = lat * u.deg, long * u.deg
altitude = 14756 * u.m

# location of obervations
obs_location = EarthLocation.from_geodetic(lat=lat, lon=long, height=altitude)

# observation time, corresponding to 2022-12-08 at 03:22:37.000969
obs_time = Time(datetime([2022,12,8,3,22,37,969]),scale='utc',location=obs_location)

# RIGHT ASCENSION AND DECLINATION at zenith for given lat, long
right_ascension_z = obs_time.sidereal_time('mean','greenwich').to('deg')
declination_z = lat

>> right_ascension_z = <Longitude 127.52856085 deg>, declination_z = <Quantity 40. deg>

Now, the calculated RA and DEC should translate back to an altitude of 90 deg (pointing zenith) for that given observing time and observing location – or at least so I thought. Going backwards:

aa = AltAz(location=obs_location, obstime=obs_time)
coord_zenith = SkyCoord(ra=right_ascension_z,dec=declination_z,frame='icrs',unit=u.deg)
altaz_coord = coord_zenith.transform_to(aa)

print(altaz_coord)
>> <SkyCoord (AltAz: obstime=2022-12-08 03:22:37.000969, location=(-2452005.67593768, 4246998.41117133, 4087470.54616891) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
(317.79729856, 7.01695887)>

The altitude coordinates using AltAz should be 90 deg. What am I missing?

2

I see on Stack Overflow that you have already discovered for yourself that you should not specify 'greenwich' when calling sidereal_time(). That’s because the local sidereal time is the same as the right ascension of the local meridian and zenith, so just as the local sidereal time is different between Greenwich and your observer at 120 deg longitude, Greenwich has a different right ascension at zenith than your observer.

The remaining small discrepancy is due to differences in the definition of ICRS, AltAz, and the output of sidereal_time(), namely in the handling of precession, nutation, and stellar aberration.

Hope that helps!

1 Like