经纬度和平面坐标互转

经纬度和平面坐标互转

Administrator 6 2025-07-03
import pyproj

def convert_coordinates(x, y):

    # 从 x 中提取带号和实际 x 坐标

    start_x = int(str(x)[:2])

    end_x = float(str(x)[2:])

    

    # 计算中央经线

    central_meridian = start_x * 3

    

    # 创建从高斯-克吕格到 WGS84 的转换

    gk_to_wgs84 = pyproj.Transformer.from_proj(

        pyproj.Proj(proj='tmerc', ellps='WGS84', lon_0=central_meridian, x_0=500000, y_0=0),

        pyproj.Proj(proj='longlat', ellps='WGS84', datum='WGS84')

    )

    

    # 执行转换

    lon, lat = gk_to_wgs84.transform(end_x, y)

    

    return lon, lat

def convert_coordinates_reverse(lon, lat, zone=None):

    """

    将经纬度坐标转换回高斯-克吕格平面坐标

    

    参数:

        lon: 经度

        lat: 纬度

        zone: 可选,带号。如果不提供,将根据经度自动计算

    

    返回:

        x: 带带号的高斯-克吕格 X 坐标

        y: 高斯-克吕格 Y 坐标

    """

    # 如果未提供带号,根据经度计算带号(3度带)

    if zone is None:

        zone = int(lon // 3) + 1

    

    # 计算中央经线

    central_meridian = zone * 3

    

    # 创建从 WGS84 到高斯-克吕格的转换

    wgs84_to_gk = pyproj.Transformer.from_proj(

        pyproj.Proj(proj='longlat', ellps='WGS84', datum='WGS84'),

        pyproj.Proj(proj='tmerc', ellps='WGS84', lon_0=central_meridian, x_0=500000, y_0=0)

    )

    

    # 执行转换

    x, y = wgs84_to_gk.transform(lon, lat)

    

    # 添加带号前缀到 x 坐标

    x_with_zone = zone * 1000000 + x

    

    return x_with_zone, y

if name == "__main__":

    # 测试原始函数

    x, y = 32402655.65, 4514313.5

    lon, lat = convert_coordinates(x, y)

    print(f"原始坐标: {x}, {y}")

    print(f"转换为经纬度: {lon}, {lat}")

    # 测试反转函数

    x_reverse, y_reverse = convert_coordinates_reverse(lon, lat)

    print(f"反转回高斯-克吕格: {x_reverse}, {y_reverse}")

    print(f"与原始坐标的差异: {x - x_reverse}, {y - y_reverse}")