《Python地理空间分析指南第2版》学习笔记-5.1 距离测量

本文链接：https://blog.csdn.net/qq_32390983/article/details/123672731

第5章 Python与地理信息系统

本章主要学习Python处理矢量数据，包含以下内容：

距离测量
坐标转换
矢量数据重投影
Shapefile 文件编辑
海量数据过滤
专题地图创建
非GIS数据类型转换
地理化编码

本篇博文记录第5.1距离测量

5.1 距离测量

地理学第一定律：“任何事物都相关，相近的事物关联更紧密”。因此距离测量显得格外重要。
地球是不规则椭圆状，计算距离时有3种地球模型可供选择：

- 平面

地球被当做没有曲率的平面，只需处理直线即可，可以利用地图投影把地球模型压扁转换成二维平面，利用十进制度数来表示，十进制度数系统会被当做包含x和y的笛卡尔直角坐标系，测量距离则就可以很方便的使用欧式几何的勾股定理了。

但是求得距离的误差比较大。
在这里插入图片描述

- 球体
使用一个完美的球体来表达地球的特征，但是地球实际上是一个不规则的椭球体
在这里插入图片描述

- 椭球体

存在多种椭球体地球模型基准，其中一套基准就是一组地球形状定义的集合，也被称为大地测量坐标系统。和其他地理坐标系统类似，基准也可以对特定区域进行优化。最常用的基准是基于全球的***WGS84基准***（会优化，2004年是最新的版本），比正圆形球体模型更符合实际。

补充：大地水准面（大地水平面模型）：是指与平均海水面重合并延伸到大陆内部的水准面。

三者对比差异

5.1.1勾股定理（也叫欧氏距离公式）——平面地球模型

例如两个点之间，计算距离
在这里插入图片描述

a2+b2=c2

import math
# 两个投影点，经密西西比横轴墨卡托投影
x1 = 456456.23
y1 = 1279721.064
x2 = 576628.34
y2 = 1071740.33

x_dist = x1 - x2
y_dist = y1 - y2

dist_sq = x_dist**2 + y_dist**2

dist = math.sqrt(dist_sq)

print(dist)

240202.668047278
大该是240202m，也就是240.2km

使用角度进行测量，必须将角度转换为弧，才能在坐标系统之间计算曲面上的距离。


import math
#经纬度坐标
x1 = -90.21
y1 = 32.31
x2 = -88.95
y2 = 30.43

x_dist = math.radians(x1 - x2)
y_dist = math.radians(y1 - y2)
dist_squre = x_dist**2 + y_dist**2
dis_rad = math.sqrt(dist_squre)
# 6371251.46m 地球半径
print("弧度是：%f"% dis_rad)
print("距离是：%f"% (dis_rad*6371251.46))

弧度是：0.039500
距离是：251664.466877，大概是251km，比第一次结果大了11km
因此不同测量方式（不用的坐标系统）和地球模型的结果可能存在很大差异

5.1.2 半正矢量公式——球体地球模型

大圆距离（也叫球面距离）：指球体表面上两点之间最短距离（根据两点的经度和纬度来确定大圆上两点之间距离的计算方法）。

球面上任意两点A和B以及球心都可以确定唯一的大圆，这个大圆被称为黎曼圆，而在大圆上连接这两点的较短的一条弧的长度就是大圆距离。

若这两点和球心正好都在球的直径上，则过这三点可以有无数大圆，但两点之间的
弧长都相等，因此该大圆能够将球体等分。

最佳的十进制度数距离测量方法：半正矢量公式。
半正矢公式：根据经度纬度计算两点间距离

在这里插入图片描述

import math
x1 = -90.212452861859035
y1 = 32.316272202663704
x2 = -88.952170968942525
y2 = 30.438559624660321
# Convert angle x from degrees to radians
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)**2 * math.cos(y1_rad) * math.cos(y2_rad)
c = 2 * math.asin(math.sqrt(a))
# 6371251.46m 地球半径
# 弧长 = 弧度*半径
distance = c * 6371 # kilometers
print(distance)

240.6359762909508
半正矢公式计算结果是240.6km，与使用优化过的重投影方法获得的240.2km相比，它们之间的误差小于0.5km。
半正矢公式是最常用的距离计算公式，因为它相对来说代码量较少并且大部分情况下足够精确，它的精度范围能够控制在1m以内

目前遇到的大部分坐标点都是未经投影变换的十进制度数格式。所以，从距离
测量的角度来看，有以下几个可供选择的方案：

重投影到一个精确的笛卡儿坐标系统中并计算距离；
根据你的精度要求只使用半正矢公式；
使用精度更高的Vincenty公式。

5.1.3Vincenty——椭球体地球模型

Vincenty公式是基于椭球体地球模型的,如果你使用的是一个本地化的椭球体模型，
那么它的计算精度可以远远小于1m。
可以自定义半长轴和扁平率来匹配任何权威的椭球体模型。


import math
distance = None
x1 = -90.212452861859035
y1 = 32.316272202663704
x2 = -88.952170968942525
y2 = 30.438559624660321
#定义椭球体参数，例如NAD83
a = 6378137 #半长轴
f = 1/298.257222101 #逆扁平度
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)
    if (sinSigma == 0):
        # 重合点
        distance = 0
        break
    cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLam
    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)))
    s = b * A * (sigma - deltaSigma)
    distance = s
    print(distance)

240091.4562741062
240091.4562741062
240091.4562741062
240091.4562741062
240091.4562741062
240091.4562741062
…
Vincenty公式的计算结果大约是240.1km，这个数据和之前经投影变换之后的欧氏距离测量公式得到的结果的误差只有100m左右。

完全使用Python实现的geopy模块包含一个Vincenty公式的实现，并且它还支持地理位置编码，能够将地名转换为经纬度坐标。它的详情可以参考如下页
面：http://geopy.readthedocs.org/en/latest/。