[AstroPy] Confused by coordinate transforms ICRS/ITRS/GCRS

[Jeffrey Brent McBeth]

On Wed, Sep 02, 2015 at 03:01:18PM -0400, Erik Bray wrote:
> Hi Jeff,
> 
> On 09/02/2015 09:06 AM, Jeffrey Brent McBeth wrote:
> >
> > Is there further information I could provide, or another place I could ask?
> 
> I think often with questions like this it helps to show what output you're 
> getting (you wrote what you expected to get but not what you actually got). 
> Sometimes seeing the output (without having to run the code one's self) can help 
> in quickly seeing what's going wrong.

You are right, I should have included the output, I am sorry.

> I have very little knowledge of this code, but I can maybe offer a few hints below:

> I think this is a probably a bug, probably stemming from the fact that your 
> "earth" coordinates are ill-defined when converting to spherical coordinates 
> (which the code does).  That's just my guess but it seems like a decent one.  

If I weren't seeing the problem no matter where I am in near earth space (I'm simulating ISS and Hubble), I would definitely agree with you.  I probably should have called out a more complicated example too.

>  >>> almostEarth2.icrs.separation_3d(almostEarth3.icrs).to(u.m)
> <Distance 0.5502365049921195 m>
> 
> So off by a little more than a meter--I don't know how much error to expect but 
> that seems reasonable considering that this is converting from AU to m.  The 
> error seems to get worse the further out I go though, but still within the 
> correct order of magnitude, at least.

And this is where my "fun little experiment" breaks down.  The hubble WC3's, for example, each subtend .632x10^-6 radians.  So, I need vector directions accurate below that to stand a chance of matching up.  It looks like going anywhere near ICRS from ITRS or comparing two GCRS frames is downright dangerous.  

Which is a shame, I was hoping to be messing with aberration (hubble should be seeing ~200 pixels), and I can't get the coordinate transforms anywhere close to even that level of accuracy :(

Thank you for your help.  Unfortunately, this is pushing me toward having to code up these transforms on my own, which as we all know is fraught with peril.

Jeff

Below are the code and results if I change things to start 1000km above the ground

#+name: confusingTransforms
#+begin_src: python
import astropy
import astropy.coordinates
import astropy.time
import astropy.units as u

# For now, lets just talk about the center of the earth
# You may complain "this isn't in the sky", but I get similar results
# even when I'm well away from the earth, this keeps things simple, and since
# I'm ultimately trying for orientation transforms, doesn't matter
earth = astropy.coordinates.SkyCoord(x=0,y=0,z=1000,frame='itrs',unit=u.km)
# And a spot one meter away
# This goes wrong with larger distances too, so it isn't simply precision
almostEarth = astropy.coordinates.SkyCoord(x=0,y=0,z=2000,frame='itrs',unit=u.km)

# I would expect _most_ prints below to be 1000 kilometer
print 'ITRS/ITRS:',earth.separation_3d(almostEarth).to(u.km)
print 'ICRS/ITRS:',earth.icrs.separation_3d(almostEarth).to(u.km)
print 'ICRS/ICRS:',earth.icrs.separation_3d(almostEarth.icrs).to(u.km)

print 'GCRS/ITRS:',earth.gcrs.separation_3d(almostEarth).to(u.km)
print 'GCRS/GCRS:',earth.gcrs.separation_3d(almostEarth.gcrs).to(u.km)

now = astropy.time.Time.now()
later = now+1*u.second

# A notional ISS, I have the real numbers, but this keeps it simple
platform = astropy.coordinates.GCRS(obsgeoloc=[0,0,400]*u.km,obsgeovel=[0,0,7.5]*u.km/u.s,obstime=now)

platformLater = astropy.coordinates.GCRS(obsgeoloc=[0,0,400]*u.km,obsgeovel=[0,0,7.5]*u.km/u.s,obstime=later)

print 'GCRS(Now  )/GCRS(Now  ):',earth.transform_to(platform).separation_3d(almostEarth.transform_to(platform)).to(u.km)
print 'GCRS(Later)/GCRS(Later):',earth.transform_to(platformLater).separation_3d(almostEarth.transform_to(platformLater)).to(u.km)

# These two, due to earth motion should probably be around 1000+/-30km
print 'GCRS(Now  )/GCRS(Later):',earth.transform_to(platform).separation_3d(almostEarth.transform_to(platformLater)).to(u.km)
print 'GCRS(Later)/GCRS(Now  ):',earth.transform_to(platformLater).separation_3d(almostEarth.transform_to(platform)).to(u.km)
#+end_src

#+begin_example
ITRS/ITRS: [ 1000.] km
ICRS/ITRS: [  2.44988467e+08] km
ICRS/ICRS: [ 390.88978995] km
GCRS/ITRS: [ 1000.] km
GCRS/GCRS: [1000.] km
GCRS(Now  )/GCRS(Now  ): [ 340.79708459] km
GCRS(Later)/GCRS(Later): [ 340.79711114] km
GCRS(Now  )/GCRS(Later): [ 59479316.98264851] km
GCRS(Later)/GCRS(Now  ): [ 59479320.25267375] km
#+end_example
-- 
"The man who does not read good books has no advantage over 
 the man who cannot read them."
 -- Mark Twain