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.

Breadwinner

Seems actually to be "trilateration" en.wikipedia.org/wiki/Trilateration as "triangulation" involves angles – Breadwinner

Solution: gis.stackexchange.com/questions/66/… – Agnesagnese

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

Brecher

While doing a little more poking around this morning, I realized that the conventional way to convert lat/lon to Cartesian coordinates (called ECEF), is to make X a function of cos(lon) and Y a function of sin(lon). I had these reversed. It makes no difference to the final calculation--it just means I was using an unconventional coordinate system for the intermediate step. But I have now corrected the equations to the more conventional form. – Brecher

Incidentally, for the greatest precision, you will need to account for altitude, and for the fact that the lat/lon coordinates are based on an ellipsoidal rather than a spherical reference surface. See here: en.wikipedia.org/wiki/ECEF and here: satsleuth.com/GPS_ECEF_Datum_transformation.htm – Brecher

Looks like I reversed the cos(lat) and sin(lat) as well. It was late. Sorry. Now fixed. (I actually did the calculation right, but screwed up transcribing the equations) – Brecher

Thank you, that is like 90% of what I was looking for, very helpful, and you deserve the bounty. I should be able to get what I need based on this answer. – Breadwinner

This post helped me a lot. Thank you very much. Just a small concern. The point I'm receiving, they are in Cartesian Coordinates. I was just wondering how should I convert them to a local coordinate systems to place it in a map? Is there any formula or so? – Spur

Can you provide us a book where I can find what you did? thank you – Catty

Answer 2

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

Breadwinner

Answer 3

On the Paul Bourke Geometry pages

intersection of two circles

Sabin

This is a broken link. – Agnesagnese

Answer 4

Consider the following 9 circles Points A,B,C and distances d1, d2, d3

• Center of A, radius d1
• Center of A, radius d2
• Center of A, radius d3
• Center of B, radius d1
• Center of B, radius d2
• Center of B, radius d3
• Center of C, radius d1
• Center of C, radius d2
• Center of C, radius d3

These are your possible circles. Now we can cull these, because we know if d1 is used on A, it won't be used on B.

This makes your possible entries, where A1 means circle with center A and radius D1:

• {A1, B2, C3}
• {A1, B3, C2}
• {A2, B1, C3}
• {A2, B3, C1}
• {A3, B1, C2}
• {A3, B2, C1}

You should be able to convert the lat/long to X,Y,Z knowing the radius of the earth, and the distances from the curved distance along the earths crust to the straight distance, and from there you can see which of them intersect at a common point. Remember to allow for small margins of error due to float imperfection.

Slaughter

I already know which radius goes with which center; what I need are the formulae for calculating the intersection point given lat/lon data as inputs – Breadwinner