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

coordinate system

I'm working on modeling the earths surface in 3D, using a combination of WGS84/EGM96 and NASA SRTM orthometric height data.

To calculate a final cartesian XYZ coordinate for a point on the earth surface, I am using the following methodology:

I convert a given (geodetic) Lat/Lon to a point on the WGS84 ellipsoid using standard math. To convert that position to a point on the surface of the geoid, I offset along the ellipsoid normal from my ellipsoid surface point by whatever the geoid undulation is at that location according to EGM96.

In order to come to a final surface point of terrain, I then add the SRTM orthometric height to my geoid surface point, along the normal to the geoid (ie not along the ellipsoid normal!).

I'm not sure if I am offsetting from ellipsoid surface to geoid surface correctly by adding the height along the ellipsoid normal. If you look at the following diagram I found, they seem to be adding the geoid undulation along some arbitrary direction (ie not normal to ellipsoid or geoid, just straight up). Is the diagram wrong, or am I wrong to add the undulation height along ellipsoid normal?

Here is the image I am referring to, and as you can see, the refer to the geoid undulation offset as "N" in the diagram:

enter image description here
Source: http://kartoweb.itc.nl/geometrics/Reference%20surfaces/refsurf.html

Thanks for any help in advance!

Best Answer

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.

Figure

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.