[Proj] Re: Locate a point from distance and backwards azimuth?

[strebe]

I wote:

> yielding two solutions, depending on the sign chosen for the radical in B. Both are valid solutions.

Which is not correct. Only one of the two solutions will produce a legal value of φ₁, the other being outside the range of ±90°. Hence, choose the one within range. There is probably some analytic way of deciding which it will be, but I haven't tried to work it out.

Regards,
-- daan Strebe


On Oct 20, 2008, at 4:59:34 PM, strebe <strebe at aol.com> wrote:

Mikæl:

The usual spherical formulæ for computing a point at a known distance and azimuth are:

φ₂ = arcsin (sin φ₁ cos c + cos φ₁ sin c cos Az)
λ₂ = λ₁ + atan2 [sin c sin Az, (cos φ₁ cos c - sin φ₁ sin c cos Az)]

where the start point is subscripted with 1, the end point is subscripted with 2, and the arc distance is c. We know φ₂, λ₂, Az, and c. Therefore these reduce to

(1) n₄ = n₁ sin φ₁ + n₂ cos φ₁
(2) λ₂ = λ₁ + atan2 [n₃, (n₁ cos φ₁ - n₂ sin φ₁)]

where we have made substitutions for symbolic clarity:

n₁ = cos c
n₂ = sin c cos Az
n₃ = sin c sin Az
n₄ = sin φ₂

Since equation 1 has only φ₁ unknown, we solve with trigonometric and algebraic manipulations:

tan φ₁ = A / B

with φ₁ resolvable using the usual atan2 function in order to preserve quadrants. Here

    n₄ - n₂ B
A = _________
        n₁

    n₂ n₄ ± √(n₂² n₁² + n₁⁴ - n₁² n₄²)
B = ----------------------------------
                n₁² + n₂²


yielding two solutions, depending on the sign chosen for the radical in B. Both are valid solutions.

Thus armed with φ₁, we can easily compute λ₁ as

λ₁ = λ₂ - atan2 [n₃, (n₁ cos φ₁ - n₂ sin φ₁)]

The ellipsoidal case is much harder. Your iterative process will work, of course, though it will be terribly inefficient and likely comes with problematic regions. In principle it seems like the problem you are trying to solve is intermediate between the forward and the inverse problems. The reason it is not as simple as the forward problem is that, in the forward case, you know both the starting point and the azimuth. This wholly characterizes which geodesic you are dealing with, since a geodesic is completely characterized by its constant and any point along it:

      cos φ sin Az
K = ________________
    √(1 - e² sin² φ)

In your case, you cannot compute the geodesic's constant because you don't know both the latitude and azimuth from any known point; instead you know the azimuth from an unknown point. On the other hand, you already know the distance, the calculation of which is what makes the inverse problem so difficult.

Warning: Do NOT use Pearson's analysis of geodesics (Frederick Pearson II 1984: "Map Projection Methods", pp. 81–82, and 1990: "Map Projections: Theory and Applications", pp. 80–81). He makes an error early on in the derivation, making the final, easily computed integral nothing but a fiction. Sadly, Donald Fenna has repeated this error in his 2007 volume, "A Compendium of Map Projections, with Derivations", pp. 394–395. 

Regards,
-- daan Strebe



-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20081020/2d16bc1e/attachment.html