eci2aer cross check against AstroPy

[AshHarvey]

Thanks for porting over these MatLab codes! I was cross checking the eci2aer function with AstroPy before incorperating it into another project I am working on, but I must be missing something or these numbers are very far off...

from pymap3d.aer import eci2aer as eci2aer
from astropy.coordinates import SkyCoord, GCRS, EarthLocation, CartesianRepresentation, AltAz
from astropy import units as u
from datetime import datetime
import numpy as np

eci = np.asarray([-34320604.87943786,2010725.129018341, 2977059.1586334566]) # units in meters
lla = np.asarray([10.0,10.0, 100]) # units in degrees, degrees, meters
t = datetime(2020, 3, 15, 0, 0, 0) # time as a datetime object


sat = SkyCoord(x=eci[0]*u.m, y=eci[1]*u.m, z=eci[2]*u.m, frame=GCRS, 
               representation_type=CartesianRepresentation)
origin = EarthLocation(lat=lla[0]*u.deg, lon=lla[1]*u.deg, height=lla[2]*u.m)
sataltaz = sat.transform_to(AltAz(obstime=t,location=origin))

az, el, srange = eci2aer(eci[0], eci[1], eci[2], lla[0], lla[1], lla[2],t)

np.array(sataltaz.az)-az
#array([70.90364332])
np.array(sataltaz.alt)-el
#array([-120.90947987])
np.array(sataltaz.distance)-srange
#array([1.77086727e+11])

Is ECI ~ GCRS ? Any feedback is appreciated. If this is my error vice the projects, I will gladly delete this.

[AshHarvey]

I went back to the MatLab code and used one of their examples to check the eci2aer function...

from: https://www.mathworks.com/help/aerotbx/ug/eci2aer.html

mat_eci = np.asarray(1e08)*[-3.8454,-0.5099,-0.3255]
mat_t = datetime(1969,7,20,21,17,40)
mat_lla = np.asarray([28.4,-80.5,2.7])
mat_aer = eci2aer_degree(mat_eci[0],mat_eci[1],mat_eci[2], mat_lla[0],mat_lla[1],mat_lla[2],mat_t,useastropy=False)
mat_answer = np.asarray(1.0e+08) * (0.0000, 0.0000, 3.8401)

np.asarray(mat_aer).flatten() - mat_answer
#array([  161.96964942,    54.85339394, 20895.93954259])

The values are also significantly off...

[scivision]

Thanks for this report, apparently I haven't checked the ECI correctly.
I took a quick look and I do see a discrepancy with Matlab example from eci2ecef and ecef2eci (which are underlying the eci2aer and aer2ecef).

[scivision]

Does the latest master commit a08dd0f work for you? I didn't try your test cases, but had better than 1% accuracy vs. Matlab Aero toolbox. I believe the prior computations used had too large error due to too few factors considered.

[scivision]

This should be fixed in PyMap3d v2.4.0. Please let us know if any further issues or this issue persists.

[AshHarvey]

That didn't do anything to resolve this issue and avoiding astropy's complex class structure is most of the motivation for trying to use this package.

[scivision]

OK I didn't expose any of Astropy's internals, just plain numbers or ndarray in/out.

Another point of comparison is Julia: https://juliaspace.github.io/SatelliteToolbox.jl/stable/man/transformations/ecef_eci/index.html

just for reference: https://earth-info.nga.mil/GandG/publications/tr8350.2/tr8350.2-a/Appendix.pdf

Although I don't have the time right today, I will leave this here for when I get a moment to compare a simple non-Astropy transform with AstroPy and Julia so that we don't give transformations with unreasonably large error--which the non-Astropy transform I had before seemed to have too much error.

[scivision]

v2.4.1 of PyMap3d restores the Numpy-only, rotation-only ECI conversion option by setting option "use_astropy=False" to the eci2* or *2eci functions. The non-AstroPy conversions have several percent error since they're rotation only. The default is to use AstroPy if available for ECI.

[scivision]

Please feel free to reopen if there's any remaining issues. Thanks for your helpful feedback.

[AshHarvey]

Thanks Michael; I ended up writing some functions that use the ERFA implementation of SOFA and IAU's EOPS to perform GCRS // ITRS transformations as well as Az El computation using the ITRS frame. If they might be any help to you, you can find a draft here - https://github.com/AshHarvey/astro-trans.

[scivision]

Okay that's interesting, I'll keep that in mind for comparison with what I do currently. I was aiming for the least lines of code but that impacts numerical accuracy as well as computation time.