空间权重矩阵的那些事(八)-球面距离权重矩阵
空间权重矩阵的那些事(八)-球面距离权重矩阵
前段时间有人向我咨询了根据经纬度计算球面距离的方法,希望我出一篇文章,所以就有了这篇文章。必要文件可通过后台回复「地理经纬度」获取。
我首先想到的是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‘;
下面给出了,基于赤道半径的球面距离和欧几里得距离法下的空间权重矩阵(这两份数据也已上传到百度云,回复可获得)的部分结果:
球面距离空间权重矩阵
欧几里得距离空间权重矩阵
可见,两者存在一定差距。毕竟,地球并不是标准的球体,以赤道半径作为地球半径本身存在一定误差。所以,在没有特殊要求的情况下,还是建议使用投影后的平面坐标来计算地理距离。
更多精彩内容,扫码关注微信公众号!