3

I am a programmer. I have a point P1 (lat,lon) given by GPS (WGS84) and I want to know the location of a point P2 if I make a straight line from P1 to the other side of the earth.

How can I achieve this if it is possible ?

2
  • 2
    Are you looking for the antipodal point (which cuts through the earth)?
    – Paul
    Jul 15, 2013 at 17:14
  • Sergio, I have interpreted your "straight line" as a vertical line at its point of origin, and changed the title accordingly to draw more interest from readers :-). If this is incorrect, please modify the title and edit your question to indicate how the direction of your line actually is determined.
    – whuber
    Jul 15, 2013 at 19:28

1 Answer 1

7

Because the WGS model of the earth is a surface of rotation (about the z axis), the longitude of P2 is the longitude of P1 plus or minus (it doesn't matter) 180 degrees. Add or subtract as you wish in order to obtain a longitude within your desired range (typically from -180 to +180 or from 0 to +360).

The latitude describes the generator of the surface, an ellipse of major radius a (in the Equatorial plane) and minor radius b (along the axis of revolution). For WGS 84, a = 6,378,137.0 meters and b = 6,356,752.314245 meters (to the nearest micron :-). Specifically, the latitude is the angle made between the local "up" direction and the Equatorial plane.

A straightforward way to solve problems of this type is to switch among two representations of this ellipse. Both of them are obtained by distorting the unit circle via a stretch by a in the x direction and a stretch by b in the y direction (switching to Cartesian coordinates in which x is horizontal--along the Equatorial plane--and y is vertical). The usual parametrization of the circle, (x, y) = (cos(q), sin(q)), becomes (x, y) = (a cos(q), b sin(q)). We have to work out the local "up" vector. It is obtained by rotating the tangent vector clockwise by 90 degrees. Using Calculus to differentiate with respect to q gives the tangent vector and rotating it yields the "up" vector as (b cos(q), a sin(q)). The angle it makes has a tangent equal to "rise" over "run," or y/x = a/b * sin(q)/cos(q) = a/b * tan(q). By definition, this is the tangent of the latitude f. Solving gives

q = ArcTan(b/a * tan(f)). 

From this solution we can compute (x, y) from the latitude using the parameterization

(x, y) = (a cos(q), b sin(q)).

Figure

In this figure (with an exaggerated value of b/a = 4/5), nine points with equally spaced latitudes from -90 through +70 degrees are projected opposite their local "up" directions until they intercept the ellipse again. The dashed arrows depict the projections. Note that although the latitudes are equally spaced, the points themselves are not equally spaced. Note, too, that the projected points are not antipodal except for the projection of the South Pole (latitude = -90 degrees, at the bottom).

Next, starting at any point on the ellipse, we travel some unknown amount opposite the "up" direction until we once again intercept the ellipse, as shown in the figure. For this, we view the ellipse as a distortion of the unit circle given as the level set x^2 + y^2 = 1, whence it must have the formula (x/a)^2 + (y/b)^2 = 1. Let the unknown amount of motion be t. This gives the equations

(x, y) = (a cos(q) + b sin(q)) - t(b cos(q), a sin(q))
(x/a)^2 + (y/b)^2 = 1

There are two solutions for t, of which one (of course) is t = 0, because we began at a point on the ellipse. Because this solution is so easily found, the other solution is obtained from a linear (not quadratic) equation and so has the relatively simple expression

t = 2 a b ((a sin(q))^2 + (b cos(q))^2) / ((a^2 sin(q))^2 + (b^2 cos(q))^2)

Plugging that into the first of the preceding equations gives (x, y). Converting those back to a latitude finishes the work.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.