Trilateration using 3 latitude and longitude points, and 3 distances

Wikipedia gives a pretty thorough discussion of the algebra here: http://en.wikipedia.org/wiki/Trilateration

The first step, not really covered in the Wikipedia entry, is to convert your lat/long coordinates to Cartesian coordinates:

x0 = cos( lon0 ) * cos( lat0 ) , y0 = sin( lon0 ) * cos( lat0 ) , z0 = sin( lat0 )
x1 = cos( lon1 ) * cos( lat0 ) , y1 = sin( lon1 ) * cos( lat1 ) , z1 = sin( lat1 )
x2 = cos( lon2 ) * cos( lat0 ) , y2 = sin( lon2 ) * cos( lat2 ) , z2 = sin( lat2 )

(To keep calculations simple, I've fudged things so we are working in units of "earth radii" instead of kilometers)

For your data, I get

         p0            p1           p2
X   -0.420442596  -0.420430618  -0.42040255
Y   -0.67380418   -0.673826567  -0.673825967
Z    0.607631426   0.607614889   0.607634975

The next step, which is covered in the Wikipedia article, is to simplify the coordinates, by translating the points so p0 is at the origin, and then rotating so that p1 is on the X axis, and p2 is in the X-Y plane.

For the translation, just subtract p0 from p1 and p2:

    p0a      p1a          p2a
X   0    1.19779E-05   4.00462E-05
Y   0   -2.23864E-05  -2.17865E-05
Z   0   -1.65372E-05   3.5486E-06

The rotation isn't much harder. p1b gets (x,y) = (d,0), where d is just the distance from the origin to p1a (Pythagorean theorem)

For p2b, we need to resolve p2a into two components: one parallel to p1a (which goes on our x axis), and one perpendicular to p1a, (which goes on our y axis in the "b" coordinate system).

To do this, we need a unit vector in the direction of p1a, which is just p1a * ( 1/d ). Take the dot product of this unit vector (call it p1a_hat, if you like) with p2a, and that's the X coordinate for p2b. The Wikipedia article calls this value "I"

Now the Y coordinate is easy. The length from the origin to p2 can't change under the coordinate transformation. So calculate p2a's length using the Pythagorean theorem, then use the Pythagorean theorem "backwards" to get what the Y coordinate for p2b has to be to keep the length the same. That's the variable that Wikipedia calls "J". (Note, there's an ambiguity that I'll leave for you to figure out over whether J is positive or negative).

Now you've got the three variables d, I and J, that the Wikipedia article uses for the calculation. You can convert them back to kilometers now, by multiplying by the earth's radius. You should be able to do the rest of the calculation from here

(Incidentally, Wikipedia gives a different calculation for the coordinate transformation. I like to avoid trig where possible).


I asked this question on the newly-formed GIS Stack Exchange, and got some good answers there as well.

https://gis.stackexchange.com/questions/66/trilateration-using-3-latitude-and-longitude-points-and-3-distances

The accepted answer there has a (presumably) working solution in Python:

https://gis.stackexchange.com/questions/66/trilateration-using-3-latitude-and-longitude-points-and-3-distances/415#415


On the Paul Bourke Geometry pages

intersection of two circles