convert - Converting (geodetic pose + cartesian offset) to geodetic position?

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.