Can I determine unique pair of (lat, lon) by using position angle and distance? (spherical geometry question.)

Question

If we are on spherically curved surface, (I will put the radius as 1 to make every unit to be treated as radian.) we can select two random coordinates. Let's call those two pairs as (lon, lat) and (lon', lat'). These coordinates are based on spherical coordinate system, however, we do not take care about the radius.

we have following formula for measuring distance (D) and position angle (PA) between two coordinates.

D = arc-cos(sin(lat) sin(lat') + cos(lat) cos(lat') cos(lon'-lon))

PA = arc-tan(sin(lon'-lon)/[cos(lat) tan(lat)' - sin(lat) cos(lon'-lon)])

What I want to do is doing the inverse. Let's suppose we have a set of (lon, lat) and (D,PA). I expected we can get an unique pair of (lon',lat') if I use inverse backtracking method (fsolve in python scipy module).

However, it seems that (D,PA cannot determine unique pair of (lon',lat') pair.

For example, If we have lon = -0.578926 (rad) lat = 0.837584 (rad), D = 0.523258 (rad), PA = -0.088221 (rad). We have at least two (lon',lat') pairs.

(-0.785721,1.354690) and (-0.532593,0.315694)

What I expected at the first time is that there would be only one certain point which would have certain distance and certain angle from the reference (in case of position angle, it would be the north pole.) This makes sense when I am using the Euclidean geometry and I was expecting this would be true for the spherical geometry. However, the calculation method shows that this might not be true. Is there any misunderstanding I am making? Then, what should I use to determine unique set of (lon',lat') other than pair of (D,PA)?

Answer 1

For reference, here is what your three points would look like if we plotted them on a map of the Earth:

The red "pin" near the center of the map (at about the same latitude as France) is at your original lon and lat and the dot labeled "Avannaata" near the top of the map (in Greenland) is at your first pair of lon' and lat' values. The dot near the bottom of the map (in the Atlantic Ocean at about the latitude of the African country Mauritania) is at your second pair of lon' and lat' values.

I reproduced your $D$ and $PA$ functions in python:

from math import sin, cos, tan, acos, atan

def D(lat1, lon1, lat2, lon2):
    d = acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon2-lon1))
    return d

def PA(lat1, lon1, lat2, lon2):
    pa = atan(sin(lon2-lon1)/(cos(lat1)*tan(lat2) - sin(lat1)*cos(lon2-lon1)))
    return pa

Here are the results of calculating $D$ and $PA$ from your original lon and lat to the first pair of lon' and lat' values:

D(lat1=0.837584, lon1=-0.578926, lat2=1.35469, lon2=-0.785721)
0.5232576219102509
PA(lat1=0.837584, lon1=-0.578926, lat2=1.35469, lon2=-0.785721)
-0.08822105650081137

The $D$ and $PA$ values agree with your $D$ and $PA$ values. The distance is hard to judge from the map due to the fact that this is a Mercator projection (which distorts distances, especially in places that are near the poles, like Greenland), but the direction is clearly northward. And indeed $PA$ returned a value in radians that converts to about $-5$ degrees, more conventionally stated as a direction of $355$ degrees (the equivalent angle in the range from $0$ to $360$) which is just a little bit west of due north. So that appears to be correct.

Now let's try the formulas with the second set of lon' and lat'.

D(lat1=0.837584, lon1=-0.578926, lat2=0.315694, lon2=-0.532593)
0.5232579146536581
PA(lat1=0.837584, lon1=-0.578926, lat2=0.315694, lon2=-0.532593)
-0.08822145475051882

We get the same results. But now the direction doesn't make any sense: the destination is the southernmost of the three sets of coordinates we plotted, almost directly south of the red pin, but the $PA$ formula is telling us that to get from the red pin to the southernmost dot we should start out in the direction $-5$ degrees (or $355$ degrees), which is almost due north. In short, $PA$ is sending us in the wrong direction entirely!

This is why I say this formula for $PA$ is wrong. It gives correct answers whenever you ask it the direction to somewhere north of you, but if you ask for the direction to someplace south of you, it sends you north instead.

---

Actually, if you start at a particular point on a sphere, point yourself in a particular direction, and travel "straight" (that is, along the great-circle path in that direction) for a particular distance, you will arrive at a particular point. There is only one possible destination of that journey.

In particular if we start at the red pin and go in the direction $-0.088221$ radians for a distance $0.523258,$ we will arrive at the point labeled "Avennaata," not the point near the bottom of the map, west of Mauritania.

So what went wrong when you tried to find that point?

You literally asked fsolve to find you any pair of coordinates lon' and lat' for which $D$ and $PA$ would give you the distance $0.523258$ and direction $-0.088221$ from the lon and lat coordinates to lon' and lat'. And we can see from the results of calculating $D$ and $PA$ above, $D = 0.523258$ and $PA = -0.088221$ for both sets of lon' and lat' coordinates. Therefore both of those points are solutions to the problem you asked fsolve to solve for you. But one of those solutions is a "solution" only because the problem you asked fsolve to solve includes a broken $PA$ function that says you are going north when you're actually going south. It gave you the correct answer, where you go in a northward direction, and also the answer where you go in a southward direction which $PA$ incorrectly says is northward.

The underlying problem, by the way, is really that atan is only capable of returning angles between $-\frac\pi2$ and $\frac\pi2$ radians, which (when you convert them to directions on the Earth) range from west to east in various northward directions, never southward directions.

You can "fix" this by adding extra code to pay attention to whether the denominator in the $PA$ function, cos(lat1)*tan(lat2) - sin(lat1)*cos(lon2-lon1), is positive, zero, or negative, but it's annoying to code and to test (especially the zero case, because the value of lon2-lon1 might be either positive or negative in an eastward direction). My advice is that the much better way to fix this is to use a function that is already designed to give you the correct answer no matter which direction your destination is in, and that function is atan2.

The implementation of $PA$ using atan2 might look like this:

from math import sin, cos, atan2

def PA2(lat1, lon1, lat2, lon2):
    pa = atan2(sin(lon2-lon1)*cos(lat2), cos(lat1)*sin(lat2) - sin(lat1)*cos(lat2)*cos(lon2-lon1))
    return pa

Basically, to convert an atan formula to an atan2 formula, just change the name of the function and put a comma between the numerator and denominator instead of performing the division to get a single input value. In the implementation above I have gone one step further and multiplied both the numerator and denominator by cos(lat2) in order to avoid having tan(lat2) in the denominator, since that function will not let you find the direction to either the north or south pole (because $\tan(\theta)$ is undefined when $\theta = \pm\frac\pi2$).

The results of calling this version of the $PA$ function on your coordinates are as follows:

PA2(lat1=0.837584, lon1=-0.578926, lat2=1.35469, lon2=-0.785721)
-0.08822105650081136
PA2(lat1=0.837584, lon1=-0.578926, lat2=0.315694, lon2=-0.532593)
3.053371198839274

In the first case, where the actual correct direction is northward, we get the northward direction $-0.088221$ (about $-5$ degrees, as before), but in the second case, instead of getting $-0.088221$ again, we get $3.05337,$ which is about $175$ degrees -- going southward in a direction $180$ degrees different from $-5$ degrees. This is the correct answer.

If you let fsolve use PA2 to find the point at a given distance and direction from a given point, it would (presumably) find a unique, correct answer. However, there is a more direct way to get the answer:

\begin{align} lat' &= \arcsin(\sin(lat) \cos(d) + \cos(lat) \sin(d) \cos(pa)) \\ lon' &= lon + \operatorname{atan2}(\sin(pa) \sin(d) \cos(lat), \cos(d) - \sin(lat) \sin(lat')) \end{align}

You may sometimes want to add or subtract $2\pi$ from the longitude calculated by this formula in order to get the result into the usual range of longitudes in radians ($-\pi$ to $\pi$).

Answer 2

On a sphere, Point A to Point B is one direction and distance. But anti-podal points are not expected to have a solution. Now points on a sphere, on the same half of sphere, could be projected to rectangular grid but the projection would be a different system.

Note that ArcTan2 makes a quadrant system.

Then "Forward" or "Direct" generally do not require a quadrant system with directions of 0 to 360 degrees but "reverse" or "inverse" requires a quadrant system.