8

I'm setting up a small program to take 2 geographical coordinates from a user and then calculate the distance between them(taking into account the curvature of the earth). So I looked up wikipedia on what the formula is here.

I basically set up my python function based on that and this is what I came up with:

def geocalc(start_lat, start_long, end_lat, end_long):
    start_lat = math.radians(start_lat)
    start_long = math.radians(start_long)
    end_lat = math.radians(end_long)
    end_long = math.radians(end_long)

    d_lat = start_lat - end_lat
    d_long = start_long - end_long

    EARTH_R = 6372.8

    c = math.atan((math.sqrt( (math.cos(end_lat)*d_long)**2 +( (math.cos(start_lat)*math.sin(end_lat)) - (math.sin(start_lat)*math.cos(end_lat)*math.cos(d_long)))**2)) / ((math.sin(start_lat)*math.sin(end_lat)) + (math.cos(start_lat)*math.cos(end_lat)*math.cos(d_long))) )

    return EARTH_R*c

The problem is that the results come out really inaccurate. I'm new to python so some help or advice would be greatly appreciated!

John Machin
  • 81,303
  • 11
  • 141
  • 189
Darkphenom
  • 647
  • 2
  • 8
  • 14

4 Answers4

12

You've got 4 or 5 or 6 problems:

(1) end_lat = math.radians(end_long) should be end_lat = math.radians(end_lat)

(2) you are missing some stuff as somebody already mentioned, most probably because

(3) your code is illegible (line far too long, redundant parentheses, 17 pointless instances of "math.")

(4) you didn't notice the remark in the Wikipedia article about using atan2()

(5) You may have been swapping lat and lon when entering your coordinates

(6) delta(latitude) is computed unnecessarily; it doesn't appear in the formula

Putting it all together:

from math import radians, sqrt, sin, cos, atan2

def geocalc(lat1, lon1, lat2, lon2):
    lat1 = radians(lat1)
    lon1 = radians(lon1)
    lat2 = radians(lat2)
    lon2 = radians(lon2)

    dlon = lon1 - lon2

    EARTH_R = 6372.8

    y = sqrt(
        (cos(lat2) * sin(dlon)) ** 2
        + (cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dlon)) ** 2
        )
    x = sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(dlon)
    c = atan2(y, x)
    return EARTH_R * c



>>> geocalc(36.12, -86.67, 33.94, -118.40)
2887.2599506071115
>>> geocalc(-6.508, 55.071, -8.886, 51.622)
463.09798886300376
>>> geocalc(55.071, -6.508, 51.622, -8.886)
414.7830891822618
John Machin
  • 81,303
  • 11
  • 141
  • 189
4

You can use the geopy module which has a built-in function for distance calculations, scroll down to "Calculating Distances" in the link below: https://pypi.python.org/pypi/geopy

MERose
  • 4,048
  • 7
  • 53
  • 79
vasek1
  • 13,541
  • 11
  • 32
  • 36
  • That's a good suggestion, but your answer would be better if you could include some code on how to achieve OP's goal. – MERose Jun 08 '15 at 11:14
4

This works (print f returns 2887.26 km as per the worked example @ http://en.wikipedia.org/wiki/Great-circle_distance):

import math

def geocalc(start_lat, start_long, end_lat, end_long):

    start_lat = math.radians(start_lat)
    start_long = math.radians(start_long)
    end_lat = math.radians(end_lat)
    end_long = math.radians(end_long)

    d_lat = math.fabs(start_lat - end_lat)
    d_long = math.fabs(start_long - end_long)

    EARTH_R = 6372.8

    y = ((math.sin(start_lat)*math.sin(end_lat)) + (math.cos(start_lat)*math.cos(end_lat)*math.cos(d_long)))

    x = math.sqrt((math.cos(end_lat)*math.sin(d_long))**2 + ( (math.cos(start_lat)*math.sin(end_lat)) - (math.sin(start_lat)*math.cos(end_lat)*math.cos(d_long)))**2)

    c = math.atan(x/y)

    return EARTH_R*c

f = geocalc(36.12, -86.67, 33.94, -118.40)
print f

Notice this line in your submission: end_lat = math.radians(end_long)

sgallen
  • 2,079
  • 13
  • 10
  • You don't need `fabs()`, because `cos(-x) == cos(x)` – John Machin Jan 14 '12 at 02:38
  • -1 This can go go spectacularly wrong because you are using `atan` instead of `atan2`. Examples: going due south from 45 deg lat to -45 deg lat is 10010 km (OK) but from 45.1 deg to -45.1 deg (a slightly longer distance) produces MINUS 9988 km. A trip from north pole to south pole produces ZERO km! – John Machin Jan 14 '12 at 04:36
  • @JohnMachin regarding the fabs() comment you are correct for this particular implementation but when teaching the concept of deltas the use of absolute values is correct not simple subtraction. – sgallen Jan 14 '12 at 13:44
  • I don't understand what you mean by "teaching the concept of deltas". If I'm at `(x1, y1)` and I want to be at `(x2, y2)`, I tell my stepper motors `deltax, deltay = (x2 - x1, y2 - y1)`, no `abs()` in sight; what do you do? – John Machin Jan 14 '12 at 20:33
  • More: (2) In any case the Wikipedia article referred to them as "differences", not "deltas". Where I come from, difference has a sign. (3) Which of the many deltas mentioned in this page http://en.wikipedia.org/wiki/Delta are you talking about? – John Machin Jan 14 '12 at 22:24
3

I think you missed a math.sin(d_long) towards the beginning, should maybe be this:

 c = math.atan((math.sqrt( (math.cos(end_lat)*math.sin(d_long))**2 +( (math.cos(start_lat)*math.sin(end_lat)) - (math.sin(start_lat)*math.cos(end_lat)*math.cos(d_long)))**2)) / ((math.sin(start_lat)*math.sin(end_lat)) + (math.cos(start_lat)*math.cos(end_lat)*math.cos(d_long))) )
nmjohn
  • 1,432
  • 7
  • 16