經緯度計算距離與方位角

距離

根據經緯度計算距離是有固定公式的，按照公式寫個函數很簡單。

import math
def cal_dis(lat1, lon1,lat2, lon2):
    latitude1 = (math.pi/180)*lat1
    latitude2 = (math.pi/180)*lat2
    longitude1 = (math.pi/180)*lon1
    longitude2= (math.pi/180)*lon2
    #因此AB兩點的球面距離為:{arccos[sinb*siny+cosb*cosy*cos(a-x)]}*R
    #地球半徑
    R = 6378.137
    d = math.acos(math.sin(latitude1)*math.sin(latitude2)+ math.cos(latitude1)*math.cos(latitude2)*math.cos(longitude2-longitude1))*R
    return d


if __name__ == '__main__':
    print cal_dis(23.0,101.1,23.06,113.34)

注意，經緯度是角度，而三角函數的輸入是弧度，角度與弧度的關系  180=pi

 

有一次我在一個項目中發現 matlab 的計算距離和上面這個函數有些出入，雖然相差不大，但是確實不一樣，當時也折騰了好久，后來發現 python 有個自帶的函數計算結果和 matlab 相同。

from geopy.distance import geodesic

geodesic((30.28708,120.12802999999997), (28.7427,115.86572000000001)).m

 

方位角

def calc_azimuth(lat1, lon1, lat2, lon2):
    lat1_rad = lat1 * math.pi / 180
    lon1_rad = lon1 * math.pi / 180
    lat2_rad = lat2 * math.pi / 180
    lon2_rad = lon2 * math.pi / 180

    y = math.sin(lon2_rad - lon1_rad) * math.cos(lat2_rad)
    x = math.cos(lat1_rad) * math.sin(lat2_rad) - \
        math.sin(lat1_rad) * math.cos(lat2_rad) * math.cos(lon2_rad - lon1_rad)

    brng = math.atan2(y, x) * 180 / math.pi

    return float((brng + 360.0) % 360.0)

所謂方位角是與正北方向、順時針之間的夾角