使用米勒投影法实现经纬度和平面坐标的转换

import math

M_PI = 3.14159265358979323846
L = (6378.2 * M_PI * 2)  # 地球周长40076
W = L  # 平面展开后，x轴等于周长
H = L / 2  # y轴约等于周长一半
mill = 2.3  # 米勒投影中的一个常数，范围大约在正负2.3之间

position = [[120.881805, 29.517331],
            [120.928116, 29.509094],
            [120.80301, 29.450337],
            [120.913428, 29.502665],
            [120.047006, 29.450389],
            [120.822468, 29.280728],
            [120.126291, 29.384722],
            [120.139216, 29.519882],
            [120.792824, 29.375422],
            [120.938022, 29.368695],
            [120.915956, 29.470399],
            [120.841442, 29.427925],
            ]


# 经纬度转换坐标
def LALtoCOORD(position):
    n = len(position)
    new_position = [[0 for i in range(2)] for i in range(n)]
    for i in range(n):
        lon = position[i][0]
        lat = position[i][1]
        x = lon * M_PI / 180  # 将经度从度数转换为弧度
        y = lat * M_PI / 180  # 将纬度从度数转换为弧度
        y1 = 1.25 * math.log(math.tan(0.25 * M_PI + 0.4 * y))  # 米勒投影的转换
        # 弧度转为实际距离
        dikaerX = (W / 2) + (W / (2 * M_PI)) * x  # 笛卡尔坐标x
        dikaerY = (H / 2) - (H / (2 * mill)) * y1  # 笛卡尔坐标y
        new_position[i][0] = dikaerX
        new_position[i][1] = dikaerY
        print("第{i}个点的坐标是({x}, {y})".format(i=i + 1, x=new_position[i][0], y=new_position[i][1]))
    return new_position


# 坐标转换经纬度

def COORDtoLAL(new_position):
    n = len(new_position)
    new_LAL = [[0 for i in range(2)] for i in range(n)]
    for i in range(n):
        dikaerX = new_position[i][0]
        dikaerY = new_position[i][1]
        x = (dikaerX - (W / 2)) / (W / (2 * M_PI))
        y1 = ((H / 2) - dikaerY) / (H / (2 * mill))
        y = (math.atan(math.exp(y1 / 1.25)) - (0.25 * M_PI)) / 0.4
        lon = x / M_PI * 180
        lat = y / M_PI * 180
        new_LAL[i][0] = lon
        new_LAL[i][1] = lat

        print("第{i}个点的经纬度是({x}, {y})".format(i=i + 1, x=new_LAL[i][0], y=new_LAL[i][1]))
    return new_LAL


arr = LALtoCOORD(position)
print("-------------------------------------------")
# arr = [[33497.8285402217,7709.64897178222]]
COORDtoLAL(arr)