geocoding - Trilateration using 3 latitude and longitude points, and 3 distances

Question

There exists an unknown target location (latitude and longitude co-ordinates). I have 3 latitude and longitude co-ordinate pairs and for each pair a distance in kilometers to the target location. How can I calculate the co-ordinates of the target location?

For example, say I have the following data points

37.418436,-121.963477   0.265710701754km
37.417243,-121.961889   0.234592423446km
37.418692,-121.960194   0.0548954278262km

What I'd like is what would the guts of the function that takes that as input and returns 37.417959,-121.961954 as output look like?

I understand how to calculate the distance between two points, from http://www.movable-type.co.uk/scripts/latlong.html I understand the general principle that with three circles you get exactly one point of overlap. What I'm hazy on is the math needed to calculate that point with this input.

Answer 1

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).