2662

Calculate distance between two coordinates on a globe

Question:

I get two coordinate pairs in the form 90°0′0″N 0°0′0″E as string and want to calculate the distance between those points on a sphere with radius R=6371km.

I found two formulas on the internet <a href="http://www.movable-type.co.uk/scripts/latlong.html" rel="nofollow">here</a>, the "haversine" and the "spherical law of cosines", but they don't seem to work. For a 90° angle which should return 2*pi*R / 4, the haversine operates correct but the cosines fail and return 0. A different point with more random coordinates returns false values with both algorithms: the haversine is too high and the cosines are too low.

Is my implementation wrong or did I chose an incorrect algorithm?

How should I make these calculations (coordinate pairs to distance on globe surface) instead?

(And yes, I know that I'm not checking for N/S and E/W yet, but the tested coordinates are all in the north-eastern hemisphere.)

Here's my Python 3 code:

<pre class="lang-py prettyprint-override">import math, re R = 6371 PAT = r'(\d+)°(\d+)′(\d+)″([NSEW])' def distance(first, second): def coords_to_rads(s): return [math.radians(int(d) +int(m)/60 +int(s)/3600) \ for d, m, s, nswe in re.findall(PAT, s)] y1, x1 = coords_to_rads(first) y2, x2 = coords_to_rads(second) dx = x1 - x2 dy = y1 - y2 print("coord string:", first, "|", second) print("coord radians:", y1, x1, "|", y2, x2) print("x/y-distances:", dy, dx) a = math.sin(dx/2)**2 + math.cos(x1) * math.cos(x2) * math.sin(dy/2)**2 c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a)) haversine = R * c law_of_cosines = math.acos( math.sin(x1) * math.sin(x2) + \ math.cos(x1) * math.cos(x2) ) * R print("HS:", round(haversine, 2), "LOC:", round(law_of_cosines, 2)) return haversine #return law_of_cosines if __name__ == '__main__': def test(result, correct): print("result: ", result) print("correct:", correct) test(distance("90°0′0″N 0°0′0″E", "0°0′0″N, 0°0′0″E"), 10007.5) test(distance("51°28′48″N 0°0′0″E", "46°12′0″N, 6°9′0″E"), 739.2) test(distance("90°0′0″N 0°0′0″E", "90°0′0″S, 0°0′0″W"), 20015.1) test(distance("33°51′31″S, 151°12′51″E", "40°46′22″N 73°59′3″W"), 15990.2)

Here is some output:

coord string: 90°0′0″N 0°0′0″E | 0°0′0″N, 0°0′0″E coord radians: 1.5707963267948966 0.0 | 0.0 0.0 x/y-distances: 1.5707963267948966 0.0 HS: 10007.54 LOC: 0.0 result: 10007.543398010286 correct: 10007.5 coord string: 51°28′48″N 0°0′0″E | 46°12′0″N, 6°9′0″E coord radians: 0.8984954989266809 0.0 | 0.8063421144213803 0.10733774899765128 x/y-distances: 0.09215338450530064 -0.10733774899765128 HS: 900.57 LOC: 683.85 result: 900.5669567853056 correct: 739.2

Answer1:

It looks like you mixed up x and y in your calculation of a. You're supposed to take the cosine of latitude (y), not longitude (x).

I discovered this by changing your distance to angular_distance (i.e. don't multiply by R) and adding some additional tests:

test(angular_distance("90°0′0″N 0°0′0″E", "89°0′0″N, 0°0′0″E"), math.radians(1)) test(angular_distance("90°0′0″N 0°0′0″E", "80°0′0″N, 0°0′0″E"), math.radians(10)) test(angular_distance("90°0′0″N 0°0′0″E", "50°0′0″N, 0°0′0″E"), math.radians(40)) test(angular_distance("90°0′0″N 0°0′0″E", "50°0′0″N, 20°0′0″E"), math.radians(40))

Recommend

  • date_sun_info() for a specific calendar date
  • Variable naming and team members who speak another language
  • Shift texture on sphere Three.js
  • CSS tabs menu looks ugly when changing browser-zoom
  • Securing web server against MITM attack in Safari
  • Passing parameters in the Form constructor, winforms c#
  • Create cameraview (mask) on surfaceview in android
  • How to get phone heading for augmented reality?
  • nginx, upstream, cors fail
  • Algorithm for placing nodes on a circle considering their distance to eachother
  • Ordering of eigenvectors when calculating eigenvectors using LAPACK's ssteqr
  • Is there a greater chance to collide when comparing GUIDs based on a hash vs \"Guid.NewGuid()?
  • Singular Value Decomposition Implementation
  • Date api for managing off days in an year?
  • how to use asyncio with boost.python?
  • How can integers (keys & values) be added and retrieved from a NSDictionary
  • Background not visible in surface view
  • Implement Gauss-Jordan elimination in Haskell
  • numpy 64bit support in PTVS and numpy System.Int64 casting
  • ng-repeat not working with table but works with list
  • Javascript unload page condition
  • How to resolve dependencies from one gradle project to another gradle project in my Eclipse workspac
  • Cordova Apache wrong module path
  • Specifying virtual keyboard type for EditText in XML
  • Regex for nested values
  • How gzip file gets stored in HDFS
  • Correct implementation of List Iterator methods
  • Android Activity.onWindowFocusChanged doesn't get called from within TabHost
  • Time complexity of a program which involves multiple variables
  • one Local Olampyad Questions on Informatic in 2011
  • Exchange data b/w iOS devices using Bluetooth 4.0
  • Fetching methods from BroadcastReceiver to update UI
  • WinForms: two way TextBox problem
  • Fill an image in a square container while keeping aspect ratio
  • Join two tables and save into third-sql
  • Convert array of 8 bytes to signed long in C++
  • php design question - will a Helper help here?
  • Understanding cpu registers
  • Recursive/Hierarchical Query Using Postgres
  • Running Map reduces the dimensions of the matrices