空間權重矩陣的那些事(八)-球面距離權重矩陣
前段時間有人向我咨詢了根據經緯度計算球面距離的方法，希望我出一篇文章，所以就有了這篇文章。必要文件可通過后台回復「地理經緯度」獲取。
我首先想到的是matlab的「jplv7」工具箱里的「distance」函數，下面是具體的程序內容：
function D = distance(xc,yc) % PURPOSE: Computes the list of pairwise distances for a given set of locations (loc). % ---------------------------------------------------------- % Usage: D = distance(xc,yc) % where: xc,yc are vectors of latt-long coordinates for each location % ---------------------------------------------------------- % Returns: D = (n x n)-matrix of pairwise distances % Written by: TONY E. SMITH, 2/10/98 n = length(xc) ; %number of locations % Start procedure. X = xc ; %column vector Y = yc ; %column vector U = ones(n,1) ; %column vector XX = X*U' - U*X' ; YY = Y*U' - U*Y' ; D = (XX.^2 + YY.^2) ;
根據程序內容，距離計算公式是：
$$
D=\,\,\left( x1-x2 \right) ^2+\left( y1-y2 \right) ^2
$$
其中，(x1, y1) 是點 A 的坐標，(x2, y2) 是點 B 的坐標。這明顯不是球面距離。實際上，如果經緯度數據如果是未投影的數據（即空間坐標，單位為 °），經過上面的函數計算得到的距離仍然不是地理矩陣，因為其單位不是 mi 或者 km 等距離單位。使用該函數時需要特別注意！總之，只有平面坐標通過上面的函數計算的結果還是真實的地理距離。
通過查找資料，發現matlab本身就存在可以計算球面距離的函數，而且函數名也是「distance」。這樣，如果在matlab中已經導入了「jplv7」包的情況下再想使用matlab自帶的distance函數，就可能存在沖突。解決方法是，調正函數的優先級。具體操作是：主頁 > 設置路徑，選中jplv7的所有路徑，點擊「移至底端」、「保存」。
這樣再使用「distance」函數時，系統就會自動調用matlab本身的函數，而不是「jplv7」下的「distance」函數。
具體的計算原理：
圖中，空間坐標系 A(x1, y1), B(x2, y2)，我們需要計算的球面AOB距離為弧長。可以通過公式：半徑 × ∠AOB來完成。地球赤道半徑 6378.140 km，地球極地半徑 6356.755 km 。那么，現在面臨的問題就是∠AOB的計算了。
matlab本身自帶的「distance」函數就可以解決這個問題。這里，取赤道半徑6378.1km，下面展示計算過程：
>> D = distance(23.0, 101.1, 23.06, 113.34) D = 11.2613 >> 11.2613/180*pi*6378.1 ans = 1.2536e+003
或者
>> D = distance(23.0, 101.1, 23.06, 113.34, 6378.1) D = 1.2536e+03
可以看到。兩種方式計算的結果是一致的。
下面，我們將演示怎樣用matlab自帶的這個函數進行空間權重矩陣的構建。首先，使用ArcGis導出地理坐標系下的經緯度數據（單位為 °），排好序后以列向量的格式導入matlab。具體過程不再演示，相關經緯度數據可通過后台回復「地理經緯度」獲取，下面直接貼出代碼內容：
% 球面距離空間權重矩陣
E = zeros(33, 33);
for i = 1:32;
for j = 1 + i:33;
E(i, j) = E(i, j) + distance(X(i), Y(i), X(j), Y(j), 6378.1);
end;
end;
wd= E + E';
下面給出了，基於赤道半徑的球面距離和歐幾里得距離法下的空間權重矩陣（這兩份數據也已上傳到百度雲，回復可獲得）的部分結果：
球面距離空間權重矩陣
歐幾里得距離空間權重矩陣
可見，兩者存在一定差距。畢竟，地球並不是標准的球體，以赤道半徑作為地球半徑本身存在一定誤差。所以，在沒有特殊要求的情況下，還是建議使用投影后的平面坐標來計算地理距離。
更多精彩內容，掃碼關注微信公眾號！