计算地球表面两点之间的距离

根据经度度，计算两点之间的距离

def geodistance(lng1, lat1, lng2, lat2):
    # lng1,lat1,lng2,lat2 = (120.12802999999997,30.28708,115.86572000000001,28.7427)
    lng1, lat1, lng2, lat2 = map(radians, [float(lng1), float(lat1), float(lng2), float(lat2)])  # 经纬度转换成弧度
    dlon = lng2 - lng1
    dlat = lat2 - lat1
    a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
    distance = 2 * asin(sqrt(a)) * 6371 * 1000  # 地球平均半径，6371km
    distance = round(distance / 1000, 3)
    return distance

lng1, lat1, lng2, lat2 = (120.12802999999997, 30.28708, 115.86572000000001, 28.7427)
print(geodistance(lng1, lat1, lng2, lat2))
# 返回 446.721 千米

利用布尔莎七参数实现坐标转换-基准面转换

在坐标转换中，除了正投影和反投影的转换之外，还有不同基准面之间的转换。基准面的转换有很多种转换模型，常见的有三参数和七参数转换。三参数的转换主要是通过对x，y，z三个坐标轴进行平移操作，适用于小范围的地方坐标基准之间的转换，对精度要求不高。相对于三参数转换，七参数转换在基基础上增加对x/y/z三个坐标轴的旋转分量以及缩放比例。
我们对一个经纬度坐标进行七参数转换时，一般先将其转换到地心坐标系（以地心为原点的三维坐标系），再利用空间直角坐标系的转换公式，使用七参数对其进行转换，将转换的结果再转回经纬度坐标，即可完成椭球基准的转换。

1. 经纬度坐标转地心空间直角坐标
   转换公式：
   $$\{\begin{array}{rl}& X=(N+H)\cos B\cos L\\ & Y=(N+H)\cos B\sin L\\ & Z=(N\ast (1-c^2)+H)\sin B\end{array}$$
   其中
   $$\{\begin{array}{rl}& AN=\frac{a}{W}\\ & W=\sqrt{1-e^2{\sin }^{2}B}\\ & e^2=\frac{a^2-b^2}{a^2}\end{array}$$

2. 基准转换

七参数基准转换的原理同三维坐标系之间的坐标转换原理类似，这里给出公式如下：
$$\left[\begin{array}{c}{X}_B\\ Y_B\\ Z_B\end{array}\right]=\left[\begin{array}{c}△X_0\\ △Y_0\\ △Z_0\end{array}\right]+(1+m)\left[\begin{array}{ccc}1& {\omega }_Z& -{\omega }_Y\\ -{\omega }_Z& 1& {\omega }_X\\ {\omega }_Y& -{\omega }_X& 1\end{array}\right]\left[\begin{array}{c}{X}_A\\ Y_A\\ Z_A\end{array}\right]$$
其中，$△X$，$△Y$, $△Z$为平移分量，$m$为缩放分量，$m$为缩放系数，$\omega$为旋转矩阵。

3. 地心空间直角坐标转经纬度
   将地心空间坐标转换为经纬度坐标，公式如下：
   $$\{\begin{array}{rl}& L=arctan(\frac{Y}{X})\\ & B=arctan(\frac{Z(N+H)}{\sqrt{(X^2+Y^2)[N(1-e^2)+H]}})\\ & H=\frac{Z}{\sin B}-N(1-e^2)\end{array}$$