Earth distance formula 0

 1import numpy
 2import math
 3from scipy import *
 4
 5eq_rad     = 6378.137 #eq radius in km
 6polar_rad  = 6356.752 #polar radius in km
 7
 8def mercator_coords(geo_pt, center):
 9    '''
10    Projects the given coordinates using Mercator projection
11    with respect to `center`.
12    '''
13
14    x = geo_pt[0: , :1] - center[0]
15    y = arctanh(numpy.sin(geo_pt[0 : , 1:]*(numpy.pi/360)))
16    
17    return hstack((x,y))
18
19def distance_in_km(lon1, lat1, lon2, lat2):
20   '''
21   Given a set of geo coordinates (in degrees) it will return the distance in km
22   '''
23
24   #convert to radians
25   lon1 = lon1*2*numpy.pi/360
26   lat1 = lat1*2*numpy.pi/360
27   lon2 = lon2*2*numpy.pi/360
28   lat2 = lat2*2*numpy.pi/360
29
30   R = earthRadius((lat1+lat2)/2) #km
31
32   #haversine formula - angles in radians
33   deltaLon = numpy.abs(lon1-lon2)
34   deltaLat = numpy.abs(lat1-lat2)
35
36   dOverR = haver_sin(deltaLat) + numpy.cos(lat1)*numpy.cos(lat2)*haver_sin(deltaLon)
37
38   return R * arc_haver_sin(dOverR)
39
40def earth_radius(lat):
41   '''
42   Given a latitude in radias returns earth radius in km
43   '''
44
45   top = (eq_rad**2 * numpy.cos(lat))**2 + (polar_rad**2 * numpy.sin(lat))**2
46   bottom = (eq_rad * numpy.cos(lat))**2 + (polar_rad * numpy.sin(lat))**2
47   
48   return numpy.sqrt(top/bottom)
49
50def haver_sin(x):
51   return numpy.sin(x/2) ** 2
52
53def arc_haver_sin(x):
54   return 2*numpy.arcsin(numpy.sqrt(x))