Skip to content

Coordinate conversion from ECEF to ECI does not match pymap3d #3

Closed

Description

@fPkX6F1nGTX

Try finding the coordinates of Null Island (geodetic latitude 0, longitude 0, elevation 0) at noon on January 1, 2000 as per the following code:

[ecef_x, ecef_y, ecef_z] = geodetic2ecef(0, 0, 0)
J2K0 = datetime(2000, 1, 1, 12, 0, 'TimeZone', 'UTCLeapSeconds')
[eci_x, eci_y, eci_z] = ecef2eci(ecef_x, ecef_y, ecef_z, J2K0)

Then run the same code via pymap3d and you will see the difference:
[J2K0_whole, J2k0_fraction] = py.matlab_datenum_to_datetime(J2K0)
J2K0_for_pymap3d = J2K0_whole + J2K0_fraction
eci_vector = py.pymap3d.ecef2eci(ecef_x, ecef_y, ecef_z, J2K0_for_pymap3d )

where matlab_datenum_to_datetime explicitly transforms the date as follows:

function [whole_days, fractional_days] = matlab_datenum_to_datetime(matlab_decimal_datenum):
    PYTHON_TO_MATLAB_CORRECTION_CONSTANT = td(days=366)
    whole_days = [dt.fromordinal(int(date)) for date in matlab_decimal_datenum]
    remainder_days = [number%1 for number in matlab_decimal_datenum]
    fractional_days_uncorrected = [td(days=remainder) for remainder in remainder_days]
    fractional_days = []
    for day in range(len(fractional_days_uncorrected)):
        next_fractional_uncorrected_day = fractional_days_uncorrected[day]
        fractional_days.append(next_fractional_uncorrected_day - PYTHON_TO_MATLAB_CORRECTION_CONSTANT)
end

I cannot figure out where the error is coming from, so I resorted to calling pymap3d from MATLAB since it uses astropy. You might want to get rid of this repo and just recommend that route if you cannot fix the issue.

Activity

  1. scivision commented on Sep 2, 2020

    @scivision
    Member

    This is expected from the simple ECI transform in this repo that does not account for nutation and other factors. Putting a position on the equator of a spheroid without those higher fidelity models in Astropy results in ECI also being on the equator.
    I am assuming the Astropy transformations are non-trivial, but I didn't investigate.

    The Matlab Aerospace Blockset has ecef2eci() if needed or yes use Python + Astropy.

  2. fPkX6F1nGTX commented on Sep 3, 2020

    @fPkX6F1nGTX
    ContributorAuthor

    I edited the MATLAB code above so it has fewer bugs. I do not have access to a MATLAB license so I cannot verify what it returns, but I did run the following Python code just now, for reference:

    import datetime as dt
    import pymap3d as pm
    
    J2K0 = dt.datetime(2000, 1, 1, 12, 0)
    pymap3d.ecef2eci(6378137, 0, 0, J2K0)
    (1158174.7900579234, -6272101.938428859, -143.09802467421062)
    >>>
    

    The magnitude of the vector returned is very, very close to 6378137, the X-coordinate of Null Island in ECEF.

  3. fPkX6F1nGTX commented on Sep 3, 2020

    @fPkX6F1nGTX
    ContributorAuthor

    @scivision If you could run the MATLAB code and report what it returns for Null Island's ECI coordinates on January 1, 2000 at noon, that would be greatly appreciated. The lat/long/alt is 0, 0, 0 and ECEF is 6378137, 0, 0; I would also appreciate knowing the difference between the ecef2eci answer and the geodetic2eci answer, if any (would depend whether lla2ecef returns the coordinates I mention above, as per https://planetcalc.com/7721/ (input 0 latitude to get radius of Earth at equator).

  4. locked and limited conversation to collaborators on Aug 11, 2022
  5. converted this issue into a discussion #5 on Aug 11, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions