经纬度和平面坐标互转

经纬度和平面坐标互转

Fre_soe 52 2025-05-30
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}")