4

I'm attempting to implement Vincenty's inverse problem as described on wiki HERE

The problem is that lambda is simply not converging. The value stays the same if I try to iterate over the sequence of formulas, and I'm really not sure why. Perhaps I've just stared myself blind on an obvious problem.

It should be noted that I'm new to Python and still learning the language, so I'm not sure if it's misuse of the language that might cause the problem, or if I do have some mistakes in some of the calculations that I perform. I just can't seem to find any mistakes in the formulas.

Basically, I've written in the code in as close of a format as I could to the wiki article, and the result is this:

import math

# Length of radius at equator of the ellipsoid
a = 6378137.0

# Flattening of the ellipsoid
f = 1/298.257223563

# Length of radius at the poles of the ellipsoid
b = (1 - f) * a

# Latitude points
la1, la2 = 10, 60

# Longitude points
lo1, lo2 = 5, 150

# For the inverse problem, we calculate U1, U2 and L.
# We set the initial value of lamb = L
u1 = math.atan( (1 - f) * math.tan(la1) )
u2 = math.atan( (1 - f) * math.tan(la2) )
L = (lo2 - lo1) * 0.0174532925

lamb = L

while True:
    sinArc = math.sqrt( math.pow(math.cos(u2) * math.sin(lamb),2) + math.pow(math.cos(u1) * math.sin(u2) - math.sin(u1) * math.cos(u2) * math.cos(lamb),2) )
    cosArc = math.sin(u1) * math.sin(u2) + math.cos(u1) * math.cos(u2) * math.cos(lamb)
    arc = math.atan2(sinArc, cosArc)
    sinAzimuth = ( math.cos(u1) * math.cos(u2) * math.sin(lamb) ) // ( sinArc )
    cosAzimuthSqr = 1 - math.pow(sinAzimuth, 2)
    cosProduct = cosArc - ((2 * math.sin(u1) * math.sin(u2) ) // (cosAzimuthSqr))
    C = (f//16) * cosAzimuthSqr  * (4 + f * (4 - 3 * cosAzimuthSqr))
    lamb = L + (1 - C) * f * sinAzimuth * ( arc + C * sinArc * ( cosProduct + C * cosArc * (-1 + 2 * math.pow(cosProduct, 2))))
    print(lamb)

As mentioned the problem is that the value "lamb" (lambda) will not become smaller. I've even tried to compare my code to other implementations, but they looked just about the same.

What am I doing wrong here? :-)

Thank you all!

CodingBeagle
  • 1,888
  • 2
  • 24
  • 52
  • 1
    In your example, which value should `lamb` converge to? – Jivan Dec 31 '14 at 03:03
  • Try `from __future__ import division` – Jonathan Davies Dec 31 '14 at 03:10
  • First off, perhaps you should get rid of integer divisions (`//`) and replace them by float divisions (`/`)? – Jivan Dec 31 '14 at 03:10
  • 1
    Then, I think there's a misunderstanding in the mathematical writing you read on the page you linked to. When they say `sin α = blabla`, it doesn't mean `sinAlpha = blabla`, but it means `α = arcsin blablabla` – Jivan Dec 31 '14 at 03:12
  • Ah, sorry, I did not put that change in the code :-) lamb should converge towards the mentioned value in the wiki article, that is: "When λ has converged to the desired degree of accuracy (10^(−12) corresponds to approximately 0.06mm)". So, I'm trying to get it to converge towards that value – CodingBeagle Dec 31 '14 at 03:13
  • In that case, I've clearly misunderstood the mathematical writing! – CodingBeagle Dec 31 '14 at 03:14
  • Well, wait. What I've said is true, but in fact it doesn't matter since you don't want `α` directly. So it seems ok for this part. – Jivan Dec 31 '14 at 03:15
  • 1
    @Bitious in `L = (lo2 - lo1) * 0.0174532925`, why `* 0.017...`? The wiki article doesn't mention this. – Jivan Dec 31 '14 at 03:33
  • First of all, thank you for the help so far jivan! :) It is wonderful! I multiplied by 0.017... in order to convert the lo2, lo1 values from degrees to radian, however I'm not sure if this is needed. It was based on another example I found on the internet. – CodingBeagle Dec 31 '14 at 03:46
  • @Bitious it's ok, I made an explanation in the comments of my answer - I refer you to that – Jivan Dec 31 '14 at 03:59
  • When calculating `u1` and `u2`, shouldn't you be converting `la1` and `la2` to radians before calling `math.tan()`? – Jim Lewis Dec 31 '14 at 04:00
  • @JimLewis this is the goal of `L = (lo2 - lo1) * 0.0174532925` - `0.017...` equals `math.pi / 180` – Jivan Dec 31 '14 at 04:06
  • 2
    Yes, but I'm talking about la1 and la2 in the immediately preceding lines. – Jim Lewis Dec 31 '14 at 04:08
  • @JimLewis - you're absolutely right – Jivan Dec 31 '14 at 04:12
  • I think you might be correct @JimLewis, they should be converted to radians the same way I believe. – CodingBeagle Dec 31 '14 at 04:12

2 Answers2

4

First, you should convert you latitudes in radians too (you already do this for your longitudes):

u1 = math.atan( (1 - f) * math.tan(math.radians(la1)) )
u2 = math.atan( (1 - f) * math.tan(math.radians(la2)) )
L = math.radians((lo2 - lo1)) # better than * 0.0174532925

Once you do this and get rid of // (int divisions) and replace them by / (float divisions), lambda stops repeating the same value through your iterations and starts following this path (based on your example coordinates):

2.5325205864224847
2.5325167509030906
2.532516759118641
2.532516759101044
2.5325167591010813
2.5325167591010813
2.5325167591010813

As you seem to expect a convergence precision of 10^(−12), it seems to make the point.

You can now exit the loop (lambda having converged) and keep going until you compute the desired geodesic distance s.

Note: you can test your final value s here.

Jivan
  • 21,522
  • 15
  • 80
  • 131
  • 1
    yes it's actually correct - except you're not done yet. Lambda is not the distance between two points but an angle. This result (2.5343...) seems correct. You now need to go on (out of the loop) with the second part until you compute `s` which is the geodesic distance, which is what you want in the end. – Jivan Dec 31 '14 at 03:58
  • Thank you Jivan! :-) Yes, I'm aware that I still need to implement the second part, however I wanted to make sure the first step of the process was correct first :) Thank you a lot for the help! – CodingBeagle Dec 31 '14 at 04:00
  • 2
    @Bitious: FWIW, the `math` module provides `degrees()` and `radians()` functions to convert radians to degrees and vice versa. Using these functions instead of multiplying or dividing by pi/180 or obscure [magic numbers](http://en.wikipedia.org/wiki/Magic_number_%28programming%29) like 0.0174532925 makes code more readable. – PM 2Ring Dec 31 '14 at 05:29
  • Cheers @PM2Ring ! :) I'll use these functions instead! – CodingBeagle Dec 31 '14 at 12:51
3

Even if it is correctly implemented, Vincenty's algorithm will fail to converge for some points. (This problem was noted by Vincenty.) I give an algorithm which is guaranteed to converge in Algorithms for geodesics; there's a python implementation available here. Finally, you can find more information on the problem at the Wikipedia page, Geodesics on an ellipsoid. (The talk page has examples of pairs of points for which Vincenty, as implemented by the NGS, fails to converge.)

cffk
  • 1,839
  • 1
  • 12
  • 17
  • I've submitted a [patch](https://github.com/geopy/geopy/pull/144) to [geopy](https://pypi.python.org/pypi/geopy) to use the accurate geographiclib algorithms if they are available. – cffk Aug 20 '15 at 19:56