一起學 Python 112 : 使用 numpy 庫 計算 經緯度間的距離

之前有個應用需要計算某一個經緯度座標與台灣各鄉鎮的距離

更具體的來說，是我們隨機給出一個座標點，去計算這個座標點離台灣的哪個鄉鎮具有最短距離(在我的資料庫中我將每個鄉鎮都用一個經緯度座標來作為代表)

由於台灣的鄉鎮有大約三百多個，我當然可以用 for 迴圈一個一個去進行計算，但這會導致我們的計算時間過久，畢竟由經緯度計算兩點間的距離會用到大量的三角函數。故我們利用 numpy 庫的矩陣運算來幫我們達成

程式碼如下

import time
import sys
from math import *
import numpy as np
from DATA import *
from datetime import datetime

## ==========================================================================================================================
def np_getDistance(A , B ):# 先緯度後經度
    ra = 6378140  # radius of equator: meter
    rb = 6356755  # radius of polar: meter
    flatten = 0.003353 # Partial rate of the earth
    # change angle to radians
    
    radLatA = np.radians(A[:,0])
    radLonA = np.radians(A[:,1])
    radLatB = np.radians(B[:,0])
    radLonB = np.radians(B[:,1])
 
    pA = np.arctan(rb / ra * np.tan(radLatA))
    pB = np.arctan(rb / ra * np.tan(radLatB))
    
    x = np.arccos( np.multiply(np.sin(pA),np.sin(pB)) + np.multiply(np.multiply(np.cos(pA),np.cos(pB)),np.cos(radLonA - radLonB)))
    c1 = np.multiply((np.sin(x) - x) , np.power((np.sin(pA) + np.sin(pB)),2)) / np.power(np.cos(x / 2),2)
    c2 = np.multiply((np.sin(x) + x) , np.power((np.sin(pA) - np.sin(pB)),2)) / np.power(np.sin(x / 2),2)
    dr = flatten / 8 * (c1 - c2)
    distance = 0.001 * ra * (x + dr)
    return distance
## ==========================================================================================================================
Find_loc = np.matrix([[24.797082,120.992036]])#求所有鄉鎮至該經緯度座標的距離

starttime = datetime.now()
# 距縣市距離(全部)
disALL = np_getDistance(loltALL_center,Find_loc)
# 距離各縣市最短距離(全部)
disMinName = nameALL[int(disALL.argmin(axis=0))]
disMinDistance = floor(disALL[int(disALL.argmin(axis=0))])

endtime = datetime.now()

print("離 %s 最近, %1.1f 公里"%(disMinName,disMinDistance))
print("use %s usec"%(endtime-starttime).microseconds)

運算速度很不錯，我很難想像用 for 迴圈去算三百多次那堆 sin , cos arctan 等等的要花多久的時間

其中各鄉鎮的經緯度座標，可由此下載


參考資料
numpy库矩阵信息的获取（最大值最小值、平均值、中值、方差标准差、求和）
Python 中的几种矩阵乘法 np.dot, np.multiply, *
Python：一篇文章掌握Numpy的基本用法
使用 Python 來認識矩陣
【python】 计算向量欧氏距离的小代码 numpy
计算Python Numpy向量之间的欧氏距离【python】numpy数组（array）扩充（复制）方法repeat和tile的使用
Python中的numpy矩阵运算
[Python] 使用Numpy建立矩陣之基本觀念釐清
python 常用函數整理
numpy中的convolve的理解
numpy教程：数学函数和基本统计函数