Repository navigation
Add short worked examples to demo intended usage for each built-in coordinate frame #11701
Description
I am attempting to simulate visibility of a low-Earth-orbit satellite from a fixed ground station. The satellite elevation angles are all coming out around -17 degrees. This is not believable. As the satellite travels around the earth, the elevation angles should vary widely, and the satellite should at least occasionally be visible to the ground station (positive elevation angle). I'm assuming that I'm doing something wrong, but it is possible that there is a documentation issue here. My code follows.
# This script calculates visibility of a satellite to a single ground station.
# Section 1: Imports and constants.
from numpy import arange, empty, pi
from pdb import set_trace
import cartopy.crs as ccrs
from cartopy.feature.nightshade import Nightshade
from astropy import coordinates, time, units as u
from astropy.coordinates import AltAz, CartesianRepresentation, EarthLocation
from poliastro.bodies import Earth
from poliastro.twobody import Orbit, propagation
from poliastro.frames.fixed import ITRS
from poliastro.frames.equatorial import GCRS
RAD2DEG= 180 / pi
# Section 2: Simulate orbit.
a = 7578.1 * u.km # semi-major axis
ecc= 0 * u.one # eccentricity
inc= 63.4 * u.deg # inclination of the orbit plane
raan= 50 * u.deg # right ascension of the ascending node
argp= 270 * u.deg # argument of perigee
nu = 0 * u.deg # true anomaly at epoch
T_sim_hrs= 24
T_step_min= 3
epoch= time.Time("2021-05-05 10:43")
orbit= Orbit.from_classical(Earth, a, ecc, inc, raan, argp, nu, epoch=epoch)
# Create time steps covering the duration of the simulation:
times= arange(0, 60*T_sim_hrs+1, T_step_min) * u.min
xyz= propagation.propagate(
orbit,
time.TimeDelta(times),
method=propagation.cowell,
rtol=1e-10,
)
epoch_times = epoch + times
# Calculate Earth-fixed (non-inertial) coordinates `itrs_xyz`:
gcrs_xyz = GCRS(xyz, obstime=epoch_times, representation_type=CartesianRepresentation)
itrs_xyz = gcrs_xyz.transform_to(ITRS(obstime=epoch_times))
# Convert satellite positions from Cartesian to spherical coordinates (radians and km):
sat_lat_lon= coordinates.cartesian_to_spherical(itrs_xyz.x, itrs_xyz.y, itrs_xyz.z)
# Section 3: Calculate satellite elevation angles.
# Reference: https://learn.astropy.org/rst-tutorials/Coordinates-Transform.html#transform-frames-to-get-to-altitude-azimuth-altaz
ground_station= EarthLocation.from_geodetic(
lat = 30 * u.deg,
lon = 30 * u.deg,
height= 0 * u.km)
elevation= empty(len(sat_lat_lon[0]))
for i in range(len(sat_lat_lon[0])):
satellite = EarthLocation.from_geodetic(
lat =sat_lat_lon[1][i],
lon =sat_lat_lon[2][i],
height=sat_lat_lon[0][i])
# In the following call to `transform_to`, it is unclear why time is needed as an input.
result= satellite.get_itrs().transform_to(AltAz(obstime=time.Time('J2010'),
location=ground_station))
elevation[i]= result.alt.value
set_trace()
Activity
pllim commentedon May 5, 2021
Are you sure this is not an issue for
poliastro? cc @astrojuanluastrojuanlu commentedon May 5, 2021
Thanks @pllim for the ping. @Phillip-M-Feldman this is a super requested use case, I'll have a deeper look later this week. We need to isolate the problem - I trust Astropy more than poliastro :)
Phillip-M-Feldman commentedon May 5, 2021
ayshih commentedon May 5, 2021
First off, your transformation code can be much simpler. You should be able to go from
gcrs_xyzstraight toAltAz, e.g.:Your direct issue is that, in your later code, you are creating your
AltAzframe with the fixedobstime=time.Time('J2010'), which is different from the defaultobstime=time.Time('J2000')that comes out ofsatellite.get_itrs(). That shift inobstimeis why all of your satellite positions are being transformed incorrectly.mkbrewer commentedon May 5, 2021
That's what I thought also. I found this attempted transformation quite curious. I also noticed that in:
The height here is wrt to the WGS84 ellipsoid, not the geocenter.
mkbrewer commentedon May 5, 2021
I still can't figure out why all the elevations would come out the same. It seems to me that you'd have an array of positions that circle the Earth several times at an altitude of 7578 km above the ellipsoid and then are asking how these would be viewed from a particular location on Earth at 0 hours UT on January 1st, 2010, so the result should be a variety of elevations.
ayshih commentedon May 5, 2021
I imagine the reason is that the barycentric location of Earth is different between J2000.0 and J2010.0 (by ~900 Mm), so you're looking at an array of positions that are much farther away from Earth than their breadth.
mkbrewer commentedon May 5, 2021
I should think that the transformation to AltAz would handle the time difference through a CIRS->CIRS transform that goes back to the ICRS to take into account the difference in the Earth's position wrt the SSB.
mkbrewer commentedon May 5, 2021
Oh, it might be that
from_geodetic()is expecting degrees and is given radians. This would result in all of the locations being within 6.28 degrees of each other in longitude.ayshih commentedon May 6, 2021
It doesn't null out the difference in Earth's position. That is,
satellite.get_itrs()defaults to anobstimeof J2000.0, so corresponds to the inertial position of the satellite in the vicinity of Earth at J2000.0, but you're calculating the vector to that from an observer based on Earth at J2010.0. Here's some code to illustrate, using justground_stationas defined above:Note the elevation of ~-17 degrees, and the distance of over a million kilometers.
Phillip-M-Feldman commentedon May 6, 2021
pllim commentedon May 6, 2021
Is the issue resolved? If so, please close it. Thanks, everyone!
mkbrewer commentedon May 6, 2021
This is quite an interesting discussion. There may be another issue involved here. I will look into it.
mkbrewer commentedon May 6, 2021
OK. I understand this now. While, intuitively, the relative positions of two locations in the ITRS have no time dependence, that is not how it works in Astropy. Positions in Astropy are ultimately positions wrt the solar system barycenter at a given point in time. Therefore, specifying an ITRS position at one time and then observing it at another time, will, in fact, include the change in the position of the Earth wrt the SSB. As I suspected, there are a couple of other issues here.
As far as I can see, the output of
propagation.propagate()is an astrometric position wrt the geocenter rather than a true geocentric GCRS position. The later includes aberration applied to the vector from the geocenter to the target body. The former does not. Therefore, just doing the transform from GCRS to AltAz as suggested will incorrectly remove aberration from the geocentric position vector before adding it to the topocentric position vector from the ground station to the target body.Also, the target body is moving rather fast, so neglecting the light time delay may result in a significant error in the observed position. Unfortunately, Astropy does not yet have an astrometric frame, so the way to do this rigorously is to calculate the time delayed ICRS position of the satellite at each time of observation and then transform from ICRS to AltAz at the location of the ground station.
This can be done by calling
ground_station.get_gcrs_posvel(obstime=epoch_times)to get the astrometric position of the observer wrt the geocenter. Then for each time of observation, take the difference in the two positions and calculate the distance from the ground station to the satellite, iteratively retard the satellite orbit by the corresponding light time, and finally add the position of the Earth wrt the SSB at each time of observation to the time delayed astrometric position of the satellite.mhvk commentedon May 7, 2021
@mkbrewer - it sounds like the above is another reason to actually ensure astropy has an astrometric frame. Since this issue should be closed (I think), would you mind transferring your text to the larger discussion so that we do not lose track of it?
mkbrewer commentedon May 7, 2021
Done.
mhvk commentedon May 7, 2021
@Phillip-M-Feldman - While there does not seem to be a true bug, it clearly is sufficiently complicated that it is not always clear how to handle things. Is there a place in the documentation where a bit of extra test might have helped you? Was there a specific place you looked/searched for?
Phillip-M-Feldman commentedon May 7, 2021
ayshih commentedon May 7, 2021
"test" --> "text"
astrojuanlu commentedon May 8, 2021
I wonder where such documentation should be, though. Not against adding stuff specific to artifical satellite tracking to the Astropy docs, but perhaps poliastro would be a better place, given that it also has groundtrack plotting and orbital element handling. What are people's thoughts on this?
mhvk commentedon May 8, 2021
@astrojuanlu - I agree that poliastro might make more sense in this example.
For astropy, what I think could be really useful a set of simpler examples, aimed at specific coordinate frames. For instance, while working on the new HADec alternative to AltAz, I realized there were no simple examples of that (there is an observation planning page).
[-]elevation angles are not believable[/-][+]Add short worked examples to demo intended usage for each built-in coordinate frame[/+]adrn commentedon Sep 22, 2026
The original issue raised here is not a bug, but the final discussion points out that it would be good to have short, self-contained examples for all coordinate frames that demonstrate how they should be used for specific use cases.
One other note relevant to the above: in the meantime, poliastro is no longer maintained and has been archived.