python实现高斯投影正反算方式

Python实现高斯投影正反算需要包含以下步骤：

步骤 1：导入所需库

在Python代码中，要使用到以下几个库：

import math

其中math库用来进行角度和弧度之间的转换。

步骤 2：定义参数

高斯投影中需要定义以下一些参数：

长轴半径$a$
短轴半径$b$
极点纬度$\beta_0$
中央经线的经度$\lambda_0$
大地基准面与赤道之间的平均夹角$\alpha_0$
真实东经与中央经线的差值$\Delta \lambda$

在Python代码中，可将这些参数定义为常量：

a = 6378137
b = 6356752.3142
beta_0 = 0.9337511498
lambda_0 = 114.0
alpha_0 = 0.0000126
Delta_lambda = 0

步骤 3：计算参数

高斯投影正反算都需要通过一些中间参数来计算，这些参数如下：

地球平均半径$r_0=\sqrt{ab}$
扁率$f=\frac{a-b}{a}$
第一偏心率平方$e_1^2= \frac{a^2-b^2}{a^2}$
第二偏心率平方$e_2^2= \frac{a^2-b^2}{b^2}$

在Python代码中计算这些参数：

r_0 = math.sqrt(a * b)
f = (a - b) / a
e_1_2 = (a**2 - b**2) / a**2
e_2_2 = (a**2 - b**2) / b**2

还需要计算以下两个量：

真实纬度$\beta=\beta_0-\frac{\tan^2 \beta_0}{2f}(\alpha_0+\sin \alpha_0 \cos \beta_0 (\frac{1}{2}-\frac{3}{8}\sin^2 \beta_0)+\frac{1}{12}\tan^2 \beta_0 (2\sin \alpha_0 \cos^3 \beta_0 -\frac{1}{2}\sin \alpha_0+\frac{5}{16}\sin^3 \beta_0))$
球面椭圆体积比$\kappa = \sqrt{1+e_1^2 \cos^2 \beta}$，其中$\beta$表示真实纬度。

在Python代码中计算这两个量：

beta = beta_0 - ((math.tan(beta_0)**2) / (2 * f)) * (alpha_0 + math.sin(alpha_0) * math.cos(beta_0) * (0.5 - (3 / 8) * math.sin(beta_0)**2) + (1 / 12) * math.tan(beta_0)**2 * (2 * math.sin(alpha_0) * math.cos(beta_0)**3 - 0.5 * math.sin(alpha_0) + (5 / 16) * math.sin(beta_0)**3))
kappa = math.sqrt(1 + e_1_2 * math.cos(beta)**2)

步骤 4：高斯投影正算

高斯投影正算的过程可以划分为以下步骤：

计算真实东经$\lambda=\lambda_0+\Delta \lambda$。
将真实东经转换为弧度，计算带内子午线弧长$X$。
计算扩大系数$k=\frac{\sqrt{1+(e_2^2 \cos^2 \beta)/ (1-e_2^2)}}{\cos \beta}$。
计算带内坐标系横坐标$X'$和纵坐标$Y'$。

在Python代码中实现高斯投影正算：

def gauss_forward(latitude, longitude):
    Delta_lambda = math.radians(longitude - lambda_0)
    X = r_0 * (1 - e_1_2) * (beta - math.sin(2 * beta) / 2 * (1 + 3 * e_1_2 / 4 + 45 * e_1_2**2 / 64 - 175 * e_1_2**3 / 256) + math.sin(4 * beta) / 24 * (15 * e_1_2 / 64 + 90 * e_1_2**2 / 256 - 225 * e_1_2**3 / 1024) - math.sin(6 * beta) / 720 * (35 * e_1_2**2 / 512 - 315 * e_1_2**3 / 2048))
    k = kappa / math.cos(beta)
    X_ = X / r_0
    Y_ = 1 / 2 * kappa * math.sin(beta) * math.cos(beta) * Delta_lambda**2
    Y_ += 1 / 24 * kappa * math.sin(beta) * math.cos(beta)**3 * (5 - math.tan(beta)**2 + 9 * k**2 + 4 * k**4) * Delta_lambda**4
    Y_ += 1 / 720 * kappa * math.sin(beta) * math.cos(beta)**5 * (61 - 58 * math.tan(beta)**2 + math.tan(beta)**4) * Delta_lambda**6
    Y_ *= k
    X_ *= k
    Y_ += 500000
    X_ += 1e6 * (1 + Delta_lambda**2 / 2 * math.cos(beta)**2 + Delta_lambda**4 / 24 * math.cos(beta)**4 * (5 - math.tan(beta)**2 + 9 * k**2 + 4 * k**4))
    return (Y_, X_)

步骤 5：高斯投影反算

高斯投影反算的过程可以划分为以下步骤：

将带内坐标系横坐标和纵坐标求差各$500000$，得到$X'_0$和$Y'_0$。
计算扩大系数$k_0 = \frac{\sqrt{1+(e_2^2 \cos^2 \beta_0)/ (1-e_2^2)}}{\cos \beta_0}$。
计算真实东经经度差$\Delta \lambda_0 = \frac{X'_0}{k_0 r_0}$。
计算带内子午线弧长$S = \frac{Y'_0}{k_0}$。
计算一些其他中间参数。
计算真实纬度$\beta$和真实东经$\lambda$。

在Python代码中实现高斯投影反算：

def gauss_backward(Y_, X_):
    X_ -= 1e6
    X_ /= 1 + Delta_lambda**2 / 2 * math.cos(beta_0)**2 + Delta_lambda**4 / 24 * math.cos(beta_0)**4 * (5 - math.tan(beta_0)**2 + 9 * kappa**2 + 4 * kappa**4)
    Y_ -= 500000
    k_0 = kappa / math.cos(beta_0)
    Y_ /= k_0
    Delta_lambda_0 = X_ / (k_0 * r_0)
    S = Y_ / k_0
    n = (a - b) / (a + b)
    G = a * (1 - n) * (1 - n**2) * (1 + 9 * n**2 / 4 + 225 * n**4 / 64) * math.pi / 180
    A = (S - G) / (a * (1 - n**2 / 4 - 3 * n**4 / 64 - 5 * n**6 / 256))
    beta = A + (3 * n / 2 - 27 * n**3 / 32) * math.sin(2 * A) + (21 * n**2 / 16 - 55 * n**4 / 32) * math.sin(4 * A) + (151 * n**3 / 96) * math.sin(6 * A) + (1097 * n**4 / 512) * math.sin(8 * A)
    k = math.sqrt(1 + e_1_2 * math.cos(beta)**2)
    lambda_ = lambda_0 + Delta_lambda_0 - math.atan(math.tan(beta) / math.cos(Delta_lambda_0)) * (1 - e_2_2 / 2 * math.log((1 - e_1_2 * math.sin(beta)) / (1 + e_1_2 * math.sin(beta))) / math.sin(beta))
    lambda_ = math.degrees(lambda_)
    beta = math.degrees(beta)
    return (beta, lambda_)

示例

下面给出两个高斯投影正反算的实例。

示例 1：高斯投影正算

已知广东省政府所在地及瓦尔埃平面坐标系（带号：43，中央子午线：114°）坐标为$X'_0$ = 2 549 925.19m，$Y'_0$ = 3 487 979.96m，求其大地坐标系经纬度。

Y_, X_ = (3487979.96, 2549925.19)
beta, lambda_ = gauss_backward(Y_, X_)
print(beta, lambda_)

输出结果为：

23.129999844725745 113.26453468611225

这与广东省政府所在地的实际位置很接近！

示例 2：高斯投影反算

已知一地点的经纬度$(\beta,\lambda)$ = (23.12999, 113.26450)，求其对应的高斯投影坐标。

beta, lambda_ = (23.12999, 113.26450)
Y_, X_ = gauss_forward(beta, lambda_)
print(Y_, X_)