I am in search of an algorithm to solve the first/direct geodetic problem precisely in 3D. The input is:

  1. a geodetic (not geocentric) position (WGS84, height is above the WGS84 reference ellipsoid)
  2. heading (relative to true north)
  3. pitch and roll (relative to gravity, which isn't necessarily pointing towards the center of the earth)
  4. a 3d-offset from the pose P (defined by 1,2,3) to a point A, given in cartesian coordinates (meters).

The output should be the geodetic position (again, WGS84) of point A.

I have read Algorithm for offsetting a latitude/longitude by some amount of meters and this is how I've computed the solution in the past. It is obviously computationally efficient, but not precise enough anymore.

Currently, my idea is to

  1. convert 1. to ECEF, populating the 4th column of a 4x4 matrix with the solution.
  2. incorporate the angles/orientation into that matrix
  3. multiply the matrix with the 3d vector to A to compute the position of A in ECEF.
  4. convert this back to Geodetic WGS84 (which isn't quite straightforward).

Of course, I still would like things to be fast, but precision has a definite priority here.

My questions:

  • What do you think of this approach? Does it work, is it precise? (I fear that compensating for pitch/roll not being relative to the center of the earth would probably imply a LOT of additional complexity).

  • Is PROJ.4's pj_transform() capable of computing the conversions in steps 1 and 4? If not, what C/C++ library can you recommend? This must be a common routine in some library?!

Note to @whuber who referred to 30448! I think I'm slowly starting to resolve misunderstandings, getting an idea of what you do. Questions about your answer to that question:

  • When you say "geocentric", you mean cartesian coordinates from earth center (=ECEF), right? (I first interpreted "geocentric coordinates" just like "geodetic coordinates", except the latitude being determined using earth-center instead of ellipsoid-normal. Now I know that "geocentric latitude" != "geocentic coordinates". Whew.)

  • In the first image of 30448, you show "the geocentric frame at the North Pole". The origin of that coordinate system really is at the center of the earth, right? But if your first image started at earth's center, simply rotating wouldn't get you to the origin of the local coordinate system. So why does it start at (0,0,6370) in the first picture? You could have also chosen (6370,0,0), but then rotations would be different. Is it that each initial offset corresponds to an Euler convention (Z-Y'-Z'')?

  • If so, it then seems to me that my drafted approach is very much similar to your solution at 30448. Are there other approaches that are more precise? Given I compute using doubles, what accuracy do you expect?

  • Slightly OT: At the end, you add a local Z displacement motivated by differences against MSL - but I thought the WGS84 datum you referenced earlier is not a tidal/MSL datum, but an ellipsoidal. Assuming you sticked to WGS84, you'd add the local Z offset based on the distance to the ellipsoid, right?

@whuber, @martin f provided quick replies and patience - I read the whole day, trying to understand more.

  • If this question is not a duplicate of gis.stackexchange.com/questions/30448/…, please help us understand how it differs.
    – whuber
    Dec 29, 2013 at 18:41
  • @whuber -- I think the Q&A you cite is geo-centric, while this is geodetic. (Maybe that simplifies things.) Also -- now adding complexity -- here the local cartesian system is tilted about three geodeic axese: vertical, E and N. Is that right @benadler? Further still, "vertical" is gravitational which is different again from the ellipsoidal normal!
    – Martin F
    Dec 29, 2013 at 19:53
  • @martin f My answer there shows how to convert back and forth between geodetic and geocentric systems as well as local and global coordinate frames. Unless you are correcting for microgravity--which requires a detailed gravitational model and no longer involves mere geodetic calculations--then the ellipsoidal normal is the same as the direction of gravitation.
    – whuber
    Dec 29, 2013 at 20:57
  • On the vertical deflection -- "the angle between the true zenith (plumb line) and the line perpendicular to the surface of the reference ellipsoid ... can amount to angles of 10″ (flat areas) or 20–50″ (alpine terrain)" -- see en.wikipedia.org/wiki/Vertical_deflection
    – Martin F
    Dec 30, 2013 at 3:09
  • @whuber -- Does your answer account for non-zero roll, pitch & yaw w.r.t local geodetic N, E & Z?
    – Martin F
    Dec 30, 2013 at 3:15

1 Answer 1


Some potentially good resources:

Not the answer you're looking for? Browse other questions tagged or ask your own question.