[GIS] Local Cartesian (NED) to Geodetic (Lat,Lon,Alt) Conversion – Max Distance from Reference Point

coordinate systemcoordinatesgeodesymatlab-mapping-toolbox

Looking at this:

https://www.mathworks.com/help/map/ref/ned2geodetic.html:

[lat,lon,h] = ned2geodetic(xNorth,yEast,zDown,lat0,lon0,h0,spheroid))

The obvious question is by how much the current location on local Cartesian NED (xNorth,yEast,zDown) can be far away from the reference point (lat0,lon0,h0)?

For example if the reference point is in Australia and the current location is in the US, does those "mathematical tricks" still will be valid?

I assume that there is an intermediate step of-

Local Cartesian NED to geocentric ECEF: https://www.mathworks.com/help/map/ref/ned2ecef.html

Since they also using Geodetic reference point – probably first converting the Geodetic reference point to ECEF, and then there is the questioned "mathematical tricks" to "add" current local Cartesian NED to that reference point (already represented in ECEF) so the current local Cartesian NED will be in ECEF – Thus the "mathematical tricks" are some sort of Cartesian/Vectorial manipulations, I suspect how long they holding in the context of maps, geodesy etc?

Regarding converting Geodetic reference point into ECEF in the intermediate step mentioned above, and in the final step converting back current Cartesian point (already in ECEF) to Geodetic, they use:

Transform geocentric (ECEF) to geodetic coordinates: https://www.mathworks.com/help/map/ref/ecef2geodetic.html

Assuming this is pretty much straight forward formulas with known limitations (Characteristics of different methods to convert ECEF to LLA) So I don't see how they interfere here…

Best Answer

I don't think there are practical limitations to ned2geodetic.

The conversion from LLA to ECEF coordinates has closed form solution.
The rotation of NED vector to ECEF has closed form solution.
Vector addition in ECEF is trivial.
Conversion of ECEF coordinates to LLA as implemented in Matlab is not closed form solution. However, the algorithm employed is documented to have good accuracy and well behaved, unless the point lies close to the center of the earth (~50km).

Related Solutions

[GIS] Lat/Lon/Alt to spherical or cartesian coordinates

Note that "Lat/Lon/Alt" is just another name for spherical coordinates, and phi/theta/rho are just another name for latitude, longitude, and altitude. :) (A minor difference: altitude is usually measured from the surface of the sphere; rho is measured from the center -- to convert, just add/subtract the radius of the sphere.)

To convert phi/theta/rho to cartesian x/y/z where (0,0,0) is the center of the earth and rho is the distance from the center:

    # python
    x = math.cos(phi) * math.cos(theta) * rho
    y = math.cos(phi) * math.sin(theta) * rho
    z = math.sin(phi) * rho # z is 'up'

(Note there's some slightly arbitrary choices here in what each axis means... you might want 'y' to point at the north pole instead of 'z', for example.)

The inverse, and a lot more information, can be found at Wikipedia: http://en.wikipedia.org/wiki/Spherical_coordinate_system

Coordinate System – Converting Geodetic Lat/Lon to ECEF Cartesian and Applying Geoid Undulation

Considering how the coordinates were originally arrived at confirms that what you are doing is correct.

Look at point P1 in the figure. To determine where it is (without a GPS signal), the first step would be to identify where the up-down (vertical) direction is. Because gravity determines that, this direction must be normal to the geoid. Therefore P1 is projected downward until the geoid is reached. The distance of that projection is the orthometric height, H. Let the point of projection on the geoid be G1: it is located where the "H" and "N" vectors meet in the diagram.

However, P1 is actually located using latitude and longitude coordinates relative to the ellipsoid. Its height above the ellipsoid is h, the "ellipsoidal height" or geodetic elevation. The direction in which P1 is projected to the ellipsoid for this positional determination is the normal to the ellipsoid, not the geoid. The line along which P1 is projected will intersect the geoid at some location, but probably not exactly at G1. In the figure, the intersection lies to the left of G1 (just above the letter "h"). Let's call this G1'.

The confusing part is that the geoid's height for point G1 is determined by projecting G1--not G1'--down to the ellipsoid. Let us therefore distinguish two points on the ellipsoid: E1, the projection of G1, and E1', the projection of G1'. Thus the geodetic height h goes from E1' to P1, the ellipsoid undulation N goes from E1 to G1, and the orthometric height H from G1 to P1. Unless G1 and G1' coincide, it sure looks like h is not quite the same as H+N. How much of an error is made by assuming h=H+N?

Because the geoid is so smooth and so close to the ellipsoid (it almost always lies within 100 meters of the ellipsoid, which amounts to less than 0.0015% of the earth's radius), the difference between the ellipsoid's normal and the geoid's normal is a few seconds or less in most places, reaching an extreme of 100 seconds. A one-second (5e-6 radian) deflection translates to a positional error of five parts per million. Therefore an extreme deflection of 100 seconds would create an error of 500 parts per million. At large heights--say 100 kilometers, just to be extreme about this calculation--the position of G1 could be as much as 500e-6 * 100 Km = 50 meters away from G1'. (At terrestrial heights and more typical vertical deflections the distance would usually be less than a meter.) Accordingly, E1 and E1' would be separated by a similar distance (because G1 and G1' are only a few meters above or below the ellipsoid.)

Because the geoid undulates smoothly and E1 and E1' are so close together, its height at E1 scarcely varies from its height at E1': the difference in heights could not exceed a millimeter (and normally would be just a matter of microns). Similarly, the normals to the geoid at G1 and G1' would essentially be the same; they would differ by no more than a thousandth of an arc-second. Using one of those normals in place of the other would introduce a position error of no more than a few hundred parts per billion. Over the 100 km distance back to the original point P1, that would create a positional error of a few hundred times 1e-9 times 100 km, which would be just a few centimeters. For points on the earth's surface the positional error would be orders of magnitude less than that: millimeters or smaller.

The upshot is that because the positional errors are second order in the geoid undulations (that is, based on differences in the discrepancies between the geoid and the ellipsoid), for almost all purposes it is safe to treat G1 and G1' as coincident and h = H+N.

Best Answer

Related Solutions

[GIS] Lat/Lon/Alt to spherical or cartesian coordinates

Coordinate System – Converting Geodetic Lat/Lon to ECEF Cartesian and Applying Geoid Undulation

Related Question