# Geodetic lat/lon to ECEF cartesian: Direction to apply geoid undulation from ellipsoid surface?

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: Source: http://kartoweb.itc.nl/geometrics/Reference%20surfaces/refsurf.html

Thanks for any help in advance!

• Check out Vertical Deflection but my gut feeling is that I think it's generally going to be negligible versus the SRTM accuracy. – mkennedy Aug 29 '14 at 20:26
• Ok, I think the part that is still confusing me is just what the direction vector of N actually represents. It doesnt seem to be along the ellipsoid normal (as h is) or along the geoid normal (as H is), but some derived direction? I understand the difference is minuscule from the ellipsoid normal, but still want to know what exactly it is... – Grumple Sep 3 '14 at 20:35

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.

• Ok, thanks very much for these explanations. It does help, but I still find myself confused on one key point. Let's say I do want to be as accurate as possible, and apply H and N along their most correct directions instead of generalizing. Based on what you have described, it sounds technically wrong to apply the undulation N along the ellipsoidal normal? I can see how generalizing to just using h = N + H is applied along ellipsoid normal, but I don't see how I find the direction vector to use for N? It seems like it depends on knowing final ECEF position of P1 which I am trying to find? – Grumple Sep 3 '14 at 20:25
• That's right: it's not an easy problem. But why tackle it unless you need sub-millimeter accuracy everywhere on earth? There are two parts that are technically suspect: one is to assume that h and N are parallel (they are not when E1 and E1' differ) and the other is to assume that H and h are parallel (they are not unless the geoid's normal vectors at G1 and G1' are parallel). It is clear, though, that a single iteration will converge to an extremely precise solution (compute P1 under these assumptions, then calculate H, N, and E1 from it and use those to correct the initial estimates). – whuber Sep 3 '14 at 20:38
• Ok thanks so much, I finally feel like I at least understand the nature of that problem. Based on this I am actually doing it wrong (even if only erring by millimeters) in my original description. So, to be 100% clear, the correct generalization would be to apply H+N along the geodetic ellipsoid normal from the surface of the ellipsoid? – Grumple Sep 3 '14 at 20:48
• Given a (lat, lon) for E1', you can determine the direction of the vector N from the ellipsoid's specifications. The geoid's height then gives you the point G1. You would need to numerically differentiate the geoid model at E1' to determine the direction vector H. The orthometric height tells you the displacement from the point G1 along H needed to arrive at P1. This will be slightly incorrect to the extent E1' and E1 differ, but should be good to 1 mm. A cruder approximation is to move from location E1 by the amount H+N along the ellipsoidal normal at E1. – whuber Sep 3 '14 at 20:57
• This makes sense when trying to stick to being as accurate as possible. What I'm trying to verify just to be 100% clear, is what happens if I am satisfied with using the simpler generalization. So, for example: If I assume the case of h = H + N, I am assuming that I can just find the ellipsoid normal at lat/lon E1', and project from the surface of the ellipsoid at E1', along that ellipsoid normal by a distance of H + N, to arrive at a very good approximation of P1? I think I'm stating the obvious at this point but I want to be sure I fully understand once and for all. =P – Grumple Sep 4 '14 at 16:19

That diagram is greatly exaggerating the real but very small difference between two vertical directions: the normal to the ellipsoid, and the direction of gravity (normal to the geoid). The difference is called the deflection of the vertical. Unless you're doing very precise geodetic calculations, you can ignore it. In other words, it's safe to assume `h = N + H`. (See New coordinate system? for reference to many alternative diagrams.)

• Just noticed, @mkennedy was correct. – Martin F Aug 30 '14 at 18:20