高德火星坐标(GCJ-02)转WGS84坐标

高德火星坐标(GCJ-02)转WGS84坐标

1 转换算法

import math

def gcj02_to_wgs84(lon, lat):
    """
    高德火星坐标(GCJ-02)转WGS84坐标
    """
    a = 6378245.0  # 长半轴
    ee = 0.00669342162296594323  # 扁率
    
    def transform_lon(x, y):
        ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * math.sqrt(abs(x))
        ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(x * math.pi) + 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0
        ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0
        return ret
    
    def transform_lat(x, y):
        ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * math.sqrt(abs(x))
        ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(y * math.pi) + 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0
        ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0
        return ret

    dlat = transform_lat(lon - 105.0, lat - 35.0)
    dlon = transform_lon(lon - 105.0, lat - 35.0)
    
    radlat = lat / 180.0 * math.pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * math.pi)
    dlon = (dlon * 180.0) / (a / sqrtmagic * math.cos(radlat) * math.pi)
    
    wgs_lat = lat - dlat
    wgs_lon = lon - dlon
    
    return wgs_lon, wgs_lat

2 使用Arcpy实现点要素类坐标转换

import arcpy
import math

def gcj02_to_wgs84(lon, lat):
    """
    高德火星坐标(GCJ-02)转WGS84坐标
    """
    a = 6378245.0  # 长半轴
    ee = 0.00669342162296594323  # 扁率
    
    def transform_lon(x, y):
        ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * math.sqrt(abs(x))
        ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(x * math.pi) + 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0
        ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0
        return ret
    
    def transform_lat(x, y):
        ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * math.sqrt(abs(x))
        ret += (20.0 * math.sin(6.0 * x * math.pi) + 20.0 * math.sin(2.0 * x * math.pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(y * math.pi) + 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0
        ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0
        return ret

    dlat = transform_lat(lon - 105.0, lat - 35.0)
    dlon = transform_lon(lon - 105.0, lat - 35.0)
    
    radlat = lat / 180.0 * math.pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * math.pi)
    dlon = (dlon * 180.0) / (a / sqrtmagic * math.cos(radlat) * math.pi)
    
    wgs_lat = lat - dlat
    wgs_lon = lon - dlon
    
    return wgs_lon, wgs_lat

# 输入参数设置
input_fc = r"Point"  # 替换为你的输入要素类路径
output_fc = r"Point2"  # 替换为你的输出要素类路径
output_sr = arcpy.SpatialReference(4490)  # CGCS2000地理坐标系

ws = r"D:\admin\Documents\Projects\tj\tj.gdb"

arcpy.env.workspace = ws

# 创建输出要素类
arcpy.management.CreateFeatureclass(
    out_path=arcpy.env.workspace,
    out_name=output_fc,
    geometry_type="POINT",
    spatial_reference=output_sr
)

# 添加字段（根据输入要素字段需要调整）
arcpy.management.AddField(output_fc, "NAME", "TEXT")

# 开始编辑
with arcpy.da.Editor(arcpy.env.workspace) as edit:
    # 使用插入游标写入新要素
    with arcpy.da.InsertCursor(output_fc, ["SHAPE@", "NAME"]) as insert_cursor:
        # 使用搜索游标读取原始要素
        with arcpy.da.SearchCursor(input_fc, ["SHAPE@", "NAME"]) as cursor:
            for row in cursor:
                # 获取原始坐标（假设输入是点要素）
                original_point = row[0]
                x = original_point.firstPoint.X
                y = original_point.firstPoint.Y
                
                # 执行坐标转换
                converted_x, converted_y = gcj02_to_wgs84(x, y)
                
                # 创建新点对象
                new_point = arcpy.Point(converted_x, converted_y)
                new_geometry = arcpy.PointGeometry(new_point, output_sr)
                
                # 写入新要素
                insert_cursor.insertRow([new_geometry, row[1]])

print("坐标转换完成！")
ut_sr)
                
                # 写入新要素
                insert_cursor.insertRow([new_geometry, row[1]])

print("坐标转换完成！")