I am writing some code dealing with an aircraft locations with respect to an observer located at ground and I got a discrepancy between the distance calculated directly from geocentric Cartesian coordinates of observer and aircraft and the distance returned by AltAz in astropy. The following minimum reproducible code illustrates the problem with the aircraft 1000 m above the observer.
#Problem calculating distance to object in the sky
#worse when close to observer. Negligible if distance > 100 km
from astropy.time import Time
from astropy import units as u
from astropy.coordinates import EarthLocation, AltAz
obs_lon = 0*u.deg
obs_lat = 0*u.deg
obs_alt = 0*u.m
aircraft_lon = 0*u.deg
aircraft_lat = 0*u.deg
aircraft_alt = 1000*u.m
obs_time = Time ('2000-01-01T00:00:00.000')
obs_location = EarthLocation.from_geodetic (obs_lon,obs_lat, obs_alt)
aircraft_location = EarthLocation.from_geodetic(aircraft_lon, aircraft_lat, \
aircraft_alt)
aircraft_altaz = aircraft_location.get_gcrs(obstime=obs_time). \
transform_to(AltAz(location=obs_location, obstime=obs_time))
dist = ((aircraft_location.x - obs_location.x)**2 + \
(aircraft_location.y - obs_location.y)**2 + \
(aircraft_location.z - obs_location.z)**2)**0.5
print("distance calculated from xyz", dist) # this gives 1000 m, as expected
print("distance returned from AltAz", aircraft_altaz.distance) # 1189.4 m ?
I was expecting to obtain similar results, but they are quite different. The difference is bigger for close objects compared to distant ones. It also depend strongly on time (especially time of day rather than date, so the problem looks related to Earth's rotation). Dependence with latitude seems to be quite small. Replacing gcrs by itrs produces the same behaviour.
I finally found a solution. I post it just in case it becomes useful for somebody else in the future: apparently the problem is related to how astropy deals with stellar aberration correction. For near-Earth objects this correction should not be included, and this can be achieved using a topocentric ITRS frame. Following this link, the steps to follow are:
Code:
#Problem calculating distance to object in the sky
#worse when close to observer. Negligible if distance > 100 km
from astropy.time import Time
from astropy import units as u
from astropy.coordinates import EarthLocation, AltAz, ITRS
obs_lon = 0*u.deg
obs_lat = 0*u.deg
obs_alt = 0*u.m
aircraft_lon = 0*u.deg
aircraft_lat = 0*u.deg
aircraft_alt = 1000*u.m
obs_time = Time ('2000-01-01T00:00:00.000')
obs_location = EarthLocation.from_geodetic (obs_lon,obs_lat, obs_alt)
aircraft_location = EarthLocation.from_geodetic(aircraft_lon, aircraft_lat, \
aircraft_alt)
aircraft_altaz_old = aircraft_location.get_gcrs(obstime=obs_time). \
transform_to(AltAz(location=obs_location, obstime=obs_time))
dist = ((aircraft_location.x - obs_location.x)**2 + \
(aircraft_location.y - obs_location.y)**2 + \
(aircraft_location.z - obs_location.z)**2)**0.5
#New block of code solving the issue:
aircraft_itrs = aircraft_location.get_itrs(obstime = obs_time)
obs_itrs = obs_location.get_itrs(obstime = obs_time)
aircraft_itrs_relative = aircraft_itrs.cartesian.without_differentials() - \
obs_itrs.cartesian
aircraft_topo = ITRS(aircraft_itrs_relative, obstime = obs_time, \
location = obs_location)
aircraft_altaz = aircraft_topo.transform_to(AltAz(obstime = obs_time, \
location=obs_location))
print("dist calculated from xyz", dist) # this gives 1000 m, as expected
print("dist returned from AltAz (old)", aircraft_altaz_old.distance) # 1189.6 m
print("dist returned from AltAz (new)", aircraft_altaz.distance) # 10000 m
print("az, alt ",aircraft_altaz.az.to_value(), aircraft_altaz.alt.to_value())
#now distance, az and alt are OK