Hızlı Haversine Yaklaşımı (Python/Pandas)

Bir Panda'nın veri çerçevesindeki her satır, 2 noktadan oluşan lat/lng koordinatlarını içerir. Aşağıdaki Python kodunu kullanarak, bu 2 nokta arasındaki mesafelerin birçok (milyon) satır için hesaplanması çok uzun zaman alıyor!

2 noktanın birbirinden 50 mil uzakta olduğu ve doğruluğunun çok önemli olmadığı düşünülürse, hesaplamayı daha hızlı yapmak mümkün müdür?

from math import radians, cos, sin, asin, sqrt
def haversine(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)) 
    km = 6367 * c
    return km


for index, row in df.iterrows():
    df.loc[index, 'distance'] = haversine(row['a_longitude'], row['a_latitude'], row['b_longitude'], row['b_latitude'])
19
Vay yorumlarda öneriler benden öte! Sadece sqrt ((lat2 - lat1) ^^ 2 + (lon2 - lon1) ^^ 2) yapıp yapamayacağımı merak ediyordum, işe yarayacak mı? Bütün noktalar New York'ta.
katma yazar Nyxynyx, kaynak
@ballsdotballs Öklid mesafesi benim için işe yarayacak gibi görünüyor! sqrt ((lat2 - lat1) ^^ 2 + (lon2 - lon1) ^^ 2) sonucu millere nasıl çevrilir?
katma yazar Nyxynyx, kaynak
Yaklaşımdan daha iyi bir yaklaşım, fonksiyonun tam olarak neden çok uzun sürdüğünü anlamak için fonksiyonun profilini yapmaktır, ardından fonksiyonun yükü çok fazla olmayan bir C fonksiyonuna çevrilmesi için ctypes/Cython/numba kullanılır. Her pandada Series veri sütununu oluşturan verilerin numcode dizisini value kullanmak için arama kuralınızı değiştirmeniz gerekebilir ve ayrıca numpy.ctypeslib bir numpy dizisinden ctypes uyumlu bir diziye kolay dönüştürme için. Çok görünüyor, ama gerçekten Python'da C işlevlerine erişmek için oldukça kolay bir yol.
katma yazar ely, kaynak
DataFrame gibi ilişkisel bir yapıya kaydetmek yerine, verilerden bir k-d ağacı oluşturmayı da düşünebilirsiniz. O zaman belirli bir noktadaki komşuları elde etmek ucuz olurdu ve belki de yalnızca talep üzerine olan mesafeleri hesaplayabilirsiniz. Uygulama her zaman her çifte ihtiyaç duyuyor mu? Diğer bir seçenek ise noktaları kümelemek ve her kümenin centroid/ortalamasını bir vekil olarak kullanmak olabilir. Sonra herhangi iki nokta arasındaki mesafeye sadece küme merkezleri arasındaki mesafeyle yaklaşılır. Bunun gibi süslü bir şeyin kaba kuvvetten gerçekten daha iyi olup olmadığı spekülatif.
katma yazar ely, kaynak
Kümelenme durumunda, aynı kümedeki tüm çiftler için gerçek mesafeleri (her küme için) hesaplamak isteyebilirsiniz, böylece bu noktalar etkin bir sıfır mesafesine sahip olmaz. Bu aynı zamanda basit bir paralelleştirmeye de uygundur: her bir işlemden sorumlu olduğu endekslerin bir listesi ile birlikte verilerin bir kopyasını (veya bir kopyasını yüklemesini sağlayın) verin. Daha sonra bu işlem, bu endekslerin diğer tüm endekslere karşı ikili mesafesini hesaplar ve bir yere yazar. Yalnızca milyonlarca satırla, mütevazı donanımda bu makul olmalıdır.
katma yazar ely, kaynak
Evet, öklid yaklaşımı, yeterince küçük mesafeler için iyi çalışacaktır. Bunun için bir Uygulaması yapmanıza bile gerek yok, doğrudan veri çerçevesindeki sütunları kullanabilirsiniz.
katma yazar reptilicus, kaynak
Çoğu aday için hesaplama yapmaktan kaçınmak mümkün olabilir. Min. Ve maksimum uzunlukları ve enlemleri başlangıç ​​noktanızın 50 mil uzağında hesaplayın. Ardından adayların çoğunu ayıklamak için bu mayınları ve makbuzları kullanın.
katma yazar Steven Rumbalski, kaynak
Pypi'de bir haversine modülü var C de uygulanır
katma yazar Steven Rumbalski, kaynak
Lat1, long1, lat1, long1 aynı satırda hesaplanmış mı?
katma yazar dustin, kaynak
@Nyxynyx Sorunuzda sağladığınız işlev, büyük çember mesafesini verir. Yorumunuzdaki hesaplama öklid mesafesini verir. Dünyanın yarıçapı çok büyük olduğundan, küçük mesafeler için kesinlikle öklid versiyonuyla yaklaşık olarak tahmin edebilirsiniz.
katma yazar derricw, kaynak

7 cevap

İşte aynı işlevin vectorized numpy versiyonu:

import numpy as np

def haversine_np(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)

    All args must be of equal length.    

    """
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

    dlon = lon2 - lon1
    dlat = lat2 - lat1

    a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2

    c = 2 * np.arcsin(np.sqrt(a))
    km = 6367 * c
    return km

Girdilerin hepsi değer dizileridir ve milyonlarca noktayı anında yapabilmelidir. Şart, girdilerin ndarralar olması ama panda tablonuzun sütunlarının çalışması.

Örneğin, rastgele oluşturulmuş değerlerle:

>>> import numpy as np
>>> import pandas
>>> lon1, lon2, lat1, lat2 = np.random.randn(4, 1000000)
>>> df = pandas.DataFrame(data={'lon1':lon1,'lon2':lon2,'lat1':lat1,'lat2':lat2})
>>> km = haversine_np(df['lon1'],df['lat1'],df['lon2'],df['lat2'])

Python'da veri dizileri arasında döngü yapmak çok yavaştır. Numpy, tüm veri dizileri üzerinde çalışan, döngüden kaçınmanıza ve performansı büyük ölçüde artırmanıza olanak sağlayan işlevler sunar.

Bu, vectorization örneğidir.

44
katma
ufak boyutlu dizilerde bile çalıştığını biliyoruz;) boyut projeksiyonu ile
katma yazar LetsPlayYahtzee, kaynak
Bunun için çok teşekkür ederim. Küçük öneri: giriş biçimini netleştirmek için rastgele değerler yerine gerçek koordinatlarla gerçek dünyaya örnek kullanım ekleyin.
katma yazar Dennis Golomazov, kaynak
Bunun bir çift argüman bir Series ve bir diğeri tuple olduğunda da işe yaradığını unutmayın: haversine_np (pd.Series ([- 74.00594, -122.41942]), pd.Series ([40.71278] , 37.77493]), -87.65005, 41.85003) (New York, San Francisco) ve Chicago arasındaki mesafeyi hesaplar.
katma yazar Dennis Golomazov, kaynak
Başka bir küçük öneri: işlev bağımsız değişkenlerinin sırasını lat, lon olarak değiştirmek isteyebilirsiniz. Birçok kaynaktan, enlem önce gider; en.wikipedia.org/wiki/Horizontal_position_representation 'de.
katma yazar Dennis Golomazov, kaynak
Dizi programlama terimini bilmeniz iyi, MATLAB ile karşılaşmadı.
katma yazar Divakar, kaynak
@DennisGolomazov afaik GIS kütüphaneleri ve yazılımları genellikle lon, lat. Enlem ilk olarak esas olarak bir Google Haritalar seçeneği olarak görünmektedir. Şahsen bunu lat olarak tuhaf buldum, Lon.
katma yazar Leandro Lima, kaynak

Tamamen açıklayıcı bir örnek olarak, @ballsdotballs'dan gelen cevapta numpy sürümünü aldım ve ayrıca ctypes aracılığıyla çağrılması için bir C arkadaşı uygulaması yaptım. numpy çok iyi bir şekilde optimize edilmiş bir araç olduğundan, C kodumun bu kadar verimli olma ihtimali çok az, ancak biraz daha yakın olması gerekir. Buradaki en büyük avantaj, C tipleri ile bir örnek kullanarak, kendi kişisel C işlevlerinizi çok fazla ek yük olmadan Python'a nasıl bağlayabileceğinizi görmenize yardımcı olabilir. Bu özellikle, bu küçük parçayı Python yerine bir C kaynağına yazarak, daha küçük bir parçayı daha büyük bir hesaplama için optimize etmek istediğinizde çok güzel. Basitçe numpy kullanmak çoğu zaman sorunu çözecektir, ancak tüm numpy kodlarına gerçekten ihtiyaç duymadığınız ve kuplajı eklemek istemediğiniz durumlarda bazı kod boyunca numpy veri türlerinin kullanılmasını gerektiren durumlarda, yerleşik ctypes kitaplığına nasıl bırakılacağını bilmek ve kendin yapmak çok kolaydır.

Öncelikle haversine.c adlı C kaynak dosyamızı oluşturalım:

#include 
#include 
#include 

int haversine(size_t n, 
              double *lon1, 
              double *lat1, 
              double *lon2, 
              double *lat2,
              double *kms){

    if (   lon1 == NULL 
        || lon2 == NULL 
        || lat1 == NULL 
        || lat2 == NULL
        || kms == NULL){
        return -1;
    }

    double km, dlon, dlat;
    double iter_lon1, iter_lon2, iter_lat1, iter_lat2;

    double km_conversion = 2.0 * 6367.0; 
    double degrees2radians = 3.14159/180.0;

    int i;
    for(i=0; i < n; i++){
        iter_lon1 = lon1[i] * degrees2radians;
        iter_lat1 = lat1[i] * degrees2radians;
        iter_lon2 = lon2[i] * degrees2radians;
        iter_lat2 = lat2[i] * degrees2radians;

        dlon = iter_lon2 - iter_lon1;
        dlat = iter_lat2 - iter_lat1;

        km = pow(sin(dlat/2.0), 2.0) 
           + cos(iter_lat1) * cos(iter_lat2) * pow(sin(dlon/2.0), 2.0);

        kms[i] = km_conversion * asin(sqrt(km));
    }

    return 0;
}

// main function for testing
int main(void) {
    double lat1[2] = {16.8, 27.4};
    double lon1[2] = {8.44, 1.23};
    double lat2[2] = {33.5, 20.07};
    double lon2[2] = {14.88, 3.05};
    double kms[2]  = {0.0, 0.0};
    size_t arr_size = 2;

    int res;
    res = haversine(arr_size, lon1, lat1, lon2, lat2, kms);
    printf("%d\n", res);

    int i;
    for (i=0; i < arr_size; i++){
        printf("%3.3f, ", kms[i]);
    }
    printf("\n");
}

C kurallarına uymaya çalıştığımızı unutmayın. Veri değişkenlerini referans olarak açıkça iletmek, bir boyut değişkeni için size_t öğesini kullanmak ve haversine işlevimizin, iletilen girdilerden birini beklenen verileri içerecek şekilde değiştirerek çalışmasını beklemek çıkışta. İşlev aslında, işlevin diğer C düzeyi tüketicileri tarafından kullanılabilecek bir başarı/başarısızlık bayrağı olan bir tamsayı döndürür.

Python'un içindeki bu C'ye özgü tüm bu konuların üstesinden gelmenin bir yolunu bulmamız gerekecek.

Şimdi fonksiyonun numpy versiyonunu, bazı ithalat ve bazı test verilerini haversine.py adlı bir dosyaya koyalım:

import time
import ctypes
import numpy as np
from math import radians, cos, sin, asin, sqrt

def haversine(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(np.radians, [lon1, lat1, lon2, lat2])
    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = (np.sin(dlat/2)**2 
         + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2)
    c = 2 * np.arcsin(np.sqrt(a)) 
    km = 6367 * c
    return km

if __name__ == "__main__":
    lat1 = 50.0 * np.random.rand(1000000)
    lon1 = 50.0 * np.random.rand(1000000)
    lat2 = 50.0 * np.random.rand(1000000)
    lon2 = 50.0 * np.random.rand(1000000)

    t0 = time.time()
    r1 = haversine(lon1, lat1, lon2, lat2)
    t1 = time.time()
    print t1-t0, r1

0 ile 50 arasında rastgele seçilen lats ve lons (derece cinsinden) yapmayı seçtim, ancak bu açıklama için çok fazla önemli değil.

Yapmamız gereken bir sonraki şey, C modülümüzü Python tarafından dinamik olarak yüklenebilecek şekilde derlemektir. Bir Linux sistemi kullanıyorum (Google'da diğer sistemlere çok kolay örnekler bulabilirsiniz), bu yüzden amacım haversine.c 'ı paylaşılan bir nesneye derlemek, şöyle:

gcc -shared -o haversine.so -fPIC haversine.c -lm

Ayrıca bir yürütülebilir dosyayı derleyebilir ve C programının main işlevinin neyi gösterdiğini görmek için çalıştırabiliriz:

> gcc haversine.c -o haversine -lm
> ./haversine
0
1964.322, 835.278, 

Şimdi haversine.so paylaşılan nesnesini derlediğimize göre, Python'a yüklemek için ctypes öğesini kullanabiliriz ve bunu yapmak için dosyanın yolunu sağlamamız gerekir:

lib_path = "/path/to/haversine.so" # Obviously use your real path here.
haversine_lib = ctypes.CDLL(lib_path)

Şimdi haversine_lib.haversine , giriş ve çıkışların doğru bir şekilde yorumlandığından emin olmak için bazı manuel tip marşlamalar yapmamız gerekebilir, ancak hemen hemen bir Python işlevi gibi davranır.

numpy actually provides some nice tools for this and the one I'll use here is numpy.ctypeslib. We're going to build a pointer type that will allow us to pass around numpy.ndarrays to these ctypes-loaded functions as through they were pointers. Here's the code:

arr_1d_double = np.ctypeslib.ndpointer(dtype=np.double, 
                                       ndim=1, 
                                       flags='CONTIGUOUS')

haversine_lib.haversine.restype = ctypes.c_int
haversine_lib.haversine.argtypes = [ctypes.c_size_t,
                                    arr_1d_double, 
                                    arr_1d_double,
                                    arr_1d_double,
                                    arr_1d_double,
                                    arr_1d_double] 

haversine_lib.haversine function proxy'ye, bağımsız değişkenlerini istediğimiz türlere göre yorumlamasını söylediğimize dikkat edin.

Şimdi, Python'dan denemek için geriye kalan sadece bir boyut değişkeni yapmak ve sonuç verilerini içerecek şekilde mutasyona uğratılacak (C kodunda olduğu gibi) bir dizi, sonra arayabiliriz o:

size = len(lat1)
output = np.empty(size, dtype=np.double)
print "====="
print output
t2 = time.time()
res = haversine_lib.haversine(size, lon1, lat1, lon2, lat2, output)
t3 = time.time()
print t3 - t2, res
print type(output), output

Hepsini haversine.py __main__ bloğuna koyarak, tüm dosya şimdi şöyle görünür:

import time
import ctypes
import numpy as np
from math import radians, cos, sin, asin, sqrt

def haversine(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(np.radians, [lon1, lat1, lon2, lat2])
    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = (np.sin(dlat/2)**2 
         + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2)
    c = 2 * np.arcsin(np.sqrt(a)) 
    km = 6367 * c
    return km

if __name__ == "__main__":
    lat1 = 50.0 * np.random.rand(1000000)
    lon1 = 50.0 * np.random.rand(1000000)
    lat2 = 50.0 * np.random.rand(1000000)
    lon2 = 50.0 * np.random.rand(1000000)

    t0 = time.time()
    r1 = haversine(lon1, lat1, lon2, lat2)
    t1 = time.time()
    print t1-t0, r1

    lib_path = "/home/ely/programming/python/numpy_ctypes/haversine.so"
    haversine_lib = ctypes.CDLL(lib_path)
    arr_1d_double = np.ctypeslib.ndpointer(dtype=np.double, 
                                           ndim=1, 
                                           flags='CONTIGUOUS')

    haversine_lib.haversine.restype = ctypes.c_int
    haversine_lib.haversine.argtypes = [ctypes.c_size_t,
                                        arr_1d_double, 
                                        arr_1d_double,
                                        arr_1d_double,
                                        arr_1d_double,
                                        arr_1d_double]

    size = len(lat1)
    output = np.empty(size, dtype=np.double)
    print "====="
    print output
    t2 = time.time()
    res = haversine_lib.haversine(size, lon1, lat1, lon2, lat2, output)
    t3 = time.time()
    print t3 - t2, res
    print type(output), output

Python ve ctypes sürümlerini ayrı ayrı çalıştıracak ve zamanlandıracak ve bazı sonuçları yazdıracak şekilde çalıştırmak için sadece yapabiliriz.

python haversine.py

hangi görüntüler:

0.111340045929 [  231.53695005  3042.84915093   169.5158946  ...,  1359.2656769
  2686.87895954  3728.54788207]
=====
[  6.92017600e-310   2.97780954e-316   2.97780954e-316 ...,
   3.20676686e-001   1.31978329e-001   5.15819721e-001]
0.148446083069 0
 [  231.53675618  3042.84723579   169.51575588 ...,  1359.26453029
  2686.87709456  3728.54493339]

Beklendiği gibi, numpy sürümü biraz daha hızlı (1 milyon uzunluğundaki vektörler için 0.11 saniye) ancak hızlı ve kirli ctypes sürümümüz hiç eğik değil: saygın bir 0.148 saniye Aynı veri.

Bunu Python'da saf bir for-loop çözümü ile karşılaştıralım:

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

def slow_haversine(lon1, lat1, lon2, lat2):
    n = len(lon1)
    kms = np.empty(n, dtype=np.double)
    for i in range(n):
       lon1_v, lat1_v, lon2_v, lat2_v = map(
           radians, 
           [lon1[i], lat1[i], lon2[i], lat2[i]]
       )

       dlon = lon2_v - lon1_v 
       dlat = lat2_v - lat1_v 
       a = (sin(dlat/2)**2 
            + cos(lat1_v) * cos(lat2_v) * sin(dlon/2)**2)
       c = 2 * asin(sqrt(a)) 
       kms[i] = 6367 * c
    return kms

Bunu diğerleri ile aynı Python dosyasına koyduğumda ve aynı milyon elementli verilere zaman verdiğimde, sürekli olarak makinemde yaklaşık 2.65 saniyelik bir zaman görüyorum.

Bu nedenle, hızlı bir şekilde ctypes 'e geçerek hızı yaklaşık 18 faktörü ile yükseltiriz. Çıplak, bitişik verilere erişimden faydalanabilecek birçok hesaplama için, bundan daha da çok kazançlar elde edersiniz.

Sadece çok açık olmak gerekirse, bunu numpy kullanmaktan daha iyi bir seçenek olarak kabul etmiyorum. Bu, tam olarak numpy 'nin çözmek için yapıldığı bir problemdir ve bu nedenle her ikisi (a) numpy ctypes kodunuzu homrew > uygulamanızdaki veri türleri ve (b) kodunuzu numpy eşdeğeriyle eşlemenin kolay bir yolu yoktur, çok verimli değildir.

Ancak, C'ye bir şey yazmayı tercih ettiğiniz, ancak Python'da çağırdığınız veya numpy bağımlılığının pratik olmadığı durumlarda (gömülü bir sistemde) bu durumlar için nasıl yapılacağını bilmek hala çok yararlıdır. örneğin numpy yüklenemiyor).

11
katma

Scikit-learn'ın kullanılmasına izin verilirse, aşağıdakileri bir şans verirdim:

from sklearn.neighbors import DistanceMetric
dist = DistanceMetric.get_metric('haversine')

# example data
lat1, lon1 = 36.4256345, -5.1510261
lat2, lon2 = 40.4165, -3.7026
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

X = [[lat1, lon1],
     [lat2, lon2]]
kms = 6367
print(kms * dist.pairwise(X))
7
katma

Scikit-learn'ın kullanılmasına izin verilirse, aşağıdakileri bir şans verirdim:

from sklearn.neighbors import DistanceMetric
dist = DistanceMetric.get_metric('haversine')

# example data
lat1, lon1 = 36.4256345, -5.1510261
lat2, lon2 = 40.4165, -3.7026
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

X = [[lat1, lon1],
     [lat2, lon2]]
kms = 6367
print(kms * dist.pairwise(X))
7
katma

Bu cevapların bazıları, dünyanın yarıçapını "yuvarlıyor". Bunları diğer mesafe hesaplayıcılarına karşı kontrol ederseniz ( coğrafi gibi), bu işlevler kapalı olacaktır.

Cevabı mil cinsinden görmek istiyorsanız, aşağıdaki dönüşüm sabiti için R = 3959.87433 kodunu değiştirebilirsiniz.

Kilometre istiyorsanız, R = 6372.8 'ı kullanın.

lon1 = -103.548851
lat1 = 32.0004311
lon2 = -103.6041946
lat2 = 33.374939


def haversine(lat1, lon1, lat2, lon2):

      R = 3959.87433 # this is in miles.  For Earth radius in kilometers use 6372.8 km

      dLat = radians(lat2 - lat1)
      dLon = radians(lon2 - lon1)
      lat1 = radians(lat1)
      lat2 = radians(lat2)

      a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2
      c = 2*asin(sqrt(a))

      return R * c

print(haversine(lat1, lon1, lat2, lon2))
0
katma

Bu cevapların bazıları, dünyanın yarıçapını "yuvarlıyor". Bunları diğer mesafe hesaplayıcılarına karşı kontrol ederseniz ( coğrafi gibi), bu işlevler kapalı olacaktır.

Cevabı mil cinsinden görmek istiyorsanız, aşağıdaki dönüşüm sabiti için R = 3959.87433 kodunu değiştirebilirsiniz.

Kilometre istiyorsanız, R = 6372.8 'ı kullanın.

lon1 = -103.548851
lat1 = 32.0004311
lon2 = -103.6041946
lat2 = 33.374939


def haversine(lat1, lon1, lat2, lon2):

      R = 3959.87433 # this is in miles.  For Earth radius in kilometers use 6372.8 km

      dLat = radians(lat2 - lat1)
      dLon = radians(lon2 - lon1)
      lat1 = radians(lat1)
      lat2 = radians(lat2)

      a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2
      c = 2*asin(sqrt(a))

      return R * c

print(haversine(lat1, lon1, lat2, lon2))
0
katma

Bu cevapların bazıları, dünyanın yarıçapını "yuvarlıyor". Bunları diğer mesafe hesaplayıcılarına karşı kontrol ederseniz ( coğrafi gibi), bu işlevler kapalı olacaktır.

Cevabı mil cinsinden görmek istiyorsanız, aşağıdaki dönüşüm sabiti için R = 3959.87433 kodunu değiştirebilirsiniz.

Kilometre istiyorsanız, R = 6372.8 'ı kullanın.

lon1 = -103.548851
lat1 = 32.0004311
lon2 = -103.6041946
lat2 = 33.374939


def haversine(lat1, lon1, lat2, lon2):

      R = 3959.87433 # this is in miles.  For Earth radius in kilometers use 6372.8 km

      dLat = radians(lat2 - lat1)
      dLon = radians(lon2 - lon1)
      lat1 = radians(lat1)
      lat2 = radians(lat2)

      a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2
      c = 2*asin(sqrt(a))

      return R * c

print(haversine(lat1, lon1, lat2, lon2))
0
katma