本章節主要參考《python地理空間分析指南》第五章的內容。

一、距離測量

距離測量包括歐式距離，球面距離，以及大地線距離（橢球距離）。主要採用math庫（標準庫，無需下載）進行運算。

1.歐式距離

計算任意兩點之間的距離可以採用距離公式：

例如：計算點A(x1, y1)和點B(x2, y2)之間的歐式距離（UTM座標）

import math

x1 = 456456.23
y1 = 1279721.064 
x2 = 576628.34
y2 = 1071740.33

ditance = math.sqrt((x1-x2)**2+(y1-y2)**2)
print(distance)

# 輸出240.63

2.球面距離（半正矢公式）

利用球面距離的計算公式（Haversine公式）進行計算。

例如：計算點A(-90.212，32.316)和點B(-88.952，30.438)之間的距離（經緯度座標）

import math

x1 = -90.212
y1 = 32.316
x2 = -88.952
y2 = 30.438

# 經緯度轉換為弧度
x_dist = math.radians(x1-x2)
y_dist = math.radians(y1-y2)
y1_rad = math.radians(y1)
y2_rad = math.radians(y2)

a = math.sin(y_dist/2)**2 + math.sin(x_dist/2) * math.cos(y1_rad) * math.cos(y2_rad)
c = 2 * math.asin(math.sqrt(a))
distance = c * 6371
print(distance)

# 輸出240.63
 

3.橢球距離（Vincenty公式）

橢球距離一般採用的是Vincenty公式和NAD83(North American Datum )橢球模型。

Vincenty公式的相關準則如下：

其中，NAD83的相關引數為：

例如：計算點A(-90.212，32.316)和點B(-88.952，30.438)之間的距離（經緯度座標）

import math

x1 = -90.212
y1 = 32.316
x2 = -88.952
y2 = 30.438

# NAD83的橢球引數
a = 6378137        # 半長軸
f = 1/298.257      # 扁平度
b = abs((f*a)-a)   # 半短軸

# 經緯度轉化為弧度
L = math.radians(x2-x1)
U1 = math.atan((1-f) * math.tan(math.radians(y1)))
U2 = math.atan((1-f) * math.tan(math.radians(y2)))

sinU1 = math.sin(U1)
cosU1 = math.cos(U1)
sinU2 = math.sin(U2)
cosU2 = math.cos(U2)
lam = L

for i in range(100):
    sinLam = math.sin(lam)
    cosLam = math.cos(lam)
    sinSigma = math.sqrt((cosU2 * sinLam)**2 + (cosU1 * sinU2 - sinU1 * cosU2 * cosLam)**2)
    cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLam    

    # 判斷是否為重合點
    if sinSigma==0:
        distance = 0
        break

    sigma = math.atan2(sinSigma, cosSigma)
    sinAlpha = cosU1 * cosU2 * sinLam / sinSigma
    cosSqAlpha = 1 - sinAlpha**2
    cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha
    
    if math.isnan(cos2SigmaM):
        cos2SigmaM = 0
    
    C = f/16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
    LP = lam
    lam = L + (1 - c) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + 
            C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)))
    
    if not abs(lam - LP) > 1e-12:
        break

uSq = cosSqAlpha * (a**2 - b**2) / b**2
A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
deltaSigma = B * sinSigma * (cos2SigmaM + B/4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
            B/6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM *cos2SigmaM)))
distance = b * A * (sigma - deltaSigma)
print(distance)
    
# 輸出結果240237.67
 

二、方位計算

方位計算是以正北方向為0°方向，順時針方向角度依次增大的方位角計算。

例如：計算點A(-90.212，32.316)和點B(-88.952，30.438)之間的方位角。

from math import atan2, cos, sin, degrees

x1 = -90.212
y1 = 32.316
x2 = -88.952
y2 = 30.438

angle = atan2(cos(y1)*sin(y2) - sin(y1)*cos(y2)*cos(x2-x1), sin(x2-x1)*cos(y2))
bearing = (degrees(angle) + 360) % 360
print(bearing)

#輸出309.36

三、座標轉換

可以利用utm模組進行墨卡託投影和大地經緯度的相互轉換。utm模組的官方下載地址為：https://pypi.org/project/utm/

例如：墨卡託投影點（5377685.825，479747.045）轉大地經緯度。

import utm

y = 479747.045
x = 5377685.825

zone = 32
band = 'U'
print(utm.to_latlon(y, x, zone, band))

#輸出（48.551，8.725）

例如：大地經緯度（48.551，8.725）轉墨卡託投影座標。

import utm
utm.from_latlon(48.551, 8.725)
# 輸出（479747.045，5377691.373，32，'U'）

四、重投影