forked from michaelpollmann/spatialTreat-example
-
Notifications
You must be signed in to change notification settings - Fork 0
/
distance.py
44 lines (39 loc) · 1.45 KB
/
distance.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
import numpy as np
import numba
from math import radians, cos, sin, asin, sqrt
def distance(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
r = 6378137 # Radius of earth at equator in meters. Use 3963.191 for miles
return c * r
@numba.jit
def distance_vec(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
#Convert decimal degrees to Radians:
lon1 = np.deg2rad(lon1)
lat1 = np.deg2rad(lat1)
lon2 = np.deg2rad(lon2)
lat2 = np.deg2rad(lat2)
#Implementing Haversine Formula:
dlon = np.subtract(lon2, lon1)
dlat = np.subtract(lat2, lat1)
# haversine calculation
a = np.add(np.power(np.sin(np.divide(dlat, 2)), 2),
np.multiply(np.cos(lat1),
np.multiply(np.cos(lat2),
np.power(np.sin(np.divide(dlon, 2)), 2))))
c = np.multiply(2, np.arcsin(np.sqrt(a)))
r = 6378137 # Radius of earth at equator in meters. Use 3963.191 for miles
return c*r