编程软件免费中文版_军事要闻军事频道_权威解读当前经济热点问题_营销案例

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

1 转换算法

import mathdef 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.0ret += (20.0 * math.sin(x * math.pi) + 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0return retdef 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.0ret += (20.0 * math.sin(y * math.pi) + 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0return retdlat = transform_lat(lon - 105.0, lat - 35.0)dlon = transform_lon(lon - 105.0, lat - 35.0)radlat = lat / 180.0 * math.pimagic = math.sin(radlat)magic = 1 - ee * magic * magicsqrtmagic = 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 - dlatwgs_lon = lon - dlonreturn wgs_lon, wgs_lat

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

import arcpy
import mathdef 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.0ret += (20.0 * math.sin(x * math.pi) + 40.0 * math.sin(x / 3.0 * math.pi)) * 2.0 / 3.0ret += (150.0 * math.sin(x / 12.0 * math.pi) + 300.0 * math.sin(x / 30.0 * math.pi)) * 2.0 / 3.0return retdef 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.0ret += (20.0 * math.sin(y * math.pi) + 40.0 * math.sin(y / 3.0 * math.pi)) * 2.0 / 3.0ret += (160.0 * math.sin(y / 12.0 * math.pi) + 320 * math.sin(y * math.pi / 30.0)) * 2.0 / 3.0return retdlat = transform_lat(lon - 105.0, lat - 35.0)dlon = transform_lon(lon - 105.0, lat - 35.0)radlat = lat / 180.0 * math.pimagic = math.sin(radlat)magic = 1 - ee * magic * magicsqrtmagic = 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 - dlatwgs_lon = lon - dlonreturn 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.Xy = 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("坐标转换完成！")