How to calculate the 3D straight line distance between 2 GPS measurements?

[aydos]

This seems like a trivial question at first, but I am struggling how to get it right. If I have 2 GPS (lat, lon, height) observations P1 and P2 taken at some height from the Earth's surface, how do I calculate the straight line distance between them? I am using the picture below to illustrate the question and get into more detail.
I understand GPS devices provide heights in reference to the WGS84 ellipsoid and perpendicular to the ellipsoid's surface (lines labelled "gps h"). I also understand latitude and longitude are angles originating at the WGS84 ellipsoid centre.
Are my assumptions ok? Can anyone out there can explain how to calculate the distance between P1 and P2?

 

[anorlunda]

This sounds like homework, is it?

 

[aydos]

No.

 

[phinds]

Possibly you are misleading yourself by the gross over-simplification of your diagram which just shows an ellipse and pretends that it is an ellipsoid in a way that doesn't really reflect what is happening.

 

[aydos]

Possibly you are misleading yourself by the gross over-simplification of your diagram which just shows an ellipse and pretends that it is an ellipsoid in a way that doesn't really reflect what is happening.

I know my diagram is simple, but the question still valid, no? Would you be able to help?

 

[aydos]

A slightly more concrete example:
A GPS mounted on a climbing aircraft produces observations P1 and P2. What is the shortest distance between P1 and P2 in meters?

 

[phinds]

A slightly more concrete example:
A GPS mounted on a climbing aircraft produces observations P1 and P2. What is the shortest distance between P1 and P2 in meters?

This does not change the situation at all; your question is exactly the same. I still think you are not seeing the solution because you have not drawn the right diagram (and perhaps you don't understand spherical trigonometry). I think the problem will require a complex calculation process, but it is conceptually trivial.

 

[berkeman]

spherical trigonometry

Yeah, I was confused by the 2-D drawing for a 3-D problem.

@aydos -- shouldn't the "heights" in your diagram be in the same plane as the ellipsoid that you've drawn? If planes climb to altitude above the surface of the Earth, they don't move parallel to the ground, which is what your drawing seems to show.

 

[Baluncore]

Are my assumptions ok? Can anyone out there can explain how to calculate the distance between P1 and P2?

First you must convert the Latitude, longitude and height to an Earth centred x, y, z.
Do that for the two points using the following code to get x1, y1, z1 and x2, y2 and z2.
Take the differences to get dx = x2 - x1, dy = y2 - y1 and dz = z2 - z1.
Sum the squares of the differences, then square root to get your straight line distance.

This is the WGS84 code to convert GPS lat, lng and height, to x, y, z.
Note: that the Latitude and Longitude may need to be in radians. With the heights in metres. Radians = degrees * Pi / 180.

Double lat, lng, h    ' GPS inputs WGS84
Double x, y, z          ' outputs in metres

Const Double a = 6378137             ' Defined semi-major axis, metres
Const Double bars = .9933056200098587  ' square of b to a ratio   bb/aa
Const Double ee = 6.694379990141316D-03 'square of first eccentricity f(2-f)

Double SinLat = Sin(lat) ' precompute
Double RofC = a / Sqrt(1 - ee * SinLat * SinLat)  ' radius of curvature
Double term = (RofC + h) * Cos(lat)
x = term * Cos(lng)    ' x-axis is Atlantic Ocean positive
y = term * Sin(lng)    ' y-axis is Indian Ocean positive
z = (h + (RofC * bars)) * SinLat  ' z-axis is North Pole positive

 

[aydos]

Hi @Baluncore , thanks so much for your reply. That is exactly what I was looking for. One doubt remains for me. Its about the (RofC + h) term at line 10.
The actual height reported by the GPS is not parallel to the ellipsoid radius, rather perpendicular to the ellipsoid surface. Am I correct to assume that the blue term above does not take that into account?

 

[Baluncore]

Your application is based on GPS from a moving vehicle, so the algorithm does not need to be particularly accurate. I believe that for an aircraft covering only a few kilometre, you could find the distance traveled to better than 100 mm by assuming a spherical Earth with a fixed radius.

The source of my code is lost in a 30 year old paper filing system. I believe that it is better than 1 mm when using double precision. The radius of the Earth is large compared with the height of an aircraft so it should make less than 1 mm difference either way. Since you are subtracting two nearby conversions, any approximation error should largely cancel.

If you want to check the algorithm, here is an accurate inverse, x,y,z to lat,long,height.
https://hal.archives-ouvertes.fr/hal-01704943/document

 

[Baluncore]

Here is an explanation of the algorithm.
http://blog.sciencenet.cn/home.php?mod=attachment&id=20102

 

[aydos]

Great references, thank you that all I needed!

 

[Baluncore]

I find now there is a GPS toolbox assoociated with the Journal “GPS Solutions”.
https://www.ngs.noaa.gov/gps-toolbox/index.html

There is a good explanation of ellipse geometry on pages 467 – 470 of;
Title: Linear Algebra, Geodesy, and GPS. Author: Gilbert Strang And Kai Borre.
ISBN 0-9614088-6-3 (Hardcover). 1997.
The Matlab M-files referred to in that text seem to have changed names and migrated to;
http://kom.aau.dk/~borre/geodesy/
http://kom.aau.dk/~borre/life-l99/
See; togeod.m and topocent.m
Found via; https://www.ngs.noaa.gov/gps-toolbox/Borre.htm