Earth distance formula
0
In Brief | A set of functions for performing various geographic calculations. Require SciPy |
Language | Python |
# 's
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))
A set of functions for performing various geographic calculations. Require SciPy
Comments