Coordinate mismatch between PyEphem/SSCWeb and this library

[letmaik]

# TLE for ISS
t1 = '1 25544U 98067A   12035.27140256  .00008132  00000-0  10934-3 0  9453'
t2 = '2 25544 051.6417 073.7902 0020652 004.6179 085.0245 15.58669131757138'
tle = twoline2rv(t1, t2, wgs72)
position, _ = tle.propagate(2012,2,4,7,55,00)
print position
# (-2395.0416350887344, 4440.7514261050846, 4492.3844513630256)

When I use http://sscweb.gsfc.nasa.gov/ I get the following coordinates: -2377.74 4447.21 4495.19

These agree with the ones I get when using PyEphem for the same data.

Any ideas?

[brandon-rhodes]

I have never used the SSCWeb before, but am happy to have another source of coordinates to compare sgp4 against! Which options did you select, on which of the site’s forms, to get the numbers that you have quoted? (I tried clicking around a bit, but could not quickly get coordinates like those to appear.) If we can figure out exactly what coordinate system they have been expressed in, then I can compare them to the coordinates output by the SGP4 algorithm here in Python.

[letmaik]

So, first go to http://sscweb.gsfc.nasa.gov/cgi-bin/Locator.cgi, then select ISS, set start time to 2012 35 7:54:0, end time to 2012 35 7:56:0. Click on Output options and make a tick for GEI/J2000 in the XYZ column. Click on the Output units/formatting button and change date and time formats to a more human format. Also change distance to kilometers. Finally click Submit query and wait for output.

[brandon-rhodes]

Well, that could be an initial source of disagreement: the SGP4 algorithm outputs coordinates in TEME (True equinox / Mean equator) coordinates, and an additional transformation is required to create J2000 coordinates. I am actually about to try implementing it over in my Skyfield project. Take a look at Appendix C of the recent SGP4 paper to see how the transformation looks:

http://www.celestrak.com/publications/AIAA/2006-6753/AIAA-2006-6753-Rev2.pdf

[brandon-rhodes]

(Once I have the transform complete, I will check and see what it does for this set of coordinates that you have generated, and I will report back here on whether it helps the agreement or not!)

[letmaik]

So, is PyEphem also calculating TEME and then converting it to J2000?

[brandon-rhodes]

Well, you now lead me to an interesting question! :) How, exactly, are you getting PyEphem to give you x,y,z coordinates for an Earth that can agree with what you are getting from the SSCWeb? I had thought that the C routines were not exposing the x,y,z coordinate in a form that Python code could access.

[letmaik]

Good question! After tle.compute(date), first, I calculate the distance from the earth:

from astropy import units
# see http://stackoverflow.com/q/19426505/60982
earthRadiusPyEphem = 6378.16 * units.km
distanceEarthCenter = tle.elevation*units.m + earthRadiusPyEphem

Then I transform RA, Dec, distance to XYZ by using astropysics. I do that by pretending I have GCRS angular coordinates which I then transform to GCRS rectangular coordinates:

issPos = coords.GCRSCoordinates(coords.AngularCoordinate(tle.ra, radians=True),
                                      coords.AngularCoordinate(tle.dec, radians=True),
                                      None, None, 2000, distanceEarthCenter.to(units.pc).value)
issPosRect = issPos.convert(coords.RectangularGCRSCoordinates)
issPosRect = units.pc.to(units.km, [issPosRect.x, issPosRect.y, issPosRect.z])

In the end, this astropysics stuff is just so that I don't have to write the following in my code (which feels more low level):

x = r*cos(ra)*cos(dec)
y = r*sin(ra)*cos(dec)
z = r*sin(dec)

I just copied the above from the corresponding conversion code and haven't tried it separately without astropysics, but I guess it will lead to the same values.

[brandon-rhodes]

I cannot find a significant difference between the predictions of this sgp4 library and PyEphem; could you help me tweak the following code so that I see the same result that you are seeing, so that I can investigate the mismatch?

t1 = '1 25544U 98067A   12035.27140256  .00008132  00000-0  10934-3 0  9453'
t2 = '2 25544 051.6417 073.7902 0020652 004.6179 085.0245 15.58669131757138'

from sgp4.earth_gravity import wgs72
from sgp4.io import twoline2rv
tle = twoline2rv(t1, t2, wgs72)
print 'sgp4:', tle.propagate(2012,2,4,7,55,00)[0]

from ephem import readtle
from math import cos, sin
sat = readtle('ISS', t1, t2)
sat.compute('2012/2/4 7:55:00')
r = 6378.16 + sat.elevation / 1000.0
x = r * cos(sat.ra) * cos(sat.dec)
y = r * sin(sat.ra) * cos(sat.dec)
z = r * sin(sat.dec)
print 'PyEphem:', x, y, z

The output that I get from the two libraries is only off by a few tens of meters, which is within the several-kilometer error bar of the SGP4 theory itself:

sgp4: [-2395.041635088735, 4440.751426105086, 4492.3844513630265]
PyEphem: -2395.05093619 4440.76887896 4492.4020381

[letmaik]

Well that's interesting. I copied your test code 1-to-1 and get the following output:

sgp4: [-2395.041635088735, 4440.751426105086, 4492.3844513630265] PyEphem: -2377.75135336 4447.22285846 4495.20640728

Just for the record, I use: sgp4==1.1 pyephem==3.7.5.1

Am 07.12.2013 15:28, schrieb Brandon Rhodes:

I cannot find a significant difference between the predictions of this |sgp4| library and |PyEphem|; could you help me tweak the following code so that I see the same result that you are seeing, so that I can investigate the mismatch?

t1 = '1 25544U 98067A 12035.27140256 .00008132 00000-0 10934-3 0 9453' t2 = '2 25544 051.6417 073.7902 0020652 004.6179 085.0245 15.58669131757138'

from sgp4.earth_gravity import wgs72 from sgp4.io import twoline2rv tle = twoline2rv(t1, t2, wgs72) print 'sgp4:', tle.propagate(2012,2,4,7,55,00)[0]

from ephem import readtle from math import cos, sin sat = readtle('ISS', t1, t2) sat.compute('2012/2/4 7:55:00') r = 6378.16 + sat.elevation / 1000.0 x = r * cos(sat.ra) * cos(sat.dec) y = r * sin(sat.ra) * cos(sat.dec) z = r * sin(sat.dec) print 'PyEphem:', x, y, z |

The output that I get from the two libraries is only off by a few tens of meters, which is within the several-kilometer error bar of the SGP4 theory itself:

sgp4: [-2395.041635088735, 4440.751426105086, 4492.3844513630265] PyEphem: -2395.05093619 4440.76887896 4492.4020381

— Reply to this email directly or view it on GitHub https://github.com/brandon-rhodes/python-sgp4/issues/3#issuecomment-30055953.

[brandon-rhodes]

It turns out that my PyEphem is giving the right coordinates because of a fix that I made but have not released!

https://github.com/brandon-rhodes/pyephem/issues/14

So we now know what the difference is (if I understand the above issue's description correctly): SSCWeb must be precessing the coordinates (to J2000), while the raw SGP4 algorithm—like satellite observations themselves—is referenced to the Earth's surface and thus to something like epoch-of-date.

[letmaik]

Ok, so... do you have an idea how I could transform these coordinates to J2000?

On 09/12/13 05:47, Brandon Rhodes wrote:

It turns out that my PyEphem is giving the right coordinates because of a fix that I made but have not released!

brandon-rhodes/pyephem#14 https://github.com/brandon-rhodes/pyephem/issues/14

So we now know what the difference is (if I understand the above issue's description correctly): SSCWeb must be precessing the coordinates (to J2000), while the raw SGP4 algorithm—like satellite observations themselves—is referenced to the Earth's surface and thus to something like epoch-of-date.

— Reply to this email directly or view it on GitHub https://github.com/brandon-rhodes/python-sgp4/issues/3#issuecomment-30105550.

[letmaik]

Just found this: "We recommend converting TEME to a truly standard coordinate frame before interfacing with other external programs. The preferred approach is to rotate to PEF using Greenwich Mean Sidereal Time (GMST), and then rotate to other standard coordinate frames. Conversions are well documented from this point. [...]" (Revisiting Spacetrack Report 3, http://www.fas.org/spp/military/program/track/spacetrack3.pdf)

On 09/12/13 05:47, Brandon Rhodes wrote:

It turns out that my PyEphem is giving the right coordinates because of a fix that I made but have not released!

brandon-rhodes/pyephem#14 https://github.com/brandon-rhodes/pyephem/issues/14

So we now know what the difference is (if I understand the above issue's description correctly): SSCWeb must be precessing the coordinates (to J2000), while the raw SGP4 algorithm—like satellite observations themselves—is referenced to the Earth's surface and thus to something like epoch-of-date.

— Reply to this email directly or view it on GitHub https://github.com/brandon-rhodes/python-sgp4/issues/3#issuecomment-30105550.

[brandon-rhodes]

Yes, I have been following that paper closely as I work out how Skyfield should take bare SGP4 coordinates as produced by this library, and convert them to PEF and then to other coordinate frames. I should be making more progress tomorrow and will hopefully finish integrating Earth satellites into its GCRS coordinate system vectors; watch the commits over on that project for my progress.

And thank for you pointing out the SSCWeb site; I will try generating some results to integrate into Skyfield's tests to see if I get the same answers!

In the meantime: my suspicion is that this sgp4 library is returning correct TEME coordinates per the SGP4 algorithm, as I am seeing matches with PyEphem, with the SGP4 test suite itself, and with the TEME example in the paper you cite. Would that justify my closing this issue, since the sgp4 module is only supposed to run the SGP4 algorithm and then return an answer?

[letmaik]

Yes absolutely. The conversion from TEME to other frames should probably happen inside another more general library like astropy in the future. But in the meantime it's probably better to include it here. So yes, please close it, and thanks!

[letmaik]

Oh I just realized you're talking about Skyfield, I'll take a look at that. So ignore my comment about integrating the conversion here.

[cnwangfeng]

Dear brandon, I also met the issue of neothemachine. I do need a conversion from X,Y,Z to J2000. I try to read your codes in SKyfield. But it's a difficult job for me. May you tell me where I can find the code? I just download SGP4 new version. Quite good.

[cnwangfeng]

Brandon, I read the codes. Now I use SGP4 to obtain the vector value of ECEF. Am I right ? I then use NOVAS to transform to GCRS? However, I need the topocentric value . How can I do ?

[brandon-rhodes]

@cnwangfeng — what coordinate system are your x,y,z coordinates in?

[cnwangfeng]

I got X,Y,Z from SGP4. But I don't know the coordinate system of them? I guess it's ECEF? Am I right?

[brandon-rhodes]

Have you tried asking Skyfield to do the conversion? It uses SGP4 under the covers, but will do all of the conversions for you to either RA/dec or alt/az:

http://rhodesmill.org/skyfield/earth-satellites.html

[cnwangfeng]

Yes. I am reading the codes. I use JPL DE405 and NOVAS to calculate the position of planets in the solar system. I need the topocentric position (RA/DEC) of the target. In the following study, I have to calculate the satellite's topocentric position (RA/DEC). I just read the paper you suggested. the values outputted by SGP4 is TEME coordinate? Not J2000?

[brandon-rhodes]

Yes, SGP4 outputs TEME coordinates.

[hshe824]

Hi there, does anyone know what the underlying mathematical/physical model that is being used to calculate the position of a Satellite from reading a TLE file in the compute() function in PyEphem is? Is it SGP4? Thanks

[brandon-rhodes]

Yes, it's SGP4, but an older version of the algorithm that sometimes crashes. Here's a publication that describes the various bugs that are found in older versions like PyEphem's:

https://celestrak.com/publications/AIAA/2006-6753/AIAA-2006-6753.pdf